#include <misc.h>
#include <params.h>

module turbulence 2,5
!---------------------------------------------------------------------------------
! Module to compute mixing coefficients associated with turbulence in the 
! planetary boundary layer and elsewhere.
!
! calling sequence:
!
!    trbinti       initializes time independent coefficients
!     .
!     .
!    vdiff_eddy
!      trbintr        interface for vertical diffusion and pbl scheme and write output
!          trbintd       initializes time dependent variables
!          pblintd       initializes time dependent variables that depend pbl depth
!          vd_k_freeatm  computes mixing coefficients for free atmosphere
!          vd_k_pbl      computes mixing coefficients for pbl
!
!---------------------------Code history--------------------------------
! Standardized:      J. Rosinski, June 1992
! Reviewed:          P. Rasch, B. Boville, August 1992
! Reviewed:          P. Rasch, April 1996
! Reviewed:          B. Boville, April 1996
! rewritten:         B. Boville, May 2000
! rewritten:         B. Stevens, August 2000
!---------------------------------------------------------------------------------
  use shr_kind_mod, only: r8 => shr_kind_r8
  use ppgrid, only    : pver, pverp, pcols
  use pmgrid, only    : masterproc
  use constituents, only: pcnst, pnats
  use history, only: outfld

  implicit none
!
! PBL limits
!
  real(r8), parameter :: ustar_min = 0.01        ! min permitted value of ustar
  real(r8), parameter :: pblmaxp   = 4.e4        ! pbl max depth in pressure units
  real(r8), parameter :: zkmin     = 0.01        ! Minimum kneutral*f(ri)
!
! PBL Parameters
!
  real(r8), parameter :: onet  = 1./3. ! 1/3 power in wind gradient expression
  real(r8), parameter :: betam = 15.0  ! Constant in wind gradient expression
  real(r8), parameter :: betas =  5.0  ! Constant in surface layer gradient expression
  real(r8), parameter :: betah = 15.0  ! Constant in temperature gradient expression 
  real(r8), parameter :: fakn  =  7.2  ! Constant in turbulent prandtl number
  real(r8), parameter :: fak   =  8.5  ! Constant in surface temperature excess         
  real(r8), parameter :: ricr  =  0.3  ! Critical richardson number
  real(r8), parameter :: sffrac=  0.1  ! Surface layer fraction of boundary layer
  real(r8), parameter :: binm  = betam*sffrac       ! betam * sffrac
  real(r8), parameter :: binh  = betah*sffrac       ! betah * sffrac
!
! Pbl constants set using values from other parts of code
!
  real(r8), save :: cpair      ! Specific heat of dry air
  real(r8), save :: rair       ! Gas const for dry air
  real(r8), save :: zvir       ! rh2o/rair - 1
  real(r8), save :: g          ! Gravitational acceleration
  real(r8), save :: ml2(pverp) ! Mixing lengths squared
  real(r8), save :: vk         ! Von Karman's constant
  real(r8), save :: ccon       ! fak * sffrac * vk

  integer, save  :: npbl       ! Maximum number of levels in pbl from surface
  integer, save  :: ntop_turb  ! Top level to which turbulent vertical diffusion is applied.
  integer, save  :: nbot_turb  ! Bottom level to which turbulent vertical diff is applied.

CONTAINS
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================

 subroutine trbinti(gravx,      cpairx, rairx, zvirx, ntop_eddy,  & 1
                    nbot_eddy,  hypm,   vkx )
