MODULE CCM3PBL
  USE Constants, ONLY :     &
       cp,            &
       grav,          &
       gasr,          &
       r8,i4,i8

 USE PBL_Entrain, ONLY :  &
       PBLNASA

  USE Options, ONLY : &
      jdt        ,&
      nferr      , &!
      nfprt      , &!
      PBLEntrain , &!
      Wgh1       , &!1.0_r8/3.0_r8 ! pbl Hostlag Boville
      Wgh2       , &!1.0_r8/3.0_r8 ! pbl Mellor Yamada 2.0
      Wgh3         !1.0_r8/3.0_r8 ! pbl Mellor Yamada 2.5 

  IMPLICIT NONE
  PRIVATE
  REAL(KIND=r8), PARAMETER :: onet  = 1.0_r8/3.0_r8   ! 1/3 power in wind gradient expression  common /compbl/
  REAL(KIND=r8), PARAMETER :: fak    = 8.5_r8     ! Constant in surface temperature excess      common /compbl/    
  REAL(KIND=r8), PARAMETER :: sffrac = 0.01_r8     ! Surface layer fraction of boundary layer common /compbl/
  REAL(KIND=r8), PARAMETER :: vk     = 0.40_r8    ! Von Karman's constant common /compbl/
  REAL(KIND=r8), PARAMETER :: fakn   = 7.2_r8     ! Constant in turbulent prandtl number
  REAL(KIND=r8), PARAMETER :: ccon  = fak*sffrac*vk  ! fak * sffrac * vk  common /compbl/

  REAL(KIND=r8), PARAMETER :: betam  = 15.0_r8 ! Constant in wind gradient expression common /compbl/
  REAL(KIND=r8), PARAMETER :: betas  =  5.0_r8 ! Constant in surface layer gradient expression common /compbl/
  REAL(KIND=r8), PARAMETER :: betah  = 15.0_r8 ! Constant in temperature gradient expression  common /compbl/
  REAL(KIND=r8), PARAMETER :: ricr   =  0.1_r8  ! Critical richardson number
    !
    ! sffrac = 0.1_r8     ! Surface layer fraction of boundary layer common /compbl/
    ! 
    ! betam  = 15.0_r8 ! Constant in wind gradient expression common /compbl/
    !
  REAL(KIND=r8), PARAMETER :: binm  = betam*sffrac  ! betam * sffrac  common /compbl/	binm = betam*sffrac
  REAL(KIND=r8), PARAMETER :: binh  = betah*sffrac  ! betah * sffrac common /compbl/   binh = betah*sffrac

  REAL(KIND=r8), PARAMETER :: cpair  = 1004.64_r8 !heat capacity dry air at const pres (j/kg/kelvin)
  REAL(KIND=r8), PARAMETER :: rair   = 287.04_r8  !gas constant for dry air (j/kg/kelvin)
  REAL(KIND=r8), PARAMETER :: gravit = 9.80616_r8 !Acceleration due to gravity common/comvd/
  REAL(KIND=r8), PARAMETER :: g= gravit      ! Gravitational acceleration  common /compbl/
    !
    !-----------------------------------------------------------------------
    !
    ! Hard-wired numbers.
    ! zkmin = minimum k = kneutral*f(ri)
    !
  REAL(KIND=r8), PARAMETER    :: zkmin   = 0.01_r8 ! Minimum kneutral*f(ri) common/comvd/
  INTEGER      , PARAMETER    :: ntopfl=1  ! Top level to which vertical diffusion is applied. common/comvd/
  INTEGER                     :: npbl      ! Maximum number of levels in pbl from surface common/comvd/

  REAL(KIND=r8), ALLOCATABLE  :: ml2    (:)    ! Mixing lengths squared common/comvd/
  REAL(KIND=r8), ALLOCATABLE  :: hypm   (:)    ! reference pressures at midpoints
  REAL(KIND=r8), ALLOCATABLE  :: qmin   (:)    ! Global minimum constituent concentration
  REAL(KIND=r8), ALLOCATABLE  :: qmincg (:)    ! Min. constituent concentration counter-gradient term
  REAL(KIND=r8), ALLOCATABLE  :: con0   (:)
  REAL(KIND=r8), ALLOCATABLE  :: sigkiv (:)
  REAL(KIND=r8), ALLOCATABLE  :: a0     (:)
  REAL(KIND=r8), ALLOCATABLE  :: b0     (:)
  REAL(KIND=r8), ALLOCATABLE  :: sigr   (:)
  REAL(KIND=r8), ALLOCATABLE  :: sigriv (:)
  REAL(KIND=r8), ALLOCATABLE  :: con1   (:)
  REAL(KIND=r8), ALLOCATABLE  :: con2   (:)
  REAL(KIND=r8), ALLOCATABLE  :: t0     (:)
  REAL(KIND=r8), ALLOCATABLE  :: t1     (:)
  REAL(KIND=r8), ALLOCATABLE  :: delsig2(:)
  REAL(KIND=r8)               :: c0

  REAL(KIND=r8),    PARAMETER :: CPWV   =  1.810e3_r8 ! specific heat of water vap
  REAL(KIND=r8),    PARAMETER :: cpwvx  = CPWV      ! spec. heat of water vapor at const. pressure
  REAL(KIND=r8),    PARAMETER :: cpairx = cpair      ! specific heat of dry air
  REAL(KIND=r8),    PARAMETER :: cpvir  = cpwvx/cpairx - 1.0_r8 ! Derived constant for cp moist air common/comvd/
  REAL(KIND=r8),    PARAMETER :: gkm0   = 1.00_r8
  REAL(KIND=r8),    PARAMETER :: gkh0   = 0.10_r8
  REAL(KIND=r8),    PARAMETER :: gkm1   = 300.0_r8
  REAL(KIND=r8),    PARAMETER :: gkh1   = 300.0_r8
  INTEGER      ,    PARAMETER :: kmean  =   1
  REAL(KIND=r8),    PARAMETER :: facl   =   0.05_r8
  REAL(KIND=r8),    PARAMETER :: a1     =   0.92_r8
  REAL(KIND=r8),    PARAMETER :: a2     =   0.74_r8
  REAL(KIND=r8),    PARAMETER :: b1     =  16.6_r8
  REAL(KIND=r8),    PARAMETER :: b2     =  10.1_r8
  REAL(KIND=r8),    PARAMETER :: c1     =   0.08_r8
  REAL(KIND=r8),    PARAMETER :: deltx  =   0.0_r8
  REAL(KIND=r8),    PARAMETER :: eps    =   0.608_r8
  !  gbyr        =(grav/gasr)**2!(m/sec**2)/(J/(Kg*K))=(m/sec**2)/((Kg*(m/sec**2)*m)/(Kg*K))
    !(m/sec**2)/((Kg*(m**2/sec**2))/(Kg*K))
    !(m/sec**2)/(m**2/sec**2*K)=K**2/m**2

  REAL(KIND=r8),    PARAMETER :: gbyr= (gravit/gasr)**2!(m/sec**2)/(J/(Kg*K))=(m/sec**2)/((Kg*(m/sec**2)*m)/(Kg*K))
    !(m/sec**2)/((Kg*(m**2/sec**2))/(Kg*K))
    !(m/sec**2)/(m**2/sec**2*K)=K**2/m**2
  INTEGER     , PARAMETER     :: nitr  =   2
  REAL(KIND=r8)               :: alfa
  REAL(KIND=r8)               :: beta
  REAL(KIND=r8)               :: gama
  REAL(KIND=r8)               :: dela
  REAL(KIND=r8)               :: r1
  REAL(KIND=r8)               :: r2
  REAL(KIND=r8)               :: r3
  REAL(KIND=r8)               :: r4
  REAL(KIND=r8)               :: s1
  REAL(KIND=r8)               :: s2
  REAL(KIND=r8)               :: rfc

  REAL(KIND=r8),PARAMETER     :: VKARMAN=0.4_r8
  !---------------------------------------------------------------------
  ! --- CONSTANTS
  !---------------------------------------------------------------------

  REAL(KIND=r8) :: aa1,   aa2,   bb1,  bb2,  ccc
  REAL(KIND=r8) :: ckm1,  ckm2,  ckm3, ckm4, ckm5, ckm6, ckm7, ckm8
  REAL(KIND=r8) :: ckh1,  ckh2,  ckh3, ckh4
  REAL(KIND=r8) :: cvfq1, cvfq2, bcq
    
  REAL(KIND=r8), PARAMETER :: aa1_old =  0.78_r8
  REAL(KIND=r8), PARAMETER :: aa2_old =  0.79_r8
  REAL(KIND=r8), PARAMETER :: bb1_old = 15.0_r8
  REAL(KIND=r8), PARAMETER :: bb2_old =  8.0_r8
  REAL(KIND=r8), PARAMETER :: ccc_old =  0.056_r8
  
  REAL(KIND=r8), PARAMETER :: aa1_new =  0.92_r8
  REAL(KIND=r8), PARAMETER :: aa2_new =  0.74_r8
  REAL(KIND=r8), PARAMETER :: bb1_new = 16.0_r8
  REAL(KIND=r8), PARAMETER :: bb2_new = 10.0_r8
  REAL(KIND=r8), PARAMETER :: ccc_new =  0.08_r8
  
  REAL(KIND=r8), PARAMETER :: cc1     =  0.27_r8
  REAL(KIND=r8), PARAMETER :: t00     =  2.7248e2_r8 
  REAL(KIND=r8), PARAMETER :: small2   =  1.0e-10_r8
  REAL(KIND=r8), PARAMETER :: VONKARM = 0.40_r8     

  !---------------------------------------------------------------------
  ! --- NAMELIST
  !---------------------------------------------------------------------

  REAL(KIND=r8), PARAMETER    :: TKEmax       =   5.0_r8
  REAL(KIND=r8), PARAMETER    :: TKEmin       =   1.0e-12_r8
  REAL(KIND=r8), PARAMETER    :: el0max1      =   1.0e6_r8
  REAL(KIND=r8), PARAMETER    :: el0min1      =   0.0_r8
  REAL(KIND=r8), PARAMETER    :: alpha_land   =   0.10_r8
  REAL(KIND=r8), PARAMETER    :: alpha_sea    =   0.10_r8
  REAL(KIND=r8), PARAMETER    :: akmax        = 300.0_r8
  REAL(KIND=r8), PARAMETER    :: akmin_land   =   5.0_r8
  REAL(KIND=r8), PARAMETER    :: akmin_sea    =   0.5_r8
  INTEGER, PARAMETER :: nk_lim       =  2
  INTEGER, PARAMETER :: init_iters   = 1! 20
  LOGICAL, PARAMETER :: do_thv_stab  = .TRUE.
  LOGICAL, PARAMETER :: use_old_cons = .FALSE.
  LOGICAL, PARAMETER :: ensamble     = .TRUE.

  REAL(KIND=r8), PARAMETER    :: kcrit        =  0.01_r8



  PUBLIC :: vdinti
  PUBLIC :: vdintr
CONTAINS

  SUBROUTINE vdinti(plev, sig, delsig, sigml,pnats,pcnst)
    !----------------------------------------------------------------------- 
    ! 
    ! Purpose: 
    ! Initialization of time independent fields for vertical diffusion.
    ! 
    ! Method: 
    ! Call initialization routine for boundary layer scheme.
    ! 
    ! Author: J. Rosinski
    !
    !-----------------------------------------------------------------------
    !   use precision
    !   use pmgrid

    IMPLICIT NONE

    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    INTEGER,  INTENT(in) :: plev         ! number of vertical levels
    REAL(KIND=r8),    INTENT(IN) :: sig   (plev)
    REAL(KIND=r8),    INTENT(IN) :: delsig(plev)
    REAL(KIND=r8),    INTENT(IN) :: sigml (plev+1)
    INTEGER,OPTIONAL, INTENT(IN) :: pnats      ! number of non-advected trace species
    INTEGER,OPTIONAL, INTENT(IN) :: pcnst      ! number of constituents (including water vapor)

    REAL(KIND=r8) :: rairx        ! gas constant for dry air
    REAL(KIND=r8) :: akappa
    REAL(KIND=r8) :: gam1
    REAL(KIND=r8) :: gam2

    REAL(KIND=r8) :: ps0    
    REAL(KIND=r8) :: sigk(plev)

    !
    !---------------------------Local workspace-----------------------------
    !
    INTEGER :: k         ! vertical loop index 
    INTEGER :: m
    !
    ALLOCATE(hypm   (plev))               ;hypm  =0.0_r8
    ALLOCATE(ml2    (plev + 1))           ;ml2   =0.0_r8
    ALLOCATE(qmin   (pcnst+pnats))        ;qmin  =0.0_r8
    ALLOCATE(qmincg (pcnst+pnats))        ;qmincg=0.0_r8 
    ALLOCATE(con0   (plev))               ;con0  =0.0_r8
    ALLOCATE(sigkiv (plev))               ;sigkiv=0.0_r8
    ALLOCATE(a0     (plev))               ;a0    =0.0_r8
    ALLOCATE(b0     (plev))               ;b0    =0.0_r8
    ALLOCATE(sigr   (plev))               ;sigr  =0.0_r8
    ALLOCATE(sigriv (plev))               ;sigriv=0.0_r8
    ALLOCATE(con1   (plev))               ;con1  =0.0_r8
    ALLOCATE(con2   (plev))               ;con2  =0.0_r8
    ALLOCATE(t0     (plev))               ;t0    =0.0_r8
    ALLOCATE(t1     (plev))               ;t1    =0.0_r8
    ALLOCATE(delsig2(plev))               ;delsig2=0.0_r8
    ps0    = 1.0e5_r8
    delsig2=delsig
    gam1=1.0_r8/3.0_r8-2.0_r8*a1/b1
    gam2=(b2+6.0_r8*a1)/b1
    alfa=b1*(gam1-c1)+3.0_r8*(a2+2.0_r8*a1) 
    beta=b1*(gam1-c1)
    gama=a2/a1*(b1*(gam1+gam2)-3.0_r8*a1)
    dela=a2/a1* b1* gam1
    r1  =0.5_r8*gama/alfa
    r2  =    beta/gama
    r3  =2.0_r8*(2.0_r8*alfa*dela-gama*beta)/(gama*gama)
    r4  =r2*r2
    s1  =3.0_r8*a2* gam1
    s2  =3.0_r8*a2*(gam1+gam2)
    !     
    !     critical flux richardson number
    !     
    rfc =s1/s2
    akappa=gasr/cp

    rairx  =rair
      
    DO k = 1, plev
       sigk  (k)=sig(k)**akappa
       sigkiv(k)=1.0_r8/sigk(k)
       con0  (k)=gasr*delsig(k)/(grav*sig(k))
    END DO
    a0    (plev)=0.0_r8
    b0    (   1)=0.0_r8
    sigr  (plev)=0.0_r8
    sigriv(   1)=0.0_r8
    !gbyr        =(grav/gasr)**2!(m/sec**2)/(J/(Kg*K))=(m/sec**2)/((Kg*(m/sec**2)*m)/(Kg*K))
    !(m/sec**2)/((Kg*(m**2/sec**2))/(Kg*K))
    !(m/sec**2)/(m**2/sec**2*K)=K**2/m**2
    DO k = 1, plev-1
       con1  (   k)=grav*sigml(k+1)/(gasr*(sig(k)-sig(k+1)))
       con2  (   k)=grav*con1(k)
       con1  (   k)=con1(k)*con1(k)
       t0    (   k)=(sig(k+1)-sigml(k+1))/(sig(k+1)-sig(k))
       t1    (   k)=(sigml(k+1)-sig(k  ))/(sig(k+1)-sig(k))
       sigr  (   k)=sigk(k)*sigkiv(k+1)
       sigriv( k+1)=sigk(k+1)*sigkiv(k)
       a0(k)=gbyr*sigml(k+1)**2/(delsig(k  )*(sig(k)-sig(k+1)))
       b0(k+1)=gbyr*sigml(k+1)**2/(delsig(k+1)*(sig(k)-sig(k+1)))
    END DO

    ! hypm     reference state midpoint pressures
    ps0    = 1.0e5_r8            ! Base state surface pressure (pascals)
    DO k=plev,1,-1
       hypm(k) =  ps0*sig(k)
    END DO
    qmin(1) = 1.0e-12_r8        ! Minimum mixing ratio for moisture
    DO m=2,pcnst
       qmin(m) = 0.0_r8
    END DO
    qmincg(1) = 0.0_r8
    DO m=2,pcnst
       qmincg(m) = qmin(m)
    END DO
    c0=grav/(gasr*delsig(1))

    !
    !-----------------------------------------------------------------------
    !
    ! Hard-wired numbers.
    ! zkmin = minimum k = kneutral*f(ri)
    !
    ! zkmin = 0.01_r8
    !
    ! Set physical constants for vertical diffusion and pbl:
    !
    !
    ! Derived constants
    ! ntopfl = top level to which v-diff is applied
    ! npbl = max number of levels (from bottom) in pbl
    !
    !
    ! Limit pbl height to regions below 400 mb
    !
    npbl = 0
    DO k=plev,1,-1
       IF (hypm(k) >= 4.e4_r8) THEN!40000
          npbl = npbl + 1
       END IF
    END DO
    npbl = MAX(npbl,1)
    WRITE(nfprt,*)'VDINTI: PBL height will be limited to bottom ',npbl, &
         ' model levels. Top is ',hypm(plev + 1-npbl),' pascals'

    !ntopfl = 1
    !IF (plev.EQ.1) ntopfl = 0
    !
    ! Set the square of the mixing lengths.
    !
    ml2(1) = 0.0_r8
    DO k=2,plev
       ml2(k) = 10.0_r8**2
    END DO
    ml2(plev + 1) = 0.0_r8
    !
    ! Initialize pbl variables
    !
    !CALL pbinti()
    !
    !---------------------------------------------------------------------
    ! --- Initialize constants
    !---------------------------------------------------------------------

    IF( use_old_cons ) THEN
       aa1 = aa1_old 
       aa2 = aa2_old  
       bb1 = bb1_old 
       bb2 = bb2_old  
       ccc = ccc_old  
    ELSE
       aa1 = aa1_new 
       aa2 = aa2_new  
       bb1 = bb1_new 
       bb2 = bb2_new  
       ccc = ccc_new  
    END IF

    ckm1 = ( 1.0_r8 - 3.0_r8*ccc )*aa1
    ckm3 =  3.0_r8 * aa1*aa2*    ( bb2 - 3.0_r8*aa2 )
    ckm4 =  9.0_r8 * aa1*aa2*ccc*( bb2 + 4.0_r8*aa1 )
    ckm5 =  6.0_r8 * aa1*aa1
    ckm6 = 18.0_r8 * aa1*aa1*aa2*( bb2 - 3.0_r8*aa2 )
    ckm7 =  3.0_r8 * aa2*        ( bb2 + 7.0_r8*aa1 )
    ckm8 = 27.0_r8 * aa1*aa2*aa2*( bb2 + 4.0_r8*aa1 )
    ckm2 =  ckm3 - ckm4
    ckh1 =  aa2
    ckh2 =  6.0_r8 * aa1*aa2
    ckh3 =  3.0_r8 * aa2*( bb2 + 4.0_r8*aa1 )
    ckh4 =  2.0e-6_r8 * aa2
    cvfq1 = 5.0_r8 * cc1 / 3.0_r8
    cvfq2 = 1.0_r8 / bb1
    bcq   = 0.5_r8 * ( bb1**(2.0_r8/3.0_r8) )

    !---------------------------------------------------------------------


    RETURN
  END SUBROUTINE vdinti



  SUBROUTINE pbinti()
    !----------------------------------------------------------------------- 
    ! 
    ! Purpose: 
    ! Initialize time independent variables of pbl package.
    ! 
    ! Method: 
    ! <Describe the algorithm(s) used in the routine.> 
    ! <Also include any applicable external references.> 
    ! 
    ! Author: B. Boville
    ! 
    !-----------------------------------------------------------------------
    IMPLICIT NONE

    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    !
    ! Basic constants
    !
    !
    ! Derived constants
    !
    !ccon = fak*sffrac*vk
    !
    ! sffrac = 0.1_r8     ! Surface layer fraction of boundary layer common /compbl/
    ! 
    ! betam  = 15.0_r8 ! Constant in wind gradient expression common /compbl/
    !
    !binm = betam*sffrac
    !binh = betah*sffrac

    !
    RETURN
  END SUBROUTINE pbinti

  SUBROUTINE vdintr(&
                    plon       , &! INTENT(IN        ) plon ! number of longitudes
                    plond      , &! INTENT(IN        ) plond ! slt extended domain longitude
                    plev       , &! INTENT(IN        ) plev ! number of vertical levels
                    pcnst      , &! INTENT(IN        ) pcnst ! number of constituents (including water vapor)
                    ztodt      , &! INTENT(IN        ) ztodt ! 2 delta-t
                    colrad     , &! INTENT(IN        ) Cosino the colatitude [radian]
                    gl0        , &! INTENT(IN        ) Maximum mixing length l0 in blackerdar's formula [m]
                    bstar      , &! INTENT(IN        ) surface_bouyancy_scale m s-2
                    FRLAND     , &! INTENT(IN        ) Fraction Land [%]
                    z0         , &! INTENT(IN        ) Rougosiness [m]
                    LwCoolRate , &! INTENT(IN        ) air_temperature_tendency_due_to_longwave [K s-1]
                    LwCoolRateC, &! INTENT(IN        ) clear sky_air_temperature_tendency_lw [K s-1]
                    cldtot     , &! INTENT(IN        ) cloud fraction [%]
                    qliq       , &! INTENT(IN        ) cloud fraction [kg/kg]
                    pmidm1     , &! INTENT(IN        ) pmidm1(plond,plev) ! midpoint pressures
                    pintm1     , &! INTENT(IN        ) pintm1(plond,plev + 1)    ! interface pressures
                    psomc      , &! INTENT(IN        ) psomc(plond,plev)         ! (psm1/pmidm1)**cappa
                    thm        , &! INTENT(IN        ) thm(plond,plev) ! potential temperature midpoints
                    zm         , &! INTENT(IN        ) zm(plond,plev) ! midpoint geopotential height above sfc
                    zhalf      , &! INTENT(IN        ) zhalf(plond,plev) ! interface pressures geopotential height
                    psfcpa     , &! INTENT(IN        ) surface pressure [mb]
                    USTAR      , &! INTENT(IN        ) scale velocity turbulent [m/s]
                    TSK        , &! INTENT(IN        ) surface temperature
                    QSFC       , &! INTENT(IN        ) surface specific temperature
                    rpdel      , &! INTENT(IN        ) rpdel(plond,plev)! 1./pdel (thickness between interfaces)
                    rpdeli     , &! INTENT(IN        ) rpdeli(plond,plev)! 1./pdeli (thickness between midpoints)
                    um1        , &! INTENT(IN        ) um1(plond,plev)         ! u-wind input
                    vm1        , &! INTENT(IN        ) vm1(plond,plev)         ! v-wind input
                    tm1        , &! INTENT(IN        ) tm1(plond,plev)         ! temperature input
                    taux       , &! INTENT(IN        ) taux(plond)                   ! x surface stress (n)
                    tauy       , &! INTENT(IN        ) tauy(plond)                   ! y surface stress (n)
                    shflx      , &! INTENT(IN        ) shflx(plond)! surface sensible heat flux (w/m2)
                    cflx       , &! INTENT(IN        ) cflx(plond,pcnst)! surface constituent flux (kg/m2/s)
                    qm1        , &! INTENT(INOUT) qm1(plond,plev,pcnst)  ! initial/final constituent field
                    dtv        , &! INTENT(OUT  ) dtv(plond,plev)     ! temperature tendency (heating)
                    dqv        , &! INTENT(OUT  ) dqv(plond,plev,pcnst)
                    duv        , &! INTENT(OUT  ) duv(plond,plev)     ! u-wind tendency
                    dvv        , &! INTENT(OUT  ) dvv(plond,plev)     ! v-wind tendency
                    up1        , &! INTENT(OUT  ) up1(plond,plev)     ! u-wind after vertical diffusion
                    vp1        , &! INTENT(OUT  ) vp1(plond,plev)     ! v-wind after vertical diffusion
                    pblh       , &! INTENT(OUT  ) pblh(plond)! planetary boundary layer height
                    rino       , &! INTENT(INOUT) bulk Richardson no. from level to ref lev
                    tpert      , &! INTENT(OUT  ) tpert(plond)! convective temperature excess
                    qpert      , &! INTENT(OUT  ) qpert(plond,pcnst)! convective humidity and constituent excess
                    TKE        , &! INTENT(INOUT) Turbulent kinetic Energy [m/s]^2
                    kvh        , &! INTENT(OUT  ) Heat Coeficient Difusivity 
                    kvm        , &! INTENT(OUT  ) Momentun Coeficient Difusivity 
                    obklen     , &! INTENT(OUT  ) Heat Coeficient Difusivity 
                    phiminv    , &! INTENT(OUT  ) Momentum Stability Function 
                    phihinv    , &! INTENT(OUT  ) Heat Stability Function 
                    tstar      , &
                    wstar        )
    !-----------------------------------------------------------------------
    !
    ! interface routine for vertical diffusion and pbl scheme
    !
    ! calling sequence:
    !
    !    vdinti        initializes vertical diffustion constants
    !    pbinti        initializes pbl constants
    !     .
    !     .
    !    vdintr        interface for vertical diffusion and pbl scheme
    !      vdiff       performs vert diff and pbl
    !        pbldif    boundary layer scheme
    !        mvdiff    diffuse momentum
    !        qvdiff    diffuse constituents
    !
    !---------------------------Code history--------------------------------
    !
    ! Original version:  J. Rosinski
    ! Standardized:      J. Rosinski, June 1992
    ! Reviewed:          P. Rasch, B. Boville, August 1992
    ! Reviewed:          P. Rasch, April 1996
    ! Reviewed:          B. Boville, April 1996
    !
    !-----------------------------------------------------------------------
    !
    ! $Id: vdintr.F,v 1.1.1.1 2001/03/09 00:29:29 mirin Exp $
    !
    !-----------------------------------------------------------------------
    !!      use precision
    !!      use pmgrid
    !-----------------------------------------------------------------------
    !!#include <implicit.h>
    !------------------------------Commons----------------------------------
    !!#include <comtrcnm.h>
    !-----------------------------------------------------------------------
    !!#include <comvd.h>
    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) ::  plon       ! number of longitudes
    INTEGER, INTENT(IN   ) ::  plev       ! number of vertical levels
    INTEGER, INTENT(IN   ) ::  plond      ! slt extended domain longitude
    INTEGER, INTENT(IN   ) ::  pcnst      ! number of constituents (including water vapor)
    REAL(kind=r8), INTENT(INOUT) ::  gl0        (plond)     !
    REAL(kind=r8), INTENT(IN   ) ::  bstar      (plond)     !
    REAL(kind=r8), INTENT(IN   ) ::  FRLAND     (plond)     !
    REAL(kind=r8), INTENT(IN   ) ::  LwCoolRate (plond,plev)     !
    REAL(kind=r8), INTENT(IN   ) ::  LwCoolRateC(plond,plev)     !
    REAL(kind=r8), INTENT(IN   ) ::  cldtot     (plond,plev)     !
    REAL(kind=r8), INTENT(IN   ) ::  qliq       (plond,plev)     !
    REAL(kind=r8), INTENT(IN   ) ::  pmidm1     (plond,plev)     ! midpoint pressures
    REAL(kind=r8), INTENT(IN   ) ::  pintm1     (plond,plev + 1) ! interface pressures
    REAL(kind=r8), INTENT(IN   ) ::  psomc      (plond,plev)     ! (psm1/pmidm1)**cappa
    REAL(kind=r8), INTENT(IN   ) ::  thm        (plond,plev)     ! potential temperature midpoints
    REAL(kind=r8), INTENT(IN   ) ::  zm         (plond,plev)     ! midpoint geopotential height above sfc
    REAL(kind=r8), INTENT(IN   ) ::  zhalf      (plond,plev+1) 
    REAL(kind=r8), INTENT(IN   ) ::  psfcpa     (plond)
    REAL(kind=r8), INTENT(INOUT) ::  USTAR      (plond)
    REAL(kind=r8), INTENT(IN   ) ::  TSK        (plond)
    REAL(kind=r8), INTENT(IN   ) ::  QSFC       (plond)
    REAL(kind=r8), INTENT(IN   ) ::  z0         (plond)
    REAL(kind=r8), INTENT(IN   ) ::  rpdel      (plond,plev)     ! 1./pdel (thickness between interfaces)
    REAL(kind=r8), INTENT(IN   ) ::  rpdeli     (plond,plev)     ! 1./pdeli (thickness between midpoints)
    REAL(kind=r8), INTENT(IN   ) ::  um1        (plond,plev)     ! u-wind input
    REAL(kind=r8), INTENT(IN   ) ::  vm1        (plond,plev)     ! v-wind input
    REAL(kind=r8), INTENT(IN   ) ::  tm1        (plond,plev)     ! temperature input
    REAL(kind=r8), INTENT(IN   ) ::  taux       (plond)          ! x surface stress (n)
    REAL(kind=r8), INTENT(IN   ) ::  tauy       (plond)          ! y surface stress (n)
    REAL(kind=r8), INTENT(IN   ) ::  shflx      (plond)          ! surface sensible heat flux (w/m2)
    REAL(kind=r8), INTENT(IN   ) ::  cflx       (plond,pcnst)    ! surface constituent flux (kg/m2/s)
    REAL(kind=r8), INTENT(IN   ) ::  ztodt                       ! 2 delta-t
    REAL(kind=r8), INTENT(IN   ) ::  colrad     (plond)    
    !
    ! Input/output arguments
    !
    REAL(kind=r8), INTENT(INOUT) :: qm1(plond,plev,pcnst)  ! initial/final constituent field
    !
    ! Output arguments
    !
    REAL(kind=r8), INTENT(OUT  ) :: dtv(plond,plev)        ! temperature tendency (heating)
    REAL(kind=r8), INTENT(OUT  ) :: dqv(plond,plev,pcnst)  ! constituent diffusion tendency

    REAL(kind=r8), INTENT(OUT  ) :: duv(plond,plev)        ! u-wind tendency
    REAL(kind=r8), INTENT(OUT  ) :: dvv(plond,plev)        ! v-wind tendency
    REAL(kind=r8), INTENT(OUT  ) :: up1(plond,plev)        ! u-wind after vertical diffusion
    REAL(kind=r8), INTENT(OUT  ) :: vp1(plond,plev)        ! v-wind after vertical diffusion
    REAL(kind=r8), INTENT(OUT  ) :: pblh(plond)            ! planetary boundary layer height
    REAL(kind=r8), INTENT(OUT  ) :: tpert(plond)           ! convective temperature excess
    REAL(kind=r8), INTENT(OUT  ) :: qpert(plond,pcnst)     ! convective humidity and constituent excess
    REAL(kind=r8), INTENT(INOUT) :: TKE(plond,plev+1)
    REAL(kind=r8), INTENT(INOUT) :: rino(plond,plev)        ! bulk Richardson no. from level to ref lev
    REAL(kind=r8), INTENT(OUT  ) :: obklen(plond)           ! Obukhov length
    REAL(kind=r8), INTENT(OUT  ) :: phiminv(plond)          ! inverse phi function for momentum
    REAL(kind=r8), INTENT(OUT  ) :: phihinv(plond)          ! inverse phi function for heat 
    REAL(kind=r8), INTENT(INOUT) :: kvh(plond,plev + 1)        ! diffusion coefficient for heat
    REAL(kind=r8), INTENT(INOUT) :: kvm(plond,plev + 1)        ! diffusion coefficient for momentum
    REAL(kind=r8), INTENT(OUT  ) ::tstar(plond)    
    REAL(kind=r8), INTENT(OUT  ) ::wstar(plond)    

    !
    !---------------------------Local storage-------------------------------
    !
    INTEGER :: i,k,m               ! longitude,level,constituent indices
    REAL(kind=r8) ::  denom                  ! denominator of expression
    REAL(kind=r8) ::  qp1   (plond,plev,pcnst)  ! constituents after vdiff
    REAL(kind=r8) ::  thp   (plond,plev)        ! potential temperature after vdiff
    REAL(kind=r8) ::  cgs   (plond,plev + 1)       ! counter-gradient star (cg/flux)
    REAL(kind=r8) ::  rztodt                 ! 1./ztodt
    REAL(kind=r8) ::  br    (plond)
    REAL(kind=r8) ::  THETA1(plond) 
    REAL(kind=r8) ::  TVS   (plond) 
    REAL(kind=r8) ::  THV1  (plond) 
    REAL(kind=r8) ::  DTV2  (plond) 
    REAL(kind=r8) ::  PS1   (plond) 
    REAL(kind=r8) ::  PS    (plond) 
    REAL(kind=r8) ::  Z1    (plond)
    REAL(kind=r8) ::  wspd  (plond)
    REAL(kind=r8) ::  TV1   (plond)
    REAL(KIND=r8),PARAMETER :: con_rv               =4.6150e+2_r8 ! gas constant H2O    (J/kg/K)
    REAL(KIND=r8),PARAMETER :: con_rd               =2.8705e+2_r8 ! gas constant air    (J/kg/K)
    REAL(KIND=r8),PARAMETER :: RVRDM1               =con_rv/con_rd-1.0_r8
     
    REAL(KIND=r8) :: KmMixl (plond,plev) 
    REAL(KIND=r8) :: KhMixl (plond,plev)
    REAL(KIND=r8) :: AKH2   (1:plond,1:plev+1)
    REAL(KIND=r8) :: AKM2   (1:plond,1:plev+1)
    REAL(KIND=r8) :: el2    (plond  ,plev + 1)         ! el       (1:nCols,1:kMax+1)!       
    REAL(KIND=r8) :: el0    (plond)
    REAL(KIND=r8) :: rrho   (1:plond)
    REAL(KIND=r8) :: ustr        
    REAL(KIND=r8) :: twodt

    !
    dtv  =0.0_r8! temperature tendency (heating)
    dqv  =0.0_r8! constituent diffusion tendency
    duv  =0.0_r8      ! u-wind tendency
    dvv  =0.0_r8      ! v-wind tendency
    obklen=0.0_r8           ! Obukhov length
    phiminv =0.0_r8        ! inverse phi function for momentum
    phihinv =0.0_r8        ! inverse phi function for heat 
    tstar   =0.0_r8  
    wstar   =0.0_r8  
    cgs =0.0_r8  
    br    =0.0_r8  
    THETA1=0.0_r8  
    TVS   =0.0_r8  
    THV1  =0.0_r8  
    DTV2  =0.0_r8  
    PS1   =0.0_r8  
    PS    =0.0_r8  
    Z1    =0.0_r8  
    wspd  =0.0_r8  
    TV1   =0.0_r8  
    KmMixl=0.0_r8  
    KhMixl=0.0_r8  
    AKH2  =0.0_r8  
    AKM2  =0.0_r8  
    el2   =0.0_r8  
    el0   =0.0_r8  
    rrho  =0.0_r8  
    ustr  =0.0_r8  
    twodt=0.0_r8  
    up1    =um1
    vp1    =vm1
    thp    =tm1
    qp1    =qm1
    pblh   =0.0_r8
    ustar  =0.0_r8
    kvh    =0.0_r8
    kvm    =0.0_r8
    tpert  =0.0_r8
    qpert  =0.0_r8
    cgs    =0.0_r8
    KmMixl =0.0_r8
    KhMixl =0.0_r8
    AKH2   =0.0_r8
    AKM2   =0.0_r8
    el2    =0.0_r8
    el0    =0.0_r8
    twodt =ztodt

    !----------------------------------------------------------------------
    !**********************************************************************
    !----------------------------------------------------------------------    
    CALL MY20_TURB (       &
                   pcnst      ,&
                   plond      ,& ! INTENT(IN   ) :: plond
                   plev       ,& ! INTENT(IN   ) :: plev
                   twodt      ,& ! INTENT(IN   ) :: ztodt
                   colrad     ,& ! INTENT(IN   ) :: colrad    (plond)         
                   tm1        ,& ! INTENT(IN   ) :: tm1       (plond,plev)          ! temperature input
                   qp1        ,& ! INTENT(IN   ) :: qp1       (plond,plev,pcnst)  ! constituents after vdiff
                   um1        ,& ! INTENT(IN   ) :: um1       (plond,plev)          ! u-wind input
                   vm1        ,& ! INTENT(IN   ) :: vm1       (plond,plev)          ! v-wind input
                   gl0        ,& ! INTENT(INOUT) :: gl0       (plond)          !
                   KmMixl     ,& ! INTENT(OUT  ) :: Pbl_KmMixl(plond,plev)   
                   KhMixl      )  ! INTENT(OUT  ) :: Pbl_KhMixl(plond,plev)
    
    !----------------------------------------------------------------------
    !**********************************************************************
    !----------------------------------------------------------------------
    DO I=1,plond
       rrho(i)  = rair*tm1(i,plev)/pmidm1(i,plev)
       ustr     = SQRT(SQRT(taux(i)**2 + tauy(i)**2)*rrho(i))
       ustar(i) = MAX(ustr,0.01_r8)
    END DO
    CALL MY25_TURB(&
       plond     ,& ! nCols
       plev      ,& ! kMax
       twodt     ,& ! delt    !       delt     -  Time step in seconds
       um1       ,& ! um       (1:nCols,1:kMax)!Zonal wind components
       vm1       ,& ! vm       (1:nCols,1:kMax)!Meridional Wind components
       thm       ,& ! theta    (1:nCols,1:kMax)!theta         -  Potential temperature
       zhalf     ,& ! zhalf    (1:nCols,1:kMax+1)!zfull         -  Height at full levels
       zm        ,& ! zfull    (1:nCols,1:kMax)!         zhalf    -  Height at half levels
       pintm1    ,& ! phalf    (1:nCols,1:kMax+1)!!          phalf    -  Pressure at half levels, 
       pmidm1    ,& ! pfull    (1:nCols,1:kMax)!!          pfull    -  Pressure at full levels! , zhalf, zfull
       z0        ,& ! z0       (1:nCols)!          z0           -  Roughness length
       FRLAND    ,& ! fracland (1:nCols)!        fracland -  Fractional amount of land beneath a grid box
       ustar     ,& ! ustar    (1:nCols)!        ustar         -  OPTIONAL:friction velocity (m/sec)
       el0       ,& ! el0      (1:nCols)!       el0  -  characteristic length scale
       akm2      ,& ! akm      (1:nCols,1:kMax+1)!       akm  -  mixing coefficient for momentum
       akh2      ,& ! akh      (1:nCols,1:kMax+1)!       akh  -  mixing coefficient for heat and moisture
       el2        ,& ! el       (1:nCols,1:kMax+1)!       el   -  master length scale
       TKE        ) !TKE       (1:nCols,1:kMax+1)!
