MIXLEN.f

      SUBROUTINE mixlen(EL0,EL)
!
!     CALCULATES LAYER-AVERAGED BLACKADAR'S MIXING LENGTH, AND PBL TOP
!     AS CPBLT*(ASYMPTOTIC EL); AND THEN EL, ACCOUNT TAKEN OF STABILITY,
!     PBL TOP AND VERTICAL GRID DISTANCE RESTRICTIONS (SEE BELOW)
!
!     SET FROM EXISTING CODES BY L. LOBOCKI, JUNE 5, 1992
!       MODIFIED BY FEDOR MESINGER, OCTOBER 13, NOVEMBER 19
!       MODIFIED BY JIM TUCCILLO FOR MPI IMPLEMENTATION
!   01-10-25  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
!   02-06-19  MIKE BALDWIN - WRF VERSION
!   21-03-11  B Cui - change local arrays to dimension (im,jsta:jend)
!   21-09-30  J MENG - 2D DECOMPOSITION
!   
!
!     INPUT:
!     ------
!
!     ZINT (IM,jsta_2l:jend_2u,LP1) - ETA INTERFACES HEIGHT FIELD
!     T    (IM,jsta_2l:jend_2u,LM)  - TEMPERATURE
!     PMID (IM,jsta_2l:jend_2u,LM)  - PRESSURE IN LAYERS
!     Q2   (IM,jsta_2l:jend_2u,LM)  - TURBULENCE KINETIC ENERGY * 2
!     HGT  (IM,jsta_2l:jend_2u)     - SURFACE ELEVATION ARRAY
!     HTM  (IM,jsta_2l:jend_2u,LM)  - HEIGHT TOPOGRAPHY MASK ARRAY
!     EL0  (IM,JM)     - ARRAY OF ASYMPTOTIC VALUES FOR MIXING LENGTH
!
!     OUTPUT:
!     -------
!
!     EL   (IM,jsta_2l:jend_2u,LM) - FIELD OF RESULTING MASTER LENGTH SCALES
!
!
!     SCRATCH AREAS:
!     --------------
!
!     VKRMZ(IM,JM)
!
!     RELEVANT CONSTANTS:
!     -------------------
!
!     VON KARMAN CONSTANT:
      use vrbls3d, only: zint, pmid, t, q2
      use masks, only: lmh, htm
      use params_mod, only: epsq2, capa
      use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, im, jm, jsta_2l, jend_2u,&
              lm, lm1, spval,&
              ista, iend, ista_m, iend_m, ista_2l, iend_2u

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
!
      real,PARAMETER :: vkrm=0.4
!     CONSTANTS NEEDED FOR THE EL(BL,ST,ZI) SCHEME:
      real,PARAMETER :: frg=4.*9.8,drdrff=0.54,cpblt=10.,     &
        csh=0.23*0.5, epsn2=1.e-7
!
!     ------------------------------------------------------------------
!
      real,intent(in) :: el0(ista_2l:iend_2u,jsta_2l:jend_2u)
      real,intent(out) ::  el(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real hgt(ista:iend,jsta:jend),ape(ista_m:iend_m,jsta_m:jend_m,2)
!
      integer i,j,l
      real zl,vkrmz,ensq,q2kl,elst,ziag,elvgd

!***********************************************************************
!
!$omp  parallel do
      DO l=1,lm
        DO j=jsta,jend
        DO i=ista,iend
          el(i,j,l)=0.
        ENDDO
        ENDDO
      ENDDO
        DO j=jsta,jend
        DO i=ista,iend
          hgt(i,j)=zint(i,j,nint(lmh(i,j))+1)
        ENDDO
        ENDDO
!
!---THE AVERAGE EL SCHEME---------------------------(FM, AUGUST 19 MEMO)
!   FIRST GET EL IN THE LAYERS
!
!$omp  parallel do private(i,j,l,vkrmz,zl)
      DO l=1,lm
        DO j=jsta,jend
          DO i=ista,iend
            IF(hgt(i,j)<spval)THEN
            zl        = 0.5*(zint(i,j,l)+zint(i,j,l+1))
            vkrmz     = (zl-hgt(i,j))*vkrm
            el(i,j,l) = el0(i,j)*vkrmz/(el0(i,j)+vkrmz)
            ELSE
            el(i,j,l) = spval
            ENDIF
          ENDDO
        ENDDO
      ENDDO
!***
!***  GET NOW THE INTERFACE EL BY TWO-POINT AVERAGING OF LAYER VALUES
!***
      DO l=1,lm1
!$omp  parallel do private(i,j)
        DO j=jsta,jend
          DO i=ista,iend
            IF(hgt(i,j)<spval)THEN
            el(i,j,l) = 0.5*(el(i,j,l)+el(i,j,l+1))*htm(i,j,l+1)
            ELSE
            el(i,j,l) = spval
            ENDIF
          ENDDO
        ENDDO
      ENDDO
!
!$omp  parallel do private(i,j)
      DO j=jsta,jend
        DO i=ista,iend
          IF(hgt(i,j)<spval)THEN
          el(i,j,lm) = 0.0
          ELSE
          el(i,j,lm) = spval
          ENDIF
        ENDDO
      ENDDO
!---STABILITY, PBL TOP, AND VERTICAL GRID DISTANCE RESTRICTIONS:--------
!   COMPUTE EL STABLE AND
!   * USE THE SMALLER OF EL BLACKADAR, EL STABLE IF WITHIN PBL;
!   * USE THE SMALLEST OF EL STABLE, ELVGD, AND VKRMZ IF ABOVE PBL
!       (ASSUME PBL TOP IS AT CPBLT*EL0(K));
!$omp  parallel do private(i,j)
      DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
          ape(i,j,1) = (1.e5/pmid(i,j,1))**capa
        ENDDO
      ENDDO
!
      DO l=1,lm1
!$omp  parallel do private(i,j,elst,elvgd,ensq,q2kl,ziag)
        DO j=jsta_m,jend_m
          DO i=ista_m,iend_m
            IF(t(i,j,l)<spval)THEN
            ape(i,j,2) = (1.e5/pmid(i,j,l+1))**capa
            ensq = htm(i,j,l+1)*                                     &
                   frg*(t(i,j,l)*ape(i,j,1)-t(i,j,l+1)*ape(i,j,2))/  &
                  ((t(i,j,l)*ape(i,j,1)+t(i,j,l+1)*ape(i,j,2))*     &
                   (zint(i,j,l)-zint(i,j,l+2))+epsn2)
            ensq = amax1(ensq,epsn2)
            q2kl = amax1(epsq2,q2(i,j,l))
            elst = drdrff*sqrt(q2kl/ensq)
!WAS        ELST = DRDRFF*SQRT(Q2(I,J,L)/ENSQ)
            ziag = zint(i,j,l+1)-hgt(i,j)
!
            IF(ziag < cpblt*el0(i,j))THEN
              el(i,j,l) = amin1(el(i,j,l),elst)
            ELSE
              elvgd     = csh*(zint(i,j,l)-zint(i,j,l+2))
              el(i,j,l) = amin1(elst,elvgd,vkrm*ziag)
            ENDIF
            ape(i,j,1) = ape(i,j,2)
            ELSE
            el(i,j,l) = spval
            ENDIF
          ENDDO
        ENDDO
      ENDDO
!
      RETURN
      END