!----------------------------------------------------------------------- 
! 
! Purpose:  
! Initialize time independent variables of turbulence/pbl package.
! 
! Author: B. Boville, B. Stevens (August 2000)
! 
!-----------------------------------------------------------------------
    implicit none
   !------------------------------Arguments--------------------------------
    real(r8), intent(in) :: gravx     ! acceleration of gravity
    real(r8), intent(in) :: cpairx    ! specific heat of dry air
    real(r8), intent(in) :: rairx     ! gas constant for dry air
    real(r8), intent(in) :: zvirx     ! rh2o/rair - 1
    real(r8), intent(in) :: hypm(pver)! reference pressures at midpoints
    real(r8), intent(in) :: vkx       ! Von Karman's constant
    integer, intent(in)  :: ntop_eddy ! Top level to which eddy vert diff is applied.
    integer, intent(in)  :: nbot_eddy ! Bottom level to which eddy vert diff is applied.
   !---------------------------Local workspace-----------------------------
    integer :: k                     ! vertical loop index
   !-----------------------------------------------------------------------
   !
   ! Basic constants
   !
    cpair = cpairx
    rair  = rairx
    g     = gravx
    zvir  = zvirx
    vk    = vkx
    ccon  = fak*sffrac*vk
    ntop_turb = ntop_eddy
    nbot_turb = nbot_eddy
   !
   ! Set the square of the mixing lengths.
   !
    ml2(ntop_turb) = 0.
    do k = ntop_turb+1, nbot_turb
       ml2(k) = 30.0**2
    end do
    ml2(nbot_turb+1) = 0.
   !
   ! Limit pbl height to regions below 400 mb
   ! npbl = max number of levels (from bottom) in pbl
   !
    npbl = 0
    do k=nbot_turb,ntop_turb,-1
       if (hypm(k) >= pblmaxp) then
          npbl = npbl + 1
       end if
    end do
    npbl = max(npbl,1)

    if (masterproc) then
       write(6,*)'TRBINTI: PBL height will be limited to bottom ',npbl, &
            ' model levels. Top is ',hypm(pverp-npbl),' pascals'
    end if

    return
 end subroutine trbinti
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================

 subroutine trbintr(lchnk   ,ncol    ,                            & 1,11
                    th      ,t       ,q       ,z       ,zi      , &
                    pmid    ,u       ,v       ,taux    ,tauy    , &
                    shflx   ,cflx    ,obklen  ,ustar   ,pblh    , &
                    kvm     ,kvh     ,cgh     ,cgs     ,kqfs    , &
                    tpert   ,qpert   ,ktopbl  ,ktopblmn,cldn    , &
                    ocnfrac ,no_outfld_in )
