LCOV - code coverage report
Current view: | top level - physics/cam - hb_diff.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 91 | 171 | 53.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 4 | 6 | 66.7 %

              Line data    Source code
           1             : module hb_diff
           2             :   !---------------------------------------------------------------------------------
           3             :   ! Module to compute mixing coefficients associated with turbulence in the
           4             :   ! planetary boundary layer and elsewhere.  PBL coefficients are based on Holtslag
           5             :   ! and Boville, 1991.
           6             :   !
           7             :   ! Public interfaces:
           8             :   !    init_hb_diff     initializes time independent coefficients
           9             :   !    compute_hb_diff  computes eddy diffusivities and counter-gradient fluxes
          10             :   !
          11             :   ! Private methods:
          12             :   !       trbintd         initializes time dependent variables
          13             :   !       pblintd         initializes time dependent variables that depend pbl depth
          14             :   !       austausch_atm   computes free atmosphere exchange coefficients
          15             :   !       austausch_pbl   computes pbl exchange coefficients
          16             :   !
          17             :   !---------------------------Code history--------------------------------
          18             :   ! Standardized:      J. Rosinski, June 1992
          19             :   ! Reviewed:          P. Rasch, B. Boville, August 1992
          20             :   ! Reviewed:          P. Rasch, April 1996
          21             :   ! Reviewed:          B. Boville, April 1996
          22             :   ! rewritten:         B. Boville, May 2000
          23             :   ! rewritten:         B. Stevens, August 2000
          24             :   ! modularized:       J. McCaa, September 2004
          25             :   !---------------------------------------------------------------------------------
          26             :   use shr_kind_mod, only: r8 => shr_kind_r8
          27             :   use spmd_utils,   only: masterproc          ! output from hb_init should be eliminated
          28             :   use ppgrid,       only: pver, pverp, pcols  ! these should be passed in
          29             :   use cam_logfile,  only: iulog
          30             :
          31             :   implicit none
          32             :   private
          33             :   save
          34             :
          35             :   ! Public interfaces
          36             :   public init_hb_diff
          37             :   public compute_hb_diff
          38             :   public compute_hb_free_atm_diff
          39             :   public pblintd
          40             :   !
          41             :   ! PBL limits
          42             :   !
          43             :   real(r8), parameter :: pblmaxp   = 4.e4_r8        ! pbl max depth in pressure units
          44             :   real(r8), parameter :: zkmin     = 0.01_r8        ! Minimum kneutral*f(ri)
          45             :   !
          46             :   ! PBL Parameters
          47             :   !
          48             :   real(r8), parameter :: onet  = 1._r8/3._r8 ! 1/3 power in wind gradient expression
          49             :   real(r8), parameter :: betam = 15.0_r8  ! Constant in wind gradient expression
          50             :   real(r8), parameter :: betas =  5.0_r8  ! Constant in surface layer gradient expression
          51             :   real(r8), parameter :: betah = 15.0_r8  ! Constant in temperature gradient expression
          52             :   real(r8), parameter :: fakn  =  7.2_r8  ! Constant in turbulent prandtl number
          53             :   real(r8), parameter :: fak   =  8.5_r8  ! Constant in surface temperature excess
          54             :   real(r8), parameter :: ricr  =  0.3_r8  ! Critical richardson number
          55             :   real(r8), parameter :: sffrac=  0.1_r8  ! Surface layer fraction of boundary layer
          56             :   real(r8), parameter :: binm  = betam*sffrac       ! betam * sffrac
          57             :   real(r8), parameter :: binh  = betah*sffrac       ! betah * sffrac
          58             :
          59             :   ! Pbl constants set using values from other parts of code
          60             :
          61             :   real(r8) :: cpair      ! Specific heat of dry air
          62             :   real(r8) :: g          ! Gravitational acceleration
          63             :   real(r8) :: ml2(pverp) ! Mixing lengths squared
          64             :   real(r8) :: vk         ! Von Karman's constant
          65             :   real(r8) :: ccon       ! fak * sffrac * vk
          66             :
          67             :   integer :: npbl       ! Maximum number of levels in pbl from surface
          68             :   integer :: ntop_turb  ! Top level to which turbulent vertical diffusion is applied.
          69             :   integer :: nbot_turb  ! Bottom level to which turbulent vertical diff is applied.
          70             :
          71             : !===============================================================================
          72             : CONTAINS
          73             : !===============================================================================
          74             :
          75        1536 : subroutine init_hb_diff(gravx, cpairx, ntop_eddy, nbot_eddy, pref_mid, &
          76        1536 :                         vkx, eddy_scheme)
          77             :
          78             :    !-----------------------------------------------------------------------
          79             :    !
          80             :    ! Initialize time independent variables of turbulence/pbl package.
          81             :    !
          82             :    !-----------------------------------------------------------------------
          83             :
          84             :    !------------------------------Arguments--------------------------------
          85             :    real(r8), intent(in) :: gravx     ! acceleration of gravity
          86             :    real(r8), intent(in) :: cpairx    ! specific heat of dry air
          87             :    real(r8), intent(in) :: pref_mid(pver)! reference pressures at midpoints
          88             :    real(r8), intent(in) :: vkx       ! Von Karman's constant
          89             :    integer, intent(in)  :: ntop_eddy ! Top level to which eddy vert diff is applied.
          90             :    integer, intent(in)  :: nbot_eddy ! Bottom level to which eddy vert diff is applied.
          91             :    character(len=16),  intent(in) :: eddy_scheme
          92             :
          93             :    !---------------------------Local workspace-----------------------------
          94             :    integer :: k                     ! vertical loop index
          95             :    !-----------------------------------------------------------------------
          96             :
          97             :    ! Basic constants
          98        1536 :    cpair = cpairx
          99        1536 :    g     = gravx
         100        1536 :    vk    = vkx
         101        1536 :    ccon  = fak*sffrac*vk
         102        1536 :    ntop_turb = ntop_eddy
         103        1536 :    nbot_turb = nbot_eddy
         104             :
         105             :    ! Set the square of the mixing lengths.
         106        1536 :    ml2(ntop_turb) = 0._r8
         107      142848 :    do k = ntop_turb+1, nbot_turb
         108      141312 :       ml2(k) = 30.0_r8**2                 ! HB scheme: length scale = 30m
         109      142848 :       if  ( eddy_scheme .eq. 'HBR' ) then
         110           0 :          ml2(k) = 1.0_r8**2               ! HBR scheme: length scale = 1m
         111             :       end if
         112             :    end do
         113        1536 :    ml2(nbot_turb+1) = 0._r8
         114             :
         115             :    ! Limit pbl height to regions below 400 mb
         116             :    ! npbl = max number of levels (from bottom) in pbl
         117             :
         118        1536 :    npbl = 0
         119      144384 :    do k=nbot_turb,ntop_turb,-1
         120      144384 :       if (pref_mid(k) >= pblmaxp) then
         121       43008 :          npbl = npbl + 1
         122             :       end if
         123             :    end do
         124        1536 :    npbl = max(npbl,1)
         125             :
         126        1536 :    if (masterproc) then
         127           2 :       write(iulog,*)'INIT_HB_DIFF: PBL height will be limited to bottom ',npbl, &
         128           4 :          ' model levels. Top is ',pref_mid(pverp-npbl),' pascals'
         129             :    end if
         130             :
         131        1536 : end subroutine init_hb_diff
         132             :
         133             : !===============================================================================
         134             :
         135           0 :   subroutine compute_hb_diff(ncol                  , &
         136             :        th      ,t       ,q       ,z       ,zi      , &
         137             :        pmid    ,u       ,v       ,taux    ,tauy    , &
         138             :        shflx   ,qflx    ,obklen  ,ustar   ,pblh    , &
         139             :        kvm     ,kvh     ,kvq     ,cgh     ,cgs     , &
         140             :        tpert   ,qpert   ,cldn    ,ocnfrac ,tke     , &
         141             :        ri      , &
         142           0 :        eddy_scheme)
         143             :     !-----------------------------------------------------------------------
         144             :     !
         145             :     ! Purpose:
         146             :     !  Interface routines for calcualtion and diatnostics of turbulence related
         147             :     !  coefficients
         148             :     !
         149             :     ! Author: B. Stevens (rewrite August 2000)
         150             :     !
         151             :     !-----------------------------------------------------------------------
         152             :
         153             :     use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm
         154             :
         155             :     !------------------------------Arguments--------------------------------
         156             :     !
         157             :     ! Input arguments
         158             :     !
         159             :     integer, intent(in) :: ncol                       ! number of atmospheric columns
         160             :
         161             :     real(r8), intent(in)  :: th(pcols,pver)           ! potential temperature [K]
         162             :     real(r8), intent(in)  :: t(pcols,pver)            ! temperature (used for density)
         163             :     real(r8), intent(in)  :: q(pcols,pver)            ! specific humidity [kg/kg]
         164             :     real(r8), intent(in)  :: z(pcols,pver)            ! height above surface [m]
         165             :     real(r8), intent(in)  :: zi(pcols,pverp)          ! height above surface [m]
         166             :     real(r8), intent(in)  :: u(pcols,pver)            ! zonal velocity
         167             :     real(r8), intent(in)  :: v(pcols,pver)            ! meridional velocity
         168             :     real(r8), intent(in)  :: taux(pcols)              ! zonal stress [N/m2]
         169             :     real(r8), intent(in)  :: tauy(pcols)              ! meridional stress [N/m2]
         170             :     real(r8), intent(in)  :: shflx(pcols)             ! sensible heat flux
         171             :     real(r8), intent(in)  :: qflx(pcols)              ! water vapor flux
         172             :     real(r8), intent(in)  :: pmid(pcols,pver)         ! midpoint pressures
         173             :     real(r8), intent(in)  :: cldn(pcols,pver)         ! new cloud fraction
         174             :     real(r8), intent(in)  :: ocnfrac(pcols)           ! Land fraction
         175             :     character(len=16), intent(in) :: eddy_scheme
         176             :
         177             :     !
         178             :     ! Output arguments
         179             :     !
         180             :     real(r8), intent(out) :: kvm(pcols,pverp)         ! eddy diffusivity for momentum [m2/s]
         181             :     real(r8), intent(out) :: kvh(pcols,pverp)         ! eddy diffusivity for heat [m2/s]
         182             :     real(r8), intent(out) :: kvq(pcols,pverp)         ! eddy diffusivity for constituents [m2/s]
         183             :     real(r8), intent(out) :: cgh(pcols,pverp)         ! counter-gradient term for heat [J/kg/m]
         184             :     real(r8), intent(out) :: cgs(pcols,pverp)         ! counter-gradient star (cg/flux)
         185             :     real(r8), intent(out) :: tpert(pcols)             ! convective temperature excess
         186             :     real(r8), intent(out) :: qpert(pcols)             ! convective humidity excess
         187             :     real(r8), intent(out) :: ustar(pcols)             ! surface friction velocity [m/s]
         188             :     real(r8), intent(out) :: obklen(pcols)            ! Obukhov length
         189             :     real(r8), intent(out) :: pblh(pcols)              ! boundary-layer height [m]
         190             :     real(r8), intent(out) :: tke(pcols,pverp)         ! turbulent kinetic energy (estimated)
         191             :     real(r8), intent(out) :: ri(pcols,pver)           ! richardson number: n2/s2
         192             :     !
         193             :     !---------------------------Local workspace-----------------------------
         194             :     !
         195             :     real(r8) :: thv(pcols,pver)         ! virtual temperature
         196             :     real(r8) :: rrho(pcols)             ! 1./bottom level density
         197             :     real(r8) :: wstar(pcols)            ! convective velocity scale [m/s]
         198             :     real(r8) :: kqfs(pcols)             ! kinematic surf constituent flux (kg/m2/s)
         199             :     real(r8) :: khfs(pcols)             ! kinimatic surface heat flux
         200             :     real(r8) :: kbfs(pcols)             ! surface buoyancy flux
         201             :     real(r8) :: kvf(pcols,pverp)        ! free atmospheric eddy diffsvty [m2/s]
         202             :     real(r8) :: s2(pcols,pver)          ! shear squared
         203             :     real(r8) :: n2(pcols,pver)          ! brunt vaisaila frequency
         204             :     real(r8) :: bge(pcols)              ! buoyancy gradient enhancment
         205             :     integer  :: ktopbl(pcols)           ! index of first midpoint inside pbl
         206             :     !
         207             :     ! Initialize time dependent variables that do not depend on pbl height
         208             :     !
         209             :
         210             :     ! virtual temperature
         211           0 :     call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:))
         212             :
         213             :     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
         214             :     call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), &
         215           0 :          rrho(:ncol),ustar(:ncol))
         216             :     call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol),  &
         217             :                      shflx(:ncol),   rrho(:ncol),     ustar(:ncol), &
         218             :                      khfs(:ncol),    kqfs(:ncol),     kbfs(:ncol),  &
         219           0 :                      obklen(:ncol))
         220             :     ! Calculate s2, n2, and Richardson number.
         221             :     call trbintd(ncol    ,                            &
         222             :          thv     ,z       ,u       ,v       , &
         223           0 :          s2      ,n2      ,ri      )
         224             :     !
         225             :     ! Initialize time dependent variables that do depend on pbl height
         226             :     !
         227             :     call  pblintd(ncol    ,                            &
         228             :          thv     ,z       ,u       ,v       , &
         229             :          ustar   ,obklen  ,kbfs    ,pblh    ,wstar   , &
         230           0 :          zi      ,cldn    ,ocnfrac ,bge     )
         231             :     !
         232             :     ! Get free atmosphere exchange coefficients
         233             :     !
         234             :     call austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
         235           0 :          ml2, ri, s2, kvf)
         236             :     !
         237             :     ! Get pbl exchange coefficients
         238             :     !
         239             :     call austausch_pbl(ncol                          , &
         240             :          z       ,kvf     ,kqfs    ,khfs    ,kbfs    , &
         241             :          obklen  ,ustar   ,wstar   ,pblh    ,kvm     , &
         242             :          kvh     ,cgh     ,cgs     ,tpert   ,qpert   , &
         243           0 :          ktopbl  ,tke     ,bge     ,eddy_scheme)
         244             :     !
         245             :
         246           0 :     kvq(:ncol,:) = kvh(:ncol,:)
         247           0 :   end subroutine compute_hb_diff
         248             :
         249     1489176 :   subroutine compute_hb_free_atm_diff(ncol,          &
         250             :        th      ,t       ,q       ,z       ,          &
         251             :        pmid    ,u       ,v       ,taux    ,tauy    , &
         252             :        shflx   ,qflx    ,obklen  ,ustar   ,          &
         253             :        kvm     ,kvh     ,kvq     ,cgh     ,cgs     , &
         254             :        ri      )
         255             :     !-----------------------------------------------------------------------
         256             :     !
         257             :     ! This is a version of compute_hb_diff that only computes free
         258             :     ! atmosphere exchange (no PBL computations)
         259             :     !
         260             :     ! Author: B. Stevens (rewrite August 2000)
         261             :     !         Modified by Thomas Toniazzo and Peter H. Lauritzen (June 2023)
         262             :     !
         263             :     !-----------------------------------------------------------------------
         264             :
         265             :     use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm_free
         266             :
         267             :     !------------------------------Arguments--------------------------------
         268             :     !
         269             :     ! Input arguments
         270             :     !
         271             :     integer, intent(in) :: ncol                       ! number of atmospheric columns
         272             :
         273             :     real(r8), intent(in)  :: th(pcols,pver)           ! potential temperature [K]
         274             :     real(r8), intent(in)  :: t(pcols,pver)            ! temperature (used for density)
         275             :     real(r8), intent(in)  :: q(pcols,pver)            ! specific humidity [kg/kg]
         276             :     real(r8), intent(in)  :: z(pcols,pver)            ! height above surface [m]
         277             :     real(r8), intent(in)  :: u(pcols,pver)            ! zonal velocity
         278             :     real(r8), intent(in)  :: v(pcols,pver)            ! meridional velocity
         279             :     real(r8), intent(in)  :: taux(pcols)              ! zonal stress [N/m2]
         280             :     real(r8), intent(in)  :: tauy(pcols)              ! meridional stress [N/m2]
         281             :     real(r8), intent(in)  :: shflx(pcols)             ! sensible heat flux
         282             :     real(r8), intent(in)  :: qflx(pcols)              ! water vapor flux
         283             :     real(r8), intent(in)  :: pmid(pcols,pver)         ! midpoint pressures
         284             :     !
         285             :     ! Output arguments
         286             :     !
         287             :     real(r8), intent(out) :: kvm(pcols,pverp)         ! eddy diffusivity for momentum [m2/s]
         288             :     real(r8), intent(out) :: kvh(pcols,pverp)         ! eddy diffusivity for heat [m2/s]
         289             :     real(r8), intent(out) :: kvq(pcols,pverp)         ! eddy diffusivity for constituents [m2/s]
         290             :     real(r8), intent(out) :: cgh(pcols,pverp)         ! counter-gradient term for heat [J/kg/m]
         291             :     real(r8), intent(out) :: cgs(pcols,pverp)         ! counter-gradient star (cg/flux)
         292             :     real(r8), intent(out) :: ustar(pcols)             ! surface friction velocity [m/s]
         293             :     real(r8), intent(out) :: obklen(pcols)            ! Obukhov length
         294             :     real(r8), intent(out) :: ri(pcols,pver)           ! richardson number: n2/s2
         295             :     !
         296             :     !---------------------------Local workspace-----------------------------
         297             :     !
         298             :     real(r8) :: thv(pcols,pver)         ! virtual potential temperature
         299             :     real(r8) :: rrho(pcols)             ! 1./bottom level density
         300             :     real(r8) :: kqfs(pcols)             ! kinematic surf constituent flux (kg/m2/s)
         301             :     real(r8) :: khfs(pcols)             ! kinimatic surface heat flux
         302             :     real(r8) :: kbfs(pcols)             ! surface buoyancy flux
         303             :     real(r8) :: kvf(pcols,pverp)        ! free atmospheric eddy diffsvty [m2/s]
         304             :     real(r8) :: s2(pcols,pver)          ! shear squared
         305             :     real(r8) :: n2(pcols,pver)          ! brunt vaisaila frequency
         306             :
         307             :     ! virtual potential temperature
         308  1628350464 :     call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:))
         309             :
         310             :     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
         311             :     call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), &
         312     1489176 :          rrho(:ncol),ustar(:ncol))
         313             :     call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol),  &
         314             :                      shflx(:ncol),   rrho(:ncol),     ustar(:ncol), &
         315             :                      khfs(:ncol),    kqfs(:ncol),     kbfs(:ncol),  &
         316     1489176 :                      obklen(:ncol))
         317             :     ! Calculate s2, n2, and Richardson number.
         318             :     call trbintd(ncol    ,                            &
         319             :          thv     ,z       ,u       ,v       , &
         320     1489176 :          s2      ,n2      ,ri      )
         321             :     !
         322             :     ! Get free atmosphere exchange coefficients
         323             :     !
         324             :     call austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
         325     1489176 :          ml2, ri, s2, kvf)
         326             :
         327  2338872120 :     kvq(:ncol,:) = kvf(:ncol,:)
         328  2338872120 :     kvm(:ncol,:) = kvf(:ncol,:)
         329  2338872120 :     kvh(:ncol,:) = kvf(:ncol,:)
         330  2338872120 :     cgh(:ncol,:) = 0._r8
         331  2338872120 :     cgs(:ncol,:) = 0._r8
         332             :
         333     1489176 :   end subroutine compute_hb_free_atm_diff
         334             :
         335             :
         336             :   !
         337             :   !===============================================================================
         338     1489176 :   subroutine trbintd(ncol    ,                            &
         339             :        thv     ,z       ,u       ,v       , &
         340             :        s2      ,n2      ,ri      )
         341             :
         342             :     !-----------------------------------------------------------------------
         343             :     !
         344             :     ! Purpose:
         345             :     !  Time dependent initialization
         346             :     !
         347             :     ! Method:
         348             :     !  Diagnosis of variables that do not depend on mixing assumptions or
         349             :     !  PBL depth
         350             :     !
         351             :     ! Author: B. Stevens (extracted from pbldiff, August, 2000)
         352             :     !
         353             :     !-----------------------------------------------------------------------
         354             :     !------------------------------Arguments--------------------------------
         355             :     !
         356             :     ! Input arguments
         357             :     !
         358             :     integer, intent(in) :: ncol                      ! number of atmospheric columns
         359             :
         360             :     real(r8), intent(in)  :: thv(pcols,pver)         ! virtual temperature
         361             :     real(r8), intent(in)  :: z(pcols,pver)           ! height above surface [m]
         362             :     real(r8), intent(in)  :: u(pcols,pver)           ! windspeed x-direction [m/s]
         363             :     real(r8), intent(in)  :: v(pcols,pver)           ! windspeed y-direction [m/s]
         364             :
         365             :     !
         366             :     ! Output arguments
         367             :     !
         368             :     real(r8), intent(out) :: s2(pcols,pver)          ! shear squared
         369             :     real(r8), intent(out) :: n2(pcols,pver)          ! brunt vaisaila frequency
         370             :     real(r8), intent(out) :: ri(pcols,pver)          ! richardson number: n2/s2
         371             :     !
         372             :     !---------------------------Local workspace-----------------------------
         373             :     !
         374             :     integer  :: i                        ! longitude index
         375             :     integer  :: k                        ! level index
         376             :
         377             :     real(r8) :: vvk                      ! velocity magnitude squared
         378             :     real(r8) :: dvdz2                    ! velocity shear squared
         379             :     real(r8) :: dz                       ! delta z between midpoints
         380             :     !
         381             :     ! Compute shear squared (s2), brunt vaisaila frequency (n2) and related richardson
         382             :     ! number (ri). Use virtual temperature to compute n2.
         383             :     !
         384             :
         385   138493368 :     do k=ntop_turb,nbot_turb-1
         386  2289140568 :        do i=1,ncol
         387  2150647200 :           dvdz2   = (u(i,k)-u(i,k+1))**2 + (v(i,k)-v(i,k+1))**2
         388  2150647200 :           dvdz2   = max(dvdz2,1.e-36_r8)
         389  2150647200 :           dz      = z(i,k) - z(i,k+1)
         390  2150647200 :           s2(i,k) = dvdz2/(dz**2)
         391  2150647200 :           n2(i,k) = g*2.0_r8*( thv(i,k) - thv(i,k+1))/((thv(i,k) + thv(i,k+1))*dz)
         392  2287651392 :           ri(i,k) = n2(i,k)/s2(i,k)
         393             :        end do
         394             :     end do
         395             :
         396     1489176 :     return
         397             :   end subroutine trbintd
         398             :   !
         399             :   !===============================================================================
         400     4467528 :   subroutine pblintd(ncol    ,                            &
         401             :        thv     ,z       ,u       ,v       , &
         402             :        ustar   ,obklen  ,kbfs    ,pblh    ,wstar   , &
         403             :        zi      ,cldn    ,ocnfrac ,bge     )
         404             :     !-----------------------------------------------------------------------
         405             :     !
         406             :     ! Purpose:
         407             :     ! Diagnose standard PBL variables
         408             :     !
         409             :     ! Method:
         410             :     ! Diagnosis of PBL depth and related variables.  In this case only wstar.
         411             :     ! The PBL depth follows:
         412             :     !    Holtslag, A.A.M., and B.A. Boville, 1993:
         413             :     !    Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate
         414             :     !    Model. J. Clim., vol. 6., p. 1825--1842.
         415             :     !
         416             :     ! Updated by Holtslag and Hack to exclude the surface layer from the
         417             :     ! definition of the boundary layer Richardson number. Ri is now defined
         418             :     ! across the outer layer of the pbl (between the top of the surface
         419             :     ! layer and the pbl top) instead of the full pbl (between the surface and
         420             :     ! the pbl top). For simiplicity, the surface layer is assumed to be the
         421             :     ! region below the first model level (otherwise the boundary layer depth
         422             :     ! determination would require iteration).
         423             :     !
         424             :     ! Modified for boundary layer height diagnosis: Bert Holtslag, june 1994
         425             :     ! >>>>>>>>>  (Use ricr = 0.3 in this formulation)
         426             :     !
         427             :     ! Author: B. Stevens (extracted from pbldiff, August 2000)
         428             :     !
         429             :     !-----------------------------------------------------------------------
         430             :     !------------------------------Arguments--------------------------------
         431             :     !
         432             :     ! Input arguments
         433             :     !
         434             :     integer, intent(in) :: ncol                      ! number of atmospheric columns
         435             :
         436             :     real(r8), intent(in)  :: thv(pcols,pver)         ! virtual temperature
         437             :     real(r8), intent(in)  :: z(pcols,pver)           ! height above surface [m]
         438             :     real(r8), intent(in)  :: u(pcols,pver)           ! windspeed x-direction [m/s]
         439             :     real(r8), intent(in)  :: v(pcols,pver)           ! windspeed y-direction [m/s]
         440             :     real(r8), intent(in)  :: ustar(pcols)            ! surface friction velocity [m/s]
         441             :     real(r8), intent(in)  :: obklen(pcols)           ! Obukhov length
         442             :     real(r8), intent(in)  :: kbfs(pcols)             ! sfc kinematic buoyancy flux [m^2/s^3]
         443             :     real(r8), intent(in)  :: zi(pcols,pverp)         ! height above surface [m]
         444             :     real(r8), intent(in)  :: cldn(pcols,pver)        ! new cloud fraction
         445             :     real(r8), intent(in)  :: ocnfrac(pcols)          ! Land fraction
         446             :
         447             :     !
         448             :     ! Output arguments
         449             :     !
         450             :     real(r8), intent(out) :: wstar(pcols)            ! convective sclae velocity [m/s]
         451             :     real(r8), intent(out) :: pblh(pcols)             ! boundary-layer height [m]
         452             :     real(r8), intent(out) :: bge(pcols)              ! buoyancy gradient enhancment
         453             :     !
         454             :     !---------------------------Local parameters----------------------------
         455             :     !
         456             :     real(r8), parameter   :: tiny = 1.e-36_r8           ! lower bound for wind magnitude
         457             :     real(r8), parameter   :: fac  = 100._r8             ! ustar parameter in height diagnosis
         458             :     !
         459             :     !---------------------------Local workspace-----------------------------
         460             :     !
         461             :     integer  :: i                       ! longitude index
         462             :     integer  :: k                       ! level index
         463             :
         464             :     real(r8) :: phiminv(pcols)          ! inverse phi function for momentum
         465             :     real(r8) :: phihinv(pcols)          ! inverse phi function for heat
         466             :     real(r8) :: rino(pcols,pver)        ! bulk Richardson no. from level to ref lev
         467             :     real(r8) :: tlv(pcols)              ! ref. level pot tmp + tmp excess
         468             :     real(r8) :: vvk                     ! velocity magnitude squared
         469             :
         470             :     logical  :: unstbl(pcols)           ! pts w/unstbl pbl (positive virtual ht flx)
         471             :     logical  :: check(pcols)            ! True=>chk if Richardson no.>critcal
         472             :     logical  :: ocncldcheck(pcols)      ! True=>if ocean surface and cloud in lowest layer
         473             :     !
         474             :     ! Compute Obukhov length virtual temperature flux and various arrays for use later:
         475             :     !
         476    74597328 :     do i=1,ncol
         477    70129800 :        check(i)     = .true.
         478    70129800 :        rino(i,pver) = 0.0_r8
         479    74597328 :        pblh(i)      = z(i,pver)
         480             :     end do
         481             :     !
         482             :     !
         483             :     ! PBL height calculation:  Scan upward until the Richardson number between
         484             :     ! the first level and the current level exceeds the "critical" value.
         485             :     !
         486   125090784 :     do k=pver-1,pver-npbl+1,-1
         487  2018595384 :        do i=1,ncol
         488  2014127856 :           if (check(i)) then
         489   464947677 :              vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
         490   464947677 :              vvk = max(vvk,tiny)
         491   464947677 :              rino(i,k) = g*(thv(i,k) - thv(i,pver))*(z(i,k)-z(i,pver))/(thv(i,pver)*vvk)
         492   464947677 :              if (rino(i,k) >= ricr) then
         493    70129800 :                 pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/(rino(i,k) - rino(i,k+1)) * &
         494    70129800 :                      (z(i,k) - z(i,k+1))
         495    70129800 :                 check(i) = .false.
         496             :              end if
         497             :           end if
         498             :        end do
         499             :     end do
         500             :     !
         501             :     ! Estimate an effective surface temperature to account for surface fluctuations
         502             :     !
         503    74597328 :     do i=1,ncol
         504    70129800 :        if (check(i)) pblh(i) = z(i,pverp-npbl)
         505    70129800 :        unstbl(i) = (kbfs(i) > 0._r8)
         506    70129800 :        check(i)  = (kbfs(i) > 0._r8)
         507    74597328 :        if (check(i)) then
         508    49815090 :           phiminv(i)   = (1._r8 - binm*pblh(i)/obklen(i))**onet
         509    49815090 :           rino(i,pver) = 0.0_r8
         510    49815090 :           tlv(i)       = thv(i,pver) + kbfs(i)*fak/( ustar(i)*phiminv(i) )
         511             :        end if
         512             :     end do
         513             :     !
         514             :     ! Improve pblh estimate for unstable conditions using the convective temperature excess:
         515             :     !
         516    74597328 :     do i = 1,ncol
         517    74597328 :        bge(i) = 1.e-8_r8
         518             :     end do
         519   125090784 :     do k=pver-1,pver-npbl+1,-1
         520  2018595384 :        do i=1,ncol
         521  2014127856 :           if (check(i)) then
         522   443316541 :              vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
         523   443316541 :              vvk = max(vvk,tiny)
         524   443316541 :              rino(i,k) = g*(thv(i,k) - tlv(i))*(z(i,k)-z(i,pver))/(thv(i,pver)*vvk)
         525   443316541 :              if (rino(i,k) >= ricr) then
         526    49815090 :                 pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/(rino(i,k) - rino(i,k+1))* &
         527    49815090 :                      (z(i,k) - z(i,k+1))
         528    49815090 :                 bge(i) = 2._r8*g/(thv(i,k)+thv(i,k+1))*(thv(i,k)-thv(i,k+1))/(z(i,k)-z(i,k+1))*pblh(i)
         529    49815090 :                 if (bge(i).lt.0._r8) then
         530        2201 :                    bge(i) = 1.e-8_r8
         531             :                 endif
         532    49815090 :                 check(i) = .false.
         533             :              end if
         534             :           end if
         535             :        end do
         536             :     end do
         537             :     !
         538             :     ! PBL height must be greater than some minimum mechanical mixing depth
         539             :     ! Several investigators have proposed minimum mechanical mixing depth
         540             :     ! relationships as a function of the local friction velocity, u*.  We
         541             :     ! make use of a linear relationship of the form h = c u* where c=700.
         542             :     ! The scaling arguments that give rise to this relationship most often
         543             :     ! represent the coefficient c as some constant over the local coriolis
         544             :     ! parameter.  Here we make use of the experimental results of Koracin
         545             :     ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
         546             :     ! where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
         547             :     ! latitude value for f so that c = 0.07/f = 700.  Also, do not allow
         548             :     ! PBL to exceed some maximum (npbl) number of allowable points
         549             :     !
         550    74597328 :     do i=1,ncol
         551    70129800 :        if (check(i)) pblh(i) = z(i,pverp-npbl)
         552    70129800 :        pblh(i) = max(pblh(i),700.0_r8*ustar(i))
         553    74597328 :        wstar(i) = (max(0._r8,kbfs(i))*g*pblh(i)/thv(i,pver))**onet
         554             :     end do
         555             :     !
         556             :     ! Final requirement on PBL heightis that it must be greater than the depth
         557             :     ! of the lowest model level over ocean if there is any cloud diagnosed in
         558             :     ! the lowest model level.  This is to deal with the inadequacies of the
         559             :     ! current "dry" formulation of the boundary layer, where this test is
         560             :     ! used to identify circumstances where there is marine stratus in the
         561             :     ! lowest level, and to provide a weak ventilation of the layer to avoid
         562             :     ! a pathology in the cloud scheme (locking in low-level stratiform cloud)
         563             :     ! If over an ocean surface, and any cloud is diagnosed in the
         564             :     ! lowest level, set pblh to 50 meters higher than top interface of lowest level
         565             :     !
         566             :     !  jrm This is being applied everywhere (not just ocean)!
         567    74597328 :     do i=1,ncol
         568    70129800 :        ocncldcheck(i) = .false.
         569    70129800 :        if (cldn(i,pver).ge.0.0_r8) ocncldcheck(i) = .true.
         570    74597328 :        if (ocncldcheck(i)) pblh(i) = max(pblh(i),zi(i,pver) + 50._r8)
         571             :     end do
         572             :     !
         573     4467528 :     return
         574             :   end subroutine pblintd
         575             :   !
         576             :   !===============================================================================
         577           0 :   subroutine austausch_pbl(ncol                    , &
         578             :        z       ,kvf     ,kqfs    ,khfs    ,kbfs    , &
         579             :        obklen  ,ustar   ,wstar   ,pblh    ,kvm     , &
         580             :        kvh     ,cgh     ,cgs     ,tpert   ,qpert   , &
         581           0 :        ktopbl  ,tke     ,bge     ,eddy_scheme)
         582             :     !-----------------------------------------------------------------------
         583             :     !
         584             :     ! Purpose:
         585             :     ! Atmospheric Boundary Layer Computation
         586             :     !
         587             :     ! Method:
         588             :     ! Nonlocal scheme that determines eddy diffusivities based on a
         589             :     ! specified boundary layer height and a turbulent velocity scale;
         590             :     ! also, countergradient effects for heat and moisture, and constituents
         591             :     ! are included, along with temperature and humidity perturbations which
         592             :     ! measure the strength of convective thermals in the lower part of the
         593             :     ! atmospheric boundary layer.
         594             :     !
         595             :     ! For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
         596             :     ! Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate
         597             :     ! Model. J. Clim., vol. 6., p. 1825--1842.
         598             :     !
         599             :     ! Updated by Holtslag and Hack to exclude the surface layer from the
         600             :     ! definition of the boundary layer Richardson number. Ri is now defined
         601             :     ! across the outer layer of the pbl (between the top of the surface
         602             :     ! layer and the pbl top) instead of the full pbl (between the surface and
         603             :     ! the pbl top). For simiplicity, the surface layer is assumed to be the
         604             :     ! region below the first model level (otherwise the boundary layer depth
         605             :     ! determination would require iteration).
         606             :     !
         607             :     ! Author: B. Boville, B. Stevens (rewrite August 2000)
         608             :     !
         609             :     !------------------------------Arguments--------------------------------
         610             :     !
         611             :     ! Input arguments
         612             :     !
         613             :     integer, intent(in) :: ncol                     ! number of atmospheric columns
         614             :
         615             :     real(r8), intent(in) :: z(pcols,pver)           ! height above surface [m]
         616             :     real(r8), intent(in) :: kvf(pcols,pverp)        ! free atmospheric eddy diffsvty [m2/s]
         617             :     real(r8), intent(in) :: kqfs(pcols)             ! kinematic surf cnstituent flux (kg/m2/s)
         618             :     real(r8), intent(in) :: khfs(pcols)             ! kinimatic surface heat flux
         619             :     real(r8), intent(in) :: kbfs(pcols)             ! surface buoyancy flux
         620             :     real(r8), intent(in) :: pblh(pcols)             ! boundary-layer height [m]
         621             :     real(r8), intent(in) :: obklen(pcols)           ! Obukhov length
         622             :     real(r8), intent(in) :: ustar(pcols)            ! surface friction velocity [m/s]
         623             :     real(r8), intent(in) :: wstar(pcols)            ! convective velocity scale [m/s]
         624             :     real(r8), intent(in) :: bge(pcols)              ! buoyancy gradient enhancment
         625             :     character(len=16), intent(in) :: eddy_scheme
         626             :
         627             :     !
         628             :     ! Output arguments
         629             :     !
         630             :     real(r8), intent(out) :: kvm(pcols,pverp)       ! eddy diffusivity for momentum [m2/s]
         631             :     real(r8), intent(out) :: kvh(pcols,pverp)       ! eddy diffusivity for heat [m2/s]
         632             :     real(r8), intent(out) :: cgh(pcols,pverp)       ! counter-gradient term for heat [J/kg/m]
         633             :     real(r8), intent(out) :: cgs(pcols,pverp)       ! counter-gradient star (cg/flux)
         634             :     real(r8), intent(out) :: tpert(pcols)           ! convective temperature excess
         635             :     real(r8), intent(out) :: qpert(pcols)           ! convective humidity excess
         636             :
         637             :     integer,  intent(out) :: ktopbl(pcols)          ! index of first midpoint inside pbl
         638             :     real(r8), intent(out) :: tke(pcols,pverp)       ! turbulent kinetic energy (estimated)
         639             :     !
         640             :     !---------------------------Local workspace-----------------------------
         641             :     !
         642             :     integer  :: i                       ! longitude index
         643             :     integer  :: k                       ! level index
         644             :
         645             :     real(r8) :: phiminv(pcols)          ! inverse phi function for momentum
         646             :     real(r8) :: phihinv(pcols)          ! inverse phi function for heat
         647             :     real(r8) :: wm(pcols)               ! turbulent velocity scale for momentum
         648             :     real(r8) :: zp(pcols)               ! current level height + one level up
         649             :     real(r8) :: fak1(pcols)             ! k*ustar*pblh
         650             :     real(r8) :: fak2(pcols)             ! k*wm*pblh
         651             :     real(r8) :: fak3(pcols)             ! fakn*wstar/wm
         652             :     real(r8) :: pblk(pcols)             ! level eddy diffusivity for momentum
         653             :     real(r8) :: pr(pcols)               ! Prandtl number for eddy diffusivities
         654             :     real(r8) :: zl(pcols)               ! zmzp / Obukhov length
         655             :     real(r8) :: zh(pcols)               ! zmzp / pblh
         656             :     real(r8) :: zzh(pcols)              ! (1-(zmzp/pblh))**2
         657             :     real(r8) :: zmzp                    ! level height halfway between zm and zp
         658             :     real(r8) :: term                    ! intermediate calculation
         659             :     real(r8) :: kve                     ! diffusivity at entrainment layer in unstable cases
         660             :
         661             :     logical  :: unstbl(pcols)           ! pts w/unstbl pbl (positive virtual ht flx)
         662             :     logical  :: pblpt(pcols)            ! pts within pbl
         663             :     !
         664             :     ! Initialize height independent arrays
         665             :     !
         666             :
         667             :     !drb initialize variables for runtime error checking
         668           0 :     kvm = 0._r8
         669           0 :     kvh = 0._r8
         670           0 :     kve = 0._r8
         671           0 :     cgh = 0._r8
         672           0 :     cgs = 0._r8
         673           0 :     tpert = 0._r8
         674           0 :     qpert = 0._r8
         675           0 :     ktopbl = 0._r8
         676           0 :     tke = 0._r8
         677             :
         678           0 :     do i=1,ncol
         679           0 :        unstbl(i) = (kbfs(i) > 0._r8)
         680           0 :        pblk(i) = 0.0_r8
         681           0 :        fak1(i) = ustar(i)*pblh(i)*vk
         682           0 :        if (unstbl(i)) then
         683           0 :           phiminv(i) = (1._r8 - binm*pblh(i)/obklen(i))**onet
         684           0 :           phihinv(i) = sqrt(1._r8 - binh*pblh(i)/obklen(i))
         685           0 :           wm(i)      = ustar(i)*phiminv(i)
         686           0 :           fak2(i)    = wm(i)*pblh(i)*vk
         687           0 :           fak3(i)    = fakn*wstar(i)/wm(i)
         688           0 :           tpert(i)   = max(khfs(i)*fak/wm(i),0._r8)
         689           0 :           qpert(i)   = max(kqfs(i)*fak/wm(i),0._r8)
         690             :        else
         691           0 :           tpert(i)   = max(khfs(i)*fak/ustar(i),0._r8)
         692           0 :           qpert(i)   = max(kqfs(i)*fak/ustar(i),0._r8)
         693             :        end if
         694             :     end do
         695             :     !
         696             :     ! Initialize output arrays with free atmosphere values
         697             :     !
         698           0 :     do k=1,pverp
         699           0 :        do i=1,ncol
         700           0 :           kvm(i,k) = kvf(i,k)
         701           0 :           kvh(i,k) = kvf(i,k)
         702           0 :           cgh(i,k) = 0.0_r8
         703           0 :           cgs(i,k) = 0.0_r8
         704             :        end do
         705             :     end do
         706             :     !
         707             :     ! Main level loop to compute the diffusivities and counter-gradient terms. These terms are
         708             :     ! only calculated at points determined to be in the interior of the pbl (pblpt(i)==.true.),
         709             :     ! and then calculations are directed toward regime: stable vs unstable, surface vs outer
         710             :     ! layer.
         711             :     !
         712           0 :     do k=pver,pver-npbl+2,-1
         713           0 :        do i=1,ncol
         714           0 :           pblpt(i) = (z(i,k) < pblh(i))
         715           0 :           if (pblpt(i)) then
         716           0 :              ktopbl(i) = k
         717           0 :              zp(i)  = z(i,k-1)
         718             :              if (zkmin == 0.0_r8 .and. zp(i) > pblh(i)) zp(i) = pblh(i)
         719           0 :              zmzp    = 0.5_r8*(z(i,k) + zp(i)) ! we think this is an approximation to the interface height (where KVs are calculated)
         720           0 :              zh(i)   = zmzp/pblh(i)
         721           0 :              zl(i)   = zmzp/obklen(i)
         722           0 :              zzh(i)  = zh(i)*max(0._r8,(1._r8 - zh(i)))**2
         723           0 :              if (unstbl(i)) then
         724           0 :                 if (zh(i) < sffrac) then
         725           0 :                    term     = (1._r8 - betam*zl(i))**onet
         726           0 :                    pblk(i)  = fak1(i)*zzh(i)*term
         727           0 :                    pr(i)    = term/sqrt(1._r8 - betah*zl(i))
         728             :                 else
         729           0 :                    pblk(i)  = fak2(i)*zzh(i)
         730           0 :                    pr(i)    = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
         731           0 :                    cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
         732           0 :                    cgh(i,k) = khfs(i)*cgs(i,k)*cpair
         733             :                 end if
         734             :              else
         735           0 :                 if (zl(i) <= 1._r8) then
         736           0 :                    pblk(i) = fak1(i)*zzh(i)/(1._r8 + betas*zl(i))
         737             :                 else
         738           0 :                    pblk(i) = fak1(i)*zzh(i)/(betas + zl(i))
         739             :                 end if
         740           0 :                 pr(i)    = 1._r8
         741             :              end if
         742           0 :              kvm(i,k) = max(pblk(i),kvf(i,k))
         743           0 :              kvh(i,k) = max(pblk(i)/pr(i),kvf(i,k))
         744             :           end if
         745             :        end do
         746             :     end do
         747             :
         748             :     !
         749             :     ! Check whether last allowed midpoint is within pbl
         750             :     !
         751             :
         752           0 :     if  ( eddy_scheme .eq. 'HBR' ) then
         753             :        ! apply new diffusivity at entrainment zone
         754           0 :        do i = 1,ncol
         755           0 :           if (bge(i) > 1.e-7_r8) then
         756           0 :              k = ktopbl(i)
         757           0 :              kve = 0.2_r8*(wstar(i)**3+5._r8*ustar(i)**3)/bge(i)
         758           0 :              kvm(i,k) = kve
         759           0 :              kvh(i,k) = kve
         760             :           end if
         761             :        end do
         762             :     end if
         763             :
         764             :     ! Crude estimate of tke (tke=0 above boundary layer)
         765           0 :     do k = max(pverp-npbl,2),pverp
         766           0 :        do i = 1, ncol
         767           0 :           if (z(i,k-1) < pblh(i)) then
         768           0 :              tke(i,k) = ( kvm(i,k) / pblh(i) ) ** 2
         769             :           endif
         770             :        end do
         771             :     end do
         772           0 :     return
         773             :   end subroutine austausch_pbl
         774             :
         775             : end module hb_diff

Generated by: LCOV version 1.14