subroutine rich_no (drhodz,dudzsq,mask,ri)
c
c=======================================================================
c                                                                    ===
c  This routine computes the Richardson number from the given        ===
c  vertical density gradient and the square of the vertical shear.   ===
c                                                                    ===
c  ------                                                            ===
c  Input:                                                            ===
c  ------                                                            ===
c                                                                    ===
c     DRHODZ   Vertical density gradient.  (real array)              ===
c     DUDZSQ   Square of vertical shear.   (real array)              ===
c     MASK     Masking array.              (real array)              ===
c                                                                    ===
c  Common Blocks:                                                    ===
c                                                                    ===
c  /SCALAR/                                                          ===
c                                                                    ===
c     GRAV   acceleration due to gravity.   (real; cm/s^2)           ===
c                                                                    ===
c  -------                                                           ===
c  Output:                                                           ===
c  -------                                                           ===
c                                                                    ===
c     RI       Richardson number.          (real array)              ===
c                                                                    ===
c  Calls:  none                                                      ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <scalar.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      FLOAT
     &      cutoff
      parameter (cutoff=c1000)
c
      integer i,k,kp1
      FLOAT
     &      denom,numer
      FLOAT
     &      drhodz(imt,km),dudzsq(imt,km),mask(imt,km),ri(imt,km)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Compute Richardson number with protection against division by zero.
c-----------------------------------------------------------------------
c
      do 10 k=1,kmm1
        kp1=k+1
        do 10 i=1,imtm1
          if (mask(i,kp1).gt.c0) then
            numer = -grav*drhodz(i,k)
            denom = dudzsq(i,k)
            if (abs(numer).lt.(cutoff*denom)) then
              ri(i,k)= numer/denom
             else
              ri(i,k)= sign(cutoff,numer)
            end if
          else
            ri(i,k)=c0
          endif
  10  continue
c
# ifdef cyclic
c-----------------------------------------------------------------------
c  Set cyclic boundary conditions.
c-----------------------------------------------------------------------
c
      do 20 k=1,kmm1
        ri(1  ,k) = ri(imtm1,k)
        ri(imt,k) = ri(2    ,k)
  20  continue
c
# endif
      return
      end