!----------------------------------------------------------------------- 
! 
! Purpose: 
!  Interface routines for calcualtion and diatnostics of turbulence related
!  coefficients
!
! Author: B. Stevens (rewrite August 2000)
! 
!-----------------------------------------------------------------------
    implicit none
   !------------------------------Arguments--------------------------------
   !
   ! Input arguments
   !
    integer, intent(in) :: lchnk                      ! chunk identifier
    integer, intent(in) :: ncol                       ! number of atmospheric columns

    real(r8), intent(in)  :: th(pcols,pver)           ! potential temperature [K]
    real(r8), intent(in)  :: t(pcols,pver)            ! temperature (used for density)
    real(r8), intent(in)  :: q(pcols,pver,pcnst+pnats)! specific humidity [kg/kg]
    real(r8), intent(in)  :: z(pcols,pver)            ! height above surface [m]
    real(r8), intent(in)  :: zi(pcols,pverp)          ! height above surface [m]
    real(r8), intent(in)  :: u(pcols,pver)            ! zonal velocity
    real(r8), intent(in)  :: v(pcols,pver)            ! meridional velocity
    real(r8), intent(in)  :: taux(pcols)              ! zonal stress [N/m2]
    real(r8), intent(in)  :: tauy(pcols)              ! meridional stress [N/m2]
    real(r8), intent(in)  :: shflx(pcols)             ! sensible heat flux
    real(r8), intent(in)  :: cflx(pcols,pcnst+pnats)  ! constituent flux
    real(r8), intent(in)  :: pmid(pcols,pver)         ! midpoint pressures
    real(r8), intent(in)  :: cldn(pcols,pver)         ! new cloud fraction
    real(r8), intent(in)  :: ocnfrac(pcols)           ! Land fraction
    logical, intent(in), optional :: no_outfld_in     ! if true, don't write out

   !
   ! Output arguments
   !
    real(r8), intent(out) :: kqfs(pcols,pcnst+pnats)  ! kinematic surf constituent flux (kg/m2/s)
    real(r8), intent(out) :: kvm(pcols,pverp)         ! eddy diffusivity for momentum [m2/s]
    real(r8), intent(out) :: kvh(pcols,pverp)         ! eddy diffusivity for heat [m2/s]
    real(r8), intent(out) :: cgh(pcols,pverp)         ! counter-gradient term for heat [J/kg/m]
    real(r8), intent(out) :: cgs(pcols,pverp)         ! counter-gradient star (cg/flux)
    real(r8), intent(out) :: tpert(pcols)             ! convective temperature excess
    real(r8), intent(out) :: qpert(pcols)             ! convective humidity excess
    real(r8), intent(out) :: ustar(pcols)             ! surface friction velocity [m/s]
    real(r8), intent(out) :: obklen(pcols)            ! Obukhov length
    real(r8), intent(out) :: pblh(pcols)              ! boundary-layer height [m]

    integer, intent(out)  :: ktopbl(pcols)            ! index of first midpoint inside pbl
    integer, intent(out)  :: ktopblmn                 ! min value of ktopbl
   !
   !---------------------------Local workspace-----------------------------
   !
    real(r8) :: wstar(pcols)            ! convective velocity scale [m/s]
    real(r8) :: khfs(pcols)             ! kinimatic surface heat flux 
    real(r8) :: kbfs(pcols)             ! surface buoyancy flux 
    real(r8) :: kvf(pcols,pverp)        ! free atmospheric eddy diffsvty [m2/s]
    real(r8) :: s2(pcols,pver)          ! shear squared
    real(r8) :: n2(pcols,pver)          ! brunt vaisaila frequency
    real(r8) :: ri(pcols,pver)          ! richardson number: n2/s2
    logical  :: no_outfld               ! default = .false. ; if .true. don't write out


    no_outfld = .false.
    if ( present(no_outfld_in) ) then
       no_outfld = no_outfld_in
    endif

   !
   ! Initialize time dependent variables that do not depend on pbl height
   !
    call trbintd(lchnk   ,ncol    ,                            &
                 th      ,q       ,z       ,u       ,v       , &
                 t       ,pmid    ,cflx    ,shflx   ,taux    , &
                 tauy    ,ustar   ,obklen  ,kqfs    ,khfs    , &
                 kbfs    ,s2      ,n2      ,ri      )
   !
   ! Initialize time dependent variables that do depend on pbl height
   !
    call  pblintd(lchnk   ,ncol    ,                            &
                  th      ,q       ,z       ,u       ,v       , &
                  ustar   ,obklen  ,kbfs    ,pblh    ,wstar   , &
                  zi      ,cldn    ,ocnfrac)
   !
   ! Get free atmosphere exchange coefficients
   !
    call austausch_atm(lchnk   ,ncol    ,ri      ,s2      ,kvf     )
   ! 
   ! Get pbl exchange coefficients
   !
    call austausch_pbl(lchnk   ,ncol    ,                            &
                       z       ,kvf     ,kqfs    ,khfs    ,kbfs    , &
                       obklen  ,ustar   ,wstar   ,pblh    ,kvm     , &
                       kvh     ,cgh     ,cgs     ,tpert   ,qpert   , &
                       ktopbl  ,ktopblmn)
   !

    if ( .not. no_outfld ) then 
       call outfld ('PBLH    ',pblh ,pcols,lchnk)
       call outfld ('TPERT   ',tpert,pcols,lchnk)
       call outfld ('QPERT   ',qpert,pcols,lchnk)
       call outfld ('USTAR   ',ustar,pcols,lchnk)
       call outfld ('KVH     ',kvh,pcols,lchnk)
       call outfld ('KVM     ',kvm,pcols,lchnk)
       call outfld ('CGS     ',cgs,pcols,lchnk)
    endif

    return
 end subroutine trbintr
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================

 subroutine trbintd(lchnk   ,ncol    ,                            & 1,1
                    th      ,q       ,z       ,u       ,v       , &
                    t       ,pmid    ,cflx    ,shflx   ,taux    , &
                    tauy    ,ustar   ,obklen  ,kqfs    ,khfs    , &
                    kbfs    ,s2      ,n2      ,ri      )
