turbulence.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! 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)
    real(r8), intent(in)  :: tauy(pcols)             ! surface v stress (N)

   !
   ! 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    ,                            &
                    th      ,q       ,z       ,u       ,v       , &
                    ustar   ,obklen  ,kbfs    ,pblh    ,wstar   )
!----------------------------------------------------------------------- 
! 
! 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]
 
   !
   ! 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
   !
   ! 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

    return
 end subroutine pblintd
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!