bldriver.F90

Go to the documentation of this file.

MODULE bldriver

!This subroutine takes the model prognostic variables as its inputs and produces
!the three components of the tri-diagonal matrix used to compute BL diffusion. 
!Diffusion coefficient is computed using Louis-Lock formulation.

USE fv3jedi_lm_const_mod
use fv3jedi_lm_kinds_mod

IMPLICIT NONE

PRIVATE
PUBLIC bl_driver

!Saturation vapor pressure table ESTBLX initialization.
!To be defined for entire module:
integer, parameter :: degsubs    =  100
real(kind_real), parameter :: tmintbl    =  150.0, tmaxtbl = 333.0
integer, parameter :: tablesize  =  nint(tmaxtbl-tmintbl)*degsubs + 1

CONTAINS

subroutine bl_driver ( IM, JM, LM, DT, U, V, TH, Q, P, QIT, QLT, FRLAND, FROCEAN, VARFLT, &
                       ZPBL, CM, CT, CQ, TURBPARAMS, TURBPARAMSI, &
                       USTAR, BSTAR, &
                       AKS, BKS, CKS, AKQ, BKQ, CKQ, AKV, BKV, CKV, EKV, FKV &
                     )

IMPLICIT NONE

!INPUTS
INTEGER, INTENT(IN) :: im, jm, lm
REAL(kind_real), INTENT(IN) :: dt

REAL(kind_real), INTENT(IN), DIMENSION(:) :: turbparams
INTEGER, INTENT(IN), DIMENSION(:) :: turbparamsi

REAL(kind_real), INTENT(IN), DIMENSION(IM,JM,LM) :: u, v, th, q, qit, qlt
REAL(kind_real), INTENT(IN), DIMENSION(IM,JM,0:LM) :: p
REAL(kind_real), INTENT(IN), DIMENSION(IM,JM) :: frland, varflt, frocean

REAL(kind_real), INTENT(IN), DIMENSION(IM,JM) :: cm, cq

REAL(kind_real), INTENT(IN), DIMENSION(IM,JM)    :: ustar, bstar

!INOUTS
REAL(kind_real), INTENT(INOUT), DIMENSION(IM,JM) :: zpbl, ct

!OUTPUTS
REAL(kind_real), INTENT(OUT), DIMENSION(IM,JM,LM) :: aks, bks, cks
REAL(kind_real), INTENT(OUT), DIMENSION(IM,JM,LM) :: akq, bkq, ckq
REAL(kind_real), INTENT(OUT), DIMENSION(IM,JM,LM) :: akv, bkv, ckv, ekv, fkv

!LOCALS
INTEGER :: i,j
REAL(kind_real), DIMENSION(IM,JM,LM) :: pf, pif, qi, ql, tv, thv, zf, dmi
REAL(kind_real), DIMENSION(IM,JM,0:LM) :: pi, z
REAL(kind_real), DIMENSION(IM,JM,1:LM-1) :: rdz
REAL(kind_real), DIMENSION(IM,JM,LM) :: t
REAL(kind_real), DIMENSION(IM,JM,0:LM) :: km, kh

REAL(kind_real), DIMENSION(IM,JM,0:LM) :: ri, du
REAL(kind_real), DIMENSION(IM,JM,0:LM) :: alh_x, kmls_x, khls_x

REAL(kind_real), DIMENSION(IM,JM,LM) :: radlw !dh Note that this is really an input but we would rather not put it in the traj
                                   !so for now do not use. Later on add approximate numbers.

!Turb constants
REAL(kind_real) :: louis,lambdam,lambdam2,lambdah,lambdah2,zkmenv,zkhenv,minthick,minshear,c_b,lambda_b,akhmmax
REAL(kind_real) :: prandtlsfc,prandtlrad,beta_rad,beta_surf,khradfac,khsfcfac,tpfac_surf,entrate_surf,pceff_surf,louis_memory
INTEGER :: kpblmin, lock_on, pblht_option, radlw_dep

!QSATVP TABLE
integer, parameter :: degsubs    =  100
real(kind_real), parameter :: tmintbl    =  150.0, tmaxtbl = 333.0
integer, parameter :: tablesize  =  nint(tmaxtbl-tmintbl)*degsubs + 1
real(kind_real), dimension(TABLESIZE) :: estblx

 !Compute saturation vapour pressure lookup table
 call esinit(estblx)

 !Initialize outputs
 aks = 0.0
 bks = 0.0
 cks = 0.0
 akq = 0.0
 bkq = 0.0
 ckq = 0.0
 akv = 0.0
 bkv = 0.0
 ckv = 0.0
 ekv = 0.0
 fkv = 0.0

 !Caclulate temperature from Exner pressure.
 pf = 0.5*(p(:,:,0:lm-1) + p(:,:,1:lm))
 pif = (pf/p00)**(rgas/cp)
 t = pif*th

 !TurbParams, this goes at the level above.
 !TurbParams(1)  = 5.0          !LOUIS
 !TurbParams(2)  = 160.0        !LAMBDAM
 !TurbParams(3)  = 1.0          !LAMBDAM2
 !TurbParams(4)  = 160.0        !LAMBDAH
 !TurbParams(5)  = 1.0          !LAMBDAH2
 !TurbParams(6)  = 3000.        !ZKMENV
 !TurbParams(7)  = 3000.        !ZKHENV
 !TurbParams(8)  = 0.1          !MINTHICK
 !TurbParams(9)  = 0.0030       !MINSHEAR
 !TurbParams(10) = 2.5101471e-8 !C_B
 !TurbParams(11) = 1500.        !LAMBDA_B 
 !TurbParams(12) = 500.         !AKHMMAX
 !TurbParams(13) = 1.0          !PRANDTLSFC
 !TurbParams(14) = 0.75         !PRANDTLRAD
 !TurbParams(15) = 0.50         !BETA_RAD
 !TurbParams(16) = 0.25         !BETA_SURF
 !TurbParams(17) = 0.85         !KHRADFAC
 !TurbParams(18) = 0.45         !KHSFCFAC
 !TurbParams(19) = 20.0         !TPFAC_SURF
 !TurbParams(20) = 1.5e-3       !ENTRATE_SURF
 !TurbParams(21) = 0.5          !PCEFF_SURF
 !TurbParams(22) = -999.        !LOUIS_MEMORY

 !TurbParamsI(1) = count(PREF < 50000.)   !KPBLMIN
 !TurbParamsI(2) = 1                      !LOCK_ON
 !TurbParamsI(3) = 1                      !PBLHT_OPTION 
 !TurbParamsI(4) = 0                      !RADLW_DEP

 !Turbulence constants
 louis        = turbparams(1)
 lambdam      = turbparams(2)
 lambdam2     = turbparams(3) 
 lambdah      = turbparams(4)
 lambdah2     = turbparams(5)
 zkmenv       = turbparams(6)
 zkhenv       = turbparams(7)
 minthick     = turbparams(8)
 minshear     = turbparams(9)
 c_b          = turbparams(10)
 lambda_b     = turbparams(11)
 akhmmax      = turbparams(12)
 prandtlsfc   = turbparams(13)
 prandtlrad   = turbparams(14)
 beta_rad     = turbparams(15)
 beta_surf    = turbparams(16)
 khradfac     = turbparams(17)
 khsfcfac     = turbparams(18)
 tpfac_surf   = turbparams(19)
 entrate_surf = turbparams(20)
 pceff_surf   = turbparams(21)
 louis_memory = turbparams(22)

 kpblmin      = turbparamsi(1)
 lock_on      = turbparamsi(2)
 pblht_option = turbparamsi(3)
 radlw_dep    = turbparamsi(4)

 !Initialise arrays
 kh    = 0.0
 km    = 0.0
 ri    = 0.0
 du    = 0.0

 !If over the ocean then Ts' = 0. That is the perturbation of surface temperature
 !is negligible. If over land then it is not. But temperature transfer is quick so
 !T1'-Ts' = 0. So set CT to 0 over land and 1 over ocean. 
 !In between use fraction to take percentage of increment.
 DO i = 1,im
    DO j = 1,jm
       IF (frocean(i,j).eq.1.0) then
          ct(i,j) = frocean(i,j) * ct(i,j)
       endif
    enddo
 enddo

 !Call the subroutine to get levels, RDZ and DMI
 call preliminary( im*jm      , & !IN
                   lm         , & !IN
                   dt         , & !IN
                   t          , & !IN
                   q          , & !IN
                   p          , & !IN
                   th         , & !IN
                   qit        , & !IN
                   qlt        , & !IN
                   zf         , & !OUT
                   z          , & !OUT
                   tv         , & !OUT
                   thv        , & !OUT
                   rdz        , & !OUT
                   dmi        , & !OUT
                   pf           & !OUT
                 )

 !GET LOUIS DIFFUSIVITIES
 call louis_diff( im*jm       , & !IN  
                  lm          , & !IN
                  kh          , & !OUT
                  km          , & !OUT
                  ri          , & !OUT
                  du          , & !OUT
                  zpbl        , & !IN
                  zf          , & !IN
                  z           , & !IN
                  thv         , & !IN
                  u           , & !IN
                  v           , & !IN
                  louis       , & !IN
                  minshear    , & !IN
                  minthick    , & !IN
                  lambdam     , & !IN
                  lambdam2    , & !IN
                  lambdah     , & !IN
                  lambdah2    , & !IN
                  zkmenv      , & !IN
                  zkhenv      , & !IN
                  akhmmax     , & !IN
                  alh_x       , & !OUT
                  kmls_x      , & !OUT
                  khls_x        & !OUT
                )

 !Optional - overwrite diffusivities using the Lock et al scheme.
 CALL lock_diff( im*jm        , & !IN
                 lm           , & !IN
                 radlw        , & !IN
                 ustar        , & !IN
                 bstar        , & !IN
                 frland       , & !IN
                 t            , & !IN
                 q            , & !IN
                 qit          , & !IN
                 qlt          , & !IN
                 u            , & !IN
                 v            , & !IN
                 zf           , & !IN
                 pf           , & !IN
                 z            , & !IN
                 p            , & !IN
                 km           , & !INOUT
                 kh           , & !INOUT
                 prandtlsfc   , & !IN
                 prandtlrad   , & !IN
                 beta_surf    , & !IN
                 beta_rad     , & !IN
                 tpfac_surf   , & !IN
                 entrate_surf , & !IN
                 pceff_surf   , & !IN
                 khradfac     , & !IN
                 khsfcfac     , & !IN
                 radlw_dep    , & !IN
                 estblx         & !IN
               )
 
 call tridiag_setup( im*jm          , & !IN
                     lm             , & !IN
                     dt             , & !IN
                     kpblmin        , & !IN
                     zf             , & !IN
                     pf             , & !IN
                     rdz            , & !IN
                     dmi            , & !IN
                     p              , & !IN
                     tv             , & !IN
                     ct             , & !IN
                     cq             , & !IN
                     cm             , & !IN
                     t              , & !IN
                     q              , & !IN
                     th             , & !IN
                     u              , & !IN
                     v              , & !IN
                     kh             , & !INOUT
                     km             , & !INOUT
                     akq, aks, akv  , & !OUT
                     bkq, bks, bkv  , & !OUT
                     ckq, cks, ckv  , & !OUT
                               ekv  , & !OUT
                     zpbl             & !OUT
                   )

 call orodrag( im*jm          , & !IN
               lm             , & !IN
               dt             , & !IN
               lambda_b       , & !IN
               c_b            , & !IN
               fkv            , & !OUT
               bkv            , & !INOUT
               u              , & !IN
               v              , & !IN
               zf             , & !IN
               varflt         , & !IN
               p                & !IN
             )


