#include #include 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 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, & 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,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 , & 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 , & 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 ) !----------------------------------------------------------------------- ! ! 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 ! !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