SUBROUTINE CALRCH(EL,RICHNO)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .     
! SUBPROGRAM:    CALRCH      COMPUTES GRD RCH NUMBER
!   PRGRMMR: TREADON         ORG: W/NP2      DATE: 93-10-11
!     
! ABSTRACT:  
!   THIS ROUTINE COMPUTES THE GRADIENT RICHARDSON NUMBER
!   AS CODED IN ETA MODEL SUBROUTINE PROFQ2.F.
!   FIX TO AVOID UNREASONABLY SMALL ANEMOMETER LEVEL WINDS.
!   .     
!     
! PROGRAM HISTORY LOG:
!   93-10-11  RUSS TREADON
!   98-06-17  T BLACK - CONVERSION FROM 1-D TO 2-D
!   00-01-04  JIM TUCCILLO - MPI VERSION
!   01-10-22  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
!   02-01-15  MIKE BALDWIN - WRF VERSION
!   05-02-25  H CHUANG - ADD COMPUTATION FOR NMM E GRID
!   05-07-07  BINBIN ZHOU - ADD RSM FOR A GRID  
!   
! USAGE:    CALL CALRCH(EL,RICHNO)
!   INPUT ARGUMENT LIST:
!     EL      - MIXING LENGTH SCALE.
!
!   OUTPUT ARGUMENT LIST: 
!     RICHNO  - GRADIENT RICHARDSON NUMBER.
!     
!   OUTPUT FILES:
!     NONE
!     
!   SUBPROGRAMS CALLED:
!     UTILITIES:
!       NONE
!     LIBRARY:
!       COMMON   - 
!                  CTLBLK
!     
!   ATTRIBUTES:
!     LANGUAGE: FORTRAN
!     MACHINE : CRAY C-90
!$$$  
!
      use vrbls3d, only: pmid, q, t, uh, vh, zmid, q2
      use masks, only: vtm
      use params_mod, only: h10e5, capa, d608,h1, epsq2, g, beta
      use ctlblk_mod, only: jsta, jend, spval, lm1, jsta_m, jend_m, &
          im, jsta_2l, jend_2u, lm
      use gridspec_mod, only: gridtype
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!     
!     DECLARE VARIABLES.
!     
      REAL,intent(in)::  EL(IM,jsta_2l:jend_2u,LM)
      REAL,intent(inout)::  RICHNO(IM,jsta_2l:jend_2u,LM)
!
      REAL, ALLOCATABLE :: THV(:,:,:)
      integer I,J,L,IW,IE
      real APE,UHKL,ULKL,VHKL,VLKL,WNDSL,WNDSLP,DZKL,RDZKL,Q2KL,QROOT, &
           ELKL, ELKLSQ,DTHVKL,DUKL,DVKL,RI,CT,CS
!
!     
!*************************************************************************
!     START CALRCH HERE.
!     
      ALLOCATE ( THV(IM,JSTA_2L:JEND_2U,LM) )
!     INITIALIZE ARRAYS.
!     
      DO L = 1,LM
        DO J=JSTA,JEND
        DO I=1,IM
          RICHNO(I,J,L)=SPVAL
        ENDDO
        ENDDO
      ENDDO
!
!     COMPUTE VIRTUAL POTENTIAL TEMPERATURE.
!
      DO L=LM,1,-1
        DO J=JSTA,JEND
        DO I=1,IM
          APE=(H10E5/PMID(I,J,L))**CAPA
          THV(I,J,L)=(Q(I,J,L)*D608+H1)*T(I,J,L)*APE
        ENDDO
        ENDDO
      ENDDO
!
!     COMPUTE GRADIENT RICHARDSON NUMBER AS CODED IN ETA MODEL
!     SUBROUTINE PROFQ2.F.  OUTER LOOP OVER THE VERTICAL. 
!     INTTER LOOP OVER THE HORIZONTAL.
!
      DO L = 1,LM1
!
        if(GRIDTYPE /= 'A')THEN  
          call exch(VTM(1,jsta_2l,L))
          call exch(UH(1,jsta_2l,L))
	  call exch(VH(1,jsta_2l,L))
          call exch(VTM(1,jsta_2l,L+1))
          call exch(UH(1,jsta_2l,L+1))
          call exch(VH(1,jsta_2l,L+1))
	end if  
	         
        DO J=JSTA_M,JEND_M
        DO I=2,IM-1