end subroutine bl_driver


subroutine preliminary( IRUN, LM, DT, T, QV, PHALF, TH, QIT, QLT, &
                        ZFULL, ZHALF, TV, PV, RDZ, DMI, PFULL &
                      )

 !Subroutine to set up priliminary

 IMPLICIT NONE

 !INTPUTS
 integer, intent(IN)                    :: IRUN, LM
 real(kind_real), intent(IN)                       :: DT
 real(kind_real), intent(IN), dimension(IRUN,LM)   :: T, QV, TH, QIT, QLT
 real(kind_real), intent(IN), dimension(IRUN,LM+1) :: PHALF

 !OUTPUTS
 real(kind_real), intent(OUT), dimension(IRUN,  LM  ) :: ZFULL, TV, PV, DMI, PFULL
 real(kind_real), intent(OUT), dimension(IRUN,  LM+1) :: ZHALF
 real(kind_real), intent(OUT), dimension(IRUN,1:LM-1) :: RDZ

 !LOCALS
 integer :: I,L
 real(kind_real) :: PKE_atL, PKE_atLp1, TVE
 real(kind_real) :: QL_tot, QI_tot 

 do i = 1, irun

    !Compute the edge heights using Arakawa-Suarez hydrostatic equation
    zhalf(i,lm+1) = 0.0
    do l = lm, 1, -1
       pke_atlp1  = (phalf(i,l+1)/p00)**kappa
       pke_atl    = (phalf(i,l  )/p00)**kappa

       zhalf(i,l) = zhalf(i,l+1) + (cp/grav)*th(i,l)*(pke_atlp1-pke_atl)
    end do

    !Layer height, pressure, and virtual temperatures
    do l = 1, lm
       ql_tot = qit(i,l)
       qi_tot = qlt(i,l)

       zfull(i,l) = 0.5*(zhalf(i,l)+zhalf(i,l+1))
       pfull(i,l) = 0.5*(phalf(i,l)+phalf(i,l+1))

       tv(i,l)  = t(i,l) *( 1.0 + vireps *qv(i,l) - ql_tot - qi_tot ) 
       pv(i,l) = tv(i,l)*(th(i,l)/t(i,l))
    end do

    do l = 1, lm-1
       tve = (tv(i,l) + tv(i,l+1))*0.5
  
       ! Miscellaneous factors
       rdz(i,l) = phalf(i,l+1) / ( rgas * tve )
       rdz(i,l) = rdz(i,l) / (zfull(i,l)-zfull(i,l+1))
    end do

    do l = 1, lm
       dmi(i,l) = (grav*dt)/(phalf(i,l+1)-phalf(i,l))
    end do

    !Running 1-2-1 smooth of bottom 5 levels of Virtual Pot. Temp.
    pv(i,lm  ) = pv(i,lm-1)*0.25 + pv(i,lm  )*0.75
    pv(i,lm-1) = pv(i,lm-2)*0.25 + pv(i,lm-1)*0.50 + pv(i,lm  )*0.25 
    pv(i,lm-2) = pv(i,lm-3)*0.25 + pv(i,lm-2)*0.50 + pv(i,lm-1)*0.25 
    pv(i,lm-3) = pv(i,lm-4)*0.25 + pv(i,lm-3)*0.50 + pv(i,lm-2)*0.25 
    pv(i,lm-4) = pv(i,lm-5)*0.25 + pv(i,lm-4)*0.50 + pv(i,lm-3)*0.25 
    pv(i,lm-5) = pv(i,lm-6)*0.25 + pv(i,lm-5)*0.50 + pv(i,lm-4)*0.25 

 end do 

end subroutine preliminary