!----------------------------------------------------------------------- 
! 
! Purpose: 
!  Time dependent initialization
! 
! Method: 
!  Diagnosis of variables that do not depend on mixing assumptions or
!  PBL depth
!
! Author: B. Stevens (extracted from pbldiff, August, 2000)
! 
!-----------------------------------------------------------------------
    implicit none
   !------------------------------Arguments--------------------------------
   !
   ! Input arguments
   !
    integer, intent(in) :: lchnk                     ! chunk identifier
    integer, intent(in) :: ncol                      ! number of atmospheric columns

    real(r8), intent(in)  :: th(pcols,pver)          ! potential temperature [K]
    real(r8), intent(in)  :: q(pcols,pver,pcnst+pnats)     ! specific humidity [kg/kg]
    real(r8), intent(in)  :: z(pcols,pver)           ! height above surface [m]
    real(r8), intent(in)  :: u(pcols,pver)           ! windspeed x-direction [m/s]
    real(r8), intent(in)  :: v(pcols,pver)           ! windspeed y-direction [m/s]
    real(r8), intent(in)  :: t(pcols,pver)           ! temperature (used for density)
    real(r8), intent(in)  :: pmid(pcols,pver)        ! midpoint pressures
    real(r8), intent(in)  :: cflx(pcols,pcnst+pnats) ! surface constituent flux (kg/m2/s)
    real(r8), intent(in)  :: shflx(pcols)            ! surface heat flux (W/m2)
    real(r8), intent(in)  :: taux(pcols)             ! surface u stress [N/m2]
    real(r8), intent(in)  :: tauy(pcols)             ! surface v stress [N/m2]

   !
   ! Output arguments
   !
    real(r8), intent(out) :: ustar(pcols)            ! surface friction velocity [m/s]
    real(r8), intent(out) :: obklen(pcols)           ! Obukhov length
    real(r8), intent(out) :: khfs(pcols)             ! surface kinematic heat flux [mK/s]
    real(r8), intent(out) :: kbfs(pcols)             ! sfc kinematic buoyancy flux [m^2/s^3]
    real(r8), intent(out) :: kqfs(pcols,pcnst+pnats) ! sfc kinematic constituent flux [m/s]
    real(r8), intent(out) :: s2(pcols,pver)          ! shear squared
    real(r8), intent(out) :: n2(pcols,pver)          ! brunt vaisaila frequency
    real(r8), intent(out) :: ri(pcols,pver)          ! richardson number: n2/s2
   !
   !---------------------------Local workspace-----------------------------
   !
    integer  :: i                        ! longitude index
    integer  :: k                        ! level index
    integer  :: m                        ! constituent index

    real(r8) :: thvsrf(pcols)            ! sfc (bottom) level virtual temperature
    real(r8) :: thv(pcols,pver)          ! bulk Richardson no. from level to ref lev
    real(r8) :: rrho(pcols)              ! 1./bottom level density (temporary)
    real(r8) :: vvk                      ! velocity magnitude squared
    real(r8) :: dvdz2                    ! velocity shear squared
    real(r8) :: dz                       ! delta z between midpoints
   !
   ! Compute ustar, and kinematic surface fluxes from surface energy fluxes
   !
    do i=1,ncol
       rrho(i)   = rair*t(i,pver)/pmid(i,pver)
       ustar(i)  = max(sqrt(sqrt(taux(i)**2 + tauy(i)**2)*rrho(i)),ustar_min)
       khfs(i)   = shflx(i)*rrho(i)/cpair
    end do
    do m=1,pcnst+pnats
       do i=1,ncol
          kqfs(i,m)= cflx(i,m)*rrho(i)
       end do
    end do
   !
   ! Compute Obukhov length virtual temperature flux and various arrays for use later:
   !
    do i=1,ncol
       kbfs(i)      = khfs(i) + 0.61*th(i,pver)*kqfs(i,1)
       thvsrf(i)    = th(i,pver)*(1.0 + 0.61*q(i,pver,1))
       obklen(i)    = -thvsrf(i)*ustar(i)**3/(g*vk*(kbfs(i) + sign(1.e-10_r8,kbfs(i))))
    end do
   !
   ! Compute shear squared (s2), brunt vaisaila frequency (n2) and related richardson
   ! number (ri).  For the n2 calcualtion use the dry theta_v derived from virtem
   !
    call virtem(ncol, pcols, pver, th      ,q(1,1,1),0.61_r8 ,thv)
    do k=ntop_turb,nbot_turb-1
       do i=1,ncol
          dvdz2   = (u(i,k)-u(i,k+1))**2 + (v(i,k)-v(i,k+1))**2
          dvdz2   = max(dvdz2,1.e-36_r8)
          dz      = z(i,k) - z(i,k+1)
          s2(i,k) = dvdz2/(dz**2)
          n2(i,k) = g*2.0*( thv(i,k) - thv(i,k+1))/((thv(i,k) + thv(i,k+1))*dz)
          ri(i,k) = n2(i,k)/s2(i,k)
       end do
    end do

    return
 end subroutine trbintd
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================

 subroutine pblintd(lchnk   ,ncol    ,                            & 1
                    th      ,q       ,z       ,u       ,v       , &
                    ustar   ,obklen  ,kbfs    ,pblh    ,wstar   , &
                    zi      ,cldn    ,ocnfrac)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Diagnose standard PBL variables