!
          IF(GRIDTYPE == 'A')THEN
           UHKL=UH(I,J,L)
           ULKL=UH(I,J,L+1)
           VHKL=VH(I,J,L)
           VLKL=VH(I,J,L+1)
          ELSE IF(GRIDTYPE == 'E')THEN
           IE=I+MOD(J+1,2) 
           IW=I+MOD(J+1,2)-1
!
!         WE NEED (U,V) WINDS AT A MASS POINT.  FOUR POINT
!         AVERAGE (U,V) WINDS TO MASS POINT.  NORMALIZE FOUR
!         POINT AVERAGE BY THE ACTUAL NUMBER OF (U,V) WINDS
!         USED IN THE AVERAGING.  VTM=1 IF WIND POINT IS
!         ABOVE GROUND.  VTM=0 IF BELOW GROUND.
!
           WNDSL=VTM(I,J-1,L)+VTM(IW,J,L)+VTM(IE,J,L)+VTM(I,J+1,L)
           WNDSLP=VTM(I,J-1,L+1)+VTM(IW,J,L+1)+                       &
                 VTM(IE,J,L+1)+VTM(I,J+1,L+1)
           IF(WNDSL.EQ.0..OR.WNDSLP.EQ.0.)GO TO 10
           UHKL=(UH(I,J-1,L)+UH(IW,J,L)+UH(IE,J,L)+UH(I,J+1,L))/WNDSL
           ULKL=(UH(I,J-1,L+1)+UH(IW,J,L+1)+UH(IE,J,L+1)+             &
                UH(I,J+1,L+1))/WNDSLP
           VHKL=(VH(I,J-1,L)+VH(IW,J,L)+VH(IE,J,L)+VH(I,J+1,L))/WNDSL
           VLKL=(VH(I,J-1,L+1)+VH(IW,J,L+1)+VH(IE,J,L+1)+             &
                 VH(I,J+1,L+1))/WNDSLP
	  ELSE IF(GRIDTYPE == 'B')THEN
           IE=I 
           IW=I-1
	   UHKL=(UH(IW,J-1,L)+UH(IW,J,L)+UH(IE,J-1,L)+UH(I,J,L))/4.0
           ULKL=(UH(IW,J-1,L+1)+UH(IW,J,L+1)+UH(IE,J-1,L+1)+           &
                UH(I,J,L+1))/4.0
           VHKL=(VH(IW,J-1,L)+VH(IW,J,L)+VH(IE,J-1,L)+VH(I,J,L))/4.0
           VLKL=(VH(IW,J-1,L+1)+VH(IW,J,L+1)+VH(IE,J-1,L+1)+           &
                 VH(I,J,L+1))/4.0	 
          END IF

          DZKL=ZMID(I,J,L)-ZMID(I,J,L+1)
          RDZKL=1./DZKL
          Q2KL=AMAX1(Q2(I,J,L),0.00001)
          QROOT=SQRT(Q2KL)
          ELKL=EL(I,J,L)
          ELKL=AMAX1(ELKL,EPSQ2)
          ELKLSQ=ELKL*ELKL
          DTHVKL=THV(I,J,L)-THV(I,J,L+1)
          DUKL=(UHKL-ULKL)
          DVKL=(VHKL-VLKL)
          CS=(DUKL*RDZKL)**2+(DVKL*RDZKL)**2
!     
!         COMPUTE GRADIENT RICHARDSON NUMBER.
!     
          IF(CS.LE.1.E-8)THEN
!
!         WIND SHEAR IS VANISHINGLY SO SET RICHARDSON
!         NUMBER TO POST PROCESSOR SPECIAL VALUE.
!
            RICHNO(I,J,L)=SPVAL
!
          ELSE
!
!         WIND SHEAR LARGE ENOUGH TO USE RICHARDSON NUMBER.
!
            CT=-1.*G*BETA*DTHVKL*RDZKL
            RI=-CT/CS
            RICHNO(I,J,L)=RI
          ENDIF
!
 10      CONTINUE
        ENDDO
        ENDDO
      ENDDO
!     
      DEALLOCATE (THV)
!     END OF ROUTINE.
!     
      RETURN
      END