subroutine louis_diff( IRUN,LM,KH,KM,RI,DU,ZPBL,ZZ,ZE,PV,UU,VV,LOUIS,MINSHEAR,MINTHICK,&
                       LAMBDAM,LAMBDAM2,LAMBDAH,LAMBDAH2,ZKMENV,ZKHENV,AKHMMAX,ALH_DIAG,KMLS_DIAG,KHLS_DIAG)

 !Subroutine to compute eddy diffussivity using the Louis (1979) method.
 !Simplifications are made compared with the NL stability functions f(Ri). 
 !Gradient is very steep for near neutral boundary layer conditions.

 IMPLICIT NONE

 !INPUTS
 integer, intent(IN   ) :: IRUN,LM
 real(kind_real), intent(IN   ) :: LOUIS         ! Louis scheme parameters (usually 5).
 real(kind_real), intent(IN   ) :: MINSHEAR      ! Min shear allowed in Ri calculation (s-1).
 real(kind_real), intent(IN   ) :: MINTHICK      ! Min layer thickness (m).
 real(kind_real), intent(IN   ) :: AKHMMAX       ! Maximum allowe diffusivity (m+2 s-1).
 real(kind_real), intent(IN   ) :: LAMBDAM       ! Blackadar(1962) length scale parameter for momentum (m).
 real(kind_real), intent(IN   ) :: LAMBDAM2      ! Second Blackadar parameter for momentum (m).
 real(kind_real), intent(IN   ) :: LAMBDAH       ! Blackadar(1962) length scale parameter for heat (m).
 real(kind_real), intent(IN   ) :: LAMBDAH2      ! Second Blackadar parameter for heat (m).
 real(kind_real), intent(IN   ) :: ZKMENV        ! Transition height for Blackadar param for momentum (m)
 real(kind_real), intent(IN   ) :: ZKHENV        ! Transition height for Blackadar param for heat     (m)
 real(kind_real), intent(IN   ) :: ZZ(IRUN,LM)   ! Height of layer center above the surface (m).
 real(kind_real), intent(IN   ) :: PV(IRUN,LM)   ! Virtual potential temperature at layer center (K).
 real(kind_real), intent(IN   ) :: UU(IRUN,LM)   ! Eastward velocity at layer center (m s-1).
 real(kind_real), intent(IN   ) :: VV(IRUN,LM)   ! Northward velocity at layer center (m s-1).
 real(kind_real), intent(IN   ) :: ZE(IRUN,LM+1) ! Height of layer base above the surface (m).
 real(kind_real), intent(IN   ) :: ZPBL(IRUN)    ! PBL Depth (m)

 !OUTPUTS
 !These are 1:LM+1 here but 0:LM in the GC - MAYBE JUST SORT THIS OUT?
 !Old code only passed in 1:LM-1 from GC which is 2:LM here.
 real(kind_real), intent(  OUT) :: KM(IRUN,LM+1) ! Momentum diffusivity at base of each layer (m+2 s-1).
 real(kind_real), intent(  OUT) :: KH(IRUN,LM+1) ! Heat diffusivity at base of each layer  (m+2 s-1).
 real(kind_real), intent(  OUT) :: RI(IRUN,LM+1) ! Richardson number
 real(kind_real), intent(  OUT) :: DU(IRUN,LM+1) ! Magnitude of wind shear (s-1).
 real(kind_real), intent(  OUT) :: ALH_DIAG(IRUN,LM+1)  ! Blackadar Length Scale diagnostic (m) [Optional] 
 real(kind_real), intent(  OUT) :: KMLS_DIAG(IRUN,LM+1) ! Momentum diffusivity at base of each layer (m+2 s-1).
 real(kind_real), intent(  OUT) :: KHLS_DIAG(IRUN,LM+1) ! Heat diffusivity at base of each layer  (m+2 s-1).

 !LOCALS
 integer :: L,Levs,i,j,k
 real(kind_real)    :: ALH, ALM, DZ, DT, TM, PS, LAMBDAM_X, LAMBDAH_X
 real(kind_real)    :: pbllocal, alhfac,almfac

 almfac = 1.2
 alhfac = 1.2

 do i = 1, irun !Loop of columns

    alh_diag(i,   1) = 0.0
    alh_diag(i,lm+1) = 0.0

    kh(i,   1) = 0.0
    kh(i,lm+1) = 0.0
    km(i,   1) = 0.0
    km(i,lm+1) = 0.0
    du(i,   1) = 0.0
    du(i,lm+1) = 0.0
    ri(i,   1) = 0.0
    ri(i,lm+1) = 0.0
    kmls_diag(i,   1) = 0.0
    kmls_diag(i,lm+1) = 0.0
    khls_diag(i,   1) = 0.0
    khls_diag(i,lm+1) = 0.0

    pbllocal = zpbl(i)

    if ( pbllocal .LE. zz(i,lm) ) then
       pbllocal = zz(i,lm)
    endif

    do k = 2, lm !This is equivalent to 1->LM-1 but we move all 0:72 elements
       kh(i,k) = 0.0
       km(i,k) = 0.0
       du(i,k) = 0.0
       ri(i,k) = 0.0

       dz      = (zz(i,k-1) - zz(i,k))
       tm      = (pv(i,k-1) + pv(i,k))*0.5
       dt      = (pv(i,k-1) - pv(i,k))
       du(i,k) = (uu(i,k-1) - uu(i,k))**2 + (vv(i,k-1) - vv(i,k))**2

       !Limits distance between layer centers and vertical shear at edges.
       dz      = max(dz, minthick)
       du(i,k) = sqrt(du(i,k))/dz

       !Calc Richardson
       ri(i,k) = grav*(dt/dz)/(tm*( max(du(i,k), minshear)**2))

       !Blackadar (1962) length scale
       lambdam_x = max( 0.1 * pbllocal * exp( -(ze(i,k) / zkmenv )**2 ) , lambdam2 )
       lambdah_x = max( 0.1 * pbllocal * exp( -(ze(i,k) / zkhenv )**2 ) , lambdah2 )

       alm = almfac * ( karman*ze(i,k)/( 1.0 + karman*(ze(i,k)/lambdam_x) ) )**2
       alh = alhfac * ( karman*ze(i,k)/( 1.0 + karman*(ze(i,k)/lambdah_x) ) )**2

       alh_diag(i,k) = sqrt( alh )

       !Unstable convectve boundary layer
       if ( ri(i,k)  < 0.0 ) then
          ps = ( (zz(i,k-1)/zz(i,k))**(1./3.) - 1.0 ) ** 3
          ps = alh*sqrt( ps/(ze(i,k)*(dz**3)) )
          ps = ri(i,k)/(1.0 + (3.0*louis*louis)*ps*sqrt(-ri(i,k)))
      
          kh(i,k) = 1.0 - (louis*3.0)*ps
          km(i,k) = 1.0 - (louis*2.0)*ps
       end if

       !Stable boundary layer
       if ( ri(i,k) >= 0.0 ) then
          ps      = sqrt(1.0 +  louis     *ri(i,k)   )
         
          kh(i,k) = 1.0 / (1.0 + (louis*3.0)*ri(i,k)*ps)
          km(i,k) = ps  / (ps  + (louis*2.0)*ri(i,k)   )
       end if

       !Dimensionalize Kz and limit diffusvity
       alm = du(i,k)*alm
       alh = du(i,k)*alh

       km(i,k) = min(km(i,k)*alm, akhmmax)
       kh(i,k) = min(kh(i,k)*alh, akhmmax)
       kmls_diag(i,k) = km(i,k)
       khls_diag(i,k) = kh(i,k)

    end do !Loop over levels

 end do !Loop over columns.

end subroutine louis_diff