! 
! Method: 
! Diagnosis of PBL depth and related variables.  In this case only wstar.
! The PBL depth follows:
!    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).
!
! Modified for boundary layer height diagnosis: Bert Holtslag, june 1994
! >>>>>>>>>  (Use ricr = 0.3 in this formulation)
! 
! Author: B. Stevens (extracted from pbldiff, August 2000)
! 
!-----------------------------------------------------------------------
    implicit none
   !------------------------------Arguments--------------------------------
   !
   ! Input arguments
   !
    integer, intent(in) :: lchnk                     ! chunk identifier
    integer, intent(in) :: ncol                      ! number of atmospheric columns

    real(r8), intent(in)  :: th(pcols,pver)          ! potential temperature [K]
    real(r8), intent(in)  :: q(pcols,pver,pcnst+pnats)     ! specific humidity [kg/kg]
    real(r8), intent(in)  :: z(pcols,pver)           ! height above surface [m]
    real(r8), intent(in)  :: u(pcols,pver)           ! windspeed x-direction [m/s]
    real(r8), intent(in)  :: v(pcols,pver)           ! windspeed y-direction [m/s]
    real(r8), intent(in)  :: ustar(pcols)            ! surface friction velocity [m/s]
    real(r8), intent(in)  :: obklen(pcols)           ! Obukhov length
    real(r8), intent(in)  :: kbfs(pcols)             ! sfc kinematic buoyancy flux [m^2/s^3]
    real(r8), intent(in)  :: zi(pcols,pverp)         ! height above surface [m]
    real(r8), intent(in)  :: cldn(pcols,pver)        ! new cloud fraction
    real(r8), intent(in)  :: ocnfrac(pcols)          ! Land fraction
 
   !
   ! Output arguments
   !
    real(r8), intent(out) :: wstar(pcols)            ! convective sclae velocity [m/s]
    real(r8), intent(out) :: pblh(pcols)             ! boundary-layer height [m]
   !
   !---------------------------Local parameters----------------------------
   !
    real(r8), parameter   :: tiny = 1.e-36           ! lower bound for wind magnitude
    real(r8), parameter   :: fac  = 100.             ! ustar parameter in height diagnosis 
   !
   !---------------------------Local workspace-----------------------------
   !
    integer  :: i                       ! longitude index
    integer  :: k                       ! level index
  
    real(r8) :: thvref(pcols)           ! reference level virtual temperature
    real(r8) :: phiminv(pcols)          ! inverse phi function for momentum
    real(r8) :: phihinv(pcols)          ! inverse phi function for heat
    real(r8) :: rino(pcols,pver)        ! bulk Richardson no. from level to ref lev
    real(r8) :: tlv(pcols)              ! ref. level pot tmp + tmp excess
    real(r8) :: tkv                     ! model level potential temperature
    real(r8) :: vvk                     ! velocity magnitude squared
  
    logical  :: unstbl(pcols)           ! pts w/unstbl pbl (positive virtual ht flx)
    logical  :: check(pcols)            ! True=>chk if Richardson no.>critcal
    logical  :: ocncldcheck(pcols)      ! True=>if ocean surface and cloud in lowest layer
   !
   ! Compute Obukhov length virtual temperature flux and various arrays for use later:
   !
    do i=1,ncol
       check(i)     = .true.
       rino(i,pver) = 0.0
       thvref(i)    = th(i,pver)*(1.0 + 0.61*q(i,pver,1))
       pblh(i)      = z(i,pver)
    end do
   !
   !
   ! PBL height calculation:  Scan upward until the Richardson number between
   ! the first level and the current level exceeds the "critical" value.
   !
    do k=pver-1,pver-npbl+1,-1
       do i=1,ncol
          if (check(i)) then
             vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
             vvk = max(vvk,tiny)
             tkv = th(i,k)*(1. + 0.61*q(i,k,1))
             rino(i,k) = g*(tkv - thvref(i))*(z(i,k)-z(i,pver))/(thvref(i)*vvk)
             if (rino(i,k) >= ricr) then
                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
   !
   ! Estimate an effective surface temperature to account for surface fluctuations
   !
    do i=1,ncol
       if (check(i)) pblh(i) = z(i,pverp-npbl)
       unstbl(i) = (kbfs(i) > 0.)
       check(i)  = (kbfs(i) > 0.)
       if (check(i)) then
          phiminv(i)   = (1. - binm*pblh(i)/obklen(i))**onet
          rino(i,pver) = 0.0
          tlv(i)       = thvref(i) + kbfs(i)*fak/( ustar(i)*phiminv(i) )
       end if
    end do
   !
   ! Improve pblh estimate for unstable conditions using the convective temperature excess:
   !
    do k=pver-1,pver-npbl+1,-1
       do i=1,ncol
          if (check(i)) then
             vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
             vvk = max(vvk,tiny)
             tkv = th(i,k)*(1. + 0.61*q(i,k,1))
             rino(i,k) = g*(tkv - tlv(i))*(z(i,k)-z(i,pver))/(thvref(i)*vvk)
             if (rino(i,k) >= ricr) then
                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
   !
   ! 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.  Also, do not allow 
   ! PBL to exceed some maximum (npbl) number of allowable points
   !
    do i=1,ncol
       if (check(i)) pblh(i) = z(i,pverp-npbl)
       pblh(i) = max(pblh(i),700.0*ustar(i))
       wstar(i) = (max(0._r8,kbfs(i))*g*pblh(i)/thvref(i))**onet
    end do
   !
   ! Final requirement on PBL heightis that it must be greater than the depth
   ! of the lowest model level over ocean if there is any cloud diagnosed in 
   ! the lowest model level.  This is to deal with the inadequacies of the 
   ! current "dry" formulation of the boundary layer, where this test is 
   ! used to identify circumstances where there is marine stratus in the 
   ! lowest level, and to provide a weak ventilation of the layer to avoid
   ! a pathology in the cloud scheme (locking in low-level stratiform cloud)
   ! If over an ocean surface, and any cloud is diagnosed in the 
   ! lowest level, set pblh to 50 meters higher than top interface of lowest level
   !
   !  jrm This is being applied everywhere (not just ocean)!
    do i=1,ncol
       ocncldcheck(i) = .false.
       if (cldn(i,pver).ge.0.0) ocncldcheck(i) = .true.
       if (ocncldcheck(i)) pblh(i) = max(pblh(i),zi(i,pver) + 50.)
    end do
   !
    return
 end subroutine pblintd
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================

  subroutine austausch_atm(lchnk   ,ncol    ,ri      ,s2      ,kvf     ) 1
