turbulence.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>
module turbulence
!---------------------------------------------------------------------------------
! 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 precision, only : r8
  use ppgrid, only    : pver, pverp, pcols
  use pmgrid, only    : masterproc
  use tracers, 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,  &
                    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    ,                            &
                    th      ,t       ,q       ,z       ,zi      , &
                    pmid    ,u       ,v       ,taux    ,tauy    , &
                    shflx   ,cflx    ,obklen  ,ustar   ,pblh    , &
                    kvm     ,kvh     ,cgh     ,cgs     ,kqfs    , &
                    tpert   ,qpert   ,ktopbl  ,ktopblmn)
!----------------------------------------------------------------------- 
! 
! 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
    real(r8), intent(in)  :: tauy(pcols)              ! meridional stress
    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

   !
   ! 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
   !
   ! 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   )
   !
   ! 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)
   !
    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)

    return
 end subroutine trbintr
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===============================================================================
 subroutine trbintd(lchnk   ,ncol    ,                            &
                    th      ,q       ,z       ,u       ,v       , &
                    t       ,pmid    ,cflx    ,shflx   ,taux    , &
                    tauy    ,ustar   ,obklen  ,kqfs    ,khfs    , &
                    kbfs    ,s2      ,n2      ,ri      )
!----------------------------------------------------------------------- 
! 
! Purpose: 
!  Time dependent initialization
!