subroutine tridiag_setup(IRUN,LM,DT,KPBLMIN,ZFULL,PFULL,RDZ,DMI,PHALF,&
                         TV,CT,CQ,CU,T,Q,TH,U,V,DIFF_T,DIFF_M,&
                         AKQ,AKS,AKV,BKQ,BKS,BKV,CKQ,CKS,CKV,EKV,ZPBL)

 IMPLICIT NONE

 !INPUTS
 integer, intent(IN)                          :: IRUN,LM,KPBLMIN
 real(kind_real),    intent(IN)                          :: DT
 real(kind_real),    intent(IN),    dimension(IRUN,LM)   :: T,Q,TH,U,V,ZFULL,PFULL,DMI,TV
 real(kind_real),    intent(IN),    dimension(IRUN,LM-1) :: RDZ
 real(kind_real),    intent(IN),    dimension(IRUN,LM+1) :: PHALF
 real(kind_real),    intent(IN),    dimension(IRUN  )    :: CT, CQ, CU

 !INOUTS
 real(kind_real),    intent(INOUT), dimension(IRUN,LM+1) :: DIFF_T, DIFF_M

 !OUTPUTS
 real(kind_real),    intent(OUT),   dimension(IRUN,LM)   :: AKQ, AKS, AKV
 real(kind_real),    intent(OUT),   dimension(IRUN,LM)   :: BKQ, BKS, BKV
 real(kind_real),    intent(OUT),   dimension(IRUN,LM)   :: CKQ, CKS, CKV
 real(kind_real),    intent(OUT),   dimension(IRUN,LM)   ::           EKV
 real(kind_real),    intent(OUT),   dimension(IRUN)      :: ZPBL

 !LOCALS
 real(kind_real), dimension(LM+1) :: temparray
 integer :: I,L,locmax
 real(kind_real)    :: maxkh,minlval,thetavs,thetavh,uv2h
 real(kind_real), parameter :: tcri_crit = 0.25 !PBL-top diagnostic
 
 real(kind_real),             dimension(LM ) :: tcrib

 do i = 1, irun

    !Update zpbl
    zpbl(i) = 10.e15
    do l=lm,2,-1
       if ( (diff_t(i,l) < 2.) .and. (diff_t(i,l+1) >= 2.) .and. (zpbl(i) == 10.e15 ) ) then
          zpbl(i) = zfull(i,l) 
       end if 
    end do
    if (  zpbl(i) .eq. 10.e15 ) then
       zpbl(i) = zfull(i,lm)
    endif
    zpbl(i) = min(zpbl(i),zfull(i,kpblmin))

    !Second difference coefficients for scalars; RDZ is RHO/DZ, DMI is (G DT)/DP
    do l = 1, lm
       if (l <= lm-1) then
          cks(i,l) = -diff_t(i,l+1) * rdz(i,l)
          if (l == 1) then
             aks(i,l) = 0.0
          endif
          aks(i,l+1) = cks(i,l) * dmi(i,l+1)
          cks(i,l) = cks(i,l) * dmi(i,l)
       end if

       if (l == lm) then
          cks(i,l) = -ct(i) * dmi(i,l)
       endif

       !Fill DIFF_T at level LM+1 with CT * RDZ for diagnostic output
       if (l == lm) then
          diff_t(i,l+1) = ct(i) * (phalf(i,l+1)/(rgas * tv(i,l))) / zfull(i,l)
       endif

       !Water vapor can differ at the surface
       ckq(i,l) = cks(i,l)
       akq(i,l) = aks(i,l)
       if (l == lm) then
          ckq(i,l) = -cq(i) * dmi(i,l)
       endif

       !Second difference coefficients for winds
       !EKV is saved to use in the frictional heating calc.
       if (l <= lm-1) then
          ekv(i,l) = -diff_m(i,l+1) * rdz(i,l)
          if (l == 1) then
             akv(i,l) = 0.0
          endif
          akv(i,l+1) = ekv(i,l) * dmi(i,l+1)
          ckv(i,l) = ekv(i,l) * dmi(i,l)
          ekv(i,l) = -grav * ekv(i,l)
       end if

       if (l == lm) then
          ckv(i,l) = -  cu(i) * dmi(i,l)
          ekv(i,l) =  grav *  cu(i)
       end if

       !Fill DIFF_M at level LM+1 with CU * RDZ for diagnostic output
       if (l == lm) then
          diff_m(i,l+1) = cu(i) * (phalf(i,l+1)/(rgas * tv(i,l))) / zfull(i,l)
       endif

       !Setup the tridiagonal matrix
       bks(i,l) = 1.00 - (aks(i,l)+cks(i,l))
       bkq(i,l) = 1.00 - (akq(i,l)+ckq(i,l))
       bkv(i,l) = 1.00 - (akv(i,l)+ckv(i,l))

    end do

 end do 
 

end subroutine tridiag_setup


subroutine orodrag(IRUN,LM,DT,LAMBDA_B,C_B,FKV,BKV,U,V,ZFULL,VARFLT,PHALF)

 !Subroutine to compute orograpic drag, Beljaars (2003).

 IMPLICIT NONE

 !INPUTS
 integer, intent(IN)                        :: IRUN, LM
 real(kind_real),    intent(IN)                        :: DT, LAMBDA_B, C_B
 real(kind_real),    intent(IN), dimension(IRUN)       :: VARFLT 
 real(kind_real),    intent(IN), dimension(IRUN,LM)    :: U, V, ZFULL
 real(kind_real),    intent(IN), dimension(IRUN,LM+1)  :: PHALF

 !INOUTS
 real(kind_real),    intent(INOUT), dimension(IRUN,LM) :: BKV

 !OUTPUTS
 real(kind_real),    intent(OUT), dimension(IRUN,LM) :: FKV
      
 !LOCALS
 integer :: I, L
 real(kind_real)    :: FKV_temp

 do i = 1, irun
    do l=lm,1,-1

       fkv(i,l) = 0.0

       if (zfull(i,l) < 4.0*lambda_b) then
          fkv_temp = zfull(i,l)*(1.0/lambda_b)
          fkv_temp = varflt(i) * exp(-fkv_temp*sqrt(fkv_temp))*(fkv_temp**(-1.2))
          fkv_temp = (c_b/lambda_b)*min( sqrt(u(i,l)**2+v(i,l)**2),5.0 )*fkv_temp

          bkv(i,l) = bkv(i,l) + dt*fkv_temp
          fkv(i,l) = fkv_temp * (phalf(i,l+1)-phalf(i,l))
       end if
    end do

 end do 

end subroutine orodrag