!----------------------------------------------------------------------- 
! 
! Purpose: 
!  Computes exchange coefficients for free turbulent flows. 
! 
! Method: 
!
! 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)
! 
! Author: B. Stevens (rewrite, August 2000)
! 
!-----------------------------------------------------------------------
    implicit none
   !------------------------------Arguments--------------------------------
   !
   ! Input arguments
   !
    integer, intent(in) :: lchnk                    ! chunk identifier
    integer, intent(in) :: ncol                     ! number of atmospheric columns

    real(r8), intent(in)  ::  s2(pcols,pver)        ! shear squared
    real(r8), intent(in)  ::  ri(pcols,pver)        ! richardson no
   !
   ! Output arguments
   !
    real(r8), intent(out) :: kvf(pcols,pverp)       ! coefficient for heat and tracers
   !
   !---------------------------Local workspace-----------------------------
   !
    real(r8) :: fofri                  ! f(ri)
    real(r8) :: kvn                    ! neutral Kv

    integer  :: i                      ! longitude index
    integer  :: k                      ! vertical index
   !
   !-----------------------------------------------------------------------
   !
   ! The surface diffusivity is always zero
   !
    kvf(:ncol,pverp) = 0.0
   !
   ! Set the vertical diffusion coefficient above the top diffusion level
   ! Note that nbot_turb != pver is not supported
   !
    kvf(:ncol,1:ntop_turb) = 0.0
   !
   ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm. 
   !
    do k = ntop_turb, nbot_turb-1
       do i=1,ncol
          if (ri(i,k) < 0.0) then
             fofri = sqrt(max(1. - 18.*ri(i,k),0._r8))
          else 
             fofri = 1.0/(1.0 + 10.0*ri(i,k)*(1.0 + 8.0*ri(i,k)))    
          end if
          kvn = ml2(k)*sqrt(s2(i,k))
          kvf(i,k+1) = max(zkmin,kvn*fofri)
       end do
    end do

    return
 end subroutine austausch_atm
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================

 subroutine austausch_pbl(lchnk   ,ncol    ,                            & 1
                          z       ,kvf     ,kqfs    ,khfs    ,kbfs    , &
                          obklen  ,ustar   ,wstar   ,pblh    ,kvm     , &
                          kvh     ,cgh     ,cgs     ,tpert   ,qpert   , &
                          ktopbl  ,ktopblmn)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Atmospheric Boundary Layer Computation
