CLMAX.f

      SUBROUTINE clmax(EL0,SQZ,SQ,RQ2L,RQ2H)
!
!     CALCULATES THE FREE-ATMOSPHERE ASYMPTOTE OF THE TURBULENCE LENGTH
!     SCALE (L-INF IN THE BLACKADAR's FORMULA) FROM THE DISTRIBUTION
!     OF THE TURBULENT ENERGY (see MY82)
!
!     EXTRACTED FROM EXISTING CODE BY L. LOBOCKI, JULY 28, 1992
!   01-10-22  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
!   02-06-19  MIKE BALDWIN - WRF VERSION
!   21-07-26  W Meng - Restrict computation from undefined grids
!
!     INPUT:
!     ------
!
!     PINT (IM,jsta_2l:jend_2u,LP1)  - PRESSURE ON INTERFACES 
!     HTM  (IM,jsta_2l:jend_2u,LM)  - HEIGHT TOPOGRAPHY MASK ARRAY
!     Q2   (IM,jsta_2l:jend_2u,LM)  - TWICE THE TURBULENT ENERGY FIELD
!     ZINT (IM,jsta_2l:jend_2u,LP1) - ETA INTERFACES HEIGHT FIELD
!     SM   (IM,jsta_2l:jend_2u)     - SEA MASK
!     HGT  (IM,jsta_2l:jend_2u)     - SURFACE ELEVATION ARRAY
!     LMH  (IM,jsta_2l:jend_2u)     - TOPOGRAPHY INDEXES ARRAY
!
!     OUTPUT:
!     -------
!
!     EL0 (IM,JM)      - ARRAY OF RESULTING ASYMPTOTIC MIXING LENGTHS
!
!
!     SCRATCH AREAS:
!     --------------
!
!     SQZ(IM,JM),SQ(IM,JM),RQ2L(IM,JM),RQ2H(IM,JM)
!
!
!     RELEVANT CONSTANTS:
!     -------------------
!
!     PROPORTIONALITY CONSTANT BETWEEN ASYMPTOTIC MIXING LENGTH AND THE
!     S.D. OF Q DISTRIBUTION, FOR LAND AND SEA AREAS, CORRESPONDINGLY:

      use vrbls3d, only: zint, q2, pint
!      use vrbls2d, only:
      use masks, only: lmh, sm
      use params_mod, only: epsq2
      use ctlblk_mod, only: jsta, jend, lm, im, spval
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
      real,PARAMETER :: alphal=0.2, alphas=0.2
!
!     ASYMPTOTIC MIXING LENGTH LIMITATIONS:
      real,PARAMETER :: el0m=300.0, elmin=11.0
!
!     MINIMAL VALUE OF TURBULENT ENERGY:
!      real,PARAMETER :: EPSQ2=0.2
!
!     ------------------------------------------------------------------
!
      real,dimension(IM,jsta:jend),intent(inout) ::  sqz,sq,rq2l,rq2h,el0
      real,dimension(IM,jsta:jend)               ::   hgt
!jw
      integer i,j,l
      real dp, rq2m
!     ------------------------------------------------------------------
!
!
!$omp  parallel do
      DO j=jsta,jend
        DO i=1,im
          sqz(i,j)  = 0.0
          sq(i,j)   = 0.0
          rq2h(i,j) = 0.0
          hgt(i,j)  = zint(i,j,nint(lmh(i,j)))
        ENDDO
      ENDDO
!
      DO l=1,lm
        DO j=jsta,jend
          DO i=1,im
            IF(q2(i,j,l) <= epsq2) THEN
              rq2l(i,j) = 0.0
            ELSE
              rq2l(i,j) = sqrt(q2(i,j,l))
            ENDIF
!
!         -----------------------------------------------------------------
!         THIS PART OF THE CODE IS LEFT FOR TESTING OTHER PARAMETERIZATION
!         SCHEMES 
!
!         IF (L>=LMH(I,J)) GOTO 215
!         RQ2L(I,J)=SQRT(Q2(I,J,L))
!         IF(Q2(I,J,L)<0.0)THEN
!           write(3,*)'NEGATIVE Q2 AT (I,J,L)=(',I,',',J,',',L,'): ',
!                     Q2(I,J,L)
!           STOP
!         ENDIF
!         -----------------------------------------------------------------
!
            dp = pint(i,j,l+1) - pint(i,j,l)
!***
!***      SUM OF Q2 AT BOTH LOWER & UPPER SURFACES:
!***
            rq2m = rq2h(i,j) + rq2l(i,j)
!***
!***      INTEGRAL OF Q*Z OVER DP
!***
            sqz(i,j) = ((zint(i,j,l)+zint(i,j,l+1))*0.5-hgt(i,j))*rq2m*dp    &
     &               +  sqz(i,j)
!***
!***      INTEGRAL OF Q OVER DP:
!***
            sq(i,j)   = rq2m*dp + sq(i,j)
            rq2h(i,j) = rq2l(i,j)
          ENDDO
        ENDDO
!215    CONTINUE
      ENDDO
!***
!***    CLIPPING & APPLYING DIFFERENT VALUES OF THE PROPORTIONALITY 
!***    CONSTANT ALPHA FOR THE LAND AND SEA AREA:
!***
!$omp  parallel do
      DO j=jsta,jend
        DO i=1,im
          IF(hgt(i,j)<spval)THEN
          el0(i,j)= max(min(                                           &
     &       ((sm(i,j)*alphas+(1.0-sm(i,j))*alphal)*sqz(i,j)           &
     &       /(sq(i,j)+epsq2)),el0m),elmin)
          ELSE
          el0(i,j)= spval
          ENDIF
        ENDDO
      ENDDO
!
      RETURN
      END