subroutine lock_diff( ncol,nlev,tdtlw_in_dev,u_star_dev,b_star_dev,frland_dev,t_dev,qv_dev, &
                      qit_dev,qlt_dev,u_dev,v_dev,zfull_dev,pfull_dev,zhalf_dev,phalf_dev, &
                      !INOUTS
                      diff_m_dev,diff_t_dev, &
                      !Constants
                      prandtlsfc_const,prandtlrad_const,beta_surf_const,beta_rad_const, &
                      tpfac_sfc_const,entrate_sfc_const,pceff_sfc_const,khradfac_const,khsfcfac_const,radlw_dep,&
                      ESTBLX )

 IMPLICIT NONE

 !INPUTS
 integer, intent(in)                               :: ncol,nlev    
 real(kind_real),    intent(in)                               :: ESTBLX(:) 
 real(kind_real),    intent(in),    dimension(ncol)           :: u_star_dev,b_star_dev,frland_dev
 real(kind_real),    intent(in),    dimension(ncol,nlev)      :: tdtlw_in_dev,t_dev,qv_dev,qit_dev,qlt_dev
 real(kind_real),    intent(in),    dimension(ncol,nlev)      :: u_dev,v_dev,zfull_dev,pfull_dev
 real(kind_real),    intent(in),    dimension(ncol,1:nlev+1)  :: zhalf_dev, phalf_dev ! 0:72 in GC, 1:73 here.
 real(kind_real),    intent(in)                               :: prandtlsfc_const,prandtlrad_const,beta_surf_const,beta_rad_const
 real(kind_real),    intent(in)                               :: khradfac_const,tpfac_sfc_const,entrate_sfc_const
 real(kind_real),    intent(in)                               :: pceff_sfc_const,khsfcfac_const

 !INOUTS
 real(kind_real),    intent(inout), dimension(ncol,1:nlev+1)  :: diff_m_dev,diff_t_dev

 !LOCALS
 real(kind_real), parameter :: akmax = 1.e4, zcldtopmax = 3.e3, ramp = 20.
 integer :: i,j,k,ibot,itmp,ipbl,kmax,kcldtop,kcldbot,kcldtop2,radlw_dep
 logical :: used,keeplook,do_jump_exit
 real(kind_real) :: maxradf, tmpradf,stab
 real(kind_real) :: maxqdot, tmpqdot,wentrmax
 real(kind_real) :: maxinv, qlcrit,ql00,qlm1,Abuoy,Ashear
 real(kind_real) :: wentr_tmp,hlf,vbulkshr,vbulk_scale
 real(kind_real) :: k_entr_tmp,tmpjump,critjump,radperturb,buoypert
 real(kind_real) :: tmp1, tmp2, slmixture
 real(kind_real) :: vsurf3, vshear3,vrad3, vbr3,dsiems
 real(kind_real) :: dslvcptmp,ztmp
 real(kind_real) :: zradtop
 real(kind_real) :: vrad,svpcp,chis,vbrv
 real(kind_real) :: vsurf, qs_dummy, ptmp
 real(kind_real) :: wentr_rad,wentr_pbl,wentr_brv
 real(kind_real) :: convpbl, radpbl
 real(kind_real), dimension(nlev) :: slv,density,qc,slvcp,hleff,radfq,pblfq,k_m_troen,k_t_troen,k_rad,dqs
 
 !Moved to local as not needed in outputs, but still used in calculations
 real(kind_real), dimension(ncol,1:nlev+1)  :: k_m_entr_dev,k_t_entr_dev
 real(kind_real), dimension(ncol)  :: zsml_dev,zradml_dev,zcloud_dev,zradbase_dev

 qlcrit     = 1.0e-6
 abuoy      = 0.23
 ashear     = 25.0
 wentrmax   = 0.05

 ibot = nlev

 do i = 1, ncol

    zradml_dev(i)     = 0.0
    zcloud_dev(i)     = 0.0
    zradbase_dev(i)   = 0.0
    zsml_dev(i)       = 0.0 ! note that this must be zero, indicates stable surface layer is used in gravity wave drag scheme

    zradtop    = 0.0
    convpbl    = 0.0
    wentr_pbl  = 0.0
    vsurf      = 0.0
    radpbl     = 0.0
    svpcp      = 0.0
    vrad       = 0.0
    wentr_rad  = 0.0

    do k = 1, nlev
       pblfq(k)      = 0.0
       k_t_troen(k)  = 0.0
       k_m_troen(k)  = 0.0
       radfq(k)      = 0.0
       k_rad(k)      = 0.0
       k_t_entr_dev(i,k)   = 0.0
       k_m_entr_dev(i,k)   = 0.0

    end do

    !For GPU we must zero out even the LM+1'th position

    k_t_entr_dev(i,nlev+1)   = 0.0
    k_m_entr_dev(i,nlev+1)   = 0.0

    !set up specific humidities and static energies and compute airdensity
    do k = 1, nlev
       if ( t_dev(i,k) <= tice-ramp ) then
          hleff(k) = alhs
       else if ( (t_dev(i,k) > tice-ramp) .and. (t_dev(i,k) < tice) ) then
          hleff(k) =  ( (t_dev(i,k)-tice+ramp)*alhl +  (tice -t_dev(i,k)    )*alhs   ) / ramp
       else
          hleff(k) = alhl
       end if

       !Compute sat specific humidity and deriv, Qs not actually needed
       !dqs(k) = dqsat(t_dev(i,k), pfull_dev(i,k), qsat=qs(k), pascals=.true. )
       !dqs(k) = dqsat(t_dev(i,k), pfull_dev(i,k), qsat=qs_dummy, pascals=.true. )

       ptmp = pfull_dev(i,k)*0.01
       call dqsat_sub_sca(dqs(k),qs_dummy,t_dev(i,k),ptmp,estblx)

       !Compute total cloud condensate - qc. These are grid box mean values.
       !qc(k) = 0.0 !(qlls_dev(i,k) + qils_dev(i,k))
       qc(k) = (qit_dev(i,k) + qlt_dev(i,k))

       !Compute liquid static energy.
       slv(k) = cp*t_dev(i,k)*(1+vireps*qv_dev(i,k)-qc(k)) + grav*zfull_dev(i,k) - hleff(k)*qc(k)

       density(k) = pfull_dev(i,k)/(rgas*(t_dev(i,k) *(1.+vireps*qv_dev(i,k)-qc(k))))              

    end do

    !Reset indices
    ipbl    = -1
    kcldtop = -1

    !Surface driven convective layers (if b_star_dev > 0., that is,upward surface buoyancy flux)
    if (b_star_dev(i) .gt. 0.) then

       !Find depth of surface driven mixing by raising a parcel from the surface with some excess buoyancy to its level 
       !of neutral buoyancy. Note the use slv as the density variable permits one to goes through phase changes to find parcel top
       call mpbl_depth(i,ncol,nlev,tpfac_sfc_const,entrate_sfc_const,pceff_sfc_const, &
                       t_dev,qv_dev,u_dev,v_dev,zfull_dev,pfull_dev,b_star_dev,u_star_dev,ipbl,zsml_dev,&
                       estblx)
   
       !Define velocity scales vsurf and vshear
       vsurf3   = u_star_dev(i)*b_star_dev(i)*zsml_dev(i) !cubed
       vshear3  = ashear*u_star_dev(i)*u_star_dev(i)*u_star_dev(i)
   
       vsurf = vsurf3**(1./3.)
   
       !Define velocity scale vbulkshr      
       vbulkshr   =  sqrt( ( u_dev(i,ipbl)-u_dev(i,ibot) )**2 + ( v_dev(i,ipbl)-v_dev(i,ibot) )**2  ) /  zsml_dev(i)
       vbulk_scale = 3.0 / 1000.   ! Non-local (PBL-deep) shear scale

       !Following Lock et al. 2000, limit height of surface well mixed layer if interior stable interface is
       !found. An interior stable interface is diagnosed if the slope between 2 full levels is greater than critjump
       critjump = 2.0
       if (ipbl .lt. ibot) then 
          do k = ibot, ipbl+1, -1
             tmpjump =(slv(k-1)-slv(k))/cp 
             if (tmpjump .gt. critjump) then
                ipbl = k
                zsml_dev(i) = zhalf_dev(i,ipbl)
                exit
             end if
          enddo
       end if

       !Compute entrainment rate
       tmp1 = grav*max(0.1,(slv(ipbl-1)-slv(ipbl))/ cp)/(slv(ipbl)/cp)
       tmp2 = ((vsurf3+vshear3)**(2./3.)) / zsml_dev(i)

       wentr_tmp = min( wentrmax,  max(0., (beta_surf_const * (vsurf3 + vshear3)/zsml_dev(i))/(tmp1+tmp2) ) )

       !Fudgey adjustment of entrainment to reduce it for shallow boundary layers, increase for deep ones
       if ( zsml_dev(i) .lt. 1600. ) then 
          wentr_tmp = wentr_tmp * ( zsml_dev(i) / 800. )
       else
          wentr_tmp = 2.*wentr_tmp
       endif

       !More fudgey adjustment of entrainment.
       k_entr_tmp = wentr_tmp*(zfull_dev(i,ipbl-1)-zfull_dev(i,ipbl))  
       k_entr_tmp = min( k_entr_tmp, akmax )

       do k = ipbl, ibot
          pblfq(k)  = 1.
       end do
       convpbl           = 1.
       wentr_pbl         = wentr_tmp
       k_t_troen(ipbl)   = k_entr_tmp
       k_m_troen(ipbl)   = k_entr_tmp
       k_t_entr_dev(i,ipbl) = k_t_entr_dev(i,ipbl) + k_entr_tmp
       k_m_entr_dev(i,ipbl) = k_m_entr_dev(i,ipbl) + k_entr_tmp

       !Compute diffusion coefficients in the interior of the PBL

       if (ipbl .lt. ibot) then
   
          call diffusivity_pbl2(i, ncol, nlev,zsml_dev,khsfcfac_const, k_entr_tmp,vsurf,&
                                frland_dev,zhalf_dev,k_m_troen,k_t_troen)

          do k = ipbl+1, ibot
             k_t_entr_dev(i,k) =  k_t_entr_dev(i,k) + k_t_troen(k)
             k_m_entr_dev(i,k) =  k_m_entr_dev(i,k) + k_m_troen(k)*prandtlsfc_const
          end do

       end if

    end if 


    !NEGATIVELY BUOYANT PLUMES DRIVEN BY LW RADIATIVE COOLING AND/OR BUOYANCY REVERSAL
    !This part is done only if a level kcldtop can be found below zcldtopmax

    kmax = ibot+1
    do k = 1, ibot
       if ( zhalf_dev(i,k) < zcldtopmax) then
          kmax = k
          exit
       end if
    end do

    !Find cloud top and bottom using GRID BOX MEAN or IN-CLOUD value of qc. Decision occurs where qc is calculated
    kcldtop  = ibot+1
    do k = ibot,kmax,-1
       qlm1 = qc(k-1)  ! qc one level UP
       ql00 = qc(k)
       stab = slv(k-1) - slv(k) 
       if ( ( ql00  .ge. qlcrit ) .and. ( qlm1 .lt. qlcrit) .and. (stab .gt. 0.) ) then
           kcldtop  = k   
           exit
       end if
    end do

    !if (kcldtop .ge. ibot+1) then 
    if (kcldtop .lt. ibot+1) then 
       !Only do this is above is true, else skip to the end

       kcldtop2=min( kcldtop+1,nlev)
       ! Look one level further down in case first guess is a thin diffusive veil
       if ( (qc(kcldtop) .lt. 10*qlcrit ) .and. (qc(kcldtop2) .ge. 10*qc(kcldtop) ) ) then
          kcldtop=kcldtop2
       endif

       kcldbot  = ibot+1
       do k = ibot,kcldtop,-1
          qlm1 = qc(k-1)  ! qc one level UP
          ql00 = qc(k)
          if ( ( ql00  .lt. qlcrit ) .and. ( qlm1 .ge. qlcrit) ) then
             kcldbot  = k   
             exit
          end if
       end do

       if (radlw_dep .eq. 1) then

          !if (kcldtop .eq. kcldbot) then 
          if (kcldtop .ne. kcldbot) then 
             !Only do this is above is true, else skip to the end

             !With diffusion of ql, qi "cloud top" found via these quantities may be above radiation max
             kcldtop2=min( kcldtop+2,nlev)
             maxradf = maxval( -1.*tdtlw_in_dev(i,kcldtop:kcldtop2) )
             maxradf = maxradf*cp*( (phalf_dev(i,kcldtop+1)-phalf_dev(i,kcldtop)) / grav )
             maxradf = max( maxradf , 0. ) ! do not consider cloud tops that are heating

             !Calculate optimal mixing fraction (chis) for buoyancy reversal. Use effective heat of evap/subl. Ignore diffs across cldtop
             hlf = hleff(kcldtop)

             tmp1 = ( slv(kcldtop-1)  -  hlf*qc(kcldtop-1) ) - (slv(kcldtop) - hlf*qc(kcldtop) )
             tmp1 = dqs(kcldtop)*tmp1/cp

             tmp2 = ( qv_dev(i,kcldtop-1)   +  qc(kcldtop-1) ) - ( qv_dev(i,kcldtop)  +  qc(kcldtop)   )  
   
             chis = -qc(kcldtop)*( 1 + hlf * dqs(kcldtop) / cp )

             if ( ( tmp2 - tmp1 ) >= 0.0 ) then
                chis = 0.
             else
                chis = chis / ( tmp2 - tmp1 ) 
             endif

             if ( chis .gt. 1.0 ) then
                chis=1.0
             endif
            
             slmixture = ( 1.0-chis ) * ( slv(kcldtop) - hlf*qc(kcldtop)   )   +  chis  * ( slv(kcldtop-1)  -  hlf*qc(kcldtop-1) )

             !Compute temperature of parcel at cloud top, svpcp.
             svpcp = slmixture /cp

             buoypert   = ( slmixture - slv(kcldtop) )/cp

             !Calculate my best guess at the LCs' D parameter attributed to Siems et al.
             stab       = slv(kcldtop-1) - slv(kcldtop)
             if (stab .eq. 0.) then 
                dsiems     =  ( slv(kcldtop) - slmixture ) ! / 1.0  ! arbitrary, needs to be re-thought 
             else
                dsiems     =  ( slv(kcldtop) - slmixture ) / stab
             endif
             dsiems     =  min( dsiems, 10. )
             dsiems     =  max( dsiems,  0. )
             zradtop = zhalf_dev(i,kcldtop)

             !Find depth of radiatively driven convection 
             !Expose radperturb and other funny business outside of radml_depth
             radperturb   = min( maxradf/100. , 0.3 ) ! dim. argument based on 100m deep cloud over 1000s
             do_jump_exit = .true.
             critjump     = 0.3
             svpcp = svpcp - radperturb

             slvcp = slv/cp

             if (kcldtop .lt. ibot) then 
                call radml_depth(i,ncol,nlev,kcldtop,ibot,svpcp,zradtop,critjump,do_jump_exit,&
                                 slvcp,zfull_dev,zhalf_dev,zradbase_dev,zradml_dev )      
             else
                zradbase_dev(i) = 0.0
                zradml_dev(i)   = zradtop  
             end if

             zcloud_dev(i) = zhalf_dev(i,kcldtop) - zhalf_dev(i,kcldbot)

             !if (zradml_dev(i) .le. 0.0 ) then 
             if (zradml_dev(i) .gt. 0.0 ) then 
                !Only do this is above is true, else skip to the end
   
                !Compute radiation driven scale
                vrad3 = grav*zradml_dev(i)*maxradf/density(kcldtop)/slv(kcldtop)   

                !tmp1 here should be the buoyancy jump at cloud top SAK has it w/ resp to parcel property - svpcp. Im not sure about that.
                tmp1 = grav*max(0.1,((slv(kcldtop-1)/cp)-svpcp))/(slv(kcldtop)/cp)

                !Straightforward buoyancy jump across cloud top
                tmp1 = grav*max( 0.1, ( slv(kcldtop-1)-slv(kcldtop) )/cp ) / ( slv(kcldtop) /cp )

                !Compute buoy rev driven scale
                vbr3  = ( max( tmp1*zcloud_dev(i), 0.)**3 )
                vbr3  = abuoy*(chis**2)*max(dsiems,0.)*sqrt( vbr3 ) 

                !Adjust velocity scales to prevent jumps near zradtop=zcldtopmax
                if ( zradtop .gt. zcldtopmax-500. ) then 
                   vrad3 = vrad3*(zcldtopmax - zradtop)/500.  
                   vbr3  = vbr3 *(zcldtopmax - zradtop)/500.  
                endif
                vrad3=max( vrad3, 0. ) ! these really should not be needed
                vbr3 =max( vbr3,  0. )

                vrad = vrad3 ** (1./3.)    
                vbrv = vbr3**(1./3.)
      
                tmp2 = (  vrad**2 + vbrv**2  ) / zradml_dev(i)
                wentr_rad = min(wentrmax,beta_rad_const*(vrad3+vbr3)/zradml_dev(i)/(tmp1+tmp2))
   
                wentr_brv =  beta_rad_const*vbr3/zradml_dev(i)/(tmp1+tmp2)

                !Fudgey adjustment of entrainment to reduce it for shallow boundary layers, increase for deep ones
                if ( zradtop .lt. 500. ) then
                   wentr_rad = 0.00
                endif
                if (( zradtop .gt. 500.) .and. (zradtop .le. 800. )) then
                   wentr_rad = wentr_rad * ( zradtop-500.) / 300.
                endif

                if ( zradtop .lt. 2400. ) then 
                   wentr_rad = wentr_rad * ( zradtop / 800. )
                else
                   wentr_rad = 3.*wentr_rad
                endif

                k_entr_tmp = min( akmax, wentr_rad*(zfull_dev(i,kcldtop-1)-zfull_dev(i,kcldtop)) )

                radfq(kcldtop)          = 1.
                radpbl                  = 1.
                k_rad(kcldtop)          = k_entr_tmp
                k_t_entr_dev(i,kcldtop) = k_t_entr_dev(i,kcldtop) + k_entr_tmp
                k_m_entr_dev(i,kcldtop) = k_m_entr_dev(i,kcldtop) + k_entr_tmp
   
                !Handle case of radiatively driven top being the same top as surface driven top
                if (ipbl .eq. kcldtop .and. ipbl .gt. 0) then

                   tmp2 = ((vbr3+vrad3+vsurf3+vshear3)**(2./3.)) / zradml_dev(i)
      
                   wentr_rad = min( wentrmax,  max(0., ((beta_surf_const *(vsurf3 + vshear3) + &
                                             beta_rad_const*(vrad3+vbr3) )/ zradml_dev(i))/(tmp1+tmp2) ) )
                   wentr_pbl       = wentr_rad

                   k_entr_tmp = min( akmax, wentr_rad*(zfull_dev(i,kcldtop-1)-zfull_dev(i,kcldtop)) )

                   pblfq(ipbl)             = 1.
                   radfq(kcldtop)          = 1.
                   radpbl                  = 1.
                   k_rad(kcldtop)          = k_entr_tmp
                   k_t_troen(ipbl)         = k_entr_tmp
                   k_m_troen(ipbl)         = k_entr_tmp
                   k_t_entr_dev(i,kcldtop) = k_entr_tmp
                   k_m_entr_dev(i,kcldtop) = k_entr_tmp
   
                end if

                !If there are any interior layers to calculate diffusivity
                if ( kcldtop .lt. ibot ) then   
   
                   do k = kcldtop+1,ibot
   
                      ztmp = max(0.,(zhalf_dev(i,k)-zradbase_dev(i))/zradml_dev(i) )
   
                      if (ztmp.gt.0.) then
   
                         radfq(k) = 1.
                         k_entr_tmp = khradfac_const*karman*( vrad+vbrv )*ztmp*zradml_dev(i)*ztmp*((1.-ztmp)**0.5)
                         k_entr_tmp = min( k_entr_tmp, akmax )
                         k_rad(k) = k_entr_tmp
                         k_t_entr_dev(i,k) = k_t_entr_dev(i,k) + k_entr_tmp
                         k_m_entr_dev(i,k) = k_m_entr_dev(i,k) + k_entr_tmp*prandtlrad_const
   
                      end if
                   enddo
      
                end if

                !Handle special case of zradbase_dev < zsml_dev, in this case there should be no entrainment from the surface.
                if (zradbase_dev(i) .lt. zsml_dev(i) .and. convpbl .eq. 1. .and. ipbl .gt. kcldtop) then
                   wentr_pbl            = 0.
                   pblfq(ipbl)          = 0.
                   k_t_entr_dev(i,ipbl) = k_t_entr_dev(i,ipbl) - k_t_troen(ipbl)
                   k_m_entr_dev(i,ipbl) = k_m_entr_dev(i,ipbl) - k_m_troen(ipbl)          
                   k_t_troen(ipbl)      = 0.
                   k_m_troen(ipbl)      = 0. 
                end if

             endif

          endif !kcldtop .ne. kcldbot
 
       endif !radlw_dep    
 
    endif

    !Modify diffusivity coefficients using MAX( A , B )        
    do k = 2, ibot     
       diff_t_dev(i,k)   = max( k_t_entr_dev(i,k), diff_t_dev(i,k) )
       diff_m_dev(i,k)   = max( k_m_entr_dev(i,k), diff_m_dev(i,k) )
       k_t_entr_dev(i,k) = max( k_t_entr_dev(i,k),             0.0 )
       k_m_entr_dev(i,k) = max( k_m_entr_dev(i,k),             0.0 )
    enddo


 end do

