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
!
!     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, jm
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      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,JM),intent(inout) ::  SQZ,SQ,RQ2L,RQ2H,EL0
      real,dimension(IM,JM) ::   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 220 L=1,LM
!$omp  parallel do
!$omp& private(dp,rq2m)
      DO J=JSTA,JEND
      DO I=1,IM
        IF(Q2(I,J,L).LE.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.GE.LMH(I,J)) GOTO 215
!         RQ2L(I,J)=SQRT(Q2(I,J,L))
!         IF(Q2(I,J,L).LT.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
 220  CONTINUE
!***
!***    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
        EL0(I,J)= AMAX1(AMIN1(                                           &
     &     ((SM(I,J)*ALPHAS+(1.0-SM(I,J))*ALPHAL)*SQZ(I,J)               &
     &     /(SQ(I,J)+EPSQ2)),EL0M),ELMIN)
      ENDDO
      ENDDO
!
      RETURN
      END