turbulence.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!===============================================================================
  subroutine austausch_atm(lchnk   ,ncol    ,ri      ,s2      ,kvf     )
!----------------------------------------------------------------------- 
! 
! 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    ,                            &
                          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
!
    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