! 
! Method: 
! Nonlocal scheme that determines eddy diffusivities based on a
! specified 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).
!
! Author: B. Boville, B. Stevens (rewrite August 2000)
! 
!-----------------------------------------------------------------------
    implicit none
!------------------------------Arguments--------------------------------
!
! Input arguments
!
    integer, intent(in) :: lchnk                    ! chunk identifier
    integer, intent(in) :: ncol                     ! number of atmospheric columns

    real(r8), intent(in) :: z(pcols,pver)           ! height above surface [m]
    real(r8), intent(in) :: kvf(pcols,pverp)        ! free atmospheric eddy diffsvty [m2/s]
    real(r8), intent(in) :: kqfs(pcols,pcnst+pnats) ! kinematic surf cnstituent flux (kg/m2/s)
    real(r8), intent(in) :: khfs(pcols)             ! kinimatic surface heat flux 
    real(r8), intent(in) :: kbfs(pcols)             ! surface buoyancy flux 
    real(r8), intent(in) :: pblh(pcols)             ! boundary-layer height [m]
    real(r8), intent(in) :: obklen(pcols)           ! Obukhov length
    real(r8), intent(in) :: ustar(pcols)            ! surface friction velocity [m/s]
    real(r8), intent(in) :: wstar(pcols)            ! convective velocity scale [m/s]
!
! Output arguments
!
    real(r8), intent(out) :: kvm(pcols,pverp)       ! eddy diffusivity for momentum [m2/s]
    real(r8), intent(out) :: kvh(pcols,pverp)       ! eddy diffusivity for heat [m2/s]
    real(r8), intent(out) :: cgh(pcols,pverp)       ! counter-gradient term for heat [J/kg/m]
    real(r8), intent(out) :: cgs(pcols,pverp)       ! counter-gradient star (cg/flux)
    real(r8), intent(out) :: tpert(pcols)           ! convective temperature excess
    real(r8), intent(out) :: qpert(pcols)           ! convective humidity excess

    integer, intent(out) :: ktopbl(pcols)           ! index of first midpoint inside pbl
    integer, intent(out) :: ktopblmn                ! min value of ktopbl
