SUBROUTINE mixlen(lmh, lpbl, hpbl, u, v, t, q, q2, ape, z, gm, gh,el)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !
  !   level 2.5 mixing length 
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: lp1=lm+1, lm1=lm-1
  INTEGER, DIMENSION(im,jm    ), INTENT(IN ) :: lmh 
  INTEGER, DIMENSION(im,jm    ), INTENT(OUT) :: lpbl
  REAL,    DIMENSION(im,jm    ), INTENT(OUT) :: hpbl
  REAL,    DIMENSION(0:im+1,0:jm+1,lm ), INTENT(IN ) :: t
  REAL,    DIMENSION(0:im+1,0:jm+1,lm ), INTENT(IN ) :: q
  REAL,    DIMENSION(im,jm,lp1), INTENT(IN ) :: u
  REAL,    DIMENSION(im,jm,lp1), INTENT(IN ) :: v
  REAL,    DIMENSION(im,jm,lm ), INTENT(IN ) :: q2
  REAL,    DIMENSION(im,jm,lm ), INTENT(IN ) :: ape 
  REAL,    DIMENSION(0:im+1,0:jm+1,lp1), INTENT(IN ) :: z
  REAL,    DIMENSION(im,jm,lm1), INTENT(OUT) :: gm, gh, el 

  REAL, PARAMETER :: epsq2=0.2
  REAL, PARAMETER :: epsl=0.32
  REAL, PARAMETER :: epsru=1.e-12
  REAL, PARAMETER :: epsrs=1.e-12
  REAL, PARAMETER :: fh=1.01
  REAL, PARAMETER :: alph=.20
  REAL, PARAMETER :: vkrm=.40
  REAL, PARAMETER :: elfc=0.23*0.5
  REAL, PARAMETER :: el0max=1000.
  REAL, PARAMETER :: g=9.8
  REAL, PARAMETER :: beta=1./270.
  REAL, PARAMETER :: btg=beta*g
  REAL, PARAMETER :: a1=0.659888514560862645
  REAL, PARAMETER :: a2=0.6574209922667784586
  REAL, PARAMETER :: b1=11.87799326209552761
  REAL, PARAMETER :: b2=7.226971804046074028
  REAL, PARAMETER :: c1=0.000830955950095854396

  !*** coefficients of the terms in the denominator

  REAL, PARAMETER :: adnm=18.*a1*a1*a2*(b2-3.*a2)*btg
  REAL, PARAMETER :: adnh= 9.*a1*a2*a2*(12.*a1+3.*b2)*btg*btg
  REAL, PARAMETER :: bdnm= 6.*a1*a1
  REAL, PARAMETER :: bdnh= 3.*a2*(7.*a1+b2)*btg

  !*** free term in the equilibrium equation for (l/q)**2

  REAL, PARAMETER :: aeqm= 3.*a1*a2*b1*(3.*a2+3.*b2*c1+18.*a1*c1-b2)*btg &
       +18.*a1*a1*a2*(b2-3.*a2)*btg
  REAL, PARAMETER :: aeqh= 9.*a1*a2*a2*b1*btg*btg &
       +9.*a1*a2*a2*(12.*a1+3.*b2)*btg*btg

  !*** forbidden turbulence area

  REAL, PARAMETER :: requ=-aeqh/aeqm*1.02
  REAL, PARAMETER :: epsgh=1.e-9
  REAL, PARAMETER :: epsgm=requ*epsgh

  !*** near isotropy for shear turbulence, ww/q2 lower limit

  REAL, PARAMETER :: ubryl=(18.*requ*a1*a1*a2*b2*c1*btg&
       +9.*a1*a2*a2*b2*btg*btg)/(requ*adnm+adnh)
  REAL, PARAMETER :: ubry=(1.+epsrs)*ubryl
  REAL, PARAMETER :: ubry3=3.*ubry  
  REAL, PARAMETER :: aubm=54.*a1*a1*a2*b2*c1*btg -adnm*ubry3 
  REAL, PARAMETER :: aubh=27.*a1*a2*a2*b2*btg*btg-adnh*ubry3 
  REAL, PARAMETER :: bubm=18.*a1*a1*c1           -bdnm*ubry3 
  REAL, PARAMETER :: bubh=(9.*a1*a2+3.*a2*b2)*btg-bdnh*ubry3 
  REAL, PARAMETER :: cubr=1.- ubry3
  REAL, PARAMETER :: rcubr=1./cubr

  INTEGER i, j, l
  REAL    qdzl, rdz, gml, ghl, aubr, bubr
  REAL    qol2st, qol2un, eloq2x, aden, bden, vkrmz
  REAL,    DIMENSION(im,jm,lm1) :: elm
  REAL,    DIMENSION(im,jm,lm ) :: q1, thv
  REAL,    DIMENSION(im,jm)     :: szq, sq, el0



  !*** locate the height of the pbl

  DO j = 3, jm-2
     DO i = 1, im
        lpbl(i,j)=1
     END DO
  END DO
  DO l=2,lm-1
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lmh(i,j)-1) THEN
              IF(q2(i,j,l).LE.epsq2*fh)THEN
                 lpbl(i,j)=l
              END IF
           ENDIF
        END DO
     END DO
  END DO

  !*** find the height of the pbl

  DO j = 3, jm-2
     DO i = 1, im
        hpbl(i,j)=z(i,j,lpbl(i,j))-z(i,j,lmh(i,j)+1)
     END DO
  END DO

  DO l=1,lm
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lmh(i,j)) THEN
              thv(i,j,l)=(0.608*q(i,j,l)+1.)*t(i,j,l)*ape(i,j,l)
           END IF
        END DO
     END DO
  END DO

  q1 = 0.0
  DO j = 3, jm-2
     DO i = 1, im
        DO l=lpbl(i,j),lmh(i,j)
           q1(i,j,l)=SQRT(q2(i,j,l))
        END DO
     END DO
  END DO

  DO j = 3, jm-2
     DO i = 1, im
        szq(i,j)=0.
        sq(i,j) =0.
     END DO
  END DO
  DO l=1,lm-1
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lmh(i,j)-1) THEN
              qdzl=(q1(i,j,l)+q1(i,j,l+1))*(z(i,j,l+1)-z(i,j,l+2))
              szq(i,j)=(z(i,j,l+1)+z(i,j,l+2)-z(i,j,lmh(i,j)+1)-z(i,j,lmh(i,j)+1))*&
                   qdzl+szq(i,j)
              sq(i,j)=qdzl+sq(i,j)
              rdz=2./(z(i,j,l)-z(i,j,l+2))
              gml=((u(i,j,l)-u(i,j,l+1))**2+(v(i,j,l)-v(i,j,l+1))**2)*rdz*rdz
              gm(i,j,l)=amax1(gml,epsgm)
              ghl=(thv(i,j,l)-thv(i,j,l+1))*rdz
              IF (ABS(ghl).LE.epsgh) ghl=epsgh
              gh(i,j,l)=ghl
           END IF
        END DO
     END DO
  END DO

  !*** computation of asymptotic l in blackadar formula

  DO j = 3, jm-2
     DO i = 1, im
        el0(i,j)=amin1(alph*szq(i,j)*0.5/sq(i,j),el0max)
     END DO
  END DO

  DO l=1,lm-1
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lmh(i,j)-1) THEN
              IF(gh(i,j,l).GE.epsgh)THEN
                 IF(gm(i,j,l)/gh(i,j,l).LE.requ)THEN
                    elm(i,j,l)=epsl
                 ELSE
                    aubr=(aubm*gm(i,j,l)+aubh*gh(i,j,l))*gh(i,j,l)
                    bubr= bubm*gm(i,j,l)+bubh*gh(i,j,l)
                    qol2st=(-0.5*bubr+SQRT(bubr*bubr*0.25-aubr*cubr))*rcubr
                    eloq2x=1./qol2st
                    elm(i,j,l)=amax1(SQRT(eloq2x*q2(i,j,l)),epsl)
                 ENDIF
              ELSE
                 aden=(adnm*gm(i,j,l)+adnh*gh(i,j,l))*gh(i,j,l)
                 bden= bdnm*gm(i,j,l)+bdnh*gh(i,j,l)
                 qol2un= -0.5*bden+SQRT(bden*bden*0.25-aden)
                 eloq2x=1./((1.+epsru)*qol2un)
                 elm(i,j,l)=amax1(SQRT(eloq2x*q2(i,j,l)),epsl)
              ENDIF
           END IF
        END DO
     END DO
  END DO

  DO l=1,lm-1
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lpbl(i,j)-1) THEN
              el(i,j,l)=amin1((z(i,j,l)-z(i,j,l+2))*elfc,elm(i,j,l))
           ELSE IF (l <= lmh(i,j)-1) THEN
              vkrmz=(z(i,j,l+1)-z(i,j,lmh(i,j)+1))*vkrm
              el(i,j,l)=amin1(vkrmz/(vkrmz/el0(i,j)+1.),elm(i,j,l))
           END IF
        END DO
     END DO
  END DO

END SUBROUTINE mixlen