end subroutine lock_diff


!======================================================================= 
subroutine mpbl_depth(i,ncol,nlev,tpfac, entrate, pceff, t, q, u, v, z, p, b_star, u_star , ipbl, ztop, ESTBLX)

 !Subroutine to calculate pbl depth

 IMPLICIT NONE

 !INPUTS
 integer, intent(in)                       :: i, nlev, ncol
 real(kind_real),    intent(in)                       :: ESTBLX(:)
 real(kind_real),    intent(in), dimension(ncol,nlev) :: t, z, q, p, u, v
 real(kind_real),    intent(in), dimension(ncol)      :: b_star, u_star
 real(kind_real),    intent(in)                       :: tpfac, entrate, pceff

 !OUTPUTS
 integer, intent(out)                       :: ipbl
 real(kind_real),    intent(out),dimension(ncol)       :: ztop

 !LOCALS
 real(kind_real) :: tep,z1,z2,t1,t2,qp,pp,qsp,dqp,dqsp,u1,v1,u2,v2,du
 real(kind_real) :: ws,k_t_ref,entfr, tpfac_x, entrate_x,vscale,ptmp
 integer :: k

 !Calculate surface parcel properties
 tep  = t(i,nlev)
 qp   = q(i,nlev)

 !Wind dependence of plume character. 
 vscale    = 0.25 / 100. ! change of .25 m s-1 in 100 m

 !tpfac scales up bstar by inv. ratio of heat-bubble area to stagnant area
 tep  = tep * (1.+ tpfac * b_star(i)/grav)

 !Search for level where this is exceeded              
 t1   = t(i,nlev)
 v1   = v(i,nlev)
 u1   = u(i,nlev)
 z1   = z(i,nlev)
 ztop(i) = z1

 do k = nlev-1 , 2, -1
    z2 = z(i,k)
    t2 = t(i,k)
    u2 = u(i,k)
    v2 = v(i,k)
    pp = p(i,k)

    du = sqrt( ( u2 - u1 )**2 + ( v2 - v1 )**2 ) / (z2-z1)
    du = min(du,1.0e-8)

    entrate_x = entrate * ( 1.0 + du / vscale )

    entfr= min( entrate_x*(z2-z1), 0.99 )

    qp  = qp  + entfr*(q(i,k)-qp)

    !Dry adiabatic ascent through one layer. Static energy conserved. 
    tep = tep - grav*( z2-z1 )/cp

    !Environmental air entrained
    tep = tep + entfr*(t(i,k)-tep)

    !dqsp = dqsat(tep , pp , qsat=qsp,  pascals=.true. )
    ptmp = pp*0.01
    call dqsat_sub_sca(dqsp,qsp,tep,ptmp,estblx)

    dqp = max( qp - qsp, 0. )/(1.+(alhl/cp)*dqsp )
    qp  = qp - dqp
    !"Precipitation efficiency" basically means fraction of condensation heating that gets applied to parcel
    tep = tep  + pceff * alhl * dqp/cp  

    !If parcel temperature (tep) colder than env (t2) OR if entrainment too big, declare this the PBL top
    if ( (t2 .ge. tep) .or. ( entfr .ge. 0.9899 ) ) then
       ztop(i) = 0.5*(z2+z1)
       ipbl = k+1
       exit
    end if

    z1 = z2
    t1 = t2
    u1 = u2
    v1 = v2

 enddo