!
!---------------------------Local workspace-----------------------------
!
    integer  :: i                       ! longitude index
    integer  :: k                       ! level index

    real(r8) :: phiminv(pcols)          ! inverse phi function for momentum
    real(r8) :: phihinv(pcols)          ! inverse phi function for heat
    real(r8) :: wm(pcols)               ! turbulent velocity scale for momentum
    real(r8) :: zp(pcols)               ! current level height + one level up 
    real(r8) :: fak1(pcols)             ! k*ustar*pblh     
    real(r8) :: fak2(pcols)             ! k*wm*pblh
    real(r8) :: fak3(pcols)             ! fakn*wstar/wm
    real(r8) :: pblk(pcols)             ! level eddy diffusivity for momentum
    real(r8) :: pr(pcols)               ! Prandtl number for eddy diffusivities
    real(r8) :: zl(pcols)               ! zmzp / Obukhov length
    real(r8) :: zh(pcols)               ! zmzp / pblh
    real(r8) :: zzh(pcols)              ! (1-(zmzp/pblh))**2
    real(r8) :: zmzp                    ! level height halfway between zm and zp
    real(r8) :: term                    ! intermediate calculation

    logical  :: unstbl(pcols)           ! pts w/unstbl pbl (positive virtual ht flx)
    logical  :: pblpt(pcols)            ! pts within pbl
!
! Initialize height independent arrays
!

!drb initialize variables for runtime error checking
      kvm = 0.	
      kvh = 0.
      cgh = 0.
      cgs = 0.
      tpert = 0.
      qpert = 0.
      ktopbl = 0.
      ktopblmn = 0.

    do i=1,ncol
       unstbl(i) = (kbfs(i) > 0.)
       pblk(i) = 0.0
       fak1(i) = ustar(i)*pblh(i)*vk
       if (unstbl(i)) then
          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
          wm(i)      = ustar(i)*phiminv(i)
          fak2(i)    = wm(i)*pblh(i)*vk
          fak3(i)    = fakn*wstar(i)/wm(i)
          tpert(i)   = max(khfs(i)*fak/wm(i),0._r8)
          qpert(i)   = max(kqfs(i,1)*fak/wm(i),0._r8)
       else
          tpert(i)   = max(khfs(i)*fak/ustar(i),0._r8)
          qpert(i)   = max(kqfs(i,1)*fak/ustar(i),0._r8)
       end if
    end do
!
! Initialize output arrays with free atmosphere values
!
    do k=1,pverp
       do i=1,ncol
          kvm(i,k) = kvf(i,k)
          kvh(i,k) = kvf(i,k)
          cgh(i,k) = 0.0
          cgs(i,k) = 0.0
       end do
    end do
!
! Main level loop to compute the diffusivities and counter-gradient terms. These terms are 
! only calculated at points determined to be in the interior of the pbl (pblpt(i)==.true.),
! and then calculations are directed toward regime: stable vs unstable, surface vs outer 
! layer.
!
    do k=pver,pver-npbl+2,-1
       do i=1,ncol
          pblpt(i) = (z(i,k) < pblh(i))
          if (pblpt(i)) then
             ktopbl(i) = k
             zp(i)  = z(i,k-1)
             if (zkmin == 0.0 .and. zp(i) > pblh(i)) zp(i) = pblh(i)
             zmzp    = 0.5*(z(i,k) + zp(i))
             zh(i)   = zmzp/pblh(i)
             zl(i)   = zmzp/obklen(i)
             zzh(i)  = zh(i)*max(0._r8,(1. - zh(i)))**2
             if (unstbl(i)) then
                if (zh(i) < sffrac) then
                   term     = (1. - betam*zl(i))**onet
                   pblk(i)  = fak1(i)*zzh(i)*term
                   pr(i)    = term/sqrt(1. - betah*zl(i))
                else
                   pblk(i)  = fak2(i)*zzh(i)
                   pr(i)    = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
                   cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
                   cgh(i,k) = khfs(i)*cgs(i,k)*cpair
                end if
             else
                if (zl(i) <= 1.) then
                   pblk(i) = fak1(i)*zzh(i)/(1. + betas*zl(i))
                else
                   pblk(i) = fak1(i)*zzh(i)/(betas + zl(i))
                end if
                pr(i)    = 1.
             end if
             kvm(i,k) = max(pblk(i),kvf(i,k))
             kvh(i,k) = max(pblk(i)/pr(i),kvf(i,k))
          end if
       end do
    end do
!
! Check whether last allowed midpoint is within pbl, determine ktopblmn
!
    ktopblmn = pver
    k = pver-npbl+1
    do i = 1, ncol
       if (z(i,k) < pblh(i)) ktopbl(i) = k
       ktopblmn = min(ktopblmn, ktopbl(i))
    end do

    return
 end subroutine austausch_pbl
END MODULE turbulence