!    !----------------------------------------------------------------------
!    !**********************************************************************
!    !----------------------------------------------------------------------

    DO i=1,plond
       wspd  (i) = SQRT((um1(i,plev)**2) + (vm1(i,plev)**2))
       wspd  (i) = MAX(wspd  (i),0.25_r8)
       THETA1(i) = tsk   (i) * psomc(i,plev)       
       TVS   (i) = tsk   (i) * (1.0_r8 + RVRDM1 * QSFC(i))
       THV1  (i) = THETA1(i) * (1.0_r8 + RVRDM1 * QSFC(i))
       TV1   (i) = tm1 (i,plev) * (1.0_r8 + RVRDM1 * QSFC(i))
       PS1   (i) = pmidm1(i,plev)
       PS    (i) = psfcpa(i)
       DTV2   (i) = THV1(i) - TVS(i)
       Z1    (i) = -con_rd * TV1(i) * LOG(PS1(i)/PS(i)) / grav
       br    (i) = grav * DTV2(i) * Z1(i) / (0.5_r8 * (THV1(i) + TVS(i))&
     &          * wspd(i) * wspd(i))
       br(i) = MAX(br(i),-5000.0_r8)
    END DO                  
    !
    !-----------------------------------------------------------------------
    !
    ! Call vertical diffusion code. No setup work is required.
    !
    CALL vdiff(bstar      , &! INTENT(IN        ) :: 
               FRLAND     , &! INTENT(IN        ) :: 
               LwCoolRate , &! INTENT(IN        ) :: 
               LwCoolRateC, &! INTENT(IN        ) :: 
               cldtot     , &! INTENT(IN        ) :: 
               qliq       , &! INTENT(IN        ) :: 
               um1        , &! INTENT(IN        ) :: um1(plond,plev)            ! u wind input
               vm1        , &! INTENT(IN        ) :: vm1(plond,plev)            ! v wind input
               tm1        , &! INTENT(IN        ) :: tm1(plond,plev)            ! temperature input
               qm1        , &! INTENT(IN        ) :: qm1(plond,plev,pcnst)  ! moisture and trace constituent input
               pmidm1  , &! INTENT(IN        ) :: pmidm1(plond,plev)     ! midpoint pressures
               pintm1  , &! INTENT(IN        ) :: pintm1(plond,plev + 1)    ! interface pressures
               rpdel   , &! INTENT(IN        ) :: rpdel(plond,plev)      ! 1./pdel  (thickness bet interfaces)
               rpdeli  , &! INTENT(IN        ) :: rpdeli(plond,plev)     ! 1./pdeli (thickness bet midpoints)
               ztodt   , &! INTENT(IN        ) :: ztodt                    ! 2 delta-t
               thm     , &! INTENT(IN        ) :: thm(plond,plev)            ! potential temperature
               zm      , &! INTENT(IN        ) :: zm(plond,plev)            ! midpoint geoptl height above sfc
               taux    , &! INTENT(IN        ) :: taux(plond)            ! x surface stress (n)
               tauy    , &! INTENT(IN        ) :: tauy(plond)            ! y surface stress (n)
               shflx   , &! INTENT(IN        ) :: shflx(plond)            ! surface sensible heat flux (w/m2)
               cflx    , &! INTENT(IN        ) :: cflx(plond,pcnst)      ! surface constituent flux (kg/m2/s)
               up1     , &! INTENT(OUT  ) :: up1(plond,plev)            ! u-wind after vertical diffusion
               vp1     , &! INTENT(OUT  ) :: vp1(plond,plev)            ! v-wind after vertical diffusion
               thp     , &! INTENT(OUT  ) :: thp(plond,plev)            ! pot temp after vert. diffusion
               qp1     , &! INTENT(OUT  ) :: qp1(plond,plev,pcnst)  ! moist, tracers after vert. diff
               pblh    , &! INTENT(OUT  ) :: pblh(plond)            ! planetary boundary layer height
               ustar   , &! INTENT(OUT  ) :: ustar(plond)            ! surface friction velocity
               kvh     , &! INTENT(OUT  ) :: kvh(plond,plev + 1)       ! coefficient for heat and tracers
               kvm     , &! INTENT(OUT  ) :: kvm(plond,plev + 1)       ! coefficient for momentum
               tpert   , &! INTENT(OUT  ) :: tpert(plond)            ! convective temperature excess
               qpert   , &! INTENT(OUT  ) :: qpert(plond)            ! convective humidity excess
               cgs     , &! INTENT(OUT  ) :: cgs(plond,plev + 1)       ! counter-grad star (cg/flux)
               plon    , &! INTENT(IN        ) :: plon        ! number of longitudes
               plev    , &! INTENT(IN        ) :: plev        ! number of vertical levels
               plond   , &! INTENT(IN        ) :: plond        ! slt extended domain longitude
               pcnst   , &! INTENT(IN        ) :: pcnst        ! number of constituents (including water vapor)
               fakn    , &! INTENT(IN        ) :: fakn    ! Constant in turbulent prandtl number
               ricr    , & ! INTENT(IN        ) :: ricr    ! Critical richardson number
               tsk     , &
               qsfc    , &
               KmMixl  , &
               KhMixl  , &
               AKH2     , &! akh      (1:nCols,1:kMax+1)!  mixing coefficient for heat and moisture 
               AKM2     , &! akm      (1:nCols,1:kMax+1)!  mixing coefficient for momentum
               rino    , &        ! bulk Richardson no. from level to ref lev
               psomc   , &
               obklen  , &
               phiminv , &
               phihinv , &
               tstar   , &
               wstar        ) 
    !----------------------------------------------------------------------
    !**********************************************************************
    !----------------------------------------------------------------------
 
    !
    ! Convert the diffused fields back to diffusion tendencies.
    ! Add the diffusion tendencies to the cummulative physics tendencies,
    ! except for constituents. The diffused values of the constituents
    ! replace the input values.
    !
    rztodt = 1.0_r8/ztodt
    DO k=ntopfl,plev
       DO i=1,plon
          duv(i,k) = (up1(i,k)*SIN( colrad(i)) - um1(i,k)*SIN( colrad(i)))*rztodt
          dvv(i,k) = (vp1(i,k)*SIN( colrad(i)) - vm1(i,k)*SIN( colrad(i)))*rztodt

          denom    = cpair*(1.0_r8 + cpvir*qm1(i,k,1))
          dtv(i,k) = (thp(i,k) - thm(i,k))*rztodt

       END DO
       DO m=1,pcnst
          DO i=1,plon
             dqv(i,k,m) = (qp1(i,k,m) - qm1(i,k,m))*rztodt
          END DO
       END DO
    END DO    
    RETURN
  END SUBROUTINE vdintr

  !
  !----------------------------------------------------------------------
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !---------------------------------------------------------------------

  SUBROUTINE vdiff(bstar , &! INTENT(IN   ) :: 
                   FRLAND, &! INTENT(IN   ) :: 
                   LwCoolRate  , &! INTENT(IN   ) :: 
                   LwCoolRateC , &! INTENT(IN   ) :: 
                   cldtot      , &! INTENT(IN   ) :: 
                   qliq    , &! INTENT(IN   ) :: 
                   um1     , &! INTENT(IN   ) :: um1(plond,plev)        ! u wind input
                   vm1     , &! INTENT(IN   ) :: vm1(plond,plev)        ! v wind input
                   tm1     , &! INTENT(IN   ) :: tm1(plond,plev)        ! temperature input
                   qm1     , &! INTENT(IN   ) :: qm1(plond,plev,pcnst)  ! moisture and trace constituent input
                   pmidm1  , &! INTENT(IN   ) :: pmidm1(plond,plev)     ! midpoint pressures
                   pintm1  , &! INTENT(IN   ) :: pintm1(plond,plev + 1)    ! interface pressures
                   rpdel   , &! INTENT(IN   ) :: rpdel(plond,plev)      ! 1./pdel  (thickness bet interfaces)
                   rpdeli  , &! INTENT(IN   ) :: rpdeli(plond,plev)     ! 1./pdeli (thickness bet midpoints)
                   ztodt   , &! INTENT(IN   ) :: ztodt                  ! 2 delta-t
                   thm     , &! INTENT(IN   ) :: thm(plond,plev)        ! potential temperature
                   zm      , &! INTENT(IN   ) :: zm(plond,plev)         ! midpoint geoptl height above sfc
                   taux           , &! INTENT(IN   ) :: taux(plond)            ! x surface stress (n)
                   tauy           , &! INTENT(IN   ) :: tauy(plond)            ! y surface stress (n)
                   shflx   , &! INTENT(IN   ) :: shflx(plond)           ! surface sensible heat flux (w/m2)
                   cflx    , &! INTENT(IN   ) :: cflx(plond,pcnst)      ! surface constituent flux (kg/m2/s)
                   up1     , &! INTENT(OUT  ) :: up1(plond,plev)        ! u-wind after vertical diffusion
                   vp1     , &! INTENT(OUT  ) :: vp1(plond,plev)        ! v-wind after vertical diffusion
                   thp     , &! INTENT(OUT  ) :: thp(plond,plev)        ! pot temp after vert. diffusion
                   qp1     , &! INTENT(OUT  ) :: qp1(plond,plev,pcnst)  ! moist, tracers after vert. diff
                   pblh    , &! INTENT(OUT  ) :: pblh(plond)            ! planetary boundary layer height
                   ustar   , &! INTENT(OUT  ) :: ustar(plond)           ! surface friction velocity
                   kvh     , &! INTENT(OUT  ) :: kvh(plond,plev + 1)       ! coefficient for heat and tracers
                   kvm           , &! INTENT(OUT  ) :: kvm(plond,plev + 1)       ! coefficient for momentum
                   tpert   , &! INTENT(OUT  ) :: tpert(plond)           ! convective temperature excess
                   qpert   , &! INTENT(OUT  ) :: qpert(plond)           ! convective humidity excess
                   cgs     , &! INTENT(OUT  ) :: cgs(plond,plev + 1)       ! counter-grad star (cg/flux)
                   plon    , &! INTENT(IN   ) :: plon       ! number of longitudes
                   plev    , &! INTENT(IN   ) :: plev       ! number of vertical levels
                   plond   , &! INTENT(IN   ) :: plond      ! slt extended domain longitude
                   pcnst   , &! INTENT(IN   ) :: pcnst      ! number of constituents (including water vapor)
                   fakn    , &! INTENT(IN   ) :: fakn    ! Constant in turbulent prandtl number
                   ricr    , &! INTENT(IN   ) :: ricr    ! Critical richardson number
                   tsk     , &
                   qsfc    , &
                   KmMixl  , &
                   KhMixl  , &
                   AKH     , &
                   AKM     , &
                   rino    , &
                   psomc   , &
                   obklen  , &
                   phiminv , &
                   phihinv , &
                   tstar   , &
                   wstar   )
    !-----------------------------------------------------------------------
    !
    ! Driver routine to compute vertical diffusion of momentum,
    ! moisture, trace constituents and potential temperature.
    !
    ! Free atmosphere diffusivities are computed first; then modified
    ! by the boundary layer scheme; then passed to individual
    ! parameterizations mvdiff, qvdiff
    !
    ! The free atmosphere diffusivities are based on standard mixing length 
    ! forms for the neutral diffusivity multiplied by functns of Richardson 
    ! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for 
    ! momentum, potential temperature, and constitutents.
    ! The stable Richardson num function (Ri>0) is taken from Holtslag and 
    ! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri))
    ! The unstable Richardson number function (Ri<0) is taken from  CCM1.
    ! f = sqrt(1 - 18*Ri)
    !
    !---------------------------Code history--------------------------------
    !
    ! Original version:  CCM1
    ! Standardized:      J. Rosinski, June 1992
    ! Reviewed:          P. Rasch, B. Boville, August 1992
    ! Reviewed:          P. Rasch, March 1996
    ! Reviewed:          B. Boville, April 1996
    !
    !-----------------------------------------------------------------------
    !
    ! $Id: vdiff.F,v 1.1.1.1 2001/03/09 00:29:28 mirin Exp $
    !
    !-----------------------------------------------------------------------
    !      use precision
    !      use pmgrid
    !#if ( ! defined CRAY )
    !      use srchutil, only: wheneq
    !#endif
    !-----------------------------------------------------------------------
    !#include <implicit.h>
    !------------------------------Commons----------------------------------
    !#include <comvd.h>
    !-----------------------------------------------------------------------
    !#include <comqmin.h>
    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) ::  plon       ! number of longitudes
    INTEGER, INTENT(IN   ) :: plev       ! number of vertical levels
    INTEGER, INTENT(IN   ) :: plond      ! slt extended domain longitude
    INTEGER, INTENT(IN   ) :: pcnst      ! number of constituents (including water vapor)
    REAL(KIND=r8), INTENT(IN   ) :: fakn    ! Constant in turbulent prandtl number
    REAL(KIND=r8), INTENT(IN   ) :: ricr    ! Critical richardson number

    REAL(KIND=r8), INTENT(IN   ) :: bstar      (plond)      
    REAL(KIND=r8), INTENT(IN   ) :: FRLAND  (plond)      
    REAL(KIND=r8), INTENT(IN   ) :: LwCoolRate (plond,plev)        
    REAL(KIND=r8), INTENT(IN   ) :: LwCoolRateC(plond,plev)        
    REAL(KIND=r8), INTENT(IN   ) :: cldtot     (plond,plev)      
    REAL(KIND=r8), INTENT(IN   ) :: qliq       (plond,plev)      
    REAL(KIND=r8), INTENT(IN   ) :: um1(plond,plev)        ! u wind input
    REAL(KIND=r8), INTENT(IN   ) :: vm1(plond,plev)        ! v wind input
    REAL(KIND=r8), INTENT(IN   ) :: tm1(plond,plev)        ! temperature input
    REAL(KIND=r8), INTENT(IN   ) :: qm1(plond,plev,pcnst)  ! moisture and trace constituent input
    REAL(KIND=r8), INTENT(IN   ) :: pmidm1(plond,plev)     ! midpoint pressures
    REAL(KIND=r8), INTENT(IN   ) :: pintm1(plond,plev + 1)    ! interface pressures
    REAL(KIND=r8), INTENT(IN   ) :: rpdel(plond,plev)      ! 1./pdel  (thickness bet interfaces)
    REAL(KIND=r8), INTENT(IN   ) :: rpdeli(plond,plev)     ! 1./pdeli (thickness bet midpoints)
    REAL(KIND=r8), INTENT(IN   ) :: ztodt                  ! 2 delta-t
    REAL(KIND=r8), INTENT(IN   ) :: thm(plond,plev)        ! potential temperature
    REAL(KIND=r8), INTENT(IN   ) :: zm(plond,plev)         ! midpoint geoptl height above sfc
    REAL(KIND=r8), INTENT(IN   ) :: taux(plond)            ! x surface stress (n)
    REAL(KIND=r8), INTENT(IN   ) :: tauy(plond)            ! y surface stress (n)
    REAL(KIND=r8), INTENT(IN   ) :: shflx(plond)           ! surface sensible heat flux (w/m2)
    REAL(KIND=r8), INTENT(IN   ) :: cflx(plond,pcnst)      ! surface constituent flux (kg/m2/s)
    REAL(KIND=r8), INTENT(IN   ) :: tsk   (plond)     
    REAL(KIND=r8), INTENT(IN   ) :: qsfc   (plond)    
    REAL(KIND=r8), INTENT(IN   ) :: KmMixl(plond,plev) 
    REAL(KIND=r8), INTENT(IN   ) :: KhMixl(plond,plev) 
    REAL(KIND=r8), INTENT(IN   ) :: AKH(plond,plev+1) 
    REAL(KIND=r8), INTENT(IN   ) :: AKM(plond,plev+1) 
    REAL(KIND=r8), INTENT(IN   ) ::  psomc(plond,plev)      ! (psm1/pmidm1)**cappa

    !
    ! Output arguments
    !
    REAL(KIND=r8), INTENT(OUT  ) :: up1     (plond,plev)        ! u-wind after vertical diffusion
    REAL(KIND=r8), INTENT(OUT  ) :: vp1     (plond,plev)        ! v-wind after vertical diffusion
    REAL(KIND=r8), INTENT(OUT  ) :: thp     (plond,plev)        ! pot temp after vert. diffusion
    REAL(KIND=r8), INTENT(OUT  ) :: qp1     (plond,plev,pcnst)  ! moist, tracers after vert. diff
    REAL(KIND=r8), INTENT(INOUT) :: pblh    (plond)            ! planetary boundary layer height
    REAL(KIND=r8), INTENT(INOUT) :: rino    (plond,plev)        ! bulk Richardson no. from level to ref lev
    REAL(KIND=r8), INTENT(OUT  ) :: ustar   (plond)           ! surface friction velocity
    REAL(KIND=r8), INTENT(OUT  ) :: kvh     (plond,plev + 1)       ! coefficient for heat and tracers
    REAL(KIND=r8), INTENT(OUT  ) :: kvm     (plond,plev + 1)       ! coefficient for momentum
    REAL(KIND=r8), INTENT(OUT  ) :: tpert   (plond)           ! convective temperature excess
    REAL(KIND=r8), INTENT(OUT  ) :: qpert   (plond)           ! convective humidity excess
    REAL(KIND=r8), INTENT(OUT  ) :: cgs     (plond,plev + 1)       ! counter-grad star (cg/flux)
    REAL(KIND=r8), INTENT(OUT  ) :: obklen  (plond)           ! Obukhov length
    REAL(KIND=r8), INTENT(OUT  ) :: phiminv (plond)          ! inverse phi function for momentum
    REAL(KIND=r8), INTENT(OUT  ) :: phihinv (plond)          ! inverse phi function for heat 
    REAL(KIND=r8), INTENT(OUT  ) :: tstar(plond)
    REAL(KIND=r8), INTENT(OUT  ) :: wstar(plond)
    !
    !---------------------------Local workspace-----------------------------
    !
    REAL(KIND=r8) cah   (plond,plev)        ! -upper diag for heat and constituts
    REAL(KIND=r8) cam   (plond,plev)        ! -upper diagonal for momentum
    REAL(KIND=r8) cch   (plond,plev)        ! -lower diag for heat and constits
    REAL(KIND=r8) ccm   (plond,plev)        ! -lower diagonal for momentum
    REAL(KIND=r8) cgh   (plond,plev + 1)       ! countergradient term for heat
    REAL(KIND=r8) cgq   (plond,plev + 1,pcnst) ! countergrad term for constituent
    REAL(KIND=r8) dvdz2                  ! (du/dz)**2 + (dv/dz)**2
    REAL(KIND=r8) dz                     ! delta z between midpoints
    REAL(KIND=r8) fstab                  ! stable f(ri)
    REAL(KIND=r8) funst                  ! unstable f(ri)
    REAL(KIND=r8) kvf   (plond,plev + 1)       ! free atmosphere kv at interfaces
    REAL(KIND=r8) rinub                  ! richardson no=(g/theta)(dtheta/dz)/
    !                                                 (du/dz**2+dv/dz**2)
    REAL(KIND=r8) sstab                  ! static stability = g/th  * dth/dz
    REAL(KIND=r8) potbar(plond,plev + 1)    ! pintm1(k)/(.5*(tm1(k)+tm1(k-1))
    REAL(KIND=r8) tmp1  (plond)            ! temporary storage
    REAL(KIND=r8) tmp2                   ! temporary storage
    REAL(KIND=r8) rcpair                 ! 1./cpair
    REAL(KIND=r8) ztodtgor               ! ztodt*gravit/rair
    REAL(KIND=r8) gorsq                  ! (gravit/rair)**2
    REAL(KIND=r8) dubot (plond)           ! lowest layer u change from stress
    REAL(KIND=r8) dvbot (plond)           ! lowest layer v change from stress
    REAL(KIND=r8) dtbot (plond)           ! lowest layer t change from heat flx
    REAL(KIND=r8) dqbot (plond,pcnst)     ! lowest layer q change from const flx
    REAL(KIND=r8) thx   (plond,plev)        ! temperature input + counter gradient
    REAL(KIND=r8) thv   (plond,plev)        ! virtual potential temperature
    REAL(KIND=r8) qmx   (plond,plev,pcnst)  ! constituents input + counter grad
    REAL(KIND=r8) zeh   (plond,plev)        ! term in tri-diag. matrix system (T & Q)
    REAL(KIND=r8) zem   (plond,plev)        ! term in tri-diag. matrix system (momentum)
    REAL(KIND=r8) termh (plond,plev)      ! 1./(1. + cah(k) + cch(k) - cch(k)*zeh(k-1))
    REAL(KIND=r8) termm (plond,plev)      ! 1./(1. + cam(k) + ccm(k) - ccm(k)*zem(k-1))
    REAL(KIND=r8) kvn                    ! neutral Kv
    REAL(KIND=r8) ksx   (plond)             ! effective surface drag factor (x)
    REAL(KIND=r8) ksy   (plond)             ! effective surface drag factor (y)
    REAL(KIND=r8) sufac (plond)           ! lowest layer u implicit stress factor
    REAL(KIND=r8) svfac (plond)           ! lowest layer v implicit stress factor
    INTEGER indx(plond)          ! array of indices of potential q<0
    INTEGER ilogic(plond)        ! 1 => adjust vertical profile
    INTEGER nval                ! num of values which meet criteria
    INTEGER ii                  ! longitude index of found points
    INTEGER i                   ! longitude index
    INTEGER k                   ! vertical index
    INTEGER m                   ! constituent index
    INTEGER ktopbl(plond)       ! index of first midpoint inside pbl
    INTEGER ktopblmn            ! min value of ktopbl
    up1=0.0_r8;vp1  =0.0_r8; thp  =0.0_r8;  qp1=0.0_r8;  pblh =0.0_r8;  rino =0.0_r8;    ustar=0.0_r8
    kvh  =0.0_r8;    kvm  =0.0_r8;    tpert=0.0_r8;    qpert=0.0_r8;    cgs  =0.0_r8;    obklen=0.0_r8;   phiminv=0.0_r8;  phihinv=0.0_r8
    cah =0.0_r8;   cam =0.0_r8;   cch =0.0_r8;   ccm =0.0_r8;   cgh =0.0_r8;   cgq =0.0_r8;   dvdz2=0.0_r8;  dz   =0.0_r8
    fstab=0.0_r8;  funst=0.0_r8;  kvf  =0.0_r8;  rinub=0.0_r8;  sstab =0.0_r8; potbar=0.0_r8; tmp1  =0.0_r8; tmp2 =0.0_r8
    rcpair  =0.0_r8;   ztodtgor =0.0_r8;  gorsq =0.0_r8;  dubot =0.0_r8;  dvbot =0.0_r8;  dtbot =0.0_r8
    dqbot =0.0_r8;  thx   =0.0_r8;  thv   =0.0_r8;  qmx   =0.0_r8;  zeh   =0.0_r8;  zem   =0.0_r8;  termh =0.0_r8;  termm =0.0_r8
    kvn   =0.0_r8;  ksx   =0.0_r8;  ksy   =0.0_r8;  sufac =0.0_r8;  svfac =0.0_r8
    !
    !-----------------------------------------------------------------------
    !
    ! Convert the surface fluxes to lowest level tendencies.
    ! Stresses are converted to effective drag coefficients if these are >0
    !
    rcpair = 1.0_r8/cpair
    DO i=1,plond
       tmp1(i) = ztodt*gravit*rpdel(i,plev)
       IF(um1(i,plev) == 0.0_r8)THEN
          ksx(i) = -taux(i) / 0.00001_r8
       ELSE
          ksx(i) = -taux(i) / um1(i,plev)
       END IF
       IF (ksx(i) .GT. 0.0_r8) THEN
          sufac(i) = tmp1(i) * ksx(i)
          dubot(i) = 0.0_r8
       ELSE
          ksx(i)   = 0.0_r8
          sufac(i) = 0.0_r8
          dubot(i) = tmp1(i) * taux(i)
       END IF
       IF(vm1(i,plev) == 0.0_r8)THEN
          ksy(i) = -tauy(i) / 0.00001_r8
       ELSE
          ksy(i) = -tauy(i) / vm1(i,plev)
       END IF
       IF (ksy(i) .GT. 0.0_r8) THEN
          svfac(i) = tmp1(i) * ksy(i)
          dvbot(i) = 0.0_r8
       ELSE
          ksy(i)   = 0.0_r8
          svfac(i) = 0.0_r8
          dvbot(i) = tmp1(i) * tauy(i)
       END IF

       dqbot(i,1)   = cflx(i,1)*tmp1(i)
       dtbot(i)     = shflx(i)*tmp1(i)*rcpair
       kvf(i,plev + 1) = 0.0_r8
    END DO
    DO m=2,pcnst
       DO i=1,plond
          dqbot(i,m) = cflx(i,m)*tmp1(i)
       END DO
    END DO
    !
    ! Set the vertical diffusion coefficient above the top diffusion level
    !
    DO k=1,ntopfl
       DO i=1,plond
          kvf(i,k) = 0.0_r8
       END DO
    END DO
    !
    ! Compute virtual potential temperature for use in static stability 
    ! calculation.  0.61 is 1. - R(water vapor)/R(dry air).  Use 0.61 instead
    ! of a computed variable in order to obtain an identical simulation to
    ! Case 414.
    !
    CALL virtem(thm     , &! INTENT(IN   )
                qm1     , &! INTENT(IN   )
                0.61_r8 , &! INTENT(IN   )
                thv     , &! INTENT(OUT  )
                plond   , &! INTENT(IN   )
                plev    , &! INTENT(IN   )
                plon      )! INTENT(IN   )
    !
    ! Compute the free atmosphere vertical diffusion coefficients
    ! kvh = kvq = kvm. 
    !
    DO k=ntopfl,plev-1
       DO i=1,plon
          !
          ! Vertical shear squared, min value of (delta v)**2 prevents zero shear.
          !
          dvdz2 = (um1(i,k)-um1(i,k+1))**2 + (vm1(i,k)-vm1(i,k+1))**2
          !
          dvdz2 = MAX(dvdz2,1.e-36_r8)
          dz    = zm(i,k) - zm(i,k+1)
          !
          !            _             _ 2        _              _ 2
          !           |               |        |               |
          !           |U(z) - U(z+1)) |     +  |V(z) - V(z+1)) |  
          !           |_             _|        |_             _|  
          !  dV^2      
          ! ------ = -------------------------------------------
          !  dz^2                _             _ 2   
          !                     |               |
          !                     |Z(z) - Z(z+1)) |  
          !                     |_             _|  
          !
          !
          dvdz2 = dvdz2/(dz**2)
          !
          ! Static stability (use virtual potential temperature)
          !
          !                _                                              _
          !               |                                                |
          !               |         THETA(z) - THETA(z+1)                  |
          !  sstab = 2*g* |----------------------------------------------- |
          !               |  _                    _     _              _   |
          !               | |                      |   |                |  |
          !               | |THETA(z) + THETA(z+1) | * |Z(z) - Z(z+1))  |  |
          !               | |_                    _|   |_              _|  |
          !               |_                                              _|
          !
          sstab = gravit*2.0_r8*( thv(i,k) - thv(i,k+1))/ &
                                ((thv(i,k) + thv(i,k+1))*dz)
          !
          ! Richardson number, stable and unstable modifying functions
          !
          !                 _                                                                          _
          !                |                                         _               _                  |    
          !                |                                        |                 |                 |   
          !                |             THETA(z) - THETA(z+1)  *   | Z(z) - Z(z+1))  |                 |   
          !                !                                        |_               _|                 !   
          !rinub =    2*g* |--------------------------------------------------------------------------- |  
          !                |                             _                                           _  |
          !                |  _                    _    |  _              _ 2      _             _ 2  | |
          !                | |                      |   | |                |      |               |   | |
          !                | |THETA(z) + THETA(z+1) | * | | U(z) - U(z+1)) |   +  |V(z) - V(z+1)) |   | |
          !                | |_                    _|   | |_              _|      |_             _|   | |
          !                |                            |_                                           _| |
          !                |_                                                                          _|
          !
          rinub = sstab/dvdz2
          !
          !
          !                              1 
          ! fstab = ------------------------------------------------
          !             _                                 _
          !            |                 _           _    | 
          !            |                |             |   | 
          !            | 1 + 10*rinub * | 1 + 8*rinub |   | 
          !            |                |_           _|   | 
          !            |_                                _|
          !
          fstab = 1.0_r8/(1.0_r8 + 10.0_r8*rinub*(1.0_r8 + 8.0_r8*rinub))
          !
          !
          !
          funst = MAX(1.0_r8 - 18.0_r8*rinub,0.0_r8)
          !
          ! Select the appropriate function of the richardson number
          !
          !         -                                              - 0.5
          !        |                                                |
          !        |                                                |
          !        |                         1                      |
          ! fstab =|------------------------------------------------|
          !        |    _                                _          |
          !        |   |                 _           _    |         | 
          !        |   |                |             |   |         | 
          !        |   | 1 + 10*rinub * | 1 + 8*rinub |   |         | 
          !        |   |                |_           _|   |         | 
          !        |   |_                                _|         | 
          !        |_                                              _|
          !
          IF (rinub < 0.0_r8) fstab = SQRT(funst)
          !
          !            Km 
          !Km = l^2 * -----
          !            l^2
          !
          !
          ! Neutral diffusion coefficient
          ! compute mixing length (z), where z is the interface height estimated
          ! with an 8 km scale height.
          !
          !      0 > l < ml2
          !
          !
          ! Set the square of the mixing lengths.
          !
          !IF(k==1)THEN
          !   MixLgh(i,k) = 0.0_r8
          !ELSE IF(k==plev + 1)THEN
          !   MixLgh(i,k) = 0.0_r8
          !ELSE
          !   MixLgh(i,k) = MIN(MAX(100.0_r8,MixLgh(i,k)),ml2(k))
          !END IF 
          !IF(XLAND(i) >1.0_r8)THEN
          ! OCEAN
          !   kvn = MixLgh(i,k)*SQRT(dvdz2)
          !ELSE
          !   ! LAND
          kvn = ml2(k)*SQRT(dvdz2)
          !END IF
          !
          ! Full diffusion coefficient (modified by f(ri)),
          !
          kvf(i,k+1) = MAX(zkmin,kvn*fstab)
          !kvf(i,k+1)=MIN(gkh1,MAX(gkh0,kvn*fstab)) !(m/sec**2)
       END DO
    END DO
    !
    ! Determine the boundary layer kvh (=kvq), kvm, 
    ! counter gradient terms (cgh, cgq, cgs)
    ! boundary layer height (pblh) and 
    ! the perturbation temperature and moisture (tpert and qpert)
    ! The free atmosphere kv is returned above the boundary layer top.
    !
    CALL pbldif(thm     , &! INTENT(IN   ) ::  th(plond,plev)               ! potential temperature [K]
                qm1     , &! INTENT(IN   ) ::  q(plond,plev,pcnst)     ! specific humidity [kg/kg]
                zm      , &! INTENT(IN   ) ::  z(plond,plev)               ! height above surface [m]
                um1     , &! INTENT(IN   ) ::  u(plond,plev)               ! windspeed x-direction [m/s]
                vm1     , &! INTENT(IN   ) ::  v(plond,plev)               ! windspeed y-direction [m/s]
                tm1     , &! INTENT(IN   ) ::  t(plond,plev)               ! temperature (used for density)
                pmidm1  , &! INTENT(IN   ) ::  pmid(plond,plev)        ! midpoint pressures
                kvf     , &! INTENT(IN   ) ::  kvf(plond,plev + 1)     ! free atmospheric eddy diffsvty [m2/s]
                cflx    , &! INTENT(IN   ) ::  cflx(plond,pcnst)       ! surface constituent flux (kg/m2/s)
                shflx   , &! INTENT(IN   ) ::  shflx(plond)               ! surface heat flux (W/m2)
                taux    , &! INTENT(IN   ) ::  taux(plond)               ! surface u stress (N)
                tauy    , &! INTENT(IN   ) ::  tauy(plond)               ! surface v stress (N)
                ustar   , &! INTENT(OUT  ) ::  ustar(plond)               ! surface friction velocity [m/s]
                kvm     , &! INTENT(OUT  ) ::  kvm(plond,plev + 1)     ! eddy diffusivity for momentum [m2/s]
                kvh     , &! INTENT(OUT  ) ::  kvh(plond,plev + 1)     ! eddy diffusivity for heat [m2/s]
                cgh     , &! INTENT(OUT  ) ::  cgh(plond,plev + 1)     ! counter-gradient term for heat [K/m]
                cgq     , &! INTENT(OUT  ) ::  cgq(plond,plev + 1,pcnst)! counter-gradient term for constituents
                cgs     , &! INTENT(OUT  ) ::  cgs(plond,plev + 1)     ! counter-gradient star (cg/flux)
                pblh    , &! INTENT(OUT  ) ::  pblh(plond)               ! boundary-layer height [m]
                tpert   , &! INTENT(OUT  ) ::  tpert(plond)               ! convective temperature excess
                qpert   , &! INTENT(OUT  ) ::  qpert(plond)               ! convective humidity excess
                ktopbl  , &! INTENT(OUT  ) ::  ktopbl(plond)               ! index of first midpoint inside pbl
                ktopblmn, &! INTENT(OUT  ) ::  ktopblmn                ! min value of ktopbl
                plond   , &! INTENT(IN   ) :: plond         ! slt extended domain longitude
                pcnst   , &! INTENT(IN   ) :: pcnst         ! number of constituents (including water vapor)
                plon    , &! INTENT(IN   ) :: plon         ! number of longitudes
                plev    , &! INTENT(IN   ) ::  plev          ! number of vertical levels
                fakn    , &! INTENT(IN   ) :: fakn    ! Constant in turbulent prandtl number
                ricr    , &! INTENT(IN   ) :: ricr    ! Critical richardson number
                tsk     , &
                qsfc    , &
                rino    , &        ! bulk Richardson no. from level to ref lev
                psomc   , &
                obklen  , &
                phiminv , &
                phihinv , &
                tstar   , &
                wstar     )

    !
    !     Wgh1=0.333!1.0_r8/3.0_r8 ! pbl Hostlag Boville
    !     Wgh2=0.333!1.0_r8/3.0_r8 ! pbl Mellor Yamada 2.0
    !     Wgh3=0.333!1.0_r8/3.0_r8 ! pbl Mellor Yamada 2.5
    ! 
    !          where kvm(plond,plev + 1) and kvh(plond,plev + 1) are the 
    !     diffusion coefficients for momentum 
    !          and heat, respectively, l is the master turbulence length scale, 
    !          q2 is the turbulent kinetic energy (so q is the magnitude of 
    !          the turbulent wind velocity), and SM and SH are momentum 
    !          flux and heat flux stability parameters, respectively    
    !         
    IF (kmean == 1) THEN
       IF (plev >= 4) THEN
              DO k = 2, plev  ! plev + 1
                 DO i = 1, plon
                    !
                    !              k=2  ****Km(k),sl*** } -----------
                    !              k=3/2----si,ric,rf,km,kh,b,l -----------
                    !              k=1  ****Km(k),sl*** } -----------
                    !              k=1/2----si ----------------------------
                    !
                    !        Km(k-1) + 2*Km(k) + Km(k+1)
                    ! Km = -------------------------------
                    !                   4
                    !
                    kvm(i,k)=(0.25_r8*(kvm(i,k-1)+2.0_r8*kvm(i,k)+kvm(i,k+1)))
                    IF(ensamble) THEN
                       IF(k<=plev-1) kvm(i,k)=( kvm(i,k)*Wgh1 + KmMixl(i,k)*Wgh2 + AKM(i,k)*Wgh3 )
                    ELSE
                       IF(k<=plev-1) kvm(i,k)=MAX( kvm(i,k), KmMixl(i,k), AKM(i,k))
                    END IF
                    !AKH
                    !        Kh(k-1) + 2*Kh(k) + Kh(k+1)
                    ! Kh = -------------------------------
                    !                   4
                    !
                    kvh(i,k)=(0.25_r8*(kvh(i,k-1)+2.0_r8*kvh(i,k)+kvh(i,k+1)))
                    IF(ensamble) THEN
                       IF(k<=plev-1) kvh(i,k)=( kvh(i,k)*Wgh1 + KhMixl(i,k)*Wgh2 + AKH(i,k)*Wgh3 )
                    ELSE
                       IF(k<=plev-1) kvh(i,k)=MAX( kvh(i,k),KhMixl(i,k), AKH(i,k))
                    END IF
                 END DO
              END DO
       END IF
       DO i = 1, plon
              !
              !          Km(1) + Km(2)
              ! Km = ---------------
              !               2
              !
              !
              !          Kh(1) + Kh(2)
              ! Kh = ---------------
              !               2
              !
              IF(ensamble) THEN
                 kvm(i,     1  )=       ( kvm(i,     1)*Wgh1 + KmMixl(i,     1)*Wgh2 + AKM(i,     1)*Wgh3 )
                 kvh(i,     1  )=       ( kvh(i,     1)*Wgh1 + KhMixl(i,     1)*Wgh2 + AKH(i,     1)*Wgh3 )
              ELSE
                 kvm(i,     1  )=       MAX( kvm(i,     1),KmMixl(i,     1), AKM(i,     1))
                 kvh(i,     1  )=       MAX( kvh(i,     1),KhMixl(i,     1), AKH(i,     1))
              END IF
              !
              !          Km(k-1) + 2*Km(k) + Km(k+1)
              ! Km = -------------------------------
              !                     4
              !
              !kvm(i,plev  )=(0.25_r8*(kvm(i,plev+1)+2.0_r8*kvm(i,plev  )+kvm(i,plev-1)))
              kvm(i,plev+1)=(0.25_r8*(kvm(i,plev-1)+2.0_r8*kvm(i,plev+1)+kvm(i,plev  )))
              !
              !          Kh(k-1) + 2*Kh(k) + Kh(k+1)
              ! Kh = -------------------------------
              !                     4
              !
              !kvh(i,plev  )=(0.25_r8*(kvh(i,plev+1)+2.0_r8*kvh(i,plev  )+kvh(i,plev-1)))
              kvh(i,plev+1)=(0.25_r8*(kvh(i,plev-1)+2.0_r8*kvh(i,plev+1)+kvh(i,plev  )))
       END DO

    END IF

    IF(PBLEntrain) THEN
    CALL PBLNASA(&
                 plond                            , &           !INTEGER         , INTENT(IN   ) :: IM
                 plev                             , &           !INTEGER         , INTENT(IN   ) :: LM
                 pcnst                            , &           !INTEGER         , INTENT(IN   ) :: pcnst
                 ztodt                            , &           !REAL(KIND=r8), INTENT(IN   ) :: DT
                 kvm        (1:plond,1:plev + 1)  , &           !REAL(KIND=r8), INTENT(INOUT) :: kvm        (IM,LM + 1) 
                 kvh        (1:plond,1:plev + 1)  , &           !REAL(KIND=r8), INTENT(INOUT) :: kvh        (IM,LM + 1 )
                 psomc      (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: psomc      (IM,LM)        !(a
                 zm         (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: zm        (IM,LM)                  ! u wind input 
                 um1        (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: um1       (IM,LM)                  ! u wind input
                 vm1        (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: vm1       (IM,LM)                  ! v wind input
                 tm1        (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: tm1       (IM,LM)                  ! temperature input
                 qm1        (1:plond,1:plev,1:pcnst), &         !REAL(KIND=r8), INTENT(IN   ) :: qm1       (IM,LM,pcnst)  ! moisture and trace constituent input
                 pmidm1     (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: pmidm1    (IM,LM)               ! midpoint pressures
                 pintm1     (1:plond,1:plev + 1)  , &           !REAL(KIND=r8), INTENT(IN   ) :: pintm1    (IM,LM + 1)    ! interface pressures
                 LwCoolRateC(1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: LwCoolRateC     (IM,LM) !clearsky_air_temperature_tendency_lw K s-1
                 LwCoolRate (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: LwCoolRate      (IM,LM)!                 air_temperature_tendency_due_to_longwave K s-1
                 cldtot     (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: cldtot    (IM,LM)             
                 qliq       (1:plond,1:plev)      , &           !REAL(KIND=r8), INTENT(IN   ) :: qliq      (IM,LM)             
                 bstar      (1:plond)             , &           !REAL(KIND=r8), INTENT(IN   ) :: bstar     (IM)   !surface_bouyancy_scale m s-2
                 ustar      (1:plond)             , &           !REAL(KIND=r8), INTENT(IN   ) :: USTAR      (IM) !surface_velocity_scale m s-1
                 FRLAND     (1:plond)                )          !REAL(KIND=r8), INTENT(IN   ) :: FRLAND     (IM)
        END IF         
    !
    ! Add the counter grad terms to potential temp, specific humidity
    ! and other constituents in the bdry layer. Note, ktopblmn gives the 
    ! minimum vertical index of the first midpoint within the boundary layer.
    !
    ! first set values above boundary layer
    !
    DO k=1,ktopblmn-2
       DO i=1,plon
          thx(i,k)   = thm(i,k)
          qmx(i,k,1) = qm1(i,k,1)
       END DO
       DO m=2,pcnst
          DO i=1,plon
             qmx(i,k,m) = qm1(i,k,m)
          END DO
       END DO
    END DO
    DO k=2,plev
       DO i=1,plon
          potbar(i,k) = pintm1(i,k)/(0.5_r8*(tm1(i,k) + tm1(i,k-1)))
       END DO
    END DO
    DO i=1,plon
       potbar(i,plev + 1) = pintm1(i,plev + 1)/tm1(i,plev)
    END DO
    !
    ! now focus on the boundary layer
    !
    ztodtgor = ztodt*gravit/rair
    DO k=ktopblmn-1,plev
       DO i=1,plon
          tmp1(i   ) = ztodtgor*rpdel(i,k)
          thx(i,k  ) = thm(i,k  ) + tmp1(i) * (potbar(i,k+1)*kvh(i,k+1)*cgh(i,k+1  ) - potbar(i,k  )*kvh(i,k  )*cgh(i,k    ))
          qmx(i,k,1) = qm1(i,k,1) + tmp1(i) * (potbar(i,k+1)*kvh(i,k+1)*cgq(i,k+1,1) - potbar(i,k  )*kvh(i,k  )*cgq(i,k  ,1))
       END DO
       DO m=2,pcnst
          DO i=1,plon
             qmx(i,k,m) = qm1(i,k,m) + tmp1(i) * (potbar(i,k+1)*kvh(i,k+1)*cgq(i,k+1,m) - potbar(i,k  )*kvh(i,k  )*cgq(i,k  ,m))
          END DO
       END DO
    END DO
    !
    ! Check for neg q's in each constituent and put the original vertical
    ! profile back if a neg value is found. A neg value implies that the
    ! quasi-equilibrium conditions assumed for the countergradient term are
    ! strongly violated.
    ! Original code rewritten by Rosinski 7/8/91 to vectorize in longitude.
    !
    DO m=1,pcnst
       DO i=1,plon
          ilogic(i) = 0
       END DO
       DO k=ktopblmn-1,plev
          DO i=1,plon
             IF (qmx(i,k,m).LT.qmincg(m)) ilogic(i) = 1
          END DO
       END DO
       !
       ! Find long indices of those columns for which negatives were found
       !
       CALL wheneq(plon  , &! INTENT(in)
                   ilogic, &! INTENT(in)
                   1     , &! INTENT(in)
                   1     , &! INTENT(in)
                   indx  , &! INTENT(out)
                   nval    )! INTENT(out)
       !
       ! Replace those columns with original values
       !
       IF (nval.GT.0) THEN
          DO k=ktopblmn-1,plev
             DO ii=1,nval
                i=indx(ii)
                qmx(i,k,m) = qm1(i,k,m)
             END DO
          END DO
       END IF
    END DO
    !
    ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
    ! tridiagonal diffusion matrix. the diagonal elements are a combination of 
    ! ca and cc; they are not explicitly provided to the solver
    !
    gorsq = (gravit/rair)**2
    DO k=ntopfl,plev-1
       DO i=1,plon
          tmp2 = ztodt*gorsq*rpdeli(i,k)*(potbar(i,k+1)**2)
          cah(i,k  ) = kvh(i,k+1)*tmp2*rpdel(i,k  )
          cam(i,k  ) = kvm(i,k+1)*tmp2*rpdel(i,k  )
          cch(i,k+1) = kvh(i,k+1)*tmp2*rpdel(i,k+1)
          ccm(i,k+1) = kvm(i,k+1)*tmp2*rpdel(i,k+1)
       END DO
    END DO
    !
    ! The last element of the upper diagonal is zero.
    !
    DO i=1,plon
       cah(i,plev) = 0.0_r8
       cam(i,plev) = 0.0_r8
    END DO
    !
    ! Calculate e(k) for heat & momentum vertical diffusion.  This term is 
    ! required in solution of tridiagonal matrix defined by implicit diffusion eqn.
    !
    DO i=1,plon
       termh(i,ntopfl) = 1.0_r8/(1.0_r8 + cah(i,ntopfl))
       termm(i,ntopfl) = 1.0_r8/(1.0_r8 + cam(i,ntopfl))
       zeh(i,ntopfl) = cah(i,ntopfl)*termh(i,ntopfl)
       zem(i,ntopfl) = cam(i,ntopfl)*termm(i,ntopfl)
    END DO
    DO k=ntopfl+1,plev-1
       DO i=1,plon
          termh(i,k) = 1.0_r8/ &
               (1.0_r8 + cah(i,k) + cch(i,k) - cch(i,k)*zeh(i,k-1))
          termm(i,k) = 1.0_r8/ &
               (1.0_r8 + cam(i,k) + ccm(i,k) - ccm(i,k)*zem(i,k-1))
          zeh(i,k) = cah(i,k)*termh(i,k)
          zem(i,k) = cam(i,k)*termm(i,k)
       END DO
    END DO
    !
    ! Diffuse momentum
    !

    CALL mvdiff(um1     , & ! INTENT(IN  ) :: um1(plond,plev)             ! u wind input
                vm1     , & ! INTENT(IN  ) :: vm1(plond,plev)             ! v wind input
                dubot   , & ! INTENT(IN  ) :: uflx(plond)    ! sfc u flux into lowest model level
                dvbot   , & ! INTENT(IN  ) :: vflx(plond)    ! sfc v flux into lowest model level
                sufac   , & ! INTENT(IN  ) :: sufac(plond)   ! lowest layer u implicit stress factor
                svfac   , & ! INTENT(IN  ) :: svfac(plond)   ! lowest layer v implicit stress factor
                ccm     , & ! INTENT(IN  ) :: cc(plond,plev) ! -lower diag coeff. of tri-diag matrix
                zem     , & ! INTENT(IN  ) :: ze(plond,plev) ! term in tri-diag. matrix system
                termm   , & ! INTENT(IN  ) :: term(plond,plev)  ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                up1     , & ! INTENT(OUT ) :: up1(plond,plev)        ! u wind after diffusion
                vp1     , & ! INTENT(OUT ) :: vp1(plond,plev)        ! v wind after diffusion
                plond   , & ! INTENT(IN  ) :: plond         ! slt extended domain longitude
                plev    , & ! INTENT(IN  ) :: plev    ! number of vertical levels
                plon      ) ! INTENT(IN  ) :: plon    ! number of longitudes
    !+
    ! Determine the difference between the implicit stress and the
    ! externally specified stress. Apply over boundary layer depth.
    !-
   
    DO i = 1, plon
       IF (ksx(i) .GT. 0.0_r8) sufac(i) = (taux(i) + ksx(i)*up1(i,plev)) &
            * ztodt * gravit / (pintm1(i,plev + 1) - pintm1(i,ktopbl(i)))
       IF (ksy(i) .GT. 0.0_r8) svfac(i) = (tauy(i) + ksy(i)*vp1(i,plev)) &
            * ztodt * gravit / (pintm1(i,plev + 1) - pintm1(i,ktopbl(i)))
       !$$$         print 1, i, lat
       !$$$     $        , taux(i), -ksx(i)*up1(i,plev)
       !$$$     $        , sufac(i)*86400.0_r8/ztodt, ksx(i), um1(i,plev),up1(i,plev)
       !$$$    1 format (1x,i4,i4, 1P, 4E20.10_r8, 2f10.4)
    END DO
    DO k = plev, ktopblmn, -1
       DO i = 1, plon
          IF (k .GE. ktopbl(i)) THEN
             up1(i,k) = up1(i,k) + sufac(i)
             vp1(i,k) = vp1(i,k) + svfac(i)
          END IF
       END DO
    END DO
    !
    ! Diffuse constituents
    !

    CALL qvdiff(pcnst   , &! INTENT(IN   ) :: ncnst                     ! number of constituents being diffused
                qmx     , &! INTENT(IN   ) :: qm1   (plond,plev,ncnst)  ! initial constituent
                dqbot   , &! INTENT(IN   ) :: qflx  (plond,ncnst)      ! sfc q flux into lowest model level
                cch     , &! INTENT(IN   ) :: cc    (plond,plev)         ! -lower diag coeff.of tri-diag matrix
                zeh     , &! INTENT(IN   ) :: ze    (plond,plev)         ! term in tri-diag. matrix system
                termh   , &! INTENT(IN   ) :: term  (plond,plev)       ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                qp1     , &! INTENT(OUT  ) :: qp1   (plond,plev,ncnst)  ! final constituent
                plon    , &! INTENT(IN   ) :: plon   ! number of longitudes
                plond   , &! INTENT(IN   ) :: plond  ! slt extended domain longitude
                plev    , &! INTENT(IN   ) :: plev   ! number of vertical levels
                pcnst   , &! INTENT(IN   ) :: pcnst  ! number of constituents (including water vapor)
                ntopfl   ) ! INTENT(IN   ) :: ntopfl ! Top level to which vertical diffusion is applied.
    !
    ! Identify and correct constituents exceeding user defined bounds
    ! 
    CALL qneg3(plond     , & ! INTENT(IN   ) :: plond           ! slt extended domain longitude
               plev      , & ! INTENT(IN   ) :: plev           ! number of vertical levels
               pcnst     , & ! INTENT(IN   ) :: pcnst           ! number of constituents (including water vapor)
               plon      , & ! INTENT(IN   ) :: plon           ! number of longitudes
               'VDIFF   ', & ! INTENT(IN   ) :: subnam       ! name of calling routine
               qp1(1,1,1)  ) ! INTENT(INOUT) :: q(plond,plev,pcnst) ! moisture/tracer field
    !
    ! Diffuse potential temperature
    !
    CALL qvdiff(1       , & ! INTENT(IN   ) :: ncnst                      ! number of constituents being diffused
                thx     , & ! INTENT(IN   ) :: qm1   (plond,plev,ncnst)  ! initial constituent
                dtbot   , & ! INTENT(IN   ) :: qflx  (plond,ncnst)        ! sfc q flux into lowest model level
                cch     , & ! INTENT(IN   ) :: cc    (plond,plev)          ! -lower diag coeff.of tri-diag matrix
                zeh     , & ! INTENT(IN   ) :: ze    (plond,plev)          ! term in tri-diag. matrix system
                termh   , & ! INTENT(IN   ) :: term  (plond,plev)        ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                thp     , & ! INTENT(OUT  ) :: qp1   (plond,plev,ncnst)  ! final constituent
                plon    , & ! INTENT(IN   ) :: plon   ! number of longitudes
                plond   , & ! INTENT(IN   ) :: plond  ! slt extended domain longitude
                plev    , & ! INTENT(IN   ) :: plev   ! number of vertical levels
                pcnst   , & ! INTENT(IN   ) :: pcnst  ! number of constituents (including water vapor)
                ntopfl   )  ! INTENT(IN   ) :: ntopfl ! Top level to which vertical diffusion is applied.
    !
    RETURN
  END SUBROUTINE vdiff



  !===============================================================================

  SUBROUTINE wheneq (n, array, inc, TARGET, index, nval)
    !----------------------------------------------------------------------- 
    ! 
    ! Purpose: Determine indices of "array" which equal "target"
    ! 
    !-----------------------------------------------------------------------

    IMPLICIT NONE
    !
    ! Arguments
    !
    INTEGER, INTENT(in) :: n

    INTEGER, INTENT(in) :: array(n)    ! array to be searched
    INTEGER, INTENT(in) :: TARGET      ! value to compare against
    INTEGER, INTENT(in) :: inc         ! increment to move through array

    INTEGER, INTENT(out) :: nval       ! number of values meeting criteria
    INTEGER, INTENT(out) :: INDEX(n)   ! output array of indices
    !
    ! Local workspace
    !
    INTEGER :: i
    INTEGER :: ina

    ina=1
    nval=0
    INDEX=0
    IF (inc .LT. 0) ina=(-inc)*(n-1)+1
    DO i=1,n
       IF(array(ina) .EQ. TARGET) THEN
          nval=nval+1
          INDEX(nval)=i
       END IF
       ina=ina+inc
    ENDDO
    RETURN
  END SUBROUTINE wheneq

  SUBROUTINE virtem(t       ,q       ,zvir    ,tv , plond  ,plev ,plon )
    !-----------------------------------------------------------------------
    !
    ! Compute the virtual temperature.
    !
    !---------------------------Code history--------------------------------
    !
    ! Original version:  B. Boville
    ! Standardized:      J. Rosinski, June 1992
    ! Reviewed:          D. Williamson, J. Hack, August 1992
    ! Reviewed:          D. Williamson, March 1996
    !
    !-----------------------------------------------------------------------
    !
    ! $Id: virtem.F,v 1.1.1.1 2001/03/09 00:29:29 mirin Exp $
    !
    !-----------------------------------------------------------------------
    !-----------------------------------------------------------------------
    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) :: plond      ! slt extended domain longitude
    INTEGER, INTENT(IN   ) ::  plev       ! number of vertical levels
    INTEGER, INTENT(IN   ) ::  plon       ! number of longitudes

    REAL(KIND=r8), INTENT(IN   ) ::  t(plond,plev)       ! temperature
    REAL(KIND=r8), INTENT(IN   ) ::  q(plond,plev)       ! specific humidity
    REAL(KIND=r8), INTENT(IN   ) ::  zvir                ! virtual temperature constant
    !
    ! Output arguments
    !
    REAL(KIND=r8), INTENT(OUT  ) ::  tv(plond,plev)      ! virtual temperature
    !
    !---------------------------Local storage-------------------------------
    !
    INTEGER i,k              ! longitude and level indexes
    !
    tv=0.0_r8
    DO k=1,plev
       DO i=1,plon
          tv(i,k) = t(i,k)*(1.0_r8 + zvir*q(i,k))
       END DO
    END DO
    !
    RETURN
  END SUBROUTINE virtem



  SUBROUTINE mvdiff(um1     , &! INTENT(IN  ) :: um1(plond,plev)        ! u wind input
                    vm1     , &! INTENT(IN  ) :: vm1(plond,plev)        ! v wind input
                    uflx    , &! INTENT(IN  ) :: uflx(plond)        ! sfc u flux into lowest model level
                    vflx    , &! INTENT(IN  ) :: vflx(plond)        ! sfc v flux into lowest model level
                    sufac   , &! INTENT(IN  ) :: sufac(plond)        ! lowest layer u implicit stress factor
                    svfac   , &! INTENT(IN  ) :: svfac(plond)        ! lowest layer v implicit stress factor
                    cc      , &! INTENT(IN  ) :: cc(plond,plev)        ! -lower diag coeff. of tri-diag matrix
                    ze      , &! INTENT(IN  ) :: ze(plond,plev)        ! term in tri-diag. matrix system
                    term    , &! INTENT(IN  ) :: term(plond,plev)  ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                    up1     , &! INTENT(OUT ) :: up1(plond,plev)   ! u wind after diffusion
                    vp1     , &! INTENT(OUT ) :: vp1(plond,plev)   ! v wind after diffusion
                    plond   , &! INTENT(IN  ) :: plond      ! slt extended domain longitude
                    plev    , &! INTENT(IN  ) :: plev         ! number of vertical levels
                    plon      )! INTENT(IN  ) :: plon         ! number of longitudes
    !-----------------------------------------------------------------------
    !
    ! Vertical momentum diffusion with explicit surface flux.
    ! Solve the vertical diffusion equation for momentum.
    ! Procedure for solution of the implicit equation follows
    ! Richtmyer and Morton (1967,pp 198-199)
    !
    !---------------------------Code history--------------------------------
    !
    ! Original version:  CCM1
    ! Standardized:      J. Rosinski, June 1992
    ! Reviewed:          P. Rasch, B. Boville, August 1992
    ! Reviewed:          P. Rasch, April 1996
    ! Reviewed:          B. Boville, May 1996
    !
    !-----------------------------------------------------------------------
    !
    ! $Id: mvdiff.F,v 1.1.1.1 2001/03/09 00:29:27 mirin Exp $
    ! $Author: mirin $
    !
    !-----------------------------------------------------------------------
    !      use precision
    !      use pmgrid
    !-----------------------------------------------------------------------
    !!#include <implicit.h>
    !------------------------------Commons----------------------------------
    !!#include <comvd.h>
    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) :: plond      ! slt extended domain longitude
    INTEGER, INTENT(IN   ) :: plev         ! number of vertical levels
    INTEGER, INTENT(IN   ) :: plon         ! number of longitudes

    REAL(KIND=r8), INTENT(IN  ) :: um1(plond,plev)        ! u wind input
    REAL(KIND=r8), INTENT(IN  ) :: vm1(plond,plev)        ! v wind input
    REAL(KIND=r8), INTENT(IN  ) :: uflx(plond)        ! sfc u flux into lowest model level
    REAL(KIND=r8), INTENT(IN  ) :: vflx(plond)        ! sfc v flux into lowest model level
    REAL(KIND=r8), INTENT(IN  ) :: sufac(plond)        ! lowest layer u implicit stress factor
    REAL(KIND=r8), INTENT(IN  ) :: svfac(plond)        ! lowest layer v implicit stress factor
    REAL(KIND=r8), INTENT(IN  ) :: cc(plond,plev)        ! -lower diag coeff. of tri-diag matrix
    REAL(KIND=r8), INTENT(IN  ) :: ze(plond,plev)        ! term in tri-diag. matrix system
    REAL(KIND=r8), INTENT(IN  ) :: term(plond,plev)  ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
    !
    ! Output arguments
    !
    REAL(KIND=r8), INTENT(OUT ) :: up1(plond,plev)   ! u wind after diffusion
    REAL(KIND=r8), INTENT(OUT ) :: vp1(plond,plev)   ! v wind after diffusion
    !
    !---------------------------Local workspace-----------------------------
    !
    REAL(KIND=r8) zfu(plond,plev)   ! terms appearing in soln of tridiag system
    REAL(KIND=r8) zfv(plond,plev)   ! terms appearing in soln of tridiag system
    INTEGER i,k            ! longitude,vertical indices
    !
    !-----------------------------------------------------------------------
    !
    ! Calc fu(k) and fv(k). These are terms required in solution of 
    ! tridiagonal matrix defined by implicit diffusion eqn.  Note that only 
    ! levels ntopfl through plev need be solved for. No diffusion is 
    ! applied above this level.
    !
    up1=0.0_r8
    vp1=0.0_r8
    zfu=0.0_r8
    zfv=0.0_r8
    DO i=1,plon
       zfu(i,ntopfl) = um1(i,ntopfl)*term(i,ntopfl)
       zfv(i,ntopfl) = vm1(i,ntopfl)*term(i,ntopfl)
    END DO
    DO k=ntopfl+1,plev-1
       DO i=1,plon
          zfu(i,k) = (um1(i,k) + cc(i,k)*zfu(i,k-1))*term(i,k)
          zfv(i,k) = (vm1(i,k) + cc(i,k)*zfv(i,k-1))*term(i,k)
       END DO
    END DO
    !
    ! Bottom level: (includes  surface fluxes as either an explicit RHS or
    ! as an implicit stress)
    !
    DO i=1,plon
       zfu(i,plev) = (um1(i,plev) + uflx(i) + cc(i,plev)*zfu(i,plev-1)) &
            / (1.0_r8 + cc(i,plev) + sufac(i) - cc(i,plev)*ze(i,plev-1))
       zfv(i,plev) = (vm1(i,plev) + vflx(i) + cc(i,plev)*zfv(i,plev-1)) &
            / (1.0_r8 + cc(i,plev) + svfac(i) - cc(i,plev)*ze(i,plev-1))
    END DO
    !
    ! Perform back substitution
    !
    DO i=1,plon
       up1(i,plev) = zfu(i,plev)
       vp1(i,plev) = zfv(i,plev)
    END DO
    DO k=plev-1,ntopfl,-1
       DO i=1,plon
          up1(i,k) = zfu(i,k) + ze(i,k)*up1(i,k+1)
          vp1(i,k) = zfv(i,k) + ze(i,k)*vp1(i,k+1)
       END DO
    END DO
    RETURN
  END SUBROUTINE mvdiff


  SUBROUTINE qneg3(plond  , &! INTENT(IN   ) :: plond      ! slt extended domain longitude
                   plev   , &! INTENT(IN   ) :: plev       ! number of vertical levels
                   pcnst  , &! INTENT(IN   ) :: pcnst      ! number of constituents (including water vapor)
                   plon   , &! INTENT(IN   ) :: plon       ! number of longitudes
                   subnam , &! INTENT(IN   ) :: subnam       ! name of calling routine
                   q       ) ! INTENT(INOUT) :: q(plond,plev,pcnst) ! moisture/tracer field
    !-----------------------------------------------------------------------
    !
    ! Check moisture and tracers for minimum value, reset any below
    ! minimum value to minimum value and return information to allow
    ! warning message to be printed. The global average is NOT preserved.
    !
    !---------------------------Code history--------------------------------
    !
    ! Original version:  J. Rosinski
    ! Standardized:      J. Rosinski, June 1992
    ! Reviewed:          D. Williamson, August 1992, March 1996
    !
    !-----------------------------------------------------------------------
    !
    ! $Id: qneg3.F,v 1.1.1.1 2001/03/09 00:29:27 mirin Exp $
    !
    !-----------------------------------------------------------------------
    !      use precision
    !      use pmgrid
    !-----------------------------------------------------------------------
    !#include <implicit.h>
    !------------------------------Commons----------------------------------
    !#include <comqmin.h>
    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) :: plon       ! number of longitudes
    INTEGER, INTENT(IN   ) :: plond      ! slt extended domain longitude
    INTEGER, INTENT(IN   ) :: plev       ! number of vertical levels
    INTEGER, INTENT(IN   ) :: pcnst      ! number of constituents (including water vapor)

    CHARACTER*8, INTENT(IN   ) :: subnam       ! name of calling routine
    !
    ! Input/Output arguments
    !
    REAL(KIND=r8), INTENT(INOUT) :: q(plond,plev,pcnst) ! moisture/tracer field
    !
    !---------------------------Local workspace-----------------------------
    !
    INTEGER indx(plond)      ! array of indices of points < qmin
    INTEGER nval             ! number of points < qmin for 1 level
    INTEGER nvals            ! number of values found < qmin
    INTEGER i,ii,k           ! longitude, level indices
    INTEGER m                ! constituent index
    INTEGER iw,kw            ! i,k indices of worst violator

    LOGICAL found            ! true => at least 1 minimum violator found

    REAL(KIND=r8) worst               ! biggest violator
    !
    !-----------------------------------------------------------------------
    !
    DO m=1,pcnst
       DO k=1,plev
          DO i=1,plon
             q(i,k,m)=max(q(i,k,m),qmin(m))
          END DO   
       END DO          
    END DO
    RETURN
    
    DO m=1,pcnst
       nvals = 0
       found = .FALSE.
       worst = 0.0_r8
       !
       ! Test all field values for being less than minimum value. Set q = qmin
       ! for all such points. Trace offenders and identify worst one.
       !
       DO k=1,plev
          nval = 0
          DO i=1,plon
             IF (q(i,k,m) < qmin(m)) THEN
                nval = nval + 1
                indx(nval) = i
             END IF
          END DO

          IF (nval.GT.0) THEN
             found = .TRUE.
             nvals = nvals + nval
             DO ii=1,nval
                i = indx(ii)
                IF (q(i,k,m).LT.worst) THEN
                   !worst = q(i,k,m)
                   kw = k
                   iw = i
                END IF
                q(i,k,m) = qmin(m)
             END DO
          END IF
       END DO
       IF (found) THEN
          WRITE(6,9000)subnam,m,nvals,qmin(m),worst,iw,kw
       END IF
    END DO
    !
    RETURN
9000 FORMAT(' QNEG3 from ',a8,':m=',i3, &
         ' Min. mixing ratio violated at ',i4,' points.  Reset to ', &
         1p,e8.1,' Worst =',e8.1,' at i,k=',i4,i3)
  END SUBROUTINE qneg3



  SUBROUTINE pbldif(th      , &! INTENT(IN   ) ::  th(plond,plev)          ! potential temperature [K]
                    q       , &! INTENT(IN   ) ::  q(plond,plev,pcnst)     ! specific humidity [kg/kg]
                    z       , &! INTENT(IN   ) ::  z(plond,plev)           ! height above surface [m]
                    u       , &! INTENT(IN   ) ::  u(plond,plev)           ! windspeed x-direction [m/s]
                    v       , &! INTENT(IN   ) ::  v(plond,plev)           ! windspeed y-direction [m/s]
                    t       , &! INTENT(IN   ) ::  t(plond,plev)           ! temperature (used for density)
                    pmid    , &! INTENT(IN   ) ::  pmid(plond,plev)        ! midpoint pressures
                    kvf     , &! INTENT(IN   ) ::  kvf(plond,plev + 1)     ! free atmospheric eddy diffsvty [m2/s]
                    cflx    , &! INTENT(IN   ) ::  cflx(plond,pcnst)           ! surface constituent flux (kg/m2/s)
                    shflx   , &! INTENT(IN   ) ::  shflx(plond)            ! surface heat flux (W/m2)
                    taux    , &! INTENT(IN   ) ::  taux(plond)             ! surface u stress (N)
                    tauy    , &! INTENT(IN   ) ::  tauy(plond)             ! surface v stress (N)
                    ustar   , &! INTENT(OUT  ) ::  ustar(plond)            ! surface friction velocity [m/s]
                    kvm     , &! INTENT(OUT  ) ::  kvm(plond,plev + 1)     ! eddy diffusivity for momentum [m2/s]
                    kvh     , &! INTENT(OUT  ) ::  kvh(plond,plev + 1)     ! eddy diffusivity for heat [m2/s]
                    cgh     , &! INTENT(OUT  ) ::  cgh(plond,plev + 1)     ! counter-gradient term for heat [K/m]
                    cgq     , &! INTENT(OUT  ) ::  cgq(plond,plev + 1,pcnst)! counter-gradient term for constituents
                    cgs     , &! INTENT(OUT  ) ::  cgs(plond,plev + 1)     ! counter-gradient star (cg/flux)
                    pblh    , &! INTENT(OUT  ) ::  pblh(plond)             ! boundary-layer height [m]
                    tpert   , &! INTENT(OUT  ) ::  tpert(plond)            ! convective temperature excess
                    qpert   , &! INTENT(OUT  ) ::  qpert(plond)            ! convective humidity excess
                    ktopbl  , &! INTENT(OUT  ) ::  ktopbl(plond)           ! index of first midpoint inside pbl
                    ktopblmn, &! INTENT(OUT  ) ::  ktopblmn                ! min value of ktopbl
                    plond   , &! INTENT(IN   ) :: plond      ! slt extended domain longitude
                    pcnst   , &! INTENT(IN   ) :: pcnst      ! number of constituents (including water vapor)
                    plon    , &! INTENT(IN   ) :: plon       ! number of longitudes
                    plev    , &! INTENT(IN   ) ::  plev       ! number of vertical levels
                    fakn    , &! INTENT(IN   ) :: fakn    ! Constant in turbulent prandtl number
                    ricr    , & ! INTENT(IN   ) :: ricr    ! Critical richardson number
                    tsk     , &
                    qsfc    , &
                    rino    , &
                    psomc   , &
                    obklen  , &
                    phiminv , &
                    phihinv , &
                    tstar   , &
                    wstar     )
    !------------------------------------------------------------------------
    ! 
    ! Atmospheric boundary layer computation.
    !
    ! Nonlocal scheme that determines eddy diffusivities based on a
    ! diagnosed boundary layer height and a turbulent velocity scale;
    ! also, countergradient effects for heat and moisture, and constituents
    ! are included, along with temperature and humidity perturbations which 
    ! measure the strength of convective thermals in the lower part of the 
    ! atmospheric boundary layer.
    !
    ! For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
    ! Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate
    ! Model. J. Clim., vol. 6., p. 1825--1842.
    !
    ! Updated by Holtslag and Hack to exclude the surface layer from the
    ! definition of the boundary layer Richardson number. Ri is now defined
    ! across the outer layer of the pbl (between the top of the surface
    ! layer and the pbl top) instead of the full pbl (between the surface and
    ! the pbl top). For simiplicity, the surface layer is assumed to be the
    ! region below the first model level (otherwise the boundary layer depth 
    ! determination would require iteration).
    !
    !---------------------------Code history--------------------------------
    !
    ! Original version:  B. Boville
    ! Standardized:      J. Rosinski, June 1992
    ! Reviewed:          B. Boville, P. Rasch, August 1992
    ! Reviewed:          B. Boville, P. Rasch, April 1996
    !
    ! Modified for boundary layer height diagnosis: Bert Holtslag, june 1994
    ! >>>>>>>>>  (Use ricr = 0.3 in this formulation)
    !
    !-----------------------------------------------------------------------
    !
    ! $Id: pbldif.F,v 1.1.1.1 2001/03/09 00:29:27 mirin Exp $
    !
    !-----------------------------------------------------------------------
    !-----------------------------------------------------------------------
    !------------------------------Arguments--------------------------------
    !
    ! Input arguments
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) :: plond      ! slt extended domain longitude
    INTEGER, INTENT(IN   ) :: pcnst      ! number of constituents (including water vapor)
    INTEGER, INTENT(IN   ) :: plon       ! number of longitudes
    INTEGER, INTENT(IN   ) ::  plev       ! number of vertical levels
    REAL(KIND=r8), INTENT(IN   ) :: fakn    ! Constant in turbulent prandtl number
    REAL(KIND=r8), INTENT(IN   ) :: ricr    ! Critical richardson number
    REAL(KIND=r8), INTENT(IN   ) ::  th(plond,plev)          ! potential temperature [K]
    REAL(KIND=r8), INTENT(IN   ) ::  q(plond,plev,pcnst)     ! specific humidity [kg/kg]
    REAL(KIND=r8), INTENT(IN   ) ::  z(plond,plev)           ! height above surface [m]
    REAL(KIND=r8), INTENT(IN   ) ::  u(plond,plev)           ! windspeed x-direction [m/s]
    REAL(KIND=r8), INTENT(IN   ) ::  v(plond,plev)           ! windspeed y-direction [m/s]
    REAL(KIND=r8), INTENT(IN   ) ::  t(plond,plev)           ! temperature (used for density)
    REAL(KIND=r8), INTENT(IN   ) ::  pmid(plond,plev)        ! midpoint pressures
    REAL(KIND=r8), INTENT(IN   ) ::  kvf(plond,plev + 1)        ! free atmospheric eddy diffsvty [m2/s]
    REAL(KIND=r8), INTENT(IN   ) ::  cflx(plond,pcnst)       ! surface constituent flux (kg/m2/s)
    REAL(KIND=r8), INTENT(IN   ) ::  shflx(plond)            ! surface heat flux (W/m2)
    REAL(KIND=r8), INTENT(IN   ) ::  taux(plond)             ! surface u stress (N)
    REAL(KIND=r8), INTENT(IN   ) ::  tauy(plond)             ! surface v stress (N)
    REAL(KIND=r8), INTENT(IN   ) ::  tsk   (plond)    
    REAL(KIND=r8), INTENT(IN   ) ::  qsfc  (plond)    
    REAL(KIND=r8), INTENT(IN   ) ::  psomc(plond,plev)      ! (psm1/pmidm1)**cappa

    !
    ! Output arguments
    !
    REAL(KIND=r8), INTENT(OUT   ) ::  ustar  (plond)            ! surface friction velocity [m/s]
    REAL(KIND=r8), INTENT(OUT   ) ::  kvm    (plond,plev + 1)        ! eddy diffusivity for momentum [m2/s]
    REAL(KIND=r8), INTENT(OUT   ) ::  kvh    (plond,plev + 1)        ! eddy diffusivity for heat [m2/s]
    REAL(KIND=r8), INTENT(OUT   ) ::  cgh    (plond,plev + 1)        ! counter-gradient term for heat [K/m]
    REAL(KIND=r8), INTENT(OUT   ) ::  cgq    (plond,plev + 1,pcnst)  ! counter-gradient term for constituents
    REAL(KIND=r8), INTENT(OUT   ) ::  cgs    (plond,plev + 1)        ! counter-gradient star (cg/flux)
    REAL(KIND=r8), INTENT(INOUT ) ::  pblh   (plond)             ! boundary-layer height [m]
    REAL(KIND=r8), INTENT(INOUT ) ::  rino   (plond,plev)        ! bulk Richardson no. from level to ref lev

    REAL(KIND=r8), INTENT(OUT   ) ::  tpert  (plond)            ! convective temperature excess
    REAL(KIND=r8), INTENT(OUT   ) ::  qpert  (plond)            ! convective humidity excess

    INTEGER      , INTENT(OUT   ) :: ktopbl  (plond)        ! index of first midpoint inside pbl
    INTEGER      , INTENT(OUT   ) :: ktopblmn             ! min value of ktopbl
    REAL(KIND=r8), INTENT(OUT   ) :: obklen  (plond)           ! Obukhov length
    REAL(KIND=r8), INTENT(OUT   ) :: phiminv (plond)          ! inverse phi function for momentum
    REAL(KIND=r8), INTENT(OUT   ) :: phihinv (plond)          ! inverse phi function for heat 
    REAL(KIND=r8), INTENT(OUT   ) :: tstar   (plond)
    REAL(KIND=r8), INTENT(OUT   ) :: wstar   (plond)
    !
    !---------------------------Local parameters----------------------------
    !
    REAL(KIND=r8) tiny                    ! lower bound for wind magnitude
    PARAMETER (tiny=1.e-36_r8)
    !
    !---------------------------Local workspace-----------------------------
    !
    INTEGER i                    ! longitude index
    INTEGER k                    ! level index
    INTEGER m                    ! constituent index
    REAL(KIND=r8) heatv  (plond)            ! surface virtual heat flux
    REAL(KIND=r8) thvsrf (plond)           ! sfc (bottom) level virtual temperature
    REAL(KIND=r8) thvref (plond)           ! reference level virtual temperature
    REAL(KIND=r8) tkv                     ! model level potential temperature
    REAL(KIND=r8) therm  (plond)            ! thermal virtual temperature excess
    REAL(KIND=r8) wm     (plond)               ! turbulent velocity scale for momentum
    REAL(KIND=r8) vvk                     ! velocity magnitude squared
    REAL(KIND=r8) zm     (plond)               ! current level height
    REAL(KIND=r8) zp     (plond)               ! current level height + one level up
    REAL(KIND=r8) khfs   (plond)             ! surface kinematic heat flux [mK/s]
    REAL(KIND=r8) kqfs   (plond,pcnst)       ! sfc kinematic constituent flux [m/s]
    REAL(KIND=r8) zmzp                    ! level height halfway between zm and zp
    REAL(KIND=r8) tlv    (plond)              ! ref. level pot tmp + tmp excess
    REAL(KIND=r8) fak1   (plond)             ! k*ustar*pblh
    REAL(KIND=r8) fak2   (plond)             ! k*wm*pblh
    REAL(KIND=r8) fak3   (plond)             ! fakn*wstr/wm 
    REAL(KIND=r8) pblk   (plond)             ! level eddy diffusivity for momentum
    REAL(KIND=r8) pr     (plond)               ! Prandtl number for eddy diffusivities
    REAL(KIND=r8) zl     (plond)               ! zmzp / Obukhov length
    REAL(KIND=r8) zh     (plond)               ! zmzp / pblh
    REAL(KIND=r8) zzh    (plond)              ! (1-(zmzp/pblh))**2
    REAL(KIND=r8) wstr   (plond)             ! w*, convective velocity scale
    REAL(KIND=r8) rrho   (plond)             ! 1./bottom level density (temporary)
    REAL(KIND=r8) ustr                    ! unbounded ustar
    REAL(KIND=r8) term                    ! intermediate calculation
    REAL(KIND=r8) fac                     ! interpolation factor
    REAL(KIND=r8) pblmin                  ! min pbl height due to mechanical mixing

    LOGICAL unstbl (plond)        ! pts w/unstbl pbl (positive virtual ht flx)
    LOGICAL stblev (plond)        ! stable pbl with levels within pbl
    LOGICAL unslev (plond)        ! unstbl pbl with levels within pbl
    LOGICAL unssrf (plond)        ! unstb pbl w/lvls within srf pbl lyr
    LOGICAL unsout (plond)        ! unstb pbl w/lvls in outer pbl lyr
    LOGICAL check  (plond)         ! True=>chk if Richardson no.>critcal
    ustar=0.0_r8;kvm    =0.0_r8;kvh    =0.0_r8;cgh    =0.0_r8  
    cgq    =0.0_r8;cgs    =0.0_r8;tpert=0.0_r8;qpert=0.0_r8  
    obklen  =0.0_r8;phiminv =0.0_r8;    phihinv =0.0_r8  
    heatv =0.0_r8;    thvsrf=0.0_r8;    thvref=0.0_r8;    tkv   =0.0_r8;
    therm =0.0_r8;    wm    =0.0_r8;    vvk   =0.0_r8;    zm    =0.0_r8;
    zp    =0.0_r8;    khfs  =0.0_r8;    kqfs  =0.0_r8;    zmzp  =0.0_r8;
    tlv   =0.0_r8;    fak1  =0.0_r8;    fak2  =0.0_r8;    fak3  =0.0_r8;
    pblk  =0.0_r8;    pr    =0.0_r8;    zl    =0.0_r8;    zh    =0.0_r8;
    zzh   =0.0_r8;    wstr  =0.0_r8;    rrho  =0.0_r8;    ustr  =0.0_r8;
    term  =0.0_r8;    fac   =0.0_r8;    pblmin=0.0_r8;
    !
    ! Compute kinematic surface fluxes
    !
    DO i=1,plon
       !j/kg/kelvin
       !
       ! P = rho * R * T
       !
       !            P
       ! rho  = -------
       !          R * T
       !
       !           1             R * T
       !  rrho = ----- =   -------
       !          rho              P
       !
       !
       rrho(i) = rair*t(i,plev)/pmid(i,plev)

       !
       !
       ! tau = rho * ustar*2
       !
       !
       ustr = SQRT(SQRT(taux(i)**2 + tauy(i)**2)*rrho(i))
       !
       ! surface friction velocity [m/s]
       !
       ustar(i) = MAX(ustr,0.01_r8)
       !
       ! 
       ! surface kinematic heat flux [mK/s] ustar*tstar
       !
       ! H=Ho=-rho*cp*ustar*tstar  
       !
       ! surface heat flux (W/m2)  (j/kg/kelvin)
       !
       !                  H           j       kg K   m*m*m           m K
       ! ustar*tstar = ---------   = ------------------------- =   ------------
       !                -rho*cp      m*m*s   j         kg              s
       !
       khfs(i) = shflx(i)*rrho(i)/cpair
       !
       !
       !
       ! surface kinematic constituent flux [m/s]
       !
       !  J            N *m         Kg * m * m             m*m
       !------ = --------- =------------------- = ---------= hvap = 2.5104e+6! latent heat of vaporization of water (J kg-1)
       !  kg             kg          s*s Kg              s*s
       !
       ! (kg/m2/s)
       !   W            J              N *m         Kg * m * m          s*s              kg 
       ! ------ = ------- = --------- = ------------  * -------  = -------
       !  m*m           m*m*s      m*m*s          s*s*s*m*m          m*m            M*m*s
       !
       !         W           kg        
       !   =   ------ *  ------  =  LFlux / hvap
       !        m*m           J  
       !
       ! E=Eo=-rho*hl*ustar*qstar
       !
       !                  E               kg          m * m * m          m 
       ! ustar*qstar = ------------ = ----------- * ------------- = -----------
       !                -rho *hl       m * m * s       kg                s
       !
       !
       ! surface constituent flux (kg/m2/s)
       !
       kqfs(i,1) = cflx(i,1)*rrho(i)
       !
    END DO
    
    DO m=2,pcnst
       DO i=1,plon
          kqfs(i,m)= cflx(i,m)*rrho(i)
       END DO
    END DO
    !
    ! Initialize output arrays with free atmosphere values
    !
    DO k=1,plev + 1
       DO i=1,plon
          kvm(i,k)   = kvf(i,k)
          kvh(i,k)   = kvf(i,k)
          cgh(i,k)   = 0.0_r8
          cgq(i,k,1) = 0.0_r8
          cgs(i,k)   = 0.0_r8
       END DO
    END DO
    DO m=2,pcnst
       DO k=1,plev + 1
          DO i=1,plon
             cgq(i,k,m) = 0.0_r8
          END DO
       END DO
    END DO
    !
    ! Compute various arrays for use later:
    !
    DO i=1,plon
       !thvsrf(i) = th(i,plev)*(1.0_r8 + 0.61_r8*q(i,plev,1))

       thvsrf(i) = (tsk(i)*psomc(i,plev))*(1.0_r8 + 0.61_r8*qsfc(i))
       !
       !
       !  ustar*tstar = ustar*tstar + ustar*qstar*(0.61*Th)
       !
       !
       !heatv(i)  = khfs(i) + 0.61_r8*th(i,plev)*kqfs(i,1)
       heatv(i)  = khfs(i) *(1.0_r8 + 0.61_r8*qsfc(i))!+ 0.61_r8*(tsk(i)*psomc(i,plev))*kqfs(i,1)
       wm(i)     = 0.0_r8
       therm(i)  = 0.0_r8
       qpert(i)  = 0.0_r8
       tpert(i)  = 0.0_r8
       fak3(i)   = 0.0_r8  
       zh(i)     = 0.0_r8  
       !
       !! Obukhov length
       !
       !
       !                 ustar^3
       ! L = ----------------------------------
       !       -                           -
       !      |         g          Ho       |
       !      | k * -------- * ------------ |
       !      |         To       rho * cp   |
       !       -                           -
       ! L > 0 indicates stable conditions
       ! L < 0 indicates unstable conditions
       ! L --> infinito  applies to neutral conditions
       ! 
       obklen(i) = -thvsrf(i)*ustar(i)**3/ &
            (g*vk*(heatv(i) + SIGN(1.e-10_r8,heatv(i))))
       
    END DO
    !
    ! 
    ! >>>> Define first a new factor fac=100 for use in Richarson number
    !      Calculate virtual potential temperature first level
    !      and initialize pbl height to z1
    !
    fac = 1.0_r8
    !
    DO i=1,plon
       thvref(i) = th(i,plev)*(1.0_r8 + 0.61_r8*q(i,plev,1))
       pblh(i)   = z(i,plev)
       check(i)  = .TRUE.
       !
       ! Initialization of lowest level Ri number 
       ! (neglected in initial Holtslag implementation)
       !
       !      vv2(i) = fac*ustar(i)**2
       !      vv2(i) = max(vv2(i),tiny)
       !      rino(i,plev) = g*(thvsrf(i) - thvref(i))*z(i,plev)/ &
       !                    (thvref(i)*vv2(i))
       rino(i,plev) = 0.0_r8
    END DO
    !
    ! PBL height calculation:
    ! Search for level of pbl. Scan upward until the Richardson number between
    ! the first level and the current level exceeds the "critical" value.
    !
!    DO k=plev-1,plev-npbl+1,-1
    DO k=plev-1,1,-1

       DO i=1,plon
          IF (check(i)) THEN
             ! y = a + b * x
             vvk = (u(i,k) - u(i,plev))**2 + (v(i,k) - v(i,plev))**2  + fac*ustar(i)**2
             !
             !     PARAMETER (tiny=1.e-36_r8)
             !
             vvk = MAX(vvk,tiny)
             tkv = th(i,k)*(1.0_r8 + 0.61_r8*q(i,k,1))
             !
             !           (Tpv(i) - Tpvo) * ( Z (i) - Zo) 
             !  Ri = g ------------------------------------------------
             !          Tpvo * ( (DU(i)^2  +  DV(i))^2  + fac*ustar^2 )
             !
             rino(i,k) = g*(tkv - thvref(i))*(z(i,k)-z(i,plev)) &
                  /(thvref(i)*vvk)                  
             IF(rino(i,k).GE.ricr) THEN
                !               -                -
                !              |                  |
                !         Ricr | u(h)^2 +  v(h)^2 |
                !              |                  |
                !               -                -
                ! h = ----------------------------------------
                !                -    -           -                    -
                !               |   g  |   |                     |
                !               | ---- | * | Thv(h) +  Ths    |
                !               |  Ths |   |                     |
                !               -    -           -                    -
                !
                !
                !
                !                Ri(k+1)       Ric/Ri - 1
                ! h = Z(k+1)   +          -------------------- *(Z(k) - Z(k+1))
                !                Ri(k+1)      Ri(k)/Ri - 1
                !
                !

                !
                !
                !                   Ric - Ri(k+1)
                ! h = Z(k+1)   +  -------------------- *(Z(k) - Z(k+1))
                !                   Ri(k) - Ri(k+1)
                !
                !
                pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/ &
                     (rino(i,k) - rino(i,k+1))*(z(i,k) - z(i,k+1))
                check(i) = .FALSE.
             END IF
          END IF
       END DO
    END DO
    !
    ! Set pbl height to maximum value where computation exceeds number of
    ! layers allowed
    !
    DO i=1,plon
       IF (check(i)) pblh(i) = z(i,plev + 1-npbl)
    END DO
    !
    ! Improve estimate of pbl height for the unstable points.
    ! Find unstable points (virtual heat flux is positive):
    !
    DO i=1,plon
          !
          ! heatv = ustar*tstar 
          !
       IF (heatv(i) .GT. 0.0_r8) THEN
          unstbl(i) = .TRUE.
          check(i) = .TRUE.
       ELSE
          unstbl(i) = .FALSE.
          check(i) = .FALSE.
       END IF
    END DO
    !
    ! For the unstable case, compute velocity scale and the
    ! convective temperature excess:
    !
    DO i=1,plon
       IF (check(i)) THEN
          !
          ! sffrac = 0.1_r8     ! Surface layer fraction of boundary layer common /compbl/
          ! 
          ! betam  = 15.0_r8 ! Constant in wind gradient expression common /compbl/
          !
          ! binm = betam*sffrac
          ! onet = 1.0_r8/3.0_r8
          !
          !              -                                        -  1/3
          !             |                                         |
          !             |                         HPBL                 |  
          !  phiminv =  |1.0_r8 - 0.1 * 15 * ------------------- |
          !             |                           L                 |  
          !             |                                         |
          !              -                                        -
          !
          phiminv(i) = (1.0_r8 - binm*pblh(i)/obklen(i))**onet
          !IF(phiminv(i) > 10.0_r8)WRITE(nfprt,*)'1',phiminv(i), pblh(i),obklen(i)

          !
          !            -                                                      -  1/3
          !           |                                                       |
          !           |                                         HPBL               |  
          ! um =   |ustar^3 - ustar^3 * 0.1 * 15 * ------------------- |
          !           |                                           L           |  
          !           |                                                       |
          !            -                                                      -
          !
          wm(i) = ustar(i)*phiminv(i)
          therm(i) = heatv(i)*fak/wm(i)       
          !--         rino(i,plev) = -g*therm(i)*z(i,plev)/(thvref(i)*vv2(i))
          rino(i,plev) = 0.0_r8
          tlv(i) = thvref(i) + therm(i)
       END IF
    END DO
    !
    ! Improve pblh estimate for unstable conditions using the
    ! convective temperature excess:
    !
    !DO k=plev-1,plev-npbl+1,-1
    DO k=plev-1,1,-1
       DO i=1,plon
          IF (check(i)) THEN
             vvk = (u(i,k) - u(i,plev))**2 + (v(i,k) - v(i,plev))**2 + fac*ustar(i)**2
             vvk = MAX(vvk,tiny)
             tkv = th(i,k)*(1.0_r8 + 0.61_r8*q(i,k,1))
             !
             !
             !             
             ! Ri = g * -----------------------
             !
             !
             !
             rino(i,k) = g*(tkv - tlv(i))*(z(i,k)-z(i,plev)) &
                  /(thvref(i)*vvk)
             IF(rino(i,k).GE.ricr) THEN
                !
                !
                !                   Ric - Ri(k+1)
                ! h = Z(k+1)   +  -------------------- *(Z(k) - Z(k+1))
                !                   Ri(k) - Ri(k+1)
                !
                !
                pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/ &
                     (rino(i,k) - rino(i,k+1))*(z(i,k) - z(i,k+1))
                check(i) = .FALSE.
             END IF
          END IF
       END DO
    END DO
    !
    ! Points for which pblh exceeds number of pbl layers allowed;
    ! set to maximum
    !
    !
    ! PBL height must be greater than some minimum mechanical mixing depth
    ! Several investigators have proposed minimum mechanical mixing depth
    ! relationships as a function of the local friction velocity, u*.  We 
    ! make use of a linear relationship of the form h = c u* where c=700.
    ! The scaling arguments that give rise to this relationship most often 
    ! represent the coefficient c as some constant over the local coriolis
    ! parameter.  Here we make use of the experimental results of Koracin 
    ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
    ! where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
    ! latitude value for f so that c = 0.07/f = 700.
    !
    DO i=1,plon
       pblmin  = 700.0_r8*ustar(i)
       pblh(i) = MAX(pblh(i),pblmin)
    END DO
    !
    ! pblh is now available; do preparation for diffusivity calculation:
    !
    DO i=1,plon
       pblk(i) = 0.0_r8
       fak1(i) = ustar(i)*pblh(i)*vk
       !
       ! Do additional preparation for unstable cases only, set temperature
       ! and moisture perturbations depending on stability.
       !
       IF (unstbl(i)) THEN
          !phiminv(i) = (1.0_r8 - binm*pblh(i)/obklen(i))**onet
          phiminv(i) = (1.0_r8 - binm*pblh(i)/obklen(i))**onet
          !IF(phiminv(i) > 10.00_r8)WRITE(nfprt,*)'2',phiminv(i), pblh(i),obklen(i)
          
          phihinv(i) = SQRT(1.0_r8 - binh*pblh(i)/obklen(i))
          wm(i)      = ustar(i)*phiminv(i)
          fak2(i)    = wm(i)*pblh(i)*vk
          wstr(i)    = (heatv(i)*g*pblh(i)/thvref(i))**onet 
          fak3(i)    = fakn*wstr(i)/wm(i)
          tpert(i)   = MAX(khfs(i)*fak/wm(i),0.0_r8)   
          qpert(i)   = MAX(kqfs(i,1)*fak/wm(i),0.0_r8)    
       ELSE
          tpert(i)   = MAX(khfs(i)*fak/ustar(i),0.0_r8) 
          qpert(i)   = MAX(kqfs(i,1)*fak/ustar(i),0.0_r8) 
       END IF
    END DO
    !
    ! temperature scale
    !
    DO i=1,plon
       !j/kg/kelvin
       !
       ! P = rho * R * T
       !
       !            P
       ! rho  = -------
       !          R * T
       !
       !           1        R * T
       !  rrho = ----- =   -------
       !          rho         P
       !
       !
       rrho(i) = gasr*t(i,plev)/pmid(i,plev)
       !
       !          _       _  2/3             _      _ -1/3
       !         |    Ho    |               | g*h    |
       ! Tstar = | -------- |           *   |--------|          =   K
       !         |   rho*cp |               |  To    |
       !          -       -                  -      -
       !
       tstar(i) = (((MAX(shflx(i),0.0001_r8)*rrho(i)/cp)**(2.0_r8))**(1.0_r8/3.0_r8))   * &
                    (1.0_r8/(((grav*MAX(pblh(i),0.1_r8))/t(i,plev))**(1.0_r8/3.0_r8)))
       !
       !
       !          _                  _ 1/3
       !         |    Ho        g*h   |
       ! wstar = | -------- * --------|   =   K
       !         |  rho*cp      To    |
       !          -                  -
       ! 
       wstar(i) =  ((MAX(shflx(i),0.0001_r8)*rrho(i)/cp) * ((grav*MAX(pblh(i),0.1_r8))/t(i,plev)) )**(1.0_r8/3.0_r8)
    END DO
    !
    ! Main level loop to compute the diffusivities and 
    ! counter-gradient terms:
    !
    ktopbl=1
    !DO 1000 k=plev,plev-npbl+2,-1
    DO 1000 k=plev,2,-1
       !
       ! Find levels within boundary layer:
       !
       DO i=1,plon
          unslev(i) = .FALSE.! unstbl pbl with levels within pbl
          stblev(i) = .FALSE.
          zm(i) = z(i,k)
          zp(i) = z(i,k-1)
          IF (zkmin.EQ.0.0_r8 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
          IF (zm(i) .LT. pblh(i)) THEN
             ktopbl(i) = k
             zmzp = 0.5_r8*(zm(i) + zp(i))
             !
             !        0.5 *( Z(k) - Z(k-1) )
             !  zh = ------------------------
             !                 h
             !
             !
             !             
             zh(i) = zmzp/pblh(i)
             !
             !        0.5 *( Z(k) - Z(k-1) )
             !  zl = ------------------------
             !                 L
             !
             !
             !
             !
             !
             !                 ustar^3
             ! L = ----------------------------------
             !       -                           -
             !      |         g          Ho       |
             !      | k * -------- * ------------ |
             !      |         To       rho * cp   |
             !       -                           -
             !

             zl(i) = zmzp/obklen(i)
             !
             zzh(i) = 0.0_r8
             !              -           - 2
             !             |          z  |
             !zzh =  |1   -   -----|
             !             |          h  |
             !              -           -
             !
             IF (zh(i).LE.1.0_r8) zzh(i) = (1.0_r8 - zh(i))**2
             !
             ! stblev for points zm < plbh and stable and neutral
             ! unslev for points zm < plbh and unstable
             !
             IF (unstbl(i)) THEN
                unslev(i) = .TRUE. ! unstbl pbl with levels within pbl
             ELSE
                stblev(i) = .TRUE.
             END IF
          END IF
       END DO
       !
       ! Stable and neutral points; set diffusivities; counter-gradient
       ! terms zero for stable case:
       !
       DO i=1,plon
          IF (stblev(i)) THEN
             IF (zl(i).LE.1.0_r8) THEN
                !
                !               -           - 2
                !              |          z  |
                !  Kc = k*wi*z*|1   -   -----|
                !              |          h  |
                !               -           -
                !
                !  fak1(i) = ustar(i)*pblh(i)*vk
                !
                pblk(i) = fak1(i)*zh(i)*zzh(i)/(1.0_r8 + betas*zl(i))
             ELSE
                !
                !               -           - 2
                !              |          z  |
                !  Kc = k*wi*z*|1   -   -----|
                !              |          h  |
                !               -           -
                !
                pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
             END IF
             kvm(i,k) = MAX(pblk(i),kvf(i,k))
             kvh(i,k) = kvm(i,k)
          END IF
       END DO
       !
       ! unssrf, unstable within surface layer of pbl
       ! unsout, unstable within outer   layer of pbl
       !
       DO i=1,plon
          unssrf(i) = .FALSE.
          unsout(i) = .FALSE.
          IF (unslev(i)) THEN
             ! unstbl pbl with levels within pbl
             IF (zh(i).LT.sffrac) THEN
                unssrf(i) = .TRUE.
             ELSE
                unsout(i) = .TRUE.
             END IF
          END IF
       END DO
       !
       ! Unstable for surface layer; counter-gradient terms zero
       !
       DO i=1,plon
          IF (unssrf(i)) THEN
             !
             !        0.5 *( Z(k) - Z(k-1) )
             !  zl = ------------------------
             !                 L
             !
             !
             !          -                                       - 1/3
             !         |                 0.5 *( Z(k) - Z(k-1) )  |
             !  term = |1   -  betam *  ------------------------ |
             !         |                          L              |
             !          -                                        -
             !
             term = (1.0_r8 - betam*zl(i))**onet
             !
             !               -           - 2
             !              |               z  |
             !  Kc = k*wi*z*|1   -   -----|
             !              |               h  |
             !               -           -
             !
             ! fak1(i) = ustar(i)*pblh(i)*vk
             !
             !                                                               -           - 2         -                                         - 1/3
             !                                     0.5 *( Z(k) - Z(k-1) )   |          z  |        |                  0.5 *( Z(k) - Z(k-1) )  |
             !  Kc =  ustar(i)* h(i) * k *  ------------------------ *|1   -   -----| * |1   -  betam *  ------------------------ |
             !                                           h               |          h  |        |                           L                  |
             !                                                               -           -         -                                          -
             !
             !                                                      -                 - 2   -                                        - 1/3
             !                                                    |               z  |   |                         0.5 *( Z(k) - Z(k-1) )  |
             !  Kc =  ustar(i)* k * 0.5 *( Z(k) - Z(k-1) ) *|1   -   -----| * |1   -  betam *  ------------------------  |
             !                                                     |               h  |   |                                L                 |
             !                                                      -                 -     -                                        -

             pblk(i) = fak1(i)*zh(i)*zzh(i)*term

             pr(i) = term/SQRT(1.0_r8 - betah*zl(i))
          END IF
       END DO
       !
       ! Unstable for outer layer; counter-gradient terms non-zero:
       !
       DO i=1,plon
          IF (unsout(i)) THEN
             !              -           - 2
             !             |          z  |
             !zzh =  |1   -   -----|
             !             |          h  |
             !              -           -
             !
             !
             !        0.5 *( Z(k) - Z(k-1) )
             !  zh = ------------------------
             !                 h
             !
             !                  ____
             !        DCo     d w*Co
             ! ----- = -------
             !  Dt      dz
             !
             !                     -         - 
             !        DCo     d   |      dCo        |   
             ! ----- = ---* |Kc * ----- | 
             !  Dt      dz  |       dz        |   
             !               -         - 
             ! fakn   = 7.2_r8     ! Constant in turbulent prandtl number

             !       ____
             !fak3 = w*Co = fakn*wstr(i)/wm(i)
             !
             ! fak2(i)    = wm(i)*h(i)*vk
             !
             !                                                        -               - 2
             !                               0.5 *( Z(k) - Z(k-1) )   |              z  |
             ! Kc = wm  * h  * vk *  ------------------------ *|1        -   -----|
             !                                         h             |              h  |
             !                                                        -               -
             !
             !                                                     -          - 2
             !                                                     |              z  |
             ! Kc = wm  *  vk *  0.5 *( Z(k) - Z(k-1) )   * |1        -   -----|
             !                                                    |         h  |
             !                                                     -          -
             !
             !               -           - 2
             !              |               z  |
             !  Kc = k*wi*z*|1   -   -----|
             !              |               h  |
             !               -           -
             !
             pblk(i) = fak2(i)*zh(i)*zzh(i)
             !         ____
             !         w*Co
             ! Yc = d*---------
             !         Wstar*h
             !
             cgs (i,k) = fak3(i)/(pblh(i)*wm(i))
             !
             !  khfs = surface kinematic heat flux [mK/s]
             !
             !
             ! 
             ! surface kinematic heat flux [mK/s] ustar*tstar
             !
             ! H=Ho=-rho*cp*ustar*tstar  
             !
             ! surface heat flux (W/m2)  (j/kg/kelvin)
             !
             !                  H           j       kg K   m*m*m           m K
             ! ustar*tstar = ---------   = ------------------------- =   ------------
             !                -rho*cp      m*m*s   j         kg              s
             !
             ! khfs(i) = shflx(i)*rrho(i)/cpair
             !
             !
             !cgh (i,k)  = counter-gradient term for heat [K/m]
             !
             cgh (i,k) = khfs(i)*cgs(i,k)
             !
             !       Qh             z     Wstar
             ! Pr = ---- + a * k * --- * -------
             !       Qm             h     Wm 
             !
             !
             !       phiminv              fak3 
             ! Pr = ----------  + ccon * ----- 
             !       phihinv              fak  
             !
             !    ccon = fak*sffrac*vk

             ! fakn   = 7.2_r8     ! Constant in turbulent prandtl number

             !       ____
             !fak3 = w*Co = fakn*wstr(i)/wm(i)
             !
             !       phiminv                        fakn*wstr(i)
             ! Pr = ----------  + fak*sffrac*vk * -------------------- 
             !       phihinv                         fak * wm(i)
             !
             !
             !       phiminv                         wstr(i)
             ! Pr = ----------  + fakn*vk*sffrac * ----------------- 
             !       phihinv                         wm(i)
             !

             pr  (i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
             
             !
             ! kqfs = sfc kinematic constituent flux [m/s]
             !
             !
             ! surface kinematic constituent flux [m/s]
             !
             ! E=Eo=-rho*cp*ustar*qstar
             !
             !                  E            kg          m * m * m          m 
             ! ustar*qstar = --------- = ----------- * ------------- = -----------
             !                -rho        m * m * s       kg                s
             !
             ! 
             ! surface constituent flux (kg/m2/s)
             !
             ! kqfs(i,1) = cflx(i,1)*rrho(i)
             !
             ! gq (i,k,1) = counter-gradient term for constituents
             !
             cgq (i,k,1) = kqfs(i,1)*cgs(i,k)
          END IF
       END DO
       DO m=2,pcnst
          DO i=1,plon
             IF (unsout(i)) cgq(i,k,m) = kqfs(i,m)*cgs(i,k)
          END DO
       END DO
       !
       ! For all unstable layers, set diffusivities
       !
       DO i=1,plon
          IF (unslev(i)) THEN
             !
             ! unstbl pbl with levels within pbl
             !
             kvm(i,k) = MAX(pblk(i)      ,kvf(i,k))
             kvh(i,k) = MAX(pblk(i)/pr(i),kvf(i,k))
          END IF
       END DO
1000   CONTINUE           ! end of level loop
       !+
       ! Check whether last allowed midpoint is within pbl, determine ktopblmn
       !-
       ktopblmn = plev
       k = plev-npbl+1
       DO i = 1, plon
          IF (z(i,k) .LT. pblh(i)) ktopbl(i) = k
          ktopblmn = MIN(ktopblmn, ktopbl(i))
       END DO
       ktopblmn=MAX(ktopblmn,2)
       !
       RETURN
     END SUBROUTINE pbldif





     SUBROUTINE qvdiff(ncnst   , &! INTENT(IN   ) :: ncnst                  ! number of constituents being diffused
                       qm1     , &! INTENT(IN   ) :: qm1   (plond,plev,ncnst)  ! initial constituent
                       qflx    , &! INTENT(IN   ) :: qflx  (plond,ncnst)      ! sfc q flux into lowest model level
                       cc      , &! INTENT(IN   ) :: cc    (plond,plev)         ! -lower diag coeff.of tri-diag matrix
                       ze      , &! INTENT(IN   ) :: ze    (plond,plev)         ! term in tri-diag. matrix system
                       term    , &! INTENT(IN   ) :: term  (plond,plev)       ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                       qp1     , &! INTENT(OUT  ) :: qp1   (plond,plev,ncnst)  ! final constituent
                       plon    , &! INTENT(IN   ) :: plon   ! number of longitudes
                       plond   , &! INTENT(IN   ) :: plond  ! slt extended domain longitude
                       plev    , &! INTENT(IN   ) :: plev   ! number of vertical levels
                       pcnst   , &! INTENT(IN   ) :: pcnst  ! number of constituents (including water vapor)
                       ntopfl   ) ! INTENT(IN   ) :: ntopfl ! Top level to which vertical diffusion is applied.
       !-----------------------------------------------------------------------
       !
       ! Solve vertical diffusion eqtn for constituent with explicit srfc flux.
       ! Procedure for solution of the implicit equation follows Richtmyer and 
       ! Morton (1967,pp 198-199).
       !
       !---------------------------Code history--------------------------------
       !
       ! Original version:  CCM1
       ! Standardized:      J. Rosinski, June 1992
       ! Reviewed:          P. Rasch, B. Boville, August 1992
       ! Reviewed:          P. Rasch, April 1996
       ! Reviewed:          B. Boville, May 1996
       !
       !-----------------------------------------------------------------------
       !
       ! $Id: qvdiff.F,v 1.1.1.1 2001/03/09 00:29:27 mirin Exp $
       !
       !-----------------------------------------------------------------------
       !-----------------------------------------------------------------------
       !------------------------------Commons----------------------------------
       !------------------------------Arguments--------------------------------
       !
       ! Input arguments
       !
       !      integer,parameter  :: plon   = 1 ! number of longitudes
       !      integer,parameter  :: plond  = 1 ! slt extended domain longitude
       !      integer,parameter  :: plev   =28 ! number of vertical levels
       !      integer,parameter  :: pcnst  = 1 ! number of constituents (including water vapor)
       !      integer,parameter  :: ntopfl =1 ! Top level to which vertical diffusion is applied.
       ! ntopfl = top level to which v-diff is applied
       IMPLICIT NONE
       INTEGER, INTENT(IN   )  :: plon   ! number of longitudes
       INTEGER, INTENT(IN   )  :: plond  ! slt extended domain longitude
       INTEGER, INTENT(IN   )  :: plev   ! number of vertical levels
       INTEGER, INTENT(IN   )  :: pcnst  ! number of constituents (including water vapor)
       INTEGER, INTENT(IN   )  :: ntopfl ! Top level to which vertical diffusion is applied.
       ! ntopfl = top level to which v-diff is applied

       INTEGER, INTENT(IN   )  :: ncnst                  ! number of constituents being diffused

       REAL(KIND=r8), INTENT(IN   ) :: qm1   (plond,plev,ncnst)  ! initial constituent
       REAL(KIND=r8), INTENT(IN   ) :: qflx  (plond,ncnst)      ! sfc q flux into lowest model level
       REAL(KIND=r8), INTENT(IN   ) :: cc    (plond,plev)         ! -lower diag coeff.of tri-diag matrix
       REAL(KIND=r8), INTENT(INOUT) :: ze    (plond,plev)         ! term in tri-diag. matrix system
       REAL(KIND=r8), INTENT(IN   ) :: term  (plond,plev)       ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
       !
       ! Output arguments
       !
       REAL(KIND=r8), INTENT(OUT  ) :: qp1   (plond,plev,ncnst)  ! final constituent
       !
       !---------------------------Local workspace-----------------------------
       !
       REAL(KIND=r8) :: zfq   (plond,plev,pcnst)  ! terms appear in soln of tri-diag sys
       REAL(KIND=r8) :: tmp1d (plond)           ! temporary workspace (1d array)

       INTEGER i,k                 ! longitude, vertical indices
       INTEGER m                   ! constituent index


       qp1=0.0_r8
       zfq=0.0_r8   
       tmp1d =0.0_r8
       !
       !-----------------------------------------------------------------------
       !
       ! Calculate fq(k).  Terms fq(k) and e(k) are required in solution of 
       ! tridiagonal matrix defined by implicit diffusion eqn.
       ! Note that only levels ntopfl through plev need be solved for.
       ! No vertical diffusion is applied above this level
       !
       DO m=1,ncnst
          DO i=1,plon
             zfq(i,ntopfl,m) = qm1(i,ntopfl,m)*term(i,ntopfl)
          END DO
          DO k=ntopfl+1,plev-1
             DO i=1,plon
                zfq(i,k,m) = (qm1(i,k,m) + cc(i,k)*zfq(i,k-1,m))*term(i,k)
             END DO
          END DO
       END DO
       !
       ! Bottom level: (includes  surface fluxes)
       !
       DO i=1,plon
          tmp1d(i) = 1.0_r8/(1.0_r8 + cc(i,plev) - cc(i,plev)*ze(i,plev-1))
          ze(i,plev) = 0.0_r8
       END DO
       DO m=1,ncnst
          DO i=1,plon
             zfq(i,plev,m) = (qm1(i,plev,m) + qflx(i,m) + &
                  cc(i,plev)*zfq(i,plev-1,m))*tmp1d(i)
          END DO
       END DO
       !
       ! Perform back substitution
       !
       DO m=1,ncnst
          DO i=1,plon
             qp1(i,plev,m) = zfq(i,plev,m)
          END DO
          DO k=plev-1,ntopfl,-1
             DO i=1,plon
                qp1(i,k,m) = zfq(i,k,m) + ze(i,k)*qp1(i,k+1,m)
             END DO
          END DO
       END DO
       RETURN
     END SUBROUTINE qvdiff
 
     SUBROUTINE MY20_TURB (       &
                     pcnst      ,& 
                     plond      ,& ! INTENT(IN   ) :: plond
                     plev       ,& ! INTENT(IN   ) :: plev
                     ztodt      ,& ! INTENT(IN   ) :: ztodt
                     colrad     ,& ! INTENT(IN   ) :: colrad        (plond)    
                     tm1        ,& ! INTENT(IN   ) :: tm1        (plond,plev)            ! temperature input
                     qp1        ,& ! INTENT(IN   ) :: qp1        (plond,plev,pcnst)  ! constituents after vdiff
                     um1        ,& ! INTENT(IN   ) :: um1        (plond,plev)            ! u-wind input
                     vm1        ,& ! INTENT(IN   ) :: vm1        (plond,plev)            ! v-wind input
                     gl0        ,& ! INTENT(INOUT) :: gl0        (plond)     !
                     KmMixl     ,& ! INTENT(OUT  ) :: KmMixl(plond,plev)   
                     KhMixl      ) ! INTENT(OUT  ) :: KhMixl(plond,plev)
      IMPLICIT NONE
      INTEGER       , INTENT(IN   ) :: pcnst
      INTEGER       , INTENT(IN   ) :: plond
      INTEGER       , INTENT(IN   ) :: plev
      REAL(kind=r8), INTENT(IN   ) :: ztodt
      REAL(kind=r8), INTENT(IN   ) :: colrad    (plond)    
      REAL(kind=r8), INTENT(IN   ) :: tm1       (plond,plev)        ! temperature input
      REAL(kind=r8), INTENT(IN   ) :: qp1       (plond,plev,pcnst)  ! constituents after vdiff
      REAL(kind=r8), INTENT(IN   ) :: um1       (plond,plev)        ! u-wind input
      REAL(kind=r8), INTENT(IN   ) :: vm1       (plond,plev)        ! v-wind input
      REAL(kind=r8), INTENT(INOUT) :: gl0       (plond)     !
      REAL(KIND=r8), INTENT(OUT  ) :: KmMixl    (plond,plev) 
      REAL(KIND=r8), INTENT(OUT  ) :: KhMixl    (plond,plev)
  
      REAL(kind=r8) :: twodt 
      REAL(kind=r8) :: twodti
      REAL(KIND=r8) :: a   (plev)
      REAL(KIND=r8) :: b   (plev)
      REAL(KIND=r8) :: c,x,y
      REAL(KIND=r8) :: Pbl_HgtLyI(plond,plev)
      REAL(KIND=r8) :: Pbl_ATemp (plond,plev) 
      REAL(KIND=r8) :: Pbl_Stabil(plond,plev)
      REAL(KIND=r8) :: Pbl_ITemp (plond,plev)
      REAL(KIND=r8) :: Pbl_Shear (plond,plev)
      REAL(KIND=r8) :: Pbl_KM    (plond,plev) 
      REAL(KIND=r8) :: Pbl_KH    (plond,plev) 
      REAL(KIND=r8) :: Pbl_BRich (plond,plev)
      REAL(KIND=r8) :: Pbl_NRich (plond,plev)
      REAL(KIND=r8) :: Pbl_ShBar (plond,plev)
      REAL(KIND=r8) :: Pbl_SmBar (plond,plev) 
      REAL(KIND=r8) :: Pbl_Sqrtw (plond,plev)
      REAL(KIND=r8) :: Pbl_EddEner(plond,plev) 
      REAL(KIND=r8) :: Pbl_MixLgh (plond,plev)
      REAL(KIND=r8) :: MixLgh     (plond,plev+1) 
      REAL(KIND=r8) :: Pbl_PotTep (plond,plev)
 
      REAL(KIND=r8) :: Pbl_KmMixl(plond,plev)
      REAL(KIND=r8) :: Pbl_KhMixl(plond,plev)
      REAL(KIND=r8) :: gln(plond)
      REAL(KIND=r8) :: gld(plond)
      REAL(KIND=r8) :: csqiv(plond)
      REAL(KIND=r8) :: s1ms2g 
      REAL(KIND=r8) :: fac
      INTEGER       :: i
      INTEGER       :: k
      INTEGER       :: itr
      KmMixl=0.0_r8
      KhMixl=0.0_r8 

      DO i=1,plond,1 
          csqiv(i)        = 1.0e0_r8/SIN(colrad(i))**2
      END DO
      twodt =ztodt
      twodti=1.0_r8/twodt
      DO k = 1, plev
          !                 --                  --
          !    1      1    | g       si(k+1)      |     con1(k)
          !  ------ =--- * |--- * ----------------| = ----------
          !    DZ     T    | R     sl(k) -sl(k+1) |       T
          !                 --                  --
          !                -- --
          !    DA      d  |     |
          !  ------ =---- | W'A'|
          !    DT      dZ |     |
          !                -- --
          !
          !                          ----
          !    DA      d         d  |    |
          !  ------ =---- * K * --- | A  |
          !    DT      dZ        dZ |    |
          !                          ----
          !                     
          !    DA         d      dA 
          !  ------ =K * ---- * ---- 
          !    DT         dZ     dZ 
          !                     
          !
          !
          !a0(k)  =gbyr*sigml(k+1)**2/(delsig(k  )*(sig(k)-sig(k+1)))
          !
          !                grav **2
          !              --------    * si(k+1)**2
          !                gasr**2 
          !a0(k)  = ---------------------------------
          !          ((si(k)-si(k+1))) * (sl(k)-sl(k+1)))
          !
          !         --  -- 2                             -- -- 2
          !        |   g  |       si(k+1)**2            |  1  | 
          !a0(k)  =| -----| * --------------------    = | --- |
          !        |   R  |     ((si(k)-si(k+1))  )     |  dZ | 
          !         --  --                               -- --
          !                   --   --      -- -- 2            --                   -- 2
          !                  |       |    |  1  |            |      m       kg * K   |  
          a(k)=twodt*a0(k)!  |  2*Dt |  * | --- |    ==> s * |   ------- * --------  | 
          !                  |       |    |  dZ |            |     s**2       J      |  
          !                   --   --      -- --              --                   --
          !J = F*DX = kg m/s**2 *m = kg * m**2/s**2
          !                   --   --      -- -- 2            --                       -- 2
          !                  |       |    |  1  |            |    m         kg * K *s**2 |    K**2 * s
          !a(k)=twodt*a0(k)! |  2*Dt |  * | --- |    ==> s * | ------- * --------------- | = -----------
          !                  |       |    |  dZ |            |   s**2     kg * m**2      |    m**2  
          !                   --   --      -- --              --                       --
          !                   --   --      -- -- 2  
          !                  |       |    |  1  |         K**2 * s
          !a(k)=twodt*a0(k)! |  2*Dt |  * | --- |    ==> -----------
          !                  |       |    |  dZ |           m**2  
          !                   --   --      -- --    

          !
          !                   gbyr*sigml(k+1)**2
          !    b0(k) = -----------------------------------------------
          !               (((si(k+1)-si(k+1+1))  ) *(sig(k)-sig(k+1)))

          b(k)=twodt*b0(k)! s * K**2/m**2 
      END DO
      c=twodt*c0
       !     
       !     Pbl_PotTep(3)   virtual potential temperature
       !     Pbl_ATemp (7)   temperature at the interface of two adjacent layers
       !
       !             k=2  ****gu,gv,gt,gq,gyu,gyv,gtd,gqd,sl*** } delsig(2)
       !             k=3/2----si,ric,rf,km,kh,b,l -----------
       !             k=1  ****gu,gv,gt,gq,gyu,gyv,gtd,gqd,sl*** } delsig(1)
       !             k=1/2----si ----------------------------
       !
       !                      sl(k+1)-si(k+1)
       !       t0    (   k)=-------------------
       !                      sl(k+1)-sl(k  )
       !
       !                      si(k+1)-sl(k  )
       !       t1    (   k)=-------------------
       !                      sl(k+1)-sl(k  )
       !
       !
       !                      sl(k+1)-si(k+1)              si(k+1)-sl(k  )
       !   Pbl_ATemp(i,k)  =-------------------*gt(i,k) + ------------------*gt(i,k+1)
       !                      sl(k+1)-sl(k  )              sl(k+1)-sl(k  )
       !
       !
       !
       !
       DO k = 1, plev-1
          DO i = 1, plond
             Pbl_ATemp(i,k)  = t0(k)*tm1 (i,(plev+1-k)) + t1(k)*tm1(i,(plev+1-k)-1)
          END DO
       END DO
       !                      si(k+1)-sl(k  )
       !       t1    (   k)=-------------------
       !                      sl(k+1)-sl(k  )
       !
       !                          1
       !      sigkiv(k   )=--------------------
       !                      (sl(k)**akappa)
       !
       !                       --                       --  
       !                      | grav        si(k+1)       |
       !   con2  (   k)=grav* |------- * ---------------- |  
       !                      | gasr      sl(k)-sl(k+1))  |
       !                       --                       -- 
       !
       !                --                       --  2
       !               |   grav        si(k+1)     |
       !   con1  (k) = |------- * ---------------- |         
       !               |  gasr      sl(k)-sl(k+1)) |
       !                --                       -- 
       !
       ! sl.........sigma value at midpoint of                  gasr/cp
       !                                         each layer : (k=287/1005)=R/cp
       !
       !                                                                     1
       !                                             +-                   + ---
       !                                             !     k+1         k+1!  k
       !                                             !si(l)   - si(l+1)   !
       !                                     sl(l) = !--------------------!
       !                                             !(k+1) (si(l)-si(l+1)!
       !             --   --  -(R/Cp)                +-                  -+     
       !            |  P    |
       !Thetav = Tv |-------|
       !            |  P0   | 
       !             --   --
       !
       !             --   --  -(R/Cp)        --   --  -(R/Cp)
       !            |  P    |               |       |
       !sigkiv(k) = |-------|            == |sl(k)  |
       !            |  P0   |               |       | 
       !             --   --                 --   --
       !  --   --      --   -- 
       ! |  P    |    |       |
       ! |-------| =  |sl(k)  |
       ! |  P0   |    |       |
       !  --   --      --   --
       !  --          --      --           -- 
       ! |  P(k)-P(k+1) |    |               |
       ! |--------------| =  |sl(k)-sl(k+1)  |
       ! |      P0      |    |               |
       !  --          --      --           --
       DO k = 1, plev
          DO i = 1, plond
             Pbl_PotTep(i,k)=sigkiv(k)*tm1 (i,(plev+1-k))*(1.0_r8+eps*qp1(i,(plev+1-k),1))
          END DO
       END DO
       !     
       !     Pbl_Stabil(2)   stability
       !     Pbl_Shear (6)   square of vertical wind shear
       !  
       gln =1.0e-5_r8 
       DO k = 1, plev-1
          DO i = 1, plond
             !
             !
             !                 g         D (Theta)
             !Pbl_Stabil(2) =------- * -------------
             !                Theta      D (Z)
             !
             !
             ! P =rho*R*T and P = rho*g*Z
             !
             !                           P
             ! DP = rho*g*DZ and rho = ----
             !                          R*T
             !        P
             ! DP = ----*g*DZ
             !       R*T
             !
             !        R*T
             ! DZ = ------*DP
             !        g*P
             !
             !    1       g*P       1   
             !  ------ = ------ * ------
             !    DZ      R*T       DP  
             !
             !    T       g         P   
             !  ------ = ------ * ------
             !    DZ      R         DP  
             !
             !  --   --      --   -- 
             ! |  P    |    |       |
             ! |-------| =  |sl(k)  |
             ! |  P0   |    |       |
             !  --   --      --   --
             !
             !    T       g            si(k+1)   
             !  ------ = ------ * ----------------- =con1(k)
             !    DZ      R         sl(k) -sl(k+1)  
             !
             !                 --                  --
             !    1      1    | g       si(k+1)      |     con1(k)
             !  ------ =--- * |--- * ----------------| = ----------
             !    DZ     T    | R     sl(k) -sl(k+1) |       T
             !                 --                  --
             !                                  --                  --
             !                    g       1    | g       si(k+1)      |   
             ! Pbl_Stabil(2) = ------- * --- * |--- * ----------------|* D(Theta)
             !                  Theta     T    | R     sl(k) -sl(k+1) | 
             !                                  --                  --
             !
             !                   con2(k)*(Pbl_PotTep(i,k+1)-Pbl_PotTep(i,k))
             !Pbl_Stabil(i,k)=-----------------------------------------------------------------
             !                  ((t0(k)*Pbl_PotTep(i,k)+t1(k)*Pbl_PotTep(i,k+1))*Pbl_ATemp(i,k))


             Pbl_Stabil(i,k)=con2(k)*(Pbl_PotTep(i,k+1)-Pbl_PotTep(i,k)) &
                  /((t0(k)*Pbl_PotTep(i,k)+t1(k)*Pbl_PotTep(i,k+1))*Pbl_ATemp(i,k))
                  
             Pbl_ITemp(i,k)=1.0_r8/(Pbl_ATemp(i,k)*Pbl_ATemp(i,k))
             !
             !                  --  -- 2    --  -- 2
             !                 |  dU  |    |  dV  |
             !Pbl_Shear(i,k) = | -----|  + | -----| 
             !                 |  DZ  |    |  DZ  |
             !                  --  --      --  --
             !
             !                  --  -- 2    --     -- 2
             !                 |   1  |    |         |
             !Pbl_Shear(i,k) = | -----|  * | DU + DV | 
             !                 |  DZ  |    |         |
             !                  --  --      --     --
             !
             !                         --                  -- 2    --     -- 2
             !                 1      | g       si(k+1)      |    |         |          1.0e0_r8
             !Pbl_Shear(i,k) =----- * |--- * ----------------|  * | Du + Dv |   * ---------------------
             !                 T*T    | R     sl(k) -sl(k+1) |    |         |       SIN(colrad(i))**2
             !                         --                  --      --     --
             !
             ! --  -- 2              --                       -- 2
             !|   1  |     1          | g            si(k+1)         | 
             !| -----|  = ----- * |--- * ----------------| 
             !|  DZ  |     T*T    | R          sl(k) -sl(k+1) | 
             ! --  --                --                       --  
             !
             ! --  --               --                       -- 
             !|   1  |     1          | g            si(k+1)         | 
             !| -----|  = ----- * |--- * ----------------| = SQRT((Pbl_ITemp(i,k) * con1(k)))
             !|  DZ  |     T      | R          sl(k) -sl(k+1) | 
             ! --  --                --                       --  
             !
             !                    --                        -- 
             !       1           | g           si(k+1)          |                  1
             !DZ = ----- *       |--- * ----------------| = ---------------------------------
             !       T           | R         sl(k) -sl(k+1)   |   SQRT((Pbl_ITemp(i,k) * con1(k)))
             !                    --                        --  
             !
             !                  1.0e0_r8
             ! csqiv (i) = ---------------------
             !               SIN(colrad(i))**2
             !
             Pbl_Shear(i,k)=(con1(k))*Pbl_ITemp(i,k)*((um1(i,(plev+1-k))-um1(i,(plev+1-k)-1))**2+&
             (vm1(i,(plev+1-k))-vm1(i,(plev+1-k)-1))**2)
             Pbl_Shear(i,k)=MAX(gln(i),Pbl_Shear(i,k))          
          END DO
       END DO
       !     
       !     Pbl_BRich(4)        richardson number
       !     Pbl_NRich(5)   flux richardson number
       !     Pbl_NRich(8)   flux richardson number
       !     
       Pbl_KmMixl=0.0_r8 ;Pbl_Km =0.0_r8
       Pbl_KhMixl=0.0_r8 ;Pbl_Kh =0.0_r8
       DO k = 1,(plev-1)
          DO i = 1, plond
             !                g         D (Theta)
             !              ------- * -------------
             !               Theta      D (Z)
             ! Pbl_BRich = -----------------------------
             !                 --  -- 2    --  -- 2
             !                |  dU  |    |  dV  |
             !                | -----|  + | -----| 
             !                |  DZ  |    |  DZ  |
             !                 --  --      --  --
             !
             Pbl_BRich(i,k)=Pbl_Stabil(i,k)/Pbl_Shear(i,k)
             !r1 = 0.5_r8*gama/alfa
             Pbl_NRich(i,k)= r1*(Pbl_BRich(i,k)+r2 &
                  -SQRT(Pbl_BRich(i,k)*(Pbl_BRich(i,k)-r3)+r4))
             Pbl_NRich(i,k)=MIN(rfc,Pbl_NRich(i,k))
             Pbl_NRich(i,k)=Pbl_NRich(i,k)
             !
             !    Pbl_SmBar and Pbl_ShBar are momentum flux and heat flux 
             !    stability parameters, respectively
             !     
             !     Pbl_ShBar(3)   shbar
             !     Pbl_SmBar(4)   smbar
             !     
             !     eliminate negative value for s1-s2*gwrk(i,1,5):
             !     gwrk(i,1,5) is s1/s2 under some circumstances
             !     which makes this expression zero.  machine roundoff
             !     can produce an unphysical negative value in this case.
             !     it is used as sqrt argument in later loop.
             !     s1  =3.0_r8*a2* gam1
             !     s2  =3.0_r8*a2*(gam1+gam2)
             !
             ! s1-s2 = 3.0*a2* gam1 - 3.0*a2*(gam1+gam2)
             ! s1-s2 = 3.0*a2* (gam1 - (gam1+gam2))
             !     
             s1ms2g=s1-s2*Pbl_NRich(i,k)
             IF (ABS(s1ms2g) < 1.0e-10_r8) s1ms2g=0.0_r8
             !
             !                     s1ms2g
             !Pbl_ShBar(i,k)=---------------------------
             !                 (1.0_r8-Pbl_NRich(i,k))
             !
             ! a1    =   0.92_r8
             ! a2    =   0.74_r8
             ! b1    =  16.6_r8
             ! b2    =  10.1_r8
             !
             !         1.0               a1
             !gam1 = -------  - 2.0  * ------
             !         3.0               b1
             !
             !       b2          a1
             !gam2=(----) + (6*------)
             !       b1          b1
             !
             !gam2=(b2  + 6.0*a1)
             !     --------------
             !          b1
             !                 --                   --
             !                | 1.0               a1  |
             !s1 = 3.0 * a2 * |------  - 2.0  * ------|
             !                | 3.0               b1  |
             !                 --                   --
             !s2  =3.0_r8*a2*(gam1+gam2)
             !
             !                   s1-s2*Pbl_NRich(i,k)
             !Pbl_ShBar(i,k)=---------------------------
             !                 (1.0_r8-Pbl_NRich(i,k))
             !
             !                        ((gam1 - (gam1+gam2)))*Pbl_NRich(i,k)
             !Pbl_ShBar(i,k)=3.0*a2* ---------------------------------------- =Sm
             !                           (1.0_r8-Pbl_NRich(i,k))

             Pbl_ShBar(i,k)=s1ms2g/(1.0_r8-Pbl_NRich(i,k))
             !     
             !     gwrk(i,1,3)=(s1-s2*gwrk(i,1,5))/(1.0_r8-gwrk(i,1,5))
             !     end of  negative sqrt argument trap
             ! c1    =   0.08_r8
             ! alfa=b1*(gam1-c1)+3.0_r8*(a2+2.0_r8*a1)
             ! beta=b1*(gam1-c1)
             !        --                              --
             !       | a2                               |
             ! gama=-|---- *(b1*(gam1+gam2)-3.0_r8*a1)  |
             !       | a1                               |
             !        --                              --
             !       --          --
             !      | a2           |
             ! dela=|---- b1* gam1 |
             !      | a1           |
             !       --          --
             !                                     ((b1*(gam1-c1))-(b1*(gam1-c1)+3.0*(a2+2.0*a1))*Pbl_NRich(i,k))
             ! Pbl_SmBar(i,k)=Pbl_ShBar(i,k) * ------------------------------------------------------------------
             !                                     --          --     --                           --
             !                                    | a2           |   | a2                            |
             !                                    |---- b1* gam1 | - |---- *(b1*(gam1+gam2)-3.0*a1)  |*Pbl_NRich(i,k))
             !                                    | a1           |   | a1                            |
             !                                     --          --     --                           --
             !
             !                                     ((b1*(gam1-c1)) - (b1*(gam1-c1) + 3.0*a2 + 6.0*a1))*Pbl_NRich(i,k))
             ! Pbl_SmBar(i,k)=Pbl_ShBar(i,k) * --------------------------------------------------------------
             !                                     --  --     --                             --
             !                                    | a2   |   |                                 |
             !                                    |----  | * |b1* gam1 -(b1*(gam1+gam2)-3.0*a1)|*Pbl_NRich(i,k))
             !                                    | a1   |   |                                 |
             !                                     --  --     --                              --

             !   
             Pbl_SmBar(i,k)=Pbl_ShBar(i,k)*(beta-alfa*Pbl_NRich(i,k))/ &
                  (dela-gama*Pbl_NRich(i,k))
             !     
             !     
             !     Pbl_Sqrtw(5)   sqrt(w) or b/l
             !     Pbl_KmMixl(4)   km/l**2
             !     Pbl_KhMixl(3)   kh/l**2
             ! 
             !     The ratio of SH to SM is equal to the ratio of the turbulent
             !     flux Richardson number to the bulk (large scale) Richardson
             !     number. 
             !
             !     u^2 =  (1-2*gam1)q^2
             !     
             !     SQRT(B1*SM*(1 - Ri)*Shear)
             !             GH    SM
             !     Ri = - ---- =-----Rf
             !             GM    SH
             !             GH
             !     GM = - ----
             !             Ri   
             !                --                  --  
             !               | --  -- 2    --  -- 2 |
             !          l^2  ||  dU  |    |  dV  |  |
             !     GM =-----*|| -----|  + | -----|  |
             !          q^2  ||  DZ  |    |  DZ  |  |
             !               | --  --               |
             !                --                  --  
             ! 
             !                          --      -- 2 
             !             l^2         |  dTHETA  |  
             !     GH = - -----*beta*g*| ---------|  
             !             q^2         |  DZ      |  
             !                          --      --   

             !
             !          --                       --  1/2
             !         |      --  -- 2    --  -- 2 |
             !     q   | 1   |  dU  |    |  dV  |  |
             !    ----=|----*| -----|  + | -----|  |
             !     l   | GM  |  DZ  |    |  DZ  |  |
             !         |      --  --               |
             !          --                       --   
             !                                               
             !  1  
             ! ---- = B1*Sm*(1 - Rf)
             !  GM 
             !
             Pbl_Sqrtw (i,k)=SQRT(b1*Pbl_SmBar(i,k)*(1.0_r8-Pbl_NRich(i,k))*Pbl_Shear(i,k))
             !
             ! Km = l*q*Sm
             !  
             ! Km     q 
             !---- = ---*Sm 
             ! l^2    l
             !
             ! Kh = l*q*Sh
             !
             !  
             ! Kh     q 
             !---- = ---*Sh 
             ! l^2    l
             !
             Pbl_KmMixl(i,k)=Pbl_Sqrtw(i,k)*Pbl_SmBar(i,k)
             Pbl_KhMixl(i,k)=Pbl_Sqrtw(i,k)*Pbl_ShBar(i,k)
          END DO
       END DO
      !          
      !          Pbl_HgtLyI(1)   height at the layer interface 
      !
      !             R*T
      ! DZ = ------ * DP
      !             g*P
      !
      !             R*T      DP
      ! DZ = ------ * -----
      !              g       P
      !
      !              R      DP
      ! DZ = ------ * ----- * T
      !              g       P
      !
      !           
      ! DZ = con0(k)  * T
      !
      !                                                   kg*m*m    
      !                                                  --------  
      !                   R          DP           (j/kg/k)           kg*K*s*s  
      ! con0(k) = ------ * ----- =  --------   =   ------------         = m/K
      !                   g           P            m/s*s              m       
      !                                                    -----    
      !                                                     s*s          
      !                 gasr             si(k) - si(k+1)
      ! con0(k)= -------- * ------------------------
      !                 grav                  sl(k)
      !
      ! Pbl_HgtLyI = con0(k)  * T
      !
      k=1
      DO i = 1, plond
         Pbl_HgtLyI(i,k)=con0(k)* tm1 (i,plev)
      END DO
    
      DO k = 2, plev
         DO i = 1, plond
             Pbl_HgtLyI(i,k)=Pbl_HgtLyI(i,k-1)+con0(k)*tm1 (i,(plev+1-k))
         END DO
      END DO


       DO itr = 1, nitr
          !     
          !     Pbl_EddEner(2)   mixing length
          !     Pbl_EddEner(2)   b     :b**2 is eddy enegy
          !     
          !..gl0    maximum mixing length l0 in blackerdar's formula
          !         this is retained as a first guess for next time step
          !                                      k0*z
          !                                l = --------
          !                                    (1 + k0*z/l0)

          gln = 0.0_r8
          gld = 0.0_r8
          Pbl_EddEner(1:plond,1) = 0.0_r8
          DO k = 1, plev-1
             DO i = 1, plond
                !
                !                 vk0*gl0(i)*Z(i,k) 
                ! Pbl_EddEner = ------------------------------
                !                (gl0(i) + vk0*Z(i,k))
                !
                Pbl_EddEner(i,k+1)=VKARMAN*gl0(i)*Pbl_HgtLyI(i,k) / (gl0(i)+VKARMAN*Pbl_HgtLyI(i,k))
                !
                !                 vk0*gl0(i)*Z(i,k)           q          q
                ! Pbl_EddEner = ------------------------ * --------- = ---------
                !                (gl0(i) + vk0*Z(i,k))        l^2        l
                !
                Pbl_EddEner(i,k+1)=Pbl_EddEner(i,k+1)*Pbl_Sqrtw(i,k)
             END DO
          END DO
          k=1
          DO i = 1, plond
             Pbl_EddEner(i,k)= 1.0e-3_r8
          END DO
          k=1
          DO i = 1, plond
             x=0.5_r8*delsig2(k)*(Pbl_EddEner(i,k)+Pbl_EddEner(i,k+1))
             y=x*0.5_r8*Pbl_HgtLyI(i,k)
             gld(i)=gld(i)+x
             gln(i)=gln(i)+y
          END DO
          k=plev
          DO i = 1, plond
             x=0.5_r8*delsig2(k)*Pbl_EddEner(i,k)
             y=x*0.5_r8*(Pbl_HgtLyI(i,k)+Pbl_HgtLyI(i,k-1))
             gln(i)=gln(i)+y
             gld(i)=gld(i)+x
          END DO
          IF (plev > 2) THEN
             DO k = 2, plev-1
                DO i = 1, plond
                   x=0.5_r8*delsig2(k)*(Pbl_EddEner(i,k)+Pbl_EddEner(i,k+1))
                   y=x*0.5_r8*(Pbl_HgtLyI(i,k)+Pbl_HgtLyI(i,k-1))
                   gln(i)=gln(i)+y
                   gld(i)=gld(i)+x
                END DO
             END DO
          END IF
          DO i = 1, plond
             !                0.5 * del(k) * (q(i,k)+q(i,k+1)) * 0.5*(Z(i,k)+Z(i,k-1))
             !gl0(i)=facl * -----------------------------------------------------------
             !                        0.5 * del(k) * (q(i,k)+q(i,k+1))
             gl0(i)=facl*gln(i)/gld(i)
          END DO
          !     
          !     iteration that determines mixing length
          !     
       END DO


       !
       !          Pbl_MixLgh(5)   mixing length
       !        
       MixLgh(:,1)=0.0_r8
       Pbl_MixLgh=0.0_r8
       DO k = 1, plev-1
          DO i=1,plond
                 !           vk0*gl0(i)*Z(i,k)             m^2
                 ! l  = -------------------------------- = ------
                 !           gl0(i) + vk0*Z(i,k)       m
                 Pbl_MixLgh(i,k)=VKARMAN*gl0(i)*Pbl_HgtLyI(i,k)/(gl0(i)+VKARMAN*Pbl_HgtLyI(i,k))
             MixLgh(i,plev-k) = (Pbl_MixLgh(i,k))**2
          END DO
       END DO
       MixLgh(:,plev+1) = 0.0_r8
       !     
       !     Pbl_CoefKm(1)   km = l*q*Sm
       !     Pbl_CoefKh(2)   kh = l*q*Sh
       ! 
       !     where KM and KH are the diffusion coefficients for momentum 
       !     and heat, respectively, l is the master turbulence length scale, 
       !     q2 is the turbulent kinetic energy (so q is the magnitude of 
       !     the turbulent wind velocity), and SM and SH are momentum 
       !     flux and heat flux stability parameters, respectively    
       !    
       IF (kmean == 0) THEN
          DO k = 1, plev-1
             DO i = 1, plond
                !
                !             Km 
                !Km = l^2 * -----
                !            l^2
                !
                KmMixl(i,plev+1-k)=MIN(gkm1,MAX(gkm0,Pbl_MixLgh(i,k)**2*Pbl_KmMixl(i,k)))
                !
                !             Kh 
                !Kh = l^2 * -----
                !            l^2
                !
                KhMixl(i,plev+1-k)=MIN(gkh1,MAX(gkh0,Pbl_MixLgh(i,k)**2*Pbl_KhMixl(i,k)))
             END DO
          END DO
       ELSE
          DO k = 1, plev-1
             DO i = 1, plond
                !
                !             Km 
                !Km = l^2 * -----
                !            l^2
                !
                Pbl_KM(i,plev+1-k)=Pbl_MixLgh(i,k)**2*Pbl_KmMixl(i,k)
                !
                !             Kh 
                !Kh = l^2 * -----
                !            l^2
                !
                Pbl_KH(i,plev+1-k)=Pbl_MixLgh(i,k)**2*Pbl_KhMixl(i,k)
             END DO
          END DO
          fac=0.25_r8
          IF (plev >= 4) THEN
             DO k = 2, plev-2
                DO i = 1, plond
                   !
                   !             k=2  ****Km(k),sl*** } -----------
                   !             k=3/2----si,ric,rf,km,kh,b,l -----------
                   !             k=1  ****Km(k),sl*** } -----------
                   !             k=1/2----si ----------------------------
                   !
                   !       Km(k-1) + 2*Km(k) + Km(k+1)
                   ! Km = -------------------------------
                   !                  4
                   !
                   KmMixl(i,k)=fac*(Pbl_KM(i,k-1)+2.0_r8*Pbl_KM(i,k)+Pbl_KM(i,k+1))
                   !
                   !       Kh(k-1) + 2*Kh(k) + Kh(k+1)
                   ! Kh = -------------------------------
                   !                  4
                   !
                   KhMixl(i,k)=fac*(Pbl_KH(i,k-1)+2.0_r8*Pbl_KH(i,k)+Pbl_KH(i,k+1))
                END DO
             END DO
          END IF
          DO i = 1, plond
             !
             !       Km(1) + Km(2)
             ! Km = ---------------
             !            2
             !
             !
             !       Kh(1) + Kh(2)
             ! Kh = ---------------
             !            2
             !
             KmMixl(i,     1)=0.5_r8*(Pbl_KM(i,     1)+Pbl_KM(i,     2))
             KhMixl(i,     1)=0.5_r8*(Pbl_KH(i,     1)+Pbl_KH(i,     2))
             KmMixl(i,plev-1)=0.5_r8*(Pbl_KM(i,plev-1)+Pbl_KM(i,plev-2))
             KhMixl(i,plev-1)=0.5_r8*(Pbl_KH(i,plev-1)+Pbl_KH(i,plev-2))
          END DO
          DO k = 1, plev-1
             DO i = 1, plond
                KmMixl(i,k)=MIN(gkm1,MAX(gkm0,KmMixl(i,k))) !(m/sec**2)
                KhMixl(i,k)=MIN(gkh1,MAX(gkh0,KhMixl(i,k))) !(m/sec**2)
             END DO
          END DO 
       END IF

 END SUBROUTINE  MY20_TURB    
     
     
  SUBROUTINE MY25_TURB(&
       nCols    ,& ! nCols
       kMax     ,& ! kMax
       delt     ,& ! delt    !       delt     -  Time step in seconds
       um       ,& ! um       (1:nCols,1:kMax)!Zonal wind components
       vm       ,& ! vm       (1:nCols,1:kMax)!Meridional Wind components
       theta    ,& ! theta    (1:nCols,1:kMax)!theta         -  Potential temperature
       zhalf    ,& ! zhalf    (1:nCols,1:kMax+1)!zfull         -  Height at full levels
       zfull    ,& ! zfull    (1:nCols,1:kMax)!         zhalf    -  Height at half levels
       phalf    ,& ! phalf    (1:nCols,1:kMax+1)!!          phalf    -  Pressure at half levels, 
       pfull    ,& ! pfull    (1:nCols,1:kMax)!!          pfull    -  Pressure at full levels! , zhalf, zfull
       z0       ,& ! z0       (1:nCols)!          z0           -  Roughness length
       fracland ,& ! fracland (1:nCols)!        fracland -  Fractional amount of land beneath a grid box
       ustar    ,& ! ustar    (1:nCols)!        ustar         -  OPTIONAL:friction velocity (m/sec)
       el0      ,& ! el0      (1:nCols)!       el0  -  characteristic length scale
       akm      ,& ! akm      (1:nCols,1:kMax+1)!       akm  -  mixing coefficient for momentum
       akh      ,& ! akh      (1:nCols,1:kMax+1)!       akh  -  mixing coefficient for heat and moisture
       el       ,& ! el       (1:nCols,1:kMax+1)!       el   -  master length scale
       TKE      )
    IMPLICIT NONE
    INTEGER, INTENT(IN   ) :: nCols
    INTEGER, INTENT(IN   ) :: kMax
    REAL(KIND=r8)   , INTENT(in   ) :: delt    !       delt     -  Time step in seconds
    REAL(KIND=r8)   , INTENT(in   ) :: um       (1:nCols,1:kMax)!Zonal wind components
    REAL(KIND=r8)   , INTENT(in   ) :: vm       (1:nCols,1:kMax)!Meridional Wind components
    REAL(KIND=r8)   , INTENT(in   ) :: theta    (1:nCols,1:kMax)!theta    -  Potential temperature
    REAL(KIND=r8)   , INTENT(in   ) :: zfull    (1:nCols,1:kMax)!zfull    -  Height at full levels
    REAL(KIND=r8)   , INTENT(in   ) :: zhalf    (1:nCols,1:kMax+1)!       zhalf    -  Height at half levels
    REAL(KIND=r8)   , INTENT(in   ) :: phalf    (1:nCols,1:kMax+1)!!       phalf    -  Pressure at half levels, 
    REAL(KIND=r8)   , INTENT(in   ) :: pfull    (1:nCols,1:kMax)!!       pfull    -  Pressure at full levels! , zhalf, zfull
    REAL(KIND=r8)   , INTENT(in   ) :: z0       (1:nCols)   !          z0          -  Roughness length
    REAL(KIND=r8)   , INTENT(in   ) :: fracland (1:nCols)!       fracland -  Fractional amount of land beneath a grid box
    REAL(KIND=r8)   , INTENT(in   ) :: ustar    (1:nCols)!       ustar    -  OPTIONAL:friction velocity (m/sec)
    REAL(KIND=r8)   , INTENT(out  ) :: el0      (1:nCols)
    REAL(KIND=r8)   , INTENT(out  ) :: akm      (1:nCols,1:kMax+1)
    REAL(KIND=r8)   , INTENT(out  ) :: akh      (1:nCols,1:kMax+1)
    REAL(KIND=r8)   , INTENT(out  ) :: el       (1:nCols,1:kMax+1)
    REAL(KIND=r8)   , INTENT(inout) :: TKE      (1:nCols,1:kMax+1)

    !       mask     -  OPTIONAL; floating point mask (0. or 1.) designating
    !                   where data is present

    REAL(KIND=r8)    :: zsfc  (1:nCols)
    REAL(KIND=r8)    :: el2   (1:nCols,1:kMax-1)

    REAL(KIND=r8)    :: dsdzh (1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: qm2   (1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: qm3   (1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: qm4   (1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: buoync(1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: dsdz  (1:nCols,1:kMax) 
    REAL(KIND=r8)    :: xx1   (1:nCols,1:kMax) 
    REAL(KIND=r8)    :: xx2   (1:nCols,1:kMax) 
    REAL(KIND=r8)    :: qm    (1:nCols,1:kMax) 
    REAL(KIND=r8)    :: x1    (1:nCols)
    REAL(KIND=r8)    :: x2    (1:nCols)
    REAL(KIND=r8)    :: akmin (1:nCols)
    REAL(KIND=r8)    :: aaa   (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: bbb   (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: ddd   (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: ccc   (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: xxm1  (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: xxm2  (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: xxm3  (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: xxm4  (1:nCols,1:kMax-1) 
    REAL(KIND=r8)    :: xxm5  (1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: shear (1:nCols,1:kMax-1)
    REAL(KIND=r8)    :: cvfqdt
    REAL(KIND=r8)    :: dvfqdt
    INTEGER  :: i
    INTEGER  :: k
    INTEGER  :: it
    INTEGER  :: klim
    el0= 0.0_r8;akm= 0.0_r8;akh= 0.0_r8;el =0.0_r8 
    zsfc= 0.0_r8;el2= 0.0_r8;dsdzh= 0.0_r8;qm2= 0.0_r8;qm3= 0.0_r8;qm4= 0.0_r8   
    buoync= 0.0_r8;dsdz= 0.0_r8; xx1 = 0.0_r8; xx2 = 0.0_r8;qm = 0.0_r8;x1  = 0.0_r8;  
    x2  = 0.0_r8;akmin = 0.0_r8;aaa  = 0.0_r8;bbb  = 0.0_r8;ddd  = 0.0_r8;ccc  = 0.0_r8; 
    xxm1 = 0.0_r8;xxm2 = 0.0_r8;xxm3 = 0.0_r8;xxm4  = 0.0_r8;xxm5  = 0.0_r8;shear = 0.0_r8;
    cvfqdt= 0.0_r8;dvfqdt= 0.0_r8
    !====================================================================
    ! --- SURFACE HEIGHT     
    !====================================================================

    !if( PRESENT( kbot ) ) then
    !  do i = 1,nCols
    !      k = kbot(i) + 1
    !     zsfc(i) = zhalf(i,k)
    !  end do
    !else
    DO i = 1,nCols
       zsfc(i) = zhalf(i,kMax+1)
    END DO
    !endif
    !====================================================================
    ! --- D( )/DZ OPERATORS: AT FULL LEVELS & AT HALF LEVELS          
    !====================================================================
    DO k = 1,kMax
       DO i = 1,nCols
          dsdz (i,k)  = 1.0_r8 / ( zhalf(i,k+1) - zhalf(i,k) )
       END DO
    END DO
    DO k = 1,kMax-1
       DO i = 1,nCols
          dsdzh(i,k) = 1.0_r8 / ( zfull(i,k+1)  - zfull(i,k) )
       END DO
    END DO
    !====================================================================
    ! --- WIND SHEAR                 
    !====================================================================
    DO k = 1,kMax-1
       DO i = 1,nCols
          xxm1 (i,k) = dsdzh(i,k) * ( um(i,k+1) - um(i,k) )
          xxm2 (i,k) = dsdzh(i,k) * ( vm(i,k+1) - vm(i,k) )
          shear(i,k) = xxm1 (i,k) * xxm1(i,k) + xxm2(i,k) * xxm2(i,k)
       END DO
    END DO

    !====================================================================
    ! --- BUOYANCY                 
    !====================================================================
    DO k = 1,kMax-1
       DO i = 1,nCols
          xxm1(i,k) = theta(i,k+1) - theta(i,k) 
          IF( do_thv_stab ) THEN
             xxm2(i,k) = 0.5_r8*( theta(i,k+1) + theta(i,k) )
          ELSE
             xxm2(i,k) = t00
          END IF

          buoync(i,k) = grav * dsdzh(i,k) * xxm1(i,k) / xxm2(i,k)

       END DO
    END DO

    !====================================================================
    ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE
    !====================================================================
    !IF(jdt<=2)
     CALL TKE_SURF ( nCols                            , &
                     kMax                             , &
                     ustar(1:nCols)                   , &
                     TKE  (1:nCols,1:kMax+1)    )
    !  if( PRESENT( mask ) ) then
    !    do k=2,kMax
    !       do i=1,nCols
    !         if(mask(i,k) < 0.1_r8) then
    !            TKE(ism+i,jsm+j,k+1) = 0.0_r8
    !            dsdz(i,j,k  ) = 0.0_r8
    !            dsdzh(i,j,k-1) = 0.0_r8
    !            shear(i,j,k-1) = 0.0_r8
    !            buoync(i,j,k-1) = 0.0_r8
    !         endif
    !       enddo
    !      enddo
    !    enddo
    !  endif

    !====================================================================
    ! --- SET ITERATION LOOP IF INITALIZING TKE
    !====================================================================


    ! $$$$$$$$$$$$$$$$$
    DO it = 1,init_iters
       ! $$$$$$$$$$$$$$$$$

       !====================================================================
       ! --- SOME TKE STUFF
       !====================================================================

       DO k=1,kMax
          DO i=1,nCols
             xx1(i,k) = MAX(2*TKE(i,k+1),1.0e-12_r8)
             IF(xx1(i,k) > 0.0_r8) THEN
                qm(i,k) = SQRT(xx1(i,k))
             ELSE
                qm(i,k) = 0.0001_r8
             END IF
          END DO
       END DO

       DO k=1,kMax-1
          DO i=1,nCols
             qm2(i,k)  = xx1(i,k) 
             qm3(i,k)  = qm (i,k) * qm2(i,k) 
             qm4(i,k)  = qm2(i,k) * qm2(i,k) 
          END DO
       END DO

       !====================================================================
       ! --- CHARACTERISTIC LENGTH SCALE                         
       !====================================================================
       DO k=1,kMax-1
          DO i=1,nCols
!             xx1(i,k) = qm(i,k)*( pfull(i,k+1) - pfull(i,k) )
             xx1(i,k) = qm(i,k) * ( MAX(pfull(i,k+1) - pfull(i,k),0.0000001_r8 ))
          END DO
       END DO

       DO k = 1, kMax-1
          DO i=1,nCols
             xx2(i,k) = xx1(i,k)  * ( zhalf(i,k+1) - zsfc(i) )
          END DO
       END DO

       !if( PRESENT( kbot ) ) then
       !  xx1(:,:,kMax) = 0.0_r8
       !  xx2(:,:,kMax) = 0.0_r8
       !   do i = 1,nCols
       !     k = kbot(i) 
       !     xx1(i,k)  =  qm(i,k) * ( phalf(i,k+1) - pfull(i,k) )
       !     xx2(i,k)  = xx1(i,k) * z0(i)
       !   end do
       !else
       DO i = 1,nCols
!          xx1(i,kMax) =  qm(i,kMax) * ( MAX(phalf(i,kMax+1) - phalf(i,kMax),1.0_r8) )
          xx1(i,kMax) =  qm(i,kMax) * ( MAX(phalf(i,kMax+1) - phalf(i,kMax),0.00000000000001_r8) )
          xx2(i,kMax) = xx1(i,kMax) * z0(i)
       END DO
       !end if

       !  if (PRESENT(mask)) then
       !    x1 = SUM( xx1, 3, mask=mask.gt.0.1_r8 )
       !    x2 = SUM( xx2, 3, mask=mask.gt.0.1_r8 )
       !  else
       DO i = 1,nCols
          x1(i) = SUM( xx1(i,1:kMax))
          x2(i) = SUM( xx2(i,1:kMax))
       END DO

       !  endif

       !---- should never be equal to zero ----

       IF (COUNT(x1 <= 0.0_r8) > 0) THEN
          WRITE(*,*)' MY25_TURB','divid by zero, x1 <= 0.0_r8'
          STOP
       END IF

       DO i = 1,nCols
          el0(i) = x2(i) / x1(i)
          el0(i) = el0(i) * (alpha_land*fracland (i)+ alpha_sea*(1.0_r8-fracland(i)))

          el0(i) = MIN( el0(i), el0max1 )
          el0(i) = MAX( el0(i), el0min1 )
       END DO
       !====================================================================
       ! --- MASTER LENGTH SCALE 
       !====================================================================

       DO k = 1, kmax-1
          DO i = 1,nCols
             xx1(i,k)  = vonkarm * ( zhalf(i,k+1) - zsfc(i) )
          END DO
       END DO
       DO i = 1,nCols
          x1(i) = vonkarm * z0(i) 
       END DO

       !if( PRESENT( kbot ) ) then
       !    do i = 1,nCols
       !      do k = kbot(i), kMax
       !        xx1(i,j,k) = x1(i,j)
       !      end do
       !    end do
       !else
       DO i = 1,nCols
          xx1(i,kMax) = x1(i)
       END DO
       !endif 

       DO k = 1,kMax
          DO i = 1,nCols
             el(i,k+1) = xx1(i,k) / ( 1.0_r8 + xx1(i,k) / el0(i) )
          END DO
       END DO
       DO i = 1,nCols
          el(i,1)   = el0(i)
       END DO
       DO k = 1,kMax-1
          DO i = 1,nCols
             el2(i,k) = el(i,k+1) * el(i,k+1) 
          END DO
       END DO
       !====================================================================
       ! --- MIXING COEFFICIENTS                     
       !====================================================================
       DO k = 1,kMax-1    
          DO i = 1,nCols
             xxm3(i,k) = el2(i,k)* buoync(i,k)
             xxm4(i,k) = el2(i,k)* shear(i,k)
             xxm5(i,k) = el (i,k+1)*qm3(i,k)
          END DO
       END DO
       !-------------------------------------------------------------------
       ! --- MOMENTUM 
       !-------------------------------------------------------------------
       DO k = 1,kMax-1
          DO i = 1,nCols
             xxm1(i,k) = xxm5(i,k)*( ckm1*qm2(i,k) + ckm2*xxm3(i,k) )
             xxm2(i,k) = qm4(i,k) + ckm5*qm2(i,k)*xxm4(i,k) + xxm3(i,k)*( ckm6*xxm4(i,k) + ckm7*qm2(i,k) + ckm8*xxm3(i,k) )

             xxm2(i,k) = MAX( xxm2(i,k), 0.2_r8*qm4(i,k) )
             xxm2(i,k) = MAX( xxm2(i,k), small2  )

             akm(i,1)    = 0.0_r8
             akm(i,k+1)  = xxm1(i,k) / xxm2(i,k)

             akm(i,k+1) = MAX( akm(i,k+1), 0.0_r8 )

             !-------------------------------------------------------------------
             ! --- HEAT AND MOISTURE 
             !-------------------------------------------------------------------

             xxm1(i,k) = ckh1*xxm5(i,k) - ckh2*xxm4(i,k)*akm(i,k+1)
             xxm2(i,k) = qm2      (i,k) + ckh3*xxm3(i,k)

             xxm1(i,k) = MAX( xxm1(i,k), ckh4*xxm5(i,k) )
             xxm2(i,k) = MAX( xxm2(i,k),  0.4_r8*qm2(i,k)   )
             xxm2(i,k) = MAX( xxm2(i,k), small2     )

             akh(i,1)    = 0.0_r8
             akh(i,k+1) = xxm1(i,k) / xxm2(i,k)

             akh(i,k+1) = MAX( akh(i,k+1), 0.0_r8 )
             !-------------------------------------------------------------------
             ! --- BOUNDS 
             !-------------------------------------------------------------------

             ! --- UPPER BOUND
             !akm(i,k) = MIN( akm(i,k), akmax )
             !akh(i,k) = MIN( akh(i,k), akmax )
          END DO
       END DO
       DO k = 1,kMax
          DO i = 1,nCols
             !-------------------------------------------------------------------
             ! --- BOUNDS 
             !-------------------------------------------------------------------

             ! --- UPPER BOUND
             akm(i,k) = MIN( akm(i,k), akmax )
             akh(i,k) = MIN( akh(i,k), akmax )
          END DO
       END DO

       ! --- LOWER BOUND 
       !klim = kMax - nk_lim + 1
       !DO  k = klim,kMax+1
       !   DO i = 1,nCols
       !      IF(akm(i,k) < small2 )  akm(i,k) = 0.0_r8
       !      IF(akh(i,k) < small2 )  akh(i,k) = 0.0_r8
       !   END DO
       !END DO

       ! --- LOWER BOUND NEAR SURFACE
       DO i = 1,nCols
          akmin(i) = akmin_land*fracland(i) + akmin_sea*(1.0_r8-fracland(i))
       END DO
       !DO i = 1,nCols
       !      !
       !      !       Km(1) + Km(2)
       !      ! Km = ---------------
       !      !            2
       !      !
       !      !
       !      !       Kh(1) + Kh(2)
       !      ! Kh = ---------------
       !      !            2
       !      !
       !      akm(i,     1)=0.5_r8*(akm(i,     1)+akm(i,     2))
       !      akh(i,     1)=0.5_r8*(akh(i,     1)+akh(i,     2))
       !      akm(i,kMax-1)=0.5_r8*(akm(i,kMax-1)+akm(i,kMax-2))
       !      akh(i,kMax-1)=0.5_r8*(akh(i,kMax-1)+akh(i,kMax-2))
       !      akm(i,kMax  )=0.5_r8*(akm(i,kMax  )+akm(i,kMax-1))
       !      akh(i,kMax  )=0.5_r8*(akh(i,kMax  )+akh(i,kMax-1))
       !      akm(i,kMax+1)=0.5_r8*(akm(i,kMax+1)+akm(i,kMax))
       !      akh(i,kMax+1)=0.5_r8*(akh(i,kMax+1)+akh(i,kMax))
       !END DO

       !if( PRESENT( kbot ) ) then
       !      do i = 1,nCols
       !           klim = kbot(i) - nk_lim + 1
       !         do  k = klim,kbot(i)
       !             akm(i,k) = MAX( akm(i,k), akmin(i) )
       !             akh(i,k) = MAX( akh(i,k), akmin(i) )
       !         end do
       !      end do
       !else
       klim = kMax - nk_lim + 1
       DO  k = klim,kMax
          DO i = 1,nCols
             akm(i,k) = MAX( akm(i,k), akmin(i) )
             akh(i,k) = MAX( akh(i,k), akmin(i) )
          END DO
       END DO
       !endif

       !-------------------------------------------------------------------
       ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE
       !-------------------------------------------------------------------

       !if( PRESENT( mask ) ) then
       !    akm(:,:,1:kx) = akm(:,:,1:kx) * mask(:,:,1:kx) 
       !    akh(:,:,1:kx) = akh(:,:,1:kx) * mask(:,:,1:kx) 
       !endif

       !====================================================================
       ! --- PROGNOSTICATE TURBULENT KE 
       !====================================================================

       cvfqdt = cvfq1 * delt
       dvfqdt = cvfq2 * delt * 2.0_r8

       !-------------------------------------------------------------------
       ! --- PART OF LINEARIZED ENERGY DISIIPATION TERM 
       !-------------------------------------------------------------------
       DO k = 1,kMax-1
          DO i = 1,nCols
             xxm1(i,k) = dvfqdt * qm(i,k) / el(i,k+1)
          END DO
       END DO

       !-------------------------------------------------------------------
       ! --- PART OF LINEARIZED VERTICAL DIFFUSION TERM
       !-------------------------------------------------------------------
       DO k = 1,kMax
          DO i = 1,nCols

             xx1(i,k ) = el(i,k+1) * qm(i,k)
          END DO
       END DO
       DO i = 1,nCols
          xx2(i,1)    = 0.5_r8*  xx1(i,1)
       END DO
       DO k = 1,kMax-1
          DO i = 1,nCols
             xx2(i,k+1 ) = 0.5_r8*( xx1(i,k+1 ) + xx1(i,k) )
          END DO
       END DO
       DO k = 1,kMax
          DO i = 1,nCols
             xx1 (i,k )= xx2(i,k ) * dsdz(i,k )
          END DO
       END DO

       !-------------------------------------------------------------------
       ! --- IMPLICIT TIME DIFFERENCING FOR VERTICAL DIFFUSION 
       ! --- AND ENERGY DISSIPATION TERM 
       !-------------------------------------------------------------------

       DO k=1,kMax-1
          DO i=1,nCols
             aaa(i,k) = -cvfqdt * xx1(i,k+1) * dsdzh(i,k)
             ccc(i,k) = -cvfqdt * xx1(i,k  ) * dsdzh(i,k)
             bbb(i,k) =     1.0_r8 - aaa(i,k  ) -   ccc(i,k) 
             bbb(i,k) =             bbb(i,k  ) +  xxm1(i,k)
             ddd(i,k) =             TKE(i,k+1)
          END DO
       END DO

       ! correction for vertical diffusion of TKE surface boundary condition

       !if (present(kbot)) then
       !   do i = 1,nCols
       !        k = kbot(i)
       !        ddd(i,k-1) = ddd(i,k-1) - aaa(i,k-1) * TKE(i,latco,k+1)
       !   end do
       !else
       DO i = 1,nCols
          ddd(i,kMax-1) = ddd(i,kMax-1) - aaa(i,kMax-1) * TKE(i,kMax+1)
       END DO
       !endif

       ! mask out terms below ground

       !  if (present(mask)) then
       !     where (mask(:,:,2:kx) < 0.1_r8) ddd(:,:,1:kxm) = 0.0_r8
       !  endif


       CALL TRI_INVERT( &
            nCols,  & 
            kmax-1, & 
            xxm1  (1:nCols,1:kmax-1)     , & ! intent(out)  :: x (:,:)
            ddd   (1:nCols,1:kmax-1)     , & ! intent(in)   :: d (:,:)
            aaa   (1:nCols,1:kmax-1)     , & ! optional     :: a (:,:)
            bbb   (1:nCols,1:kmax-1)     , & ! optional     :: b (:,:)
            ccc   (1:nCols,1:kmax-1)       ) ! optional     :: c (:,:)


       !-------------------------------------------------------------------
       ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE
       !-------------------------------------------------------------------

       ! if( PRESENT( mask ) ) then
       !   do k=1,kxm
       !   do j=1,jx
       !   do i=1,ix
       !     if(mask(i,j,k+1) < 0.1_r8) xxm1(i,j,k) = TKE(ism+i,jsm+j,k+1)
       !   enddo
       !   enddo
       !   enddo
       ! endif

       !-------------------------------------------------------------------
       ! --- SHEAR AND BUOYANCY TERMS
       !-------------------------------------------------------------------
       DO k=1,kMax-1
          DO i=1,nCols

             xxm2(i,k) =  delt*( akm(i, k+1 )* shear (i,k)    &
                  - akh(i,k+1  )* buoync(i,k ) )
          END DO
       END DO

       !-------------------------------------------------------------------
       ! --- UPDATE TURBULENT KINETIC ENERGY
       !-------------------------------------------------------------------

       DO i=1,nCols
          TKE(i,1) = 0.0_r8
       ENDDO

       DO k=2,kMax
          DO i=1,nCols
             TKE(i,k) = xxm1(i,k-1) + xxm2(i,k-1)
          ENDDO
       ENDDO

       !====================================================================
       ! --- BOUND TURBULENT KINETIC ENERGY
       !====================================================================
       DO k=1,kMax+1
          DO i=1,nCols
             TKE(i,k) = MIN( TKE(i,k), TKEmax )
             TKE(i,k) = MAX( TKE(i,k), TKEmin )
          ENDDO
       ENDDO     

       !DO k=1,kMax+1
       !  DO i=1,nCols
             !IF(TKE(i,k) == TKEmin) TKE(i,k)=0.0_r8
       !   ENDDO
       !ENDDO     

       !  if( PRESENT( mask ) ) then
       !    do k=1,kx
       !    do j=1,jx
       !    do i=1,ix
       !      if(mask(i,j,k) < 0.1_r8) TKE(ism+i,jsm+j,k+1) = 0.0_r8
       !    enddo
       !    enddo
       !    enddo
       !  endif

       !====================================================================
       ! --- COMPUTE PBL DEPTH IF DESIRED
       !====================================================================

       !  if (present(h)) then

       !      if (.not.present(ustar)) then
       !          CALL ERROR_MESG( ' MY25_TURB',     &
       !              'cannot request pbl depth diagnostic if ustar'// &
       !              ' and bstar are not also supplied', FATAL )
       !      end if
       !if (present(kbot)) then
       !    call k_pbl_depth(ustar (1:nCols)                , &
       !                           akm   (1:nCols,1:kMax+1)       , &
       !                           akh   (1:nCols,1:kMax+1)       , &
       !                           zsfc  (1:nCols)                , &
       !                           zfull (1:nCols,1:kMax)         , &
       !                           zhalf (1:nCols,1:kMax+1)       , &
       !                           h     (1:nCols)                , &
       !                           kbot(1:nCols)                    )
       !      else
       !          call k_pbl_depth(ustar (1:nCols)                  , &
       !                           akm   (1:nCols,1:kMax+1)          , &
       !                           akh   (1:nCols,1:kMax+1)          , &
       !                           zsfc  (1:nCols)                  , &
       !                           zfull (1:nCols,1:kMax)          , &
       !                           zhalf (1:nCols,1:kMax+1)          , &
       !                           h     (1:nCols)                    )
       !      end if

       !  end if
       !====================================================================

       ! $$$$$$$$$$$$$$$$$
    END DO
    ! $$$$$$$$$$$$$$$$$

  !#######################################################################
  END SUBROUTINE MY25_TURB

  SUBROUTINE TKE_SURF ( nCols,kMax,u_star,TKE )

    !=======================================================================
    !---------------------------------------------------------------------
    ! Arguments (Intent in)
    !       u_star -  surface friction velocity (m/s)
    !       kbot   -  OPTIONAL;lowest model level index (integer);
    !                 at levels > Kbot, Mask = 0.
    !---------------------------------------------------------------------
    INTEGER , INTENT(in)   :: nCols
    INTEGER , INTENT(in)   :: kMax

    REAL(KIND=r8), INTENT(in)   :: u_star(nCols) 
    REAL(KIND=r8), INTENT(inout)  :: TKE   (nCols,kMax+1) 
    !---------------------------------------------------------------------
    !  (Intent local)
    !---------------------------------------------------------------------
    REAL(KIND=r8) :: x1(nCols)
    INTEGER  :: ix, kxp, i

    !=======================================================================

    ix  = nCols
    kxp = kMax+1
    !TKE=0.0_r8
    !---------------------------------------------------------------------
    DO i = 1,ix
       ! bb1 = 16.0_r8
       ! bcq   = 0.5 * ( 16.0_r8**(2.0/3.0) )
       ! tk2(i) = ( 0.5 * ( 16.0_r8**(2.0/3.0) )) * u_star(i) * u_star(i)
       x1(i) = MAX(bcq * u_star(i) * u_star(i),1.0e-12_r8 )
       TKE(i,kxp) = MAX(x1(i),TKE(i,kxp)) 
    END DO

    !=======================================================================
  END SUBROUTINE TKE_SURF


  SUBROUTINE tri_invert(&
       nCols   , &
       kMax    , &
       x       , &  ! intent(out)  :: x (:,:)
       d       , &  ! intent(in)   :: d (:,:)
       a       , &  ! optional          :: a (:,:)
       b       , &  ! optional          :: b (:,:)
       c         )  ! optional          :: c (:,:)

    IMPLICIT NONE
    INTEGER, INTENT(IN   ) :: nCols
    INTEGER, INTENT(IN   ) :: kMax
    REAL(KIND=r8), INTENT(out)  :: x  (nCols,kMax)
    REAL(KIND=r8), INTENT(in)   :: d  (nCols,kMax)
    REAL(KIND=r8), OPTIONAL     :: a  (nCols,kMax)
    REAL(KIND=r8), OPTIONAL     :: b  (nCols,kMax)
    REAL(KIND=r8), OPTIONAL     :: c  (nCols,kMax)

    REAL(KIND=r8)  :: f (nCols,kMax)

    REAL(KIND=r8)  :: e (nCols,kMax)
    REAL(KIND=r8)  :: g (nCols,kMax)
    REAL(KIND=r8)  :: bb(nCols)
    REAL(KIND=r8)  :: cc(nCols,kMax)

    INTEGER :: k
    INTEGER :: i
    
    f =0.0_r8;e=0.0_r8;g=0.0_r8;bb=0.0_r8;cc=0.0_r8;x=0.0_r8
    
    DO i=1,nCols
       e(i,1) = - a(i,1)/b(i,1)
       a(i,kMax) = 0.0_r8
    END DO

    DO  k= 2,kMax
       DO i=1,nCols
          g(i,k) = 1/(b(i,k)+c(i,k)*e(i,k-1))
          e(i,k) = - a(i,k)*g(i,k)
       END DO
    END DO
    DO  k= 1,kMax
       DO i=1,nCols
          cc(i,k) = c(i,k)
       END DO
    END DO
    DO i=1,nCols
       bb(i) = 1.0_r8/b(i,1)
    END DO

    ! if(.not.init_tridiagonal) error
    DO i=1,nCols
       f(i,1) =  d(i,1)*bb(i)  
    END DO
    DO k= 2, kMax
       DO i=1,nCols
          f(i,k) = (d(i,k) - cc(i,k)*f(i,k-1))*g(i,k)
       END DO
    END DO
    DO i=1,nCols
       x(i,kMax) = f(i,kMax)
    END DO

    DO k = kMax-1,1,-1
       DO i=1,nCols
          !IF(ABS(e(i,k)) < 1.0e-12_r8)PRINT*,'e(i,k)=',e(i,k)
          !IF(ABS(x(i,k)) < 1.0e-12_r8)PRINT*,'x(i,k)=',x(i,k+1)
          !IF(ABS(f(i,k)) < 1.0e-12_r8)PRINT*,'f(i,k)=',f(i,k)
          x(i,k) = e(i,k)*x(i,k+1)+f(i,k)
       END DO
    END DO

    RETURN
  END SUBROUTINE tri_invert



   END MODULE CCM3PBL