end subroutine mpbl_depth

subroutine radml_depth(i, ncol, nlev, toplev, botlev, svp, zt, critjump, do_jump_exit, t, zf, zh, zb, zml)

 !Subroutine to calculate bottom and depth of radiatively driven mixed layer

 IMPLICIT NONE

 !INPUTS
 integer, intent(in)                         :: i, toplev, botlev, ncol, nlev
 real(kind_real),    intent(in)                         :: svp, zt, critjump
 real(kind_real),    intent(in), dimension(nlev)        :: t
 real(kind_real),    intent(in), dimension(ncol,nlev)   :: zf
 real(kind_real),    intent(in), dimension(ncol,nlev+1) :: zh
 logical, intent(in)                         :: do_jump_exit

 !OUTPUTS
 real(kind_real),    intent(out), dimension(ncol)        :: zb, zml

 !LOCALS
 real(kind_real)    :: svpar,h1,h2,t1,t2,entrate,entfr
 integer :: k

 !Initialize zml
 zml(i) = 0.

 !Calculate cloud top parcel properties
 svpar  = svp
 h1   = zf(i,toplev)
 t1   = t(toplev)
 entrate = 0.2/200.

 !Search for level where parcel is warmer than env. First cut out if parcel is already warmer than cloudtop. 
 if (t1.lt.svpar) then
    zb(i)  = h1
    zml(i) = 0.00
    return !udh - get rid of this if linearising
 endif

 !If above is false keep looking
 do k = 1,botlev
    if (k > toplev) then
       h2 = zf(i,k)
       t2 = t(k)

       if (t2.lt.svpar) then
          if ( abs(t1 - t2 ) .gt. 0.2 ) then 
             zb(i) = h2 + (h1 - h2)*(svpar - t2)/(t1 - t2)
             zb(i) = max( zb(i) , 0. )  ! final protection against interp problems
          else
             zb(i) = h2
          endif
          zml(i) = zt - zb(i)
          return !dh get rid of this if linearising
       end if

       if (do_jump_exit .and. (t1-t2) .gt. critjump .and. k .gt. toplev+1) then
          zb(i) = zh(i,k)
          zml(i) = zt - zb(i)
          return !dh get rid of this if linearising
       end if

       entfr = min( entrate*(h1-h2), 1.0 )
       svpar = svpar + entfr*(t2-svpar)

       h1 = h2
       t1 = t2

    end if

 enddo

 zb(i) = 0.
 zml(i) = zt

end subroutine radml_depth

subroutine diffusivity_pbl2(i, ncol, lm, h, kfac, k_ent, vsurf, frland, zm, k_m, k_t)

 !Subroutine to return the vertical K-profile of diffusion coefficients for the surface driven convective mixed layer.

 IMPLICIT NONE

 !INPUTS
 integer, intent(in)                         :: i, lm, ncol
 real(kind_real),    intent(in), dimension(ncol)        :: h, frland
 real(kind_real),    intent(in)                         :: k_ent, vsurf, kfac 
 real(kind_real),    intent(in), dimension(ncol,1:lm+1) :: zm

 !OUTPUTS
 real(kind_real),    intent(out), dimension(lm)         :: k_m, k_t

 !LOCALS
 real(kind_real) :: EE, hin, kfacx 
 integer :: k, kk

 !Kluge. Raise KHs over land
 if ( frland(i) < 0.5 ) then 
    kfacx = kfac 
 else
    kfacx = kfac*2.0
 endif

 hin    = 0.0 ! 200.  ! "Organization" scale for plume (m).

 do k = 1, lm
    k_m(k) = 0.0
    k_t(k) = 0.0
 end do

 if ( vsurf*h(i) .gt. 0. ) then
    ee  = 1.0 - sqrt( k_ent / ( kfacx * karman * vsurf * h(i) ) )
    ee  = max( ee , 0.7 )  ! If EE is too small, then punt, as LCs
    do k = 1, lm
       if ( ( zm(i,k) .le. h(i) ) .and.  ( zm(i,k) .gt. hin )  ) then
          k_t(k) = kfacx * karman * vsurf * ( zm(i,k)-hin ) * ( 1. - ee*( (zm(i,k)-hin)/(h(i)-hin) ))**2
       end if
    end do
 endif

 do k = 1, lm
    k_m(k) = k_t(k)
 end do

