CALRCH.f

Go to the documentation of this file.

      SUBROUTINE calrch(EL,RICHNO)

!
      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,rdzkl,             &
           dthvkl,dukl,dvkl,ri,ct,cs
!     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.
!     
!$omp  parallel do
      DO l = 1,lm
        DO j=jsta,jend
        DO i=1,im
          richno(i,j,l)=spval
        ENDDO
        ENDDO
      ENDDO
!
!     COMPUTE VIRTUAL POTENTIAL TEMPERATURE.
!
!$omp  parallel do private(i,j,ape)
      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.
!
!!$omp  parallel do private(i,j,l,ie,iw,cs,ct,dthvkl,dukl,dvkl,             &
!!$omp&         rdzkl,ri,uhkl,ulkl,vhkl,vlkl,wndsl,wndslp)
      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 == 0. .OR. wndslp == 0.) cycle
              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

            rdzkl   = 1.0 / (zmid(i,j,l)-zmid(i,j,l+1))

!           Q2KL   = MAX(Q2(I,J,L),0.00001)
!           QROOT  = SQRT(Q2KL)
!           ELKL   = EL(I,J,L)
!           ELKL   = MAX(ELKL,EPSQ2)
!           ELKLSQ = ELKL*ELKL

            dthvkl = thv(i,j,l)-thv(i,j,l+1)
            dukl   = (uhkl-ulkl) * rdzkl
            dvkl   = (vhkl-vlkl) * rdzkl
            cs     = dukl*dukl + dvkl*dvkl
!     
!         COMPUTE GRADIENT RICHARDSON NUMBER.
!     
            IF(cs <= 1.e-8) THEN
!
!            WIND SHEAR IS VANISHINGLY SMALL - 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
!
          ENDDO
        ENDDO
      ENDDO           ! end of l loop
!     
      DEALLOCATE (thv)
!     END OF ROUTINE.
!     
      RETURN
      END