end subroutine diffusivity_pbl2

subroutine esinit(ESTBLX)
!Initialise look-up table

 IMPLICIT NONE

 !OUTPUT
 real(kind_real), dimension(TABLESIZE) :: ESTBLX

 !LOCALS
 real(kind_real), parameter :: ZEROC = 273.16, tmix = -20.0

 real(kind_real), dimension(TABLESIZE) :: ESTBLE, ESTBLW

 integer :: I
 real(kind_real)    :: T, DELTA_T

 delta_t = 1.0/degsubs

 do i=1,tablesize

    t = (i-1)*delta_t + tmintbl

    if(t>zeroc) then
       call qsatlqu0(estble(i),t)
    else
       call qsatice0(estble(i),t)
    end if

    call qsatlqu0(estblw(i),t)

    t = t-zeroc
    if(t>=tmix .and. t<0.0) then
       estblx(i) = ( t/tmix )*( estble(i) - estblw(i) ) + estblw(i)
    else
       estblx(i) = estble(i)
    end if

 end do

 end subroutine esinit


subroutine qsatlqu0(QS,TL)
!SUPERSATURATED AS LIQUID

 IMPLICIT NONE

 !INPUTS
 real(kind_real) :: TL!, TMAXTBL

 !OUTPUTS
 real(kind_real) :: QS

 !LOCALS
 real(kind_real), parameter :: ZEROC   = 273.16
 real(kind_real), parameter :: TMINLQU = zeroc - 40.0

 real*8,  parameter :: B6 = 6.136820929e-11*100.0
 real*8,  parameter :: B5 = 2.034080948e-8 *100.0
 real*8,  parameter :: B4 = 3.031240396e-6 *100.0
 real*8,  parameter :: B3 = 2.650648471e-4 *100.0
 real*8,  parameter :: B2 = 1.428945805e-2 *100.0
 real*8,  parameter :: B1 = 4.436518521e-1 *100.0
 real*8,  parameter :: B0 = 6.107799961e+0 *100.0

 real(8) :: TX, EX, TI, TT

 tx = tl

 if    (tx<tminlqu) then
    ti = tminlqu
 elseif(tx>tmaxtbl) then
    ti = tmaxtbl
 else
    ti = tx
 end if

 tt = ti-zeroc  !Starr polynomial fit
 ex = (tt*(tt*(tt*(tt*(tt*(tt*b6+b5)+b4)+b3)+b2)+b1)+b0)

 tl = tx
 qs = ex

 return

end subroutine qsatlqu0


subroutine qsatice0(QS,TL)
!SUPERSATURATED AS ICE

 IMPLICIT NONE   
      
 !INPUTS
 real(kind_real) :: TL

 !OUTPUTS
 real(kind_real) :: QS

 !LOCALS
 real(kind_real), parameter :: ZEROC = 273.16, tminstr = -95.0
 real(kind_real), parameter :: TMINICE = zeroc + tminstr

 real(kind_real), parameter :: TSTARR1 = -75.0, tstarr2 = -65.0, tstarr3 = -50.0,  tstarr4 = -40.0

 real*8,  parameter :: BI6= 1.838826904e-10*100.0
 real*8,  parameter :: BI5= 4.838803174e-8 *100.0
 real*8,  parameter :: BI4= 5.824720280e-6 *100.0
 real*8,  parameter :: BI3= 4.176223716e-4 *100.0
 real*8,  parameter :: BI2= 1.886013408e-2 *100.0
 real*8,  parameter :: BI1= 5.034698970e-1 *100.0
 real*8,  parameter :: BI0= 6.109177956e+0 *100.0
 real*8,  parameter :: S16= 0.516000335e-11*100.0
 real*8,  parameter :: S15= 0.276961083e-8 *100.0
 real*8,  parameter :: S14= 0.623439266e-6 *100.0
 real*8,  parameter :: S13= 0.754129933e-4 *100.0
 real*8,  parameter :: S12= 0.517609116e-2 *100.0
 real*8,  parameter :: S11= 0.191372282e+0 *100.0
 real*8,  parameter :: S10= 0.298152339e+1 *100.0
 real*8,  parameter :: S26= 0.314296723e-10*100.0
 real*8,  parameter :: S25= 0.132243858e-7 *100.0
 real*8,  parameter :: S24= 0.236279781e-5 *100.0
 real*8,  parameter :: S23= 0.230325039e-3 *100.0
 real*8,  parameter :: S22= 0.129690326e-1 *100.0
 real*8,  parameter :: S21= 0.401390832e+0 *100.0
 real*8,  parameter :: S20= 0.535098336e+1 *100.0

 real(kind_real) :: TX, TI, TT, W, EX

 tx = tl

 if (tx<tminice) then
    ti = tminice
 elseif(tx>zeroc  ) then
    ti = zeroc
 else
    ti = tx
 end if

 tt = ti - zeroc
 if (tt < tstarr1 ) then
     ex = (tt*(tt*(tt*(tt*(tt*(tt*s16+s15)+s14)+s13)+s12)+s11)+s10)
 elseif(tt >= tstarr1 .and. tt < tstarr2) then
     w = (tstarr2 - tt)/(tstarr2-tstarr1)
     ex =       w *(tt*(tt*(tt*(tt*(tt*(tt*s16+s15)+s14)+s13)+s12)+s11)+s10) &
              + (1.-w)*(tt*(tt*(tt*(tt*(tt*(tt*s26+s25)+s24)+s23)+s22)+s21)+s20)
 elseif(tt >= tstarr2 .and. tt < tstarr3) then
     ex = (tt*(tt*(tt*(tt*(tt*(tt*s26+s25)+s24)+s23)+s22)+s21)+s20)
 elseif(tt >= tstarr3 .and. tt < tstarr4) then
     w = (tstarr4 - tt)/(tstarr4-tstarr3)
     ex =       w *(tt*(tt*(tt*(tt*(tt*(tt*s26+s25)+s24)+s23)+s22)+s21)+s20) &
              + (1.-w)*(tt*(tt*(tt*(tt*(tt*(tt*bi6+bi5)+bi4)+bi3)+bi2)+bi1)+bi0)
 else
     ex = (tt*(tt*(tt*(tt*(tt*(tt*bi6+bi5)+bi4)+bi3)+bi2)+bi1)+bi0)
 endif

 qs = ex

 return
 
end subroutine qsatice0

subroutine dqsat_sub_sca(DQSi,QSSi,TEMP,PLO,ESTBLX)
!COMPUTES SATURATION VAPOUR PRESSURE QSSi AND GRADIENT w.r.t TEMPERATURE DQSi.
!INPUTS ARE TEMPERATURE AND PLO (PRESSURE AT T-LEVELS)
!VALES ARE COMPUTED FROM LOOK-UP TALBE (PIECEWISE LINEAR)

 IMPLICIT NONE

 !Inputs
 real(kind_real) :: TEMP, PLO
 real(kind_real) :: ESTBLX(:)

 !Outputs
 real(kind_real) :: DQSi, QSSi

 !Locals
 real(kind_real), parameter :: MAX_MIXING_RATIO = 1.0
 real(kind_real),    parameter :: ESFAC = h2omw/airmw

 real(kind_real) :: TL, TT, TI, DQSAT, QSAT, DQQ, QQ, PL, PP, DD
 integer :: IT

 tl = temp
 pl = plo

 pp = pl*100.0
 
 if (tl<=tmintbl) then
    ti = tmintbl
 elseif(tl>=tmaxtbl-.001) then
    ti = tmaxtbl-.001
 else
    ti = tl
 end if

 tt = (ti - tmintbl)*degsubs+1
 it = int(tt)

 dqq =  estblx(it+1) - estblx(it)
 qq  =  (tt-it)*dqq + estblx(it)

 if (pp <= qq) then
    qsat = max_mixing_ratio
    dqsat = 0.0
 else
    dd = 1.0/(pp - (1.0-esfac)*qq)
    qsat = esfac*qq*dd
    dqsat = (esfac*degsubs)*dqq*pp*(dd*dd)
 end if

 dqsi = dqsat
 qssi = qsat

end subroutine dqsat_sub_sca



END MODULE bldriver