#include <misc.h>
#include <params.h>
module zm_conv 2,6
!---------------------------------------------------------------------------------
! Purpose:
!
! Interface from Zhang-McFarlane convection scheme, includes evaporation of convective
! precip from the ZM scheme
!
! Author: Byron Boville, from code in tphysbc
!
!---------------------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use ppgrid, only: pcols, pver, pverp
use moistconvection, only: cp, grav, rgrav, rgas, limcnv
use cldwat, only: cldwat_fice, conke
use history, only: outfld, addfld, add_default, phys_decomp
use physconst, only: cpair, epsilo, gravit, latice, latvap, tmelt
implicit none
save
private ! Make default type private to the module
!
! PUBLIC: interfaces
!
public zm_convi ! ZM schemea
public zm_convr ! ZM schemea
public zm_conv_evap ! evaporation of precip from ZM schemea
!
! Private data
!
public rl, cpres, capelmt
real(r8) rl ! wg latent heat of vaporization.
real(r8) cpres ! specific heat at constant pressure in j/kg-degk.
real(r8), parameter :: capelmt = 70. ! threshold value for cape for deep convection.
real(r8) :: ke ! Tunable evaporation efficiency
real(r8) c0
real(r8) tau ! convective time scale
real(r8),parameter :: a = 21.656
real(r8),parameter :: b = 5418.
real(r8),parameter :: c1 = 6.112
real(r8),parameter :: c2 = 17.67
real(r8),parameter :: c3 = 243.5
real(r8) :: tfreez
real(r8) :: eps1
contains
subroutine zm_convi( x1, x2, x3, x4 ) 1,12
use dycore, only: dycore_is, get_resolution
use pmgrid, only: masterproc
!
! Initialization of ZM constants
!
real(r8), intent(in) :: x1 ! unused
real(r8), intent(in) :: x2 ! unused
real(r8), intent(in) :: x3 ! unused
real(r8), intent(in) :: x4 ! unused
tfreez = tmelt
eps1 = epsilo
rl = latvap
cpres = cpair
! tau=4800. were used in canadian climate center. however, in echam3 t42,
! convection is too weak, thus adjusted to 2400.
if ( dycore_is ('LR') ) then
if ( get_resolution() == '1x1.25' ) then
if(masterproc) write(6,*) 'ZM: Found LR dycore at 1x1.25'
tau = 3600.
c0 = 3.5E-3
ke = 1.0E-6
elseif ( get_resolution() == '4x5' ) then
if(masterproc) write(6,*) 'ZM: Found LR dycore at 4x5'
tau = 3600.
c0 = 3.5E-3
ke = 1.0E-6
else
if(masterproc) write(6,*) 'ZM: Found LR dycore at default resolution'
tau = 3600.
c0 = 3.5E-3
ke = 1.0E-6
endif
else
if(get_resolution() == 'T85')then
tau = 3600.
c0 = 4.E-3
ke = 1.0E-6
if(masterproc) write(6,*) 'ZM: Found spectral dycore at T85 resolution'
elseif(get_resolution() == 'T31')then
tau = 3600.
c0 = 2.E-3
ke = 3.0E-6
if(masterproc) write(6,*) 'ZM: Found spectral dycore at non-T85 resolution'
else
tau = 3600.
c0 = 3.E-3
ke = 3.0E-6
if(masterproc) write(6,*) 'ZM: Found spectral dycore at default resolution'
endif
endif
!-----------------------------------------------------------------------
call addfld ('ZMFLXPRC','kg/m2/s ',pverp, 'A','Flux of precipitation from ZM convection' ,phys_decomp)
call addfld ('ZMFLXSNW','kg/m2/s ',pverp, 'A','Flux of snow from ZM convection' ,phys_decomp)
call addfld ('ZMNTPRPD','kg/kg/s ',pver , 'A','Net precipitation production from ZM convection',phys_decomp)
call addfld ('ZMNTSNPD','kg/kg/s ',pver , 'A','Net snow production from ZM convection' ,phys_decomp)
call addfld ('ZMEIHEAT','W/kg' ,pver , 'A','Heating by ice and evaporation in ZM convection',phys_decomp)
call addfld ('HKFLXPRC','kg/m2/s ',pverp, 'A','Flux of precipitation from HK convection' ,phys_decomp)
call addfld ('HKFLXSNW','kg/m2/s ',pverp, 'A','Flux of snow from HK convection' ,phys_decomp)
call addfld ('HKNTPRPD','kg/kg/s ',pver , 'A','Net precipitation production from HK convection',phys_decomp)
call addfld ('HKNTSNPD','kg/kg/s ',pver , 'A','Net snow production from HK convection' ,phys_decomp)
call addfld ('HKEIHEAT','W/kg' ,pver , 'A','Heating by ice and evaporation in HK convection',phys_decomp)
! These fields removed 10/30/2003 per CRB decision.
! call add_default ('ZMFLXPRC', 1, ' ')
! call add_default ('ZMFLXSNW', 1, ' ')
! call add_default ('ZMNTPRPD', 1, ' ')
! call add_default ('ZMNTSNPD', 1, ' ')
! call add_default ('ZMEIHEAT', 1, ' ')
! call add_default ('HKFLXPRC', 1, ' ')
! call add_default ('HKFLXSNW', 1, ' ')
! call add_default ('HKNTPRPD', 1, ' ')
! call add_default ('HKNTSNPD', 1, ' ')
! call add_default ('HKEIHEAT', 1, ' ')
end subroutine zm_convi
subroutine zm_convr(lchnk ,ncol , & 1,5
t ,qh ,prec ,jctop ,jcbot , &
pblh ,zm ,geos ,zi ,qtnd , &
heat ,pap ,paph ,dpp , &
delt ,mcon ,cme , &
tpert ,dlf ,pflx ,zdu ,rprd , &
mu ,md ,du ,eu ,ed , &
dp ,dsubcld ,jt ,maxg ,ideep , &
lengath ,ql ,rliq )
!-----------------------------------------------------------------------
!
! Purpose:
! Main driver for zhang-mcfarlane convection scheme
!
! Method:
! performs deep convective adjustment based on mass-flux closure
! algorithm.
!
! Author:guang jun zhang, m.lazare, n.mcfarlane. CAM Contact: P. Rasch
!
! This is contributed code not fully standardized by the CAM core group.
! All variables have been typed, where most are identified in comments
! The current procedure will be reimplemented in a subsequent version
! of the CAM where it will include a more straightforward formulation
! and will make use of the standard CAM nomenclature
!
!-----------------------------------------------------------------------
use constituents, only: pcnst, pnats
!
! ************************ index of variables **********************
!
! wg * alpha array of vertical differencing used (=1. for upstream).
! w * cape convective available potential energy.
! wg * capeg gathered convective available potential energy.
! c * capelmt threshold value for cape for deep convection.
! ic * cpres specific heat at constant pressure in j/kg-degk.
! i * dpp
! ic * delt length of model time-step in seconds.
! wg * dp layer thickness in mbs (between upper/lower interface).
! wg * dqdt mixing ratio tendency at gathered points.
! wg * dsdt dry static energy ("temp") tendency at gathered points.
! wg * dudt u-wind tendency at gathered points.
! wg * dvdt v-wind tendency at gathered points.
! wg * dsubcld layer thickness in mbs between lcl and maxi.
! ic * grav acceleration due to gravity in m/sec2.
! wg * du detrainment in updraft. specified in mid-layer
! wg * ed entrainment in downdraft.
! wg * eu entrainment in updraft.
! wg * hmn moist static energy.
! wg * hsat saturated moist static energy.
! w * ideep holds position of gathered points vs longitude index.
! ic * pver number of model levels.
! wg * j0 detrainment initiation level index.
! wg * jd downdraft initiation level index.
! ic * jlatpr gaussian latitude index for printing grids (if needed).
! wg * jt top level index of deep cumulus convection.
! w * lcl base level index of deep cumulus convection.
! wg * lclg gathered values of lcl.
! w * lel index of highest theoretical convective plume.
! wg * lelg gathered values of lel.
! w * lon index of onset level for deep convection.
! w * maxi index of level with largest moist static energy.
! wg * maxg gathered values of maxi.
! wg * mb cloud base mass flux.
! wg * mc net upward (scaled by mb) cloud mass flux.
! wg * md downward cloud mass flux (positive up).
! wg * mu upward cloud mass flux (positive up). specified
! at interface
! ic * msg number of missing moisture levels at the top of model.
! w * p grid slice of ambient mid-layer pressure in mbs.
! i * pblt row of pbl top indices.
! w * pcpdh scaled surface pressure.
! w * pf grid slice of ambient interface pressure in mbs.
! wg * pg grid slice of gathered values of p.
! w * q grid slice of mixing ratio.
! wg * qd grid slice of mixing ratio in downdraft.
! wg * qg grid slice of gathered values of q.
! i/o * qh grid slice of specific humidity.
! w * qh0 grid slice of initial specific humidity.
! wg * qhat grid slice of upper interface mixing ratio.
! wg * ql grid slice of cloud liquid water.
! wg * qs grid slice of saturation mixing ratio.
! w * qstp grid slice of parcel temp. saturation mixing ratio.
! wg * qstpg grid slice of gathered values of qstp.
! wg * qu grid slice of mixing ratio in updraft.
! ic * rgas dry air gas constant.
! wg * rl latent heat of vaporization.
! w * s grid slice of scaled dry static energy (t+gz/cp).
! wg * sd grid slice of dry static energy in downdraft.
! wg * sg grid slice of gathered values of s.
! wg * shat grid slice of upper interface dry static energy.
! wg * su grid slice of dry static energy in updraft.
! i/o * t
! o * jctop row of top-of-deep-convection indices passed out.
! O * jcbot row of base of cloud indices passed out.
! wg * tg grid slice of gathered values of t.
! w * tl row of parcel temperature at lcl.
! wg * tlg grid slice of gathered values of tl.
! w * tp grid slice of parcel temperatures.
! wg * tpg grid slice of gathered values of tp.
! i/o * u grid slice of u-wind (real).
! wg * ug grid slice of gathered values of u.
! i/o * utg grid slice of u-wind tendency (real).
! i/o * v grid slice of v-wind (real).
! w * va work array re-used by called subroutines.
! wg * vg grid slice of gathered values of v.
! i/o * vtg grid slice of v-wind tendency (real).
! i * w grid slice of diagnosed large-scale vertical velocity.
! w * z grid slice of ambient mid-layer height in metres.
! w * zf grid slice of ambient interface height in metres.
! wg * zfg grid slice of gathered values of zf.
! wg * zg grid slice of gathered values of z.
!
!-----------------------------------------------------------------------
!
! multi-level i/o fields:
! i => input arrays.
! i/o => input/output arrays.
! w => work arrays.
! wg => work arrays operating only on gathered points.
! ic => input data constants.
! c => data constants pertaining to subroutine itself.
!
! input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: t(pcols,pver) ! grid slice of temperature at mid-layer.
real(r8), intent(in) :: qh(pcols,pver,pcnst+pnats) ! grid slice of specific humidity.
real(r8), intent(in) :: pap(pcols,pver)
real(r8), intent(in) :: paph(pcols,pver+1)
real(r8), intent(in) :: dpp(pcols,pver) ! local sigma half-level thickness (i.e. dshj).
real(r8), intent(in) :: zm(pcols,pver)
real(r8), intent(in) :: geos(pcols)
real(r8), intent(in) :: zi(pcols,pver+1)
real(r8), intent(in) :: pblh(pcols)
real(r8), intent(in) :: tpert(pcols)
!
! output arguments
!
real(r8), intent(out) :: qtnd(pcols,pver) ! specific humidity tendency (kg/kg/s)
real(r8), intent(out) :: heat(pcols,pver) ! heating rate (dry static energy tendency, W/kg)
real(r8), intent(out) :: mcon(pcols,pverp)
real(r8), intent(out) :: dlf(pcols,pver) ! scattrd version of the detraining cld h2o tend
real(r8), intent(out) :: pflx(pcols,pverp) ! scattered precip flux at each level
real(r8), intent(out) :: cme(pcols,pver)
real(r8), intent(out) :: zdu(pcols,pver)
real(r8), intent(out) :: rprd(pcols,pver) ! rain production rate
! move these vars from local storage to output so that convective
! transports can be done in outside of conv_cam.
real(r8), intent(out) :: mu(pcols,pver)
real(r8), intent(out) :: eu(pcols,pver)
real(r8), intent(out) :: du(pcols,pver)
real(r8), intent(out) :: md(pcols,pver)
real(r8), intent(out) :: ed(pcols,pver)
real(r8), intent(out) :: dp(pcols,pver) ! wg layer thickness in mbs (between upper/lower interface).
real(r8), intent(out) :: dsubcld(pcols) ! wg layer thickness in mbs between lcl and maxi.
real(r8), intent(out) :: jctop(pcols) ! o row of top-of-deep-convection indices passed out.
real(r8), intent(out) :: jcbot(pcols) ! o row of base of cloud indices passed out.
real(r8), intent(out) :: prec(pcols)
real(r8), intent(out) :: rliq(pcols) ! reserved liquid (not yet in cldliq) for energy integrals
real(r8) zs(pcols)
real(r8) dlg(pcols,pver) ! gathrd version of the detraining cld h2o tend
real(r8) pflxg(pcols,pverp) ! gather precip flux at each level
real(r8) cug(pcols,pver) ! gathered condensation rate
real(r8) evpg(pcols,pver) ! gathered evap rate of rain in downdraft
real(r8) mumax(pcols)
integer jt(pcols) ! wg top level index of deep cumulus convection.
integer maxg(pcols) ! wg gathered values of maxi.
integer ideep(pcols) ! w holds position of gathered points vs longitude index.
integer lengath
! diagnostic field used by chem/wetdep codes
real(r8) ql(pcols,pver) ! wg grid slice of cloud liquid water.
!
real(r8) pblt(pcols) ! i row of pbl top indices.
!
!-----------------------------------------------------------------------
!
! general work fields (local variables):
!
real(r8) q(pcols,pver) ! w grid slice of mixing ratio.
real(r8) p(pcols,pver) ! w grid slice of ambient mid-layer pressure in mbs.
real(r8) z(pcols,pver) ! w grid slice of ambient mid-layer height in metres.
real(r8) s(pcols,pver) ! w grid slice of scaled dry static energy (t+gz/cp).
real(r8) tp(pcols,pver) ! w grid slice of parcel temperatures.
real(r8) zf(pcols,pver+1) ! w grid slice of ambient interface height in metres.
real(r8) pf(pcols,pver+1) ! w grid slice of ambient interface pressure in mbs.
real(r8) qstp(pcols,pver) ! w grid slice of parcel temp. saturation mixing ratio.
real(r8) cape(pcols) ! w convective available potential energy.
real(r8) tl(pcols) ! w row of parcel temperature at lcl.
integer lcl(pcols) ! w base level index of deep cumulus convection.
integer lel(pcols) ! w index of highest theoretical convective plume.
integer lon(pcols) ! w index of onset level for deep convection.
integer maxi(pcols) ! w index of level with largest moist static energy.
integer index(pcols)
real(r8) precip
!
! gathered work fields:
!
real(r8) qg(pcols,pver) ! wg grid slice of gathered values of q.
real(r8) tg(pcols,pver) ! w grid slice of temperature at interface.
real(r8) pg(pcols,pver) ! wg grid slice of gathered values of p.
real(r8) zg(pcols,pver) ! wg grid slice of gathered values of z.
real(r8) sg(pcols,pver) ! wg grid slice of gathered values of s.
real(r8) tpg(pcols,pver) ! wg grid slice of gathered values of tp.
real(r8) zfg(pcols,pver+1) ! wg grid slice of gathered values of zf.
real(r8) qstpg(pcols,pver) ! wg grid slice of gathered values of qstp.
real(r8) ug(pcols,pver) ! wg grid slice of gathered values of u.
real(r8) vg(pcols,pver) ! wg grid slice of gathered values of v.
real(r8) cmeg(pcols,pver)
real(r8) rprdg(pcols,pver) ! wg gathered rain production rate
real(r8) capeg(pcols) ! wg gathered convective available potential energy.
real(r8) tlg(pcols) ! wg grid slice of gathered values of tl.
integer lclg(pcols) ! wg gathered values of lcl.
integer lelg(pcols)
!
! work fields arising from gathered calculations.
!
real(r8) dqdt(pcols,pver) ! wg mixing ratio tendency at gathered points.
real(r8) dsdt(pcols,pver) ! wg dry static energy ("temp") tendency at gathered points.
! real(r8) alpha(pcols,pver) ! array of vertical differencing used (=1. for upstream).
real(r8) sd(pcols,pver) ! wg grid slice of dry static energy in downdraft.
real(r8) qd(pcols,pver) ! wg grid slice of mixing ratio in downdraft.
real(r8) mc(pcols,pver) ! wg net upward (scaled by mb) cloud mass flux.
real(r8) qhat(pcols,pver) ! wg grid slice of upper interface mixing ratio.
real(r8) qu(pcols,pver) ! wg grid slice of mixing ratio in updraft.
real(r8) su(pcols,pver) ! wg grid slice of dry static energy in updraft.
real(r8) qs(pcols,pver) ! wg grid slice of saturation mixing ratio.
real(r8) shat(pcols,pver) ! wg grid slice of upper interface dry static energy.
real(r8) hmn(pcols,pver) ! wg moist static energy.
real(r8) hsat(pcols,pver) ! wg saturated moist static energy.
real(r8) qlg(pcols,pver)
real(r8) dudt(pcols,pver) ! wg u-wind tendency at gathered points.
real(r8) dvdt(pcols,pver) ! wg v-wind tendency at gathered points.
! real(r8) ud(pcols,pver)
! real(r8) vd(pcols,pver)
real(r8) mb(pcols) ! wg cloud base mass flux.
integer jlcl(pcols)
integer j0(pcols) ! wg detrainment initiation level index.
integer jd(pcols) ! wg downdraft initiation level index.
real(r8) delt ! length of model time-step in seconds.
integer i
integer ii
integer k
integer msg ! ic number of missing moisture levels at the top of model.
real(r8) qdifr
real(r8) sdifr
!
!--------------------------Data statements------------------------------
!
! Set internal variable "msg" (convection limit) to "limcnv-1"
!
msg = limcnv - 1
!
! initialize necessary arrays.
! zero out variables not used in cam
!
qtnd(:,:) = 0.
heat(:,:) = 0.
mcon(:,:) = 0.
rliq(:ncol) = 0.
!
! initialize convective tendencies
!
prec(:ncol) = 0.
do k = 1,pver
do i = 1,ncol
dqdt(i,k) = 0.
dsdt(i,k) = 0.
dudt(i,k) = 0.
dvdt(i,k) = 0.
pflx(i,k) = 0.
pflxg(i,k) = 0.
cme(i,k) = 0.
rprd(i,k) = 0.
zdu(i,k) = 0.
ql(i,k) = 0.
qlg(i,k) = 0.
dlf(i,k) = 0.
dlg(i,k) = 0.
end do
end do
do i = 1,ncol
pflx(i,pverp) = 0
pflxg(i,pverp) = 0
end do
!
do i = 1,ncol
pblt(i) = pver
dsubcld(i) = 0.
jctop(i) = pver
jcbot(i) = 1
end do
!
! calculate local pressure (mbs) and height (m) for both interface
! and mid-layer locations.
!
do i = 1,ncol
zs(i) = geos(i)*rgrav
pf(i,pver+1) = paph(i,pver+1)*0.01
zf(i,pver+1) = zi(i,pver+1) + zs(i)
end do
do k = 1,pver
do i = 1,ncol
p(i,k) = pap(i,k)*0.01
pf(i,k) = paph(i,k)*0.01
z(i,k) = zm(i,k) + zs(i)
zf(i,k) = zi(i,k) + zs(i)
end do
end do
!
do k = pver - 1,msg + 1,-1
do i = 1,ncol
if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5) pblt(i) = k
end do
end do
!
! store incoming specific humidity field for subsequent calculation
! of precipitation (through change in storage).
! define dry static energy (normalized by cp).
!
do k = 1,pver
do i = 1,ncol
q(i,k) = qh(i,k,1)
s(i,k) = t(i,k) + (grav/cpres)*z(i,k)
tp(i,k)=0.0
shat(i,k) = s(i,k)
qhat(i,k) = q(i,k)
end do
end do
do i = 1,ncol
capeg(i) = 0.
lclg(i) = 1
lelg(i) = pver
maxg(i) = 1
tlg(i) = 400.
dsubcld(i) = 0.
end do
!
! evaluate covective available potential energy (cape).
!
call buoyan(lchnk ,ncol , &
q ,t ,p ,z ,pf , &
tp ,qstp ,tl ,rl ,cape , &
pblt ,lcl ,lel ,lon ,maxi , &
rgas ,grav ,cpres ,msg , &
tpert )
!
! determine whether grid points will undergo some deep convection
! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel
! (require cape.gt. 0 and lel<lcl as minimum conditions).
!
lengath = 0
do i=1,ncol
if (cape(i) > capelmt) then
lengath = lengath + 1
index(lengath) = i
end if
end do
if (lengath.eq.0) return
do ii=1,lengath
i=index(ii)
ideep(ii)=i
end do
!
! obtain gathered arrays necessary for ensuing calculations.
!
do k = 1,pver
do i = 1,lengath
dp(i,k) = 0.01*dpp(ideep(i),k)
qg(i,k) = q(ideep(i),k)
tg(i,k) = t(ideep(i),k)
pg(i,k) = p(ideep(i),k)
zg(i,k) = z(ideep(i),k)
sg(i,k) = s(ideep(i),k)
tpg(i,k) = tp(ideep(i),k)
zfg(i,k) = zf(ideep(i),k)
qstpg(i,k) = qstp(ideep(i),k)
ug(i,k) = 0.
vg(i,k) = 0.
end do
end do
!
do i = 1,lengath
zfg(i,pver+1) = zf(ideep(i),pver+1)
end do
do i = 1,lengath
capeg(i) = cape(ideep(i))
lclg(i) = lcl(ideep(i))
lelg(i) = lel(ideep(i))
maxg(i) = maxi(ideep(i))
tlg(i) = tl(ideep(i))
end do
!
! calculate sub-cloud layer pressure "thickness" for use in
! closure and tendency routines.
!
do k = msg + 1,pver
do i = 1,lengath
if (k >= maxg(i)) then
dsubcld(i) = dsubcld(i) + dp(i,k)
end if
end do
end do
!
! define array of factors (alpha) which defines interfacial
! values, as well as interfacial values for (q,s) used in
! subsequent routines.
!
do k = msg + 2,pver
do i = 1,lengath
! alpha(i,k) = 0.5
sdifr = 0.
qdifr = 0.
if (sg(i,k) > 0. .or. sg(i,k-1) > 0.) &
sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
if (qg(i,k) > 0. .or. qg(i,k-1) > 0.) &
qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
if (sdifr > 1.E-6) then
shat(i,k) = log(sg(i,k-1)/sg(i,k))*sg(i,k-1)*sg(i,k)/(sg(i,k-1)-sg(i,k))
else
shat(i,k) = 0.5* (sg(i,k)+sg(i,k-1))
end if
if (qdifr > 1.E-6) then
qhat(i,k) = log(qg(i,k-1)/qg(i,k))*qg(i,k-1)*qg(i,k)/(qg(i,k-1)-qg(i,k))
else
qhat(i,k) = 0.5* (qg(i,k)+qg(i,k-1))
end if
end do
end do
!
! obtain cloud properties.
!
call cldprp(lchnk , &
qg ,tg ,ug ,vg ,pg , &
zg ,sg ,mu ,eu ,du , &
md ,ed ,sd ,qd ,mc , &
qu ,su ,zfg ,qs ,hmn , &
hsat ,shat ,qlg , &
cmeg ,maxg ,lelg ,jt ,jlcl , &
maxg ,j0 ,jd ,rl ,lengath , &
rgas ,grav ,cpres ,msg , &
pflxg ,evpg ,cug ,rprdg ,limcnv )
!
! convert detrainment from units of "1/m" to "1/mb".
!
do k = msg + 1,pver
do i = 1,lengath
du (i,k) = du (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
eu (i,k) = eu (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
ed (i,k) = ed (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
cug (i,k) = cug (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
end do
end do
call closure(lchnk , &
qg ,tg ,pg ,zg ,sg , &
tpg ,qs ,qu ,su ,mc , &
du ,mu ,md ,qd ,sd , &
qhat ,shat ,dp ,qstpg ,zfg , &
qlg ,dsubcld ,mb ,capeg ,tlg , &
lclg ,lelg ,jt ,maxg ,1 , &
lengath ,rgas ,grav ,cpres ,rl , &
msg ,capelmt )
!
! limit cloud base mass flux to theoretical upper bound.
!
do i=1,lengath
mumax(i) = 0
end do
do k=msg + 2,pver
do i=1,lengath
mumax(i) = max(mumax(i), mu(i,k)/dp(i,k))
end do
end do
do i=1,lengath
if (mumax(i) > 0.) then
mb(i) = min(mb(i),0.5/(delt*mumax(i)))
else
mb(i) = 0.
endif
end do
do k=msg+1,pver
do i=1,lengath
mu (i,k) = mu (i,k)*mb(i)
md (i,k) = md (i,k)*mb(i)
mc (i,k) = mc (i,k)*mb(i)
du (i,k) = du (i,k)*mb(i)
eu (i,k) = eu (i,k)*mb(i)
ed (i,k) = ed (i,k)*mb(i)
cmeg (i,k) = cmeg (i,k)*mb(i)
rprdg(i,k) = rprdg(i,k)*mb(i)
cug (i,k) = cug (i,k)*mb(i)
evpg (i,k) = evpg (i,k)*mb(i)
pflxg(i,k+1)= pflxg(i,k+1)*mb(i)*100./grav
end do
end do
!
! compute temperature and moisture changes due to convection.
!
call q1q2_pjr(lchnk , &
dqdt ,dsdt ,qg ,qs ,qu , &
su ,du ,qhat ,shat ,dp , &
mu ,md ,sd ,qd ,qlg , &
dsubcld ,jt ,maxg ,1 ,lengath , &
cpres ,rl ,msg , &
dlg ,evpg ,cug )
!
! gather back temperature and mixing ratio.
!
do k = msg + 1,pver
do i = 1,lengath
!
! q is updated to compute net precip.
!
q(ideep(i),k) = qh(ideep(i),k,1) + 2.*delt*dqdt(i,k)
qtnd(ideep(i),k) = dqdt (i,k)
cme (ideep(i),k) = cmeg (i,k)
rprd(ideep(i),k) = rprdg(i,k)
zdu (ideep(i),k) = du (i,k)
mcon(ideep(i),k) = mc (i,k)
heat(ideep(i),k) = dsdt (i,k)*cpres
dlf (ideep(i),k) = dlg (i,k)
pflx(ideep(i),k) = pflxg(i,k)
ql (ideep(i),k) = qlg (i,k)
end do
end do
!
do i = 1,lengath
jctop(ideep(i)) = jt(i)
!++bee
jcbot(ideep(i)) = maxg(i)
!--bee
pflx(ideep(i),pverp) = pflxg(i,pverp)
end do
! Compute precip by integrating change in water vapor minus detrained cloud water
do k = pver,msg + 1,-1
do i = 1,ncol
prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k,1)) - dpp(i,k)*dlf(i,k)*2*delt
end do
end do
! obtain final precipitation rate in m/s.
do i = 1,ncol
prec(i) = rgrav*max(prec(i),0._r8)/ (2.*delt)/1000.
end do
! Compute reserved liquid (not yet in cldliq) for energy integrals.
! Treat rliq as flux out bottom, to be added back later.
do k = 1, pver
do i = 1, ncol
rliq(i) = rliq(i) + dlf(i,k)*dpp(i,k)/gravit
end do
end do
rliq(:ncol) = rliq(:ncol) /1000.
return
end subroutine zm_convr
!===============================================================================
subroutine zm_conv_evap(state, ptend, prdprec, cldfrc, deltat, prec, snow, hackflag) 2,15
!-----------------------------------------------------------------------
! Compute tendencies due to evaporation of rain from ZM scheme
!--
! Compute the total precipitation and snow fluxes at the surface.
! Add in the latent heat of fusion for snow formation and melt, since it not dealt with
! in the Zhang-MacFarlane parameterization.
! Evaporate some of the precip directly into the environment using a Sundqvist type algorithm
!-----------------------------------------------------------------------
use physics_types, only: physics_state, physics_ptend
use wv_saturation, only: aqsat
use phys_grid, only: get_rlat_all_p
!------------------------------Arguments--------------------------------
type(physics_state), intent(in) :: state ! Physics state variables
type(physics_ptend), intent(inout) :: ptend ! indivdual parameterization tendencies
real(r8), intent(in ) :: prdprec(pcols,pver)! precipitation production (kg/ks/s)
real(r8), intent(in ) :: cldfrc(pcols,pver) ! cloud fraction
real(r8), intent(in ) :: deltat ! time step
logical, intent(in ) :: hackflag ! input is from Hack instead of ZM scheme
! only used to specify output variables
real(r8), intent(inout) :: prec(pcols) ! Convective-scale preciptn rate
real(r8), intent(out) :: snow(pcols) ! Convective-scale snowfall rate
!
!---------------------------Local storage-------------------------------
real(r8) :: est (pcols,pver) ! Saturation vapor pressure
real(r8) :: fice (pcols,pver) ! ice fraction in precip production
real(r8) :: fsnow_conv(pcols,pver) ! snow fraction in precip production
real(r8) :: qsat (pcols,pver) ! saturation specific humidity
real(r8) :: flxprec(pcols,pverp) ! Convective-scale flux of precip at interfaces (kg/m2/s)
real(r8) :: flxsnow(pcols,pverp) ! Convective-scale flux of snow at interfaces (kg/m2/s)
real(r8) :: ntprprd(pcols,pver) ! net precip production in layer
real(r8) :: ntsnprd(pcols,pver) ! net snow production in layer
real(r8) :: work1 ! temp variable (pjr)
real(r8) :: work2 ! temp variable (pjr)
real(r8) :: evpvint(pcols) ! vertical integral of evaporation
real(r8) :: evpprec(pcols) ! evaporation of precipitation (kg/kg/s)
real(r8) :: evpsnow(pcols) ! evaporation of snowfall (kg/kg/s)
real(r8) :: snowmlt(pcols) ! snow melt tendency in layer
real(r8) :: flxsntm(pcols) ! flux of snow into layer, after melting
real(r8) :: evplimit ! temp variable for evaporation limits
real(r8) :: rlat(pcols)
integer :: i,k ! longitude,level indices
integer :: ncol, lchnk ! number of columns and chunk index
!-----------------------------------------------------------------------
!!$ ke = conke
!!$ ke = 2.0E-6
ncol = state%ncol
lchnk = state%lchnk
! heating and specific humidity tendencies produced
ptend%name = 'zm_conv_evap'
ptend%ls = .TRUE.
ptend%lq(1) = .TRUE.
! convert input precip to kg/m2/s
prec(:ncol) = prec(:ncol)*1000.
! determine saturation vapor pressure
call aqsat (state%t ,state%pmid ,est ,qsat ,pcols , &
ncol ,pver ,1 ,pver )
! determine ice fraction in rain production (use cloud water parameterization fraction at present)
call cldwat_fice(ncol, state%t, fice, fsnow_conv)
! zero the flux integrals on the top boundary
flxprec(:ncol,1) = 0.
flxsnow(:ncol,1) = 0.
evpvint(:ncol) = 0.
do k = 1, pver
do i = 1, ncol
! Melt snow falling into layer, if necessary.
if (state%t(i,k) > tmelt) then
flxsntm(i) = 0.
snowmlt(i) = flxsnow(i,k) * gravit/ state%pdel(i,k)
else
flxsntm(i) = flxsnow(i,k)
snowmlt(i) = 0.
end if
! relative humidity depression must be > 0 for evaporation
evplimit = max(1. - state%q(i,k,1)/qsat(i,k), 0._r8)
! total evaporation depends on flux in the top of the layer
! flux prec is the net production above layer minus evaporation into environmet
evpprec(i) = ke * (1. - cldfrc(i,k)) * evplimit * sqrt(flxprec(i,k))
!**********************************************************
!!$ evpprec(i) = 0. ! turn off evaporation for now
!**********************************************************
! Don't let evaporation supersaturate layer (approx). Layer may already be saturated.
! Currently does not include heating/cooling change to qsat
evplimit = max(0._r8, (qsat(i,k)-state%q(i,k,1)) / deltat)
! Don't evaporate more than is falling into the layer - do not evaporate rain formed
! in this layer but if precip production is negative, remove from the available precip
! Negative precip production occurs because of evaporation in downdrafts.
!!$ evplimit = flxprec(i,k) * gravit / state%pdel(i,k) + min(prdprec(i,k), 0.)
evplimit = min(evplimit, flxprec(i,k) * gravit / state%pdel(i,k))
! Total evaporation cannot exceed input precipitation
evplimit = min(evplimit, (prec(i) - evpvint(i)) * gravit / state%pdel(i,k))
evpprec(i) = min(evplimit, evpprec(i))
! evaporation of snow depends on snow fraction of total precipitation in the top after melting
if (flxprec(i,k) > 0.) then
! evpsnow(i) = evpprec(i) * flxsntm(i) / flxprec(i,k)
! prevent roundoff problems
work1 = min(max(0._r8,flxsntm(i)/flxprec(i,k)),1._r8)
evpsnow(i) = evpprec(i) * work1
else
evpsnow(i) = 0.
end if
! vertically integrated evaporation
evpvint(i) = evpvint(i) + evpprec(i) * state%pdel(i,k)/gravit
! net precip production is production - evaporation
ntprprd(i,k) = prdprec(i,k) - evpprec(i)
! net snow production is precip production * ice fraction - evaporation - melting
!pjrworks ntsnprd(i,k) = prdprec(i,k)*fice(i,k) - evpsnow(i) - snowmlt(i)
!pjrwrks2 ntsnprd(i,k) = prdprec(i,k)*fsnow_conv(i,k) - evpsnow(i) - snowmlt(i)
! the small amount added to flxprec in the work1 expression has been increased from
! 1e-36 to 8.64e-11 (1e-5 mm/day). This causes the temperature based partitioning
! scheme to be used for small flxprec amounts. This is to address error growth problems.
#ifdef PERGRO
work1 = min(max(0._r8,flxsnow(i,k)/(flxprec(i,k)+8.64e-11_r8)),1._r8)
#else
if (flxprec(i,k).gt.0.) then
work1 = min(max(0._r8,flxsnow(i,k)/flxprec(i,k)),1._r8)
else
work1 = 0.
endif
#endif
work2 = max(fsnow_conv(i,k), work1)
if (snowmlt(i).gt.0.) work2 = 0.
! work2 = fsnow_conv(i,k)
ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i)
! precipitation fluxes
flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * state%pdel(i,k)/gravit
flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * state%pdel(i,k)/gravit
! protect against rounding error
flxprec(i,k+1) = max(flxprec(i,k+1), 0._r8)
flxsnow(i,k+1) = max(flxsnow(i,k+1), 0._r8)
! more protection (pjr)
! flxsnow(i,k+1) = min(flxsnow(i,k+1), flxprec(i,k+1))
! heating (cooling) and moistening due to evaporation
! - latent heat of vaporization for precip production has already been accounted for
! - snow is contained in prec
ptend%s(i,k) =-evpprec(i)*latvap + ntsnprd(i,k)*latice
ptend%q(i,k,1) = evpprec(i)
end do
end do
! set output precipitation rates (m/s)
prec(:ncol) = flxprec(:ncol,pver+1) / 1000.
snow(:ncol) = flxsnow(:ncol,pver+1) / 1000.
! record history variables
if (hackflag) then
call outfld('HKFLXPRC', flxprec, pcols, lchnk)
call outfld('HKFLXSNW', flxsnow, pcols, lchnk)
call outfld('HKNTPRPD', ntprprd, pcols, lchnk)
call outfld('HKNTSNPD', ntsnprd, pcols, lchnk)
call outfld('HKEIHEAT', ptend%s, pcols, lchnk)
else
call outfld('ZMFLXPRC', flxprec, pcols, lchnk)
call outfld('ZMFLXSNW', flxsnow, pcols, lchnk)
call outfld('ZMNTPRPD', ntprprd, pcols, lchnk)
call outfld('ZMNTSNPD', ntsnprd, pcols, lchnk)
call outfld('ZMEIHEAT', ptend%s, pcols, lchnk)
end if
!**********************************************************
!!$ ptend%s(:ncol,:) = 0. ! turn heating off
!**********************************************************
end subroutine zm_conv_evap
subroutine buoyan(lchnk ,ncol , & 1
q ,t ,p ,z ,pf , &
tp ,qstp ,tl ,rl ,cape , &
pblt ,lcl ,lel ,lon ,mx , &
rd ,grav ,cp ,msg , &
tpert )
!-----------------------------------------------------------------------
!
! Purpose:
! <Say what the routine does>
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author:
! This is contributed code not fully standardized by the CCM core group.
! The documentation has been enhanced to the degree that we are able.
! Reviewed: P. Rasch, April 1996
!
!-----------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
!
! input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: q(pcols,pver) ! spec. humidity
real(r8), intent(in) :: t(pcols,pver) ! temperature
real(r8), intent(in) :: p(pcols,pver) ! pressure
real(r8), intent(in) :: z(pcols,pver) ! height
real(r8), intent(in) :: pf(pcols,pver+1) ! pressure at interfaces
real(r8), intent(in) :: pblt(pcols) ! index of pbl depth
real(r8), intent(in) :: tpert(pcols) ! perturbation temperature by pbl processes
!
! output arguments
!
real(r8), intent(out) :: tp(pcols,pver) ! parcel temperature
real(r8), intent(out) :: qstp(pcols,pver) ! saturation mixing ratio of parcel
real(r8), intent(out) :: tl(pcols) ! parcel temperature at lcl
real(r8), intent(out) :: cape(pcols) ! convective aval. pot. energy.
integer lcl(pcols) !
integer lel(pcols) !
integer lon(pcols) ! level of onset of deep convection
integer mx(pcols) ! level of max moist static energy
!
!--------------------------Local Variables------------------------------
!
real(r8) capeten(pcols,5) ! provisional value of cape
real(r8) tv(pcols,pver) !
real(r8) tpv(pcols,pver) !
real(r8) buoy(pcols,pver)
real(r8) a1(pcols)
real(r8) a2(pcols)
real(r8) estp(pcols)
real(r8) pl(pcols)
real(r8) plexp(pcols)
real(r8) hmax(pcols)
real(r8) hmn(pcols)
real(r8) y(pcols)
logical plge600(pcols)
integer knt(pcols)
integer lelten(pcols,5)
real(r8) cp
real(r8) e
real(r8) grav
integer i
integer k
integer msg
integer n
real(r8) rd
real(r8) rl
#ifdef PERGRO
real(r8) rhd
#endif
!
!-----------------------------------------------------------------------
!
do n = 1,5
do i = 1,ncol
lelten(i,n) = pver
capeten(i,n) = 0.
end do
end do
!
do i = 1,ncol
lon(i) = pver
knt(i) = 0
lel(i) = pver
mx(i) = lon(i)
cape(i) = 0.
hmax(i) = 0.
end do
tp(:ncol,:) = t(:ncol,:)
qstp(:ncol,:) = q(:ncol,:)
!
! set "launching" level(mx) to be at maximum moist static energy.
! search for this level stops at planetary boundary layer top.
!
#ifdef PERGRO
do k = pver,msg + 1,-1
do i = 1,ncol
hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
!
! Reset max moist static energy level when relative difference exceeds 1.e-4
!
rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i))
if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4) then
hmax(i) = hmn(i)
mx(i) = k
end if
end do
end do
#else
do k = pver,msg + 1,-1
do i = 1,ncol
hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
hmax(i) = hmn(i)
mx(i) = k
end if
end do
end do
#endif
!
do i = 1,ncol
lcl(i) = mx(i)
e = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
tl(i) = 2840./ (3.5*log(t(i,mx(i)))-log(e)-4.805) + 55.
if (tl(i) < t(i,mx(i))) then
plexp(i) = (1./ (0.2854* (1.-0.28*q(i,mx(i)))))
pl(i) = p(i,mx(i))* (tl(i)/t(i,mx(i)))**plexp(i)
else
tl(i) = t(i,mx(i))
pl(i) = p(i,mx(i))
end if
end do
!
! calculate lifting condensation level (lcl).
!
do k = pver,msg + 2,-1
do i = 1,ncol
if (k <= mx(i) .and. (p(i,k) > pl(i) .and. p(i,k-1) <= pl(i))) then
lcl(i) = k - 1
end if
end do
end do
!
! if lcl is above the nominal level of non-divergence (600 mbs),
! no deep convection is permitted (ensuing calculations
! skipped and cape retains initialized value of zero).
!
do i = 1,ncol
plge600(i) = pl(i).ge.600.
end do
!
! initialize parcel properties in sub-cloud layer below lcl.
!
do k = pver,msg + 1,-1
do i=1,ncol
if (k > lcl(i) .and. k <= mx(i) .and. plge600(i)) then
tv(i,k) = t(i,k)* (1.+1.608*q(i,k))/ (1.+q(i,k))
qstp(i,k) = q(i,mx(i))
tp(i,k) = t(i,mx(i))* (p(i,k)/p(i,mx(i)))**(0.2854* (1.-0.28*q(i,mx(i))))
!
! buoyancy is increased by 0.5 k as in tiedtke
!
!-jjh tpv (i,k)=tp(i,k)*(1.+1.608*q(i,mx(i)))/
!-jjh 1 (1.+q(i,mx(i)))
tpv(i,k) = (tp(i,k)+tpert(i))*(1.+1.608*q(i,mx(i)))/ (1.+q(i,mx(i)))
buoy(i,k) = tpv(i,k) - tv(i,k) + 0.5
end if
end do
end do
!
! define parcel properties at lcl (i.e. level immediately above pl).
!
do k = pver,msg + 1,-1
do i=1,ncol
if (k == lcl(i) .and. plge600(i)) then
tv(i,k) = t(i,k)* (1.+1.608*q(i,k))/ (1.+q(i,k))
qstp(i,k) = q(i,mx(i))
tp(i,k) = tl(i)* (p(i,k)/pl(i))**(0.2854* (1.-0.28*qstp(i,k)))
! estp(i) =exp(a-b/tp(i,k))
! use of different formulas for est has about 1 g/kg difference
! in qs at t= 300k, and 0.02 g/kg at t=263k, with the formula
! above giving larger qs.
!
estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3))
qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
a1(i) = cp / rl + qstp(i,k) * (1.+ qstp(i,k) / eps1) * rl * eps1 / &
(rd * tp(i,k) ** 2)
a2(i) = .5* (qstp(i,k)* (1.+2./eps1*qstp(i,k))* &
(1.+qstp(i,k)/eps1)*eps1**2*rl*rl/ &
(rd**2*tp(i,k)**4)-qstp(i,k)* &
(1.+qstp(i,k)/eps1)*2.*eps1*rl/ &
(rd*tp(i,k)**3))
a1(i) = 1./a1(i)
a2(i) = -a2(i)*a1(i)**3
y(i) = q(i,mx(i)) - qstp(i,k)
tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2
! estp(i) =exp(a-b/tp(i,k))
estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3))
qstp(i,k) = eps1*estp(i) / (p(i,k)-estp(i))
!
! buoyancy is increased by 0.5 k in cape calculation.
! dec. 9, 1994
!-jjh tpv(i,k) =tp(i,k)*(1.+1.608*qstp(i,k))/(1.+q(i,mx(i)))
!
tpv(i,k) = (tp(i,k)+tpert(i))* (1.+1.608*qstp(i,k)) / (1.+q(i,mx(i)))
buoy(i,k) = tpv(i,k) - tv(i,k) + 0.5
end if
end do
end do
!
! main buoyancy calculation.
!
do k = pver - 1,msg + 1,-1
do i=1,ncol
if (k < lcl(i) .and. plge600(i)) then
tv(i,k) = t(i,k)* (1.+1.608*q(i,k))/ (1.+q(i,k))
qstp(i,k) = qstp(i,k+1)
tp(i,k) = tp(i,k+1)* (p(i,k)/p(i,k+1))**(0.2854* (1.-0.28*qstp(i,k)))
! estp(i) = exp(a-b/tp(i,k))
estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3))
qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
a1(i) = cp/rl + qstp(i,k)* (1.+qstp(i,k)/eps1)*rl*eps1/ (rd*tp(i,k)**2)
a2(i) = .5* (qstp(i,k)* (1.+2./eps1*qstp(i,k))* &
(1.+qstp(i,k)/eps1)*eps1**2*rl*rl/ &
(rd**2*tp(i,k)**4)-qstp(i,k)* &
(1.+qstp(i,k)/eps1)*2.*eps1*rl/ &
(rd*tp(i,k)**3))
a1(i) = 1./a1(i)
a2(i) = -a2(i)*a1(i)**3
y(i) = qstp(i,k+1) - qstp(i,k)
tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2
! estp(i) =exp(a-b/tp(i,k))
estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3))
qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
!-jjh tpv(i,k) =tp(i,k)*(1.+1.608*qstp(i,k))/
!jt (1.+q(i,mx(i)))
tpv(i,k) = (tp(i,k)+tpert(i))* (1.+1.608*qstp(i,k))/(1.+q(i,mx(i)))
buoy(i,k) = tpv(i,k) - tv(i,k) + 0.5
end if
end do
end do
!
do k = msg + 2,pver
do i = 1,ncol
if (k < lcl(i) .and. plge600(i)) then
if (buoy(i,k+1) > 0. .and. buoy(i,k) <= 0.) then
knt(i) = min(5,knt(i) + 1)
lelten(i,knt(i)) = k
end if
end if
end do
end do
!
! calculate convective available potential energy (cape).
!
do n = 1,5
do k = msg + 1,pver
do i = 1,ncol
if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k))
end if
end do
end do
end do
!
! find maximum cape from all possible tentative capes from
! one sounding,
! and use it as the final cape, april 26, 1995
!
do n = 1,5
do i = 1,ncol
if (capeten(i,n) > cape(i)) then
cape(i) = capeten(i,n)
lel(i) = lelten(i,n)
end if
end do
end do
!
! put lower bound on cape for diagnostic purposes.
!
do i = 1,ncol
cape(i) = max(cape(i), 0._r8)
end do
!
return
end subroutine buoyan
subroutine cldprp(lchnk , & 1
q ,t ,u ,v ,p , &
z ,s ,mu ,eu ,du , &
md ,ed ,sd ,qd ,mc , &
qu ,su ,zf ,qst ,hmn , &
hsat ,shat ,ql , &
cmeg ,jb ,lel ,jt ,jlcl , &
mx ,j0 ,jd ,rl ,il2g , &
rd ,grav ,cp ,msg , &
pflx ,evp ,cu ,rprd ,limcnv )
!-----------------------------------------------------------------------
!
! Purpose:
! <Say what the routine does>
!
! Method:
! may 09/91 - guang jun zhang, m.lazare, n.mcfarlane.
! original version cldprop.
!
! Author: See above, modified by P. Rasch
! This is contributed code not fully standardized by the CCM core group.
!
! this code is very much rougher than virtually anything else in the CCM
! there are debug statements left strewn about and code segments disabled
! these are to facilitate future development. We expect to release a
! cleaner code in a future release
!
! the documentation has been enhanced to the degree that we are able
!
!-----------------------------------------------------------------------
implicit none
!------------------------------------------------------------------------------
!
! Input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
real(r8), intent(in) :: q(pcols,pver) ! spec. humidity of env
real(r8), intent(in) :: t(pcols,pver) ! temp of env
real(r8), intent(in) :: p(pcols,pver) ! pressure of env
real(r8), intent(in) :: z(pcols,pver) ! height of env
real(r8), intent(in) :: s(pcols,pver) ! normalized dry static energy of env
real(r8), intent(in) :: zf(pcols,pverp) ! height of interfaces
real(r8), intent(in) :: u(pcols,pver) ! zonal velocity of env
real(r8), intent(in) :: v(pcols,pver) ! merid. velocity of env
integer, intent(in) :: jb(pcols) ! updraft base level
integer, intent(in) :: lel(pcols) ! updraft launch level
integer, intent(out) :: jt(pcols) ! updraft plume top
integer, intent(out) :: jlcl(pcols) ! updraft lifting cond level
integer, intent(in) :: mx(pcols) ! updraft base level (same is jb)
integer, intent(out) :: j0(pcols) ! level where updraft begins detraining
integer, intent(out) :: jd(pcols) ! level of downdraft
integer, intent(in) :: limcnv ! convection limiting level
integer, intent(in) :: il2g !CORE GROUP REMOVE
integer, intent(in) :: msg ! missing moisture vals (always 0)
real(r8), intent(in) :: rl ! latent heat of vap
real(r8), intent(in) :: shat(pcols,pver) ! interface values of dry stat energy
!
! output
!
real(r8), intent(out) :: rprd(pcols,pver) ! rate of production of precip at that layer
real(r8), intent(out) :: du(pcols,pver) ! detrainement rate of updraft
real(r8), intent(out) :: ed(pcols,pver) ! entrainment rate of downdraft
real(r8), intent(out) :: eu(pcols,pver) ! entrainment rate of updraft
real(r8), intent(out) :: hmn(pcols,pver) ! moist stat energy of env
real(r8), intent(out) :: hsat(pcols,pver) ! sat moist stat energy of env
real(r8), intent(out) :: mc(pcols,pver) ! net mass flux
real(r8), intent(out) :: md(pcols,pver) ! downdraft mass flux
real(r8), intent(out) :: mu(pcols,pver) ! updraft mass flux
real(r8), intent(out) :: pflx(pcols,pverp) ! precipitation flux thru layer
real(r8), intent(out) :: qd(pcols,pver) ! spec humidity of downdraft
real(r8), intent(out) :: ql(pcols,pver) ! liq water of updraft
real(r8), intent(out) :: qst(pcols,pver) ! saturation spec humidity of env.
real(r8), intent(out) :: qu(pcols,pver) ! spec hum of updraft
real(r8), intent(out) :: sd(pcols,pver) ! normalized dry stat energy of downdraft
real(r8), intent(out) :: su(pcols,pver) ! normalized dry stat energy of updraft
real(r8) rd ! gas constant for dry air
real(r8) grav ! gravity
real(r8) cp ! heat capacity of dry air
!
! Local workspace
!
real(r8) gamma(pcols,pver)
real(r8) dz(pcols,pver)
real(r8) iprm(pcols,pver)
real(r8) hu(pcols,pver)
real(r8) hd(pcols,pver)
real(r8) eps(pcols,pver)
real(r8) f(pcols,pver)
real(r8) k1(pcols,pver)
real(r8) i2(pcols,pver)
real(r8) ihat(pcols,pver)
real(r8) i3(pcols,pver)
real(r8) idag(pcols,pver)
real(r8) i4(pcols,pver)
real(r8) qsthat(pcols,pver)
real(r8) hsthat(pcols,pver)
real(r8) gamhat(pcols,pver)
real(r8) cu(pcols,pver)
real(r8) evp(pcols,pver)
real(r8) cmeg(pcols,pver)
real(r8) qds(pcols,pver)
real(r8) hmin(pcols)
real(r8) expdif(pcols)
real(r8) expnum(pcols)
real(r8) ftemp(pcols)
real(r8) eps0(pcols)
real(r8) rmue(pcols)
real(r8) zuef(pcols)
real(r8) zdef(pcols)
real(r8) epsm(pcols)
real(r8) ratmjb(pcols)
real(r8) est(pcols)
real(r8) totpcp(pcols)
real(r8) totevp(pcols)
real(r8) alfa(pcols)
real(r8) ql1
real(r8) tu
real(r8) estu
real(r8) qstu
real(r8) small
real(r8) mdt
integer khighest
integer klowest
integer kount
integer i,k
logical doit(pcols)
logical done(pcols)
!
!------------------------------------------------------------------------------
!
do i = 1,il2g
ftemp(i) = 0.
expnum(i) = 0.
expdif(i) = 0.
end do
!
!jr Change from msg+1 to 1 to prevent blowup
!
do k = 1,pver
do i = 1,il2g
dz(i,k) = zf(i,k) - zf(i,k+1)
end do
end do
!
! initialize many output and work variables to zero
!
pflx(:il2g,1) = 0
do k = 1,pver
do i = 1,il2g
k1(i,k) = 0.
i2(i,k) = 0.
i3(i,k) = 0.
i4(i,k) = 0.
mu(i,k) = 0.
f(i,k) = 0.
eps(i,k) = 0.
eu(i,k) = 0.
du(i,k) = 0.
ql(i,k) = 0.
cu(i,k) = 0.
evp(i,k) = 0.
cmeg(i,k) = 0.
qds(i,k) = q(i,k)
md(i,k) = 0.
ed(i,k) = 0.
sd(i,k) = s(i,k)
qd(i,k) = q(i,k)
mc(i,k) = 0.
qu(i,k) = q(i,k)
su(i,k) = s(i,k)
! est(i)=exp(a-b/t(i,k))
est(i) = c1*exp((c2* (t(i,k)-tfreez))/((t(i,k)-tfreez)+c3))
!++bee
if ( p(i,k)-est(i) > 0. ) then
qst(i,k) = eps1*est(i)/ (p(i,k)-est(i))
else
qst(i,k) = 1.0
end if
!--bee
gamma(i,k) = qst(i,k)*(1. + qst(i,k)/eps1)*eps1*rl/(rd*t(i,k)**2)*rl/cp
hmn(i,k) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
hsat(i,k) = cp*t(i,k) + grav*z(i,k) + rl*qst(i,k)
hu(i,k) = hmn(i,k)
hd(i,k) = hmn(i,k)
rprd(i,k) = 0.
end do
end do
!
!jr Set to zero things which make this routine blow up
!
do k=1,msg
do i=1,il2g
rprd(i,k) = 0.
end do
end do
!
! interpolate the layer values of qst, hsat and gamma to
! layer interfaces
!
do i = 1,il2g
hsthat(i,msg+1) = hsat(i,msg+1)
qsthat(i,msg+1) = qst(i,msg+1)
gamhat(i,msg+1) = gamma(i,msg+1)
totpcp(i) = 0.
totevp(i) = 0.
end do
do k = msg + 2,pver
do i = 1,il2g
if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6) then
qsthat(i,k) = log(qst(i,k-1)/qst(i,k))*qst(i,k-1)*qst(i,k)/ (qst(i,k-1)-qst(i,k))
else
qsthat(i,k) = qst(i,k)
end if
hsthat(i,k) = cp*shat(i,k) + rl*qsthat(i,k)
if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6) then
gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))*gamma(i,k-1)*gamma(i,k)/ &
(gamma(i,k-1)-gamma(i,k))
else
gamhat(i,k) = gamma(i,k)
end if
end do
end do
!
! initialize cloud top to highest plume top.
!jr changed hard-wired 4 to limcnv+1 (not to exceed pver)
!
do i = 1,il2g
jt(i) = max(lel(i),limcnv+1)
jt(i) = min(jt(i),pver)
jd(i) = pver
jlcl(i) = lel(i)
hmin(i) = 1.E6
end do
!
! find the level of minimum hsat, where detrainment starts
!
do k = msg + 1,pver
do i = 1,il2g
if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then
hmin(i) = hsat(i,k)
j0(i) = k
end if
end do
end do
do i = 1,il2g
j0(i) = min(j0(i),jb(i)-2)
j0(i) = max(j0(i),jt(i)+2)
!
! Fix from Guang Zhang to address out of bounds array reference
!
j0(i) = min(j0(i),pver)
end do
!
! Initialize certain arrays inside cloud
!
do k = msg + 1,pver
do i = 1,il2g
if (k >= jt(i) .and. k <= jb(i)) then
hu(i,k) = hmn(i,mx(i)) + cp*0.5
su(i,k) = s(i,mx(i)) + 0.5
end if
end do
end do
!
! *********************************************************
! compute taylor series for approximate eps(z) below
! *********************************************************
!
do k = pver - 1,msg + 1,-1
do i = 1,il2g
if (k < jb(i) .and. k >= jt(i)) then
k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k)
ihat(i,k) = 0.5* (k1(i,k+1)+k1(i,k))
i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k)
idag(i,k) = 0.5* (i2(i,k+1)+i2(i,k))
i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k)
iprm(i,k) = 0.5* (i3(i,k+1)+i3(i,k))
i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k)
end if
end do
end do
!
! re-initialize hmin array for ensuing calculation.
!
do i = 1,il2g
hmin(i) = 1.E6
end do
do k = msg + 1,pver
do i = 1,il2g
if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then
hmin(i) = hmn(i,k)
expdif(i) = hmn(i,mx(i)) - hmin(i)
end if
end do
end do
!
! *********************************************************
! compute approximate eps(z) using above taylor series
! *********************************************************
!
do k = msg + 2,pver
do i = 1,il2g
expnum(i) = 0.
ftemp(i) = 0.
if (k < jt(i) .or. k >= jb(i)) then
k1(i,k) = 0.
expnum(i) = 0.
else
expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + &
hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k))
end if
if ((expdif(i) > 100. .and. expnum(i) > 0.) .and. &
k1(i,k) > expnum(i)*dz(i,k)) then
ftemp(i) = expnum(i)/k1(i,k)
f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)*ftemp(i)**2 + &
(2.*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* &
ftemp(i)**3 + (-5.*k1(i,k)*i2(i,k)*i3(i,k)+ &
5.*i2(i,k)**3+k1(i,k)**2*i4(i,k))/ &
k1(i,k)**3*ftemp(i)**4
f(i,k) = max(f(i,k),0._r8)
f(i,k) = min(f(i,k),0.0002_r8)
end if
end do
end do
do i = 1,il2g
if (j0(i) < jb(i)) then
if (f(i,j0(i)) < 1.E-6 .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1
end if
end do
do k = msg + 2,pver
do i = 1,il2g
if (k >= jt(i) .and. k <= j0(i)) then
f(i,k) = max(f(i,k),f(i,k-1))
end if
end do
end do
do i = 1,il2g
eps0(i) = f(i,j0(i))
eps(i,jb(i)) = eps0(i)
end do
!
! This is set to match the Rasch and Kristjansson paper
!
do k = pver,msg + 1,-1
do i = 1,il2g
if (k >= j0(i) .and. k <= jb(i)) then
eps(i,k) = f(i,j0(i))
end if
end do
end do
do k = pver,msg + 1,-1
do i = 1,il2g
if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k)
end do
end do
!
! specify the updraft mass flux mu, entrainment eu, detrainment du
! and moist static energy hu.
! here and below mu, eu,du, md and ed are all normalized by mb
!
do i = 1,il2g
if (eps0(i) > 0.) then
mu(i,jb(i)) = 1.
eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i))
end if
end do
do k = pver,msg + 1,-1
do i = 1,il2g
if (eps0(i) > 0. .and. (k >= jt(i) .and. k < jb(i))) then
zuef(i) = zf(i,k) - zf(i,jb(i))
rmue(i) = (1./eps0(i))* (exp(eps(i,k+1)*zuef(i))-1.)/zuef(i)
mu(i,k) = (1./eps0(i))* (exp(eps(i,k )*zuef(i))-1.)/zuef(i)
eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k)
du(i,k) = (rmue(i)-mu(i,k))/dz(i,k)
end if
end do
end do
!
khighest = pverp
klowest = 1
do i=1,il2g
khighest = min(khighest,lel(i))
klowest = max(klowest,jb(i))
end do
do k = klowest-1,khighest,-1
!cdir$ ivdep
do i = 1,il2g
if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0.) then
if (mu(i,k) < 0.01) then
hu(i,k) = hu(i,jb(i))
mu(i,k) = 0.
eu(i,k) = 0.
du(i,k) = mu(i,k+1)/dz(i,k)
else
hu(i,k) = mu(i,k+1)/mu(i,k)*hu(i,k+1) + &
dz(i,k)/mu(i,k)* (eu(i,k)*hmn(i,k)- du(i,k)*hsat(i,k))
end if
end if
end do
end do
!
! reset cloud top index beginning from two layers above the
! cloud base (i.e. if cloud is only one layer thick, top is not reset
!
do i=1,il2g
doit(i) = .true.
end do
do k=klowest-2,khighest-1,-1
do i=1,il2g
if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then
if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) &
.and. mu(i,k) >= 0.02) then
if (hu(i,k)-hsthat(i,k) < -2000.) then
jt(i) = k + 1
doit(i) = .false.
else
jt(i) = k
if (eps0(i) <= 0.) doit(i) = .false.
end if
else if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.01) then
jt(i) = k + 1
doit(i) = .false.
end if
end if
end do
end do
do k = pver,msg + 1,-1
!cdir$ ivdep
do i = 1,il2g
if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0.) then
mu(i,k) = 0.
eu(i,k) = 0.
du(i,k) = 0.
hu(i,k) = hu(i,jb(i))
end if
if (k == jt(i) .and. eps0(i) > 0.) then
du(i,k) = mu(i,k+1)/dz(i,k)
eu(i,k) = 0.
mu(i,k) = 0.
end if
end do
end do
!
! specify downdraft properties (no downdrafts if jd.ge.jb).
! scale down downward mass flux profile so that net flux
! (up-down) at cloud base in not negative.
!
do i = 1,il2g
!
! in normal downdraft strength run alfa=0.2. In test4 alfa=0.1
!
alfa(i) = 0.1
jt(i) = min(jt(i),jb(i)-1)
jd(i) = max(j0(i),jt(i)+1)
jd(i) = min(jd(i),jb(i))
hd(i,jd(i)) = hmn(i,jd(i)-1)
if (jd(i) < jb(i) .and. eps0(i) > 0.) then
epsm(i) = eps0(i)
md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i)
end if
end do
do k = msg + 1,pver
do i = 1,il2g
if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0.) then
zdef(i) = zf(i,jd(i)) - zf(i,k)
md(i,k) = -alfa(i)/ (2.*eps0(i))*(exp(2.*epsm(i)*zdef(i))-1.)/zdef(i)
end if
end do
end do
do k = msg + 1,pver
do i = 1,il2g
if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0. .and. jd(i) < jb(i)) then
ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._r8)
md(i,k) = md(i,k)*ratmjb(i)
end if
end do
end do
small = 1.e-20
do k = msg + 1,pver
do i = 1,il2g
if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0.) then
ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1)
mdt = min(md(i,k),-small)
hd(i,k) = (md(i,k-1)*hd(i,k-1) - dz(i,k-1)*ed(i,k-1)*hmn(i,k-1))/mdt
end if
end do
end do
!
! calculate updraft and downdraft properties.
!
do k = msg + 2,pver
do i = 1,il2g
if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0. .and. jd(i) < jb(i)) then
qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ &
(rl*(1. + gamhat(i,k)))
end if
end do
end do
!
do i = 1,il2g
done(i) = .false.
end do
kount = 0
do k = pver,msg + 2,-1
do i = 1,il2g
if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0.) then
su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + &
dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k)
qu(i,k) = mu(i,k+1)/mu(i,k)*qu(i,k+1) + dz(i,k)/mu(i,k)* (eu(i,k)*q(i,k)- &
du(i,k)*qst(i,k))
tu = su(i,k) - grav/cp*zf(i,k)
estu = c1*exp((c2* (tu-tfreez))/ ((tu-tfreez)+c3))
qstu = eps1*estu/ ((p(i,k)+p(i,k-1))/2.-estu)
if (qu(i,k) >= qstu) then
jlcl(i) = k
kount = kount + 1
done(i) = .true.
end if
end if
end do
if (kount >= il2g) goto 690
end do
690 continue
do k = msg + 2,pver
do i = 1,il2g
if (k == jb(i) .and. eps0(i) > 0.) then
qu(i,k) = q(i,mx(i))
su(i,k) = (hu(i,k)-rl*qu(i,k))/cp
end if
if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0.) then
su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1.+gamhat(i,k)))
qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ &
(rl* (1.+gamhat(i,k)))
end if
end do
end do
! compute condensation in updraft
do k = pver,msg + 2,-1
do i = 1,il2g
if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0.) then
cu(i,k) = ((mu(i,k)*su(i,k)-mu(i,k+1)*su(i,k+1))/ &
dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp)
if (k == jt(i)) cu(i,k) = 0.
cu(i,k) = max(0._r8,cu(i,k))
end if
end do
end do
! compute condensed liquid, rain production rate
! accumulate total precipitation (condensation - detrainment of liquid)
! Note ql1 = ql(k) + rprd(k)*dz(k)/mu(k)
! The differencing is somewhat strange (e.g. du(i,k)*ql(i,k+1)) but is
! consistently applied.
! mu, ql are interface quantities
! cu, du, eu, rprd are midpoint quantites
!!$ c0 = 1.E-3
!!$ c0 = 0.25E-3
!!$ c0 = 0.5E-3
do k = pver,msg + 2,-1
do i = 1,il2g
rprd(i,k) = 0.
if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0. .and. mu(i,k) >= 0.0) then
if (mu(i,k) > 0.) then
ql1 = 1./mu(i,k)* (mu(i,k+1)*ql(i,k+1)- &
dz(i,k)*du(i,k)*ql(i,k+1)+dz(i,k)*cu(i,k))
ql(i,k) = ql1/ (1.+dz(i,k)*c0)
else
ql(i,k) = 0.
end if
totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1))
rprd(i,k) = c0*mu(i,k)*ql(i,k)
end if
end do
end do
!
do i = 1,il2g
qd(i,jd(i)) = qds(i,jd(i))
sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp
end do
!
do k = msg + 2,pver
do i = 1,il2g
if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0.) then
qd(i,k+1) = qds(i,k+1)
evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k)-md(i,k+1)*qd(i,k+1))/dz(i,k)
evp(i,k) = max(evp(i,k),0._r8)
mdt = min(md(i,k+1),-small)
sd(i,k+1) = ((rl/cp*evp(i,k)-ed(i,k)*s(i,k))*dz(i,k) + md(i,k)*sd(i,k))/mdt
totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
end if
end do
end do
do i = 1,il2g
!*guang totevp(i) = totevp(i) + md(i,jd(i))*q(i,jd(i)-1) -
totevp(i) = totevp(i) + md(i,jd(i))*qd(i,jd(i)) - md(i,jb(i))*qd(i,jb(i))
end do
!!$ if (.true.) then
if (.false.) then
do i = 1,il2g
k = jb(i)
if (eps0(i) > 0.) then
evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k))/dz(i,k)
evp(i,k) = max(evp(i,k),0._r8)
totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
end if
end do
endif
do i = 1,il2g
totpcp(i) = max(totpcp(i),0._r8)
totevp(i) = max(totevp(i),0._r8)
end do
!
do k = msg + 2,pver
do i = 1,il2g
if (totevp(i) > 0. .and. totpcp(i) > 0.) then
md(i,k) = md (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
ed(i,k) = ed (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
evp(i,k) = evp(i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
else
md(i,k) = 0.
ed(i,k) = 0.
evp(i,k) = 0.
end if
! cmeg is the cloud water condensed - rain water evaporated
! rprd is the cloud water converted to rain - (rain evaporated)
cmeg(i,k) = cu(i,k) - evp(i,k)
rprd(i,k) = rprd(i,k)-evp(i,k)
end do
end do
! compute the net precipitation flux across interfaces
pflx(:il2g,1) = 0.
do k = 2,pverp
do i = 1,il2g
pflx(i,k) = pflx(i,k-1) + rprd(i,k-1)*dz(i,k-1)
end do
end do
!
do k = msg + 1,pver
do i = 1,il2g
mc(i,k) = mu(i,k) + md(i,k)
end do
end do
!
return
end subroutine cldprp
subroutine closure(lchnk , & 1,1
q ,t ,p ,z ,s , &
tp ,qs ,qu ,su ,mc , &
du ,mu ,md ,qd ,sd , &
qhat ,shat ,dp ,qstp ,zf , &
ql ,dsubcld ,mb ,cape ,tl , &
lcl ,lel ,jt ,mx ,il1g , &
il2g ,rd ,grav ,cp ,rl , &
msg ,capelmt )
!-----------------------------------------------------------------------
!
! Purpose:
! <Say what the routine does>
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: G. Zhang and collaborators. CCM contact:P. Rasch
! This is contributed code not fully standardized by the CCM core group.
!
! this code is very much rougher than virtually anything else in the CCM
! We expect to release cleaner code in a future release
!
! the documentation has been enhanced to the degree that we are able
!
!-----------------------------------------------------------------------
use dycore, only: dycore_is, get_resolution
implicit none
!
!-----------------------------Arguments---------------------------------
!
integer, intent(in) :: lchnk ! chunk identifier
real(r8), intent(inout) :: q(pcols,pver) ! spec humidity
real(r8), intent(inout) :: t(pcols,pver) ! temperature
real(r8), intent(inout) :: p(pcols,pver) ! pressure (mb)
real(r8), intent(inout) :: mb(pcols) ! cloud base mass flux
real(r8), intent(in) :: z(pcols,pver) ! height (m)
real(r8), intent(in) :: s(pcols,pver) ! normalized dry static energy
real(r8), intent(in) :: tp(pcols,pver) ! parcel temp
real(r8), intent(in) :: qs(pcols,pver) ! sat spec humidity
real(r8), intent(in) :: qu(pcols,pver) ! updraft spec. humidity
real(r8), intent(in) :: su(pcols,pver) ! normalized dry stat energy of updraft
real(r8), intent(in) :: mc(pcols,pver) ! net convective mass flux
real(r8), intent(in) :: du(pcols,pver) ! detrainment from updraft
real(r8), intent(in) :: mu(pcols,pver) ! mass flux of updraft
real(r8), intent(in) :: md(pcols,pver) ! mass flux of downdraft
real(r8), intent(in) :: qd(pcols,pver) ! spec. humidity of downdraft
real(r8), intent(in) :: sd(pcols,pver) ! dry static energy of downdraft
real(r8), intent(in) :: qhat(pcols,pver) ! environment spec humidity at interfaces
real(r8), intent(in) :: shat(pcols,pver) ! env. normalized dry static energy at intrfcs
real(r8), intent(in) :: dp(pcols,pver) ! pressure thickness of layers
real(r8), intent(in) :: qstp(pcols,pver) ! spec humidity of parcel
real(r8), intent(in) :: zf(pcols,pver+1) ! height of interface levels
real(r8), intent(in) :: ql(pcols,pver) ! liquid water mixing ratio
real(r8), intent(in) :: cape(pcols) ! available pot. energy of column
real(r8), intent(in) :: tl(pcols)
real(r8), intent(in) :: dsubcld(pcols) ! thickness of subcloud layer
integer, intent(in) :: lcl(pcols) ! index of lcl
integer, intent(in) :: lel(pcols) ! index of launch leve
integer, intent(in) :: jt(pcols) ! top of updraft
integer, intent(in) :: mx(pcols) ! base of updraft
!
!--------------------------Local variables------------------------------
!
real(r8) dtpdt(pcols,pver)
real(r8) dqsdtp(pcols,pver)
real(r8) dtmdt(pcols,pver)
real(r8) dqmdt(pcols,pver)
real(r8) dboydt(pcols,pver)
real(r8) thetavp(pcols,pver)
real(r8) thetavm(pcols,pver)
real(r8) dtbdt(pcols),dqbdt(pcols),dtldt(pcols)
real(r8) beta
real(r8) capelmt
real(r8) cp
real(r8) dadt(pcols)
real(r8) debdt
real(r8) dltaa
real(r8) eb
real(r8) grav
integer i
integer il1g
integer il2g
integer k, kmin, kmax
integer msg
real(r8) rd
real(r8) rl
! change of subcloud layer properties due to convection is
! related to cumulus updrafts and downdrafts.
! mc(z)=f(z)*mb, mub=betau*mb, mdb=betad*mb are used
! to define betau, betad and f(z).
! note that this implies all time derivatives are in effect
! time derivatives per unit cloud-base mass flux, i.e. they
! have units of 1/mb instead of 1/sec.
!
do i = il1g,il2g
mb(i) = 0.
eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
dtbdt(i) = (1./dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ &
md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i))))
dqbdt(i) = (1./dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ &
md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i))))
debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i)
dtldt(i) = -2840.* (3.5/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
(3.5*log(t(i,mx(i)))-log(eb)-4.805)**2
end do
!
! dtmdt and dqmdt are cumulus heating and drying.
!
do k = msg + 1,pver
do i = il1g,il2g
dtmdt(i,k) = 0.
dqmdt(i,k) = 0.
end do
end do
!
do k = msg + 1,pver - 1
do i = il1g,il2g
if (k == jt(i)) then
dtmdt(i,k) = (1./dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- &
rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)))
dqmdt(i,k) = (1./dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- &
qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1)))
end if
end do
end do
!
beta = 0.
do k = msg + 1,pver - 1
do i = il1g,il2g
if (k > jt(i) .and. k < mx(i)) then
dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ &
dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1))
! dqmdt(i,k)=(mc(i,k)*(qhat(i,k)-q(i,k))
! 1 +mc(i,k+1)*(q(i,k)-qhat(i,k+1)))/dp(i,k)
! 2 +du(i,k)*(qs(i,k)-q(i,k))
! 3 +du(i,k)*(beta*ql(i,k)+(1-beta)*ql(i,k+1))
dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- &
mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* &
(qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* &
(qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + &
du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1))
end if
end do
end do
!
do k = msg + 1,pver
do i = il1g,il2g
if (k >= lel(i) .and. k <= lcl(i)) then
thetavp(i,k) = tp(i,k)* (1000./p(i,k))** (rd/cp)*(1.+1.608*qstp(i,k)-q(i,mx(i)))
thetavm(i,k) = t(i,k)* (1000./p(i,k))** (rd/cp)*(1.+0.608*q(i,k))
dqsdtp(i,k) = qstp(i,k)* (1.+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2)
!
! dtpdt is the parcel temperature change due to change of
! subcloud layer properties during convection.
!
dtpdt(i,k) = tp(i,k)/ (1.+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* &
(dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ &
tl(i)**2*dtldt(i)))
!
! dboydt is the integrand of cape change.
!
dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1./(1.+1.608*qstp(i,k)-q(i,mx(i)))* &
(1.608 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608/ &
(1.+0.608*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k)
end if
end do
end do
!
do k = msg + 1,pver
do i = il1g,il2g
if (k > lcl(i) .and. k < mx(i)) then
thetavp(i,k) = tp(i,k)* (1000./p(i,k))** (rd/cp)*(1.+0.608*q(i,mx(i)))
thetavm(i,k) = t(i,k)* (1000./p(i,k))** (rd/cp)*(1.+0.608*q(i,k))
!
! dboydt is the integrand of cape change.
!
dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608/ (1.+0.608*q(i,mx(i)))*dqbdt(i)- &
dtmdt(i,k)/t(i,k)-0.608/ (1.+0.608*q(i,k))*dqmdt(i,k))* &
grav*thetavp(i,k)/thetavm(i,k)
end if
end do
end do
!
! buoyant energy change is set to 2/3*excess cape per 3 hours
!
dadt(il1g:il2g) = 0.
kmin = minval(lel(il1g:il2g))
kmax = maxval(mx(il1g:il2g)) - 1
do k = kmin, kmax
do i = il1g,il2g
if ( k >= lel(i) .and. k <= mx(i) - 1) then
dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1))
endif
end do
end do
do i = il1g,il2g
dltaa = -1.* (cape(i)-capelmt)
if (dadt(i) /= 0.) mb(i) = max(dltaa/tau/dadt(i),0._r8)
end do
!
return
end subroutine closure
subroutine q1q2_pjr(lchnk , & 1
dqdt ,dsdt ,q ,qs ,qu , &
su ,du ,qhat ,shat ,dp , &
mu ,md ,sd ,qd ,ql , &
dsubcld ,jt ,mx ,il1g ,il2g , &
cp ,rl ,msg , &
dl ,evp ,cu )
implicit none
!-----------------------------------------------------------------------
!
! Purpose:
! <Say what the routine does>
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: phil rasch dec 19 1995
!
!-----------------------------------------------------------------------
real(r8), intent(in) :: cp
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: il1g
integer, intent(in) :: il2g
integer, intent(in) :: msg
real(r8), intent(in) :: q(pcols,pver)
real(r8), intent(in) :: qs(pcols,pver)
real(r8), intent(in) :: qu(pcols,pver)
real(r8), intent(in) :: su(pcols,pver)
real(r8), intent(in) :: du(pcols,pver)
real(r8), intent(in) :: qhat(pcols,pver)
real(r8), intent(in) :: shat(pcols,pver)
real(r8), intent(in) :: dp(pcols,pver)
real(r8), intent(in) :: mu(pcols,pver)
real(r8), intent(in) :: md(pcols,pver)
real(r8), intent(in) :: sd(pcols,pver)
real(r8), intent(in) :: qd(pcols,pver)
real(r8), intent(in) :: ql(pcols,pver)
real(r8), intent(in) :: evp(pcols,pver)
real(r8), intent(in) :: cu(pcols,pver)
real(r8), intent(in) :: dsubcld(pcols)
real(r8),intent(out) :: dqdt(pcols,pver),dsdt(pcols,pver)
real(r8),intent(out) :: dl(pcols,pver)
integer kbm
integer ktm
integer jt(pcols)
integer mx(pcols)
!
! work fields:
!
integer i
integer k
real(r8) emc
real(r8) rl
!-------------------------------------------------------------------
do k = msg + 1,pver
do i = il1g,il2g
dsdt(i,k) = 0.
dqdt(i,k) = 0.
dl(i,k) = 0.
end do
end do
!
! find the highest level top and bottom levels of convection
!
ktm = pver
kbm = pver
do i = il1g, il2g
ktm = min(ktm,jt(i))
kbm = min(kbm,mx(i))
end do
do k = ktm,pver-1
do i = il1g,il2g
emc = -cu (i,k) & ! condensation in updraft
+evp(i,k) ! evaporating rain in downdraft
dsdt(i,k) = -rl/cp*emc &
+ (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) &
-mu(i,k)* (su(i,k)-shat(i,k)) &
+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) &
-md(i,k)* (sd(i,k)-shat(i,k)) &
)/dp(i,k)
dqdt(i,k) = emc + &
(+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) &
-mu(i,k)* (qu(i,k)-qhat(i,k)) &
+md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) &
-md(i,k)* (qd(i,k)-qhat(i,k)) &
)/dp(i,k)
dl(i,k) = du(i,k)*ql(i,k+1)
end do
end do
!
do k = kbm,pver
do i = il1g,il2g
if (k == mx(i)) then
dsdt(i,k) = (1./dsubcld(i))* &
(-mu(i,k)* (su(i,k)-shat(i,k)) &
-md(i,k)* (sd(i,k)-shat(i,k)) &
)
dqdt(i,k) = (1./dsubcld(i))* &
(-mu(i,k)*(qu(i,k)-qhat(i,k)) &
-md(i,k)*(qd(i,k)-qhat(i,k)) &
)
else if (k > mx(i)) then
dsdt(i,k) = dsdt(i,k-1)
dqdt(i,k) = dqdt(i,k-1)
end if
end do
end do
!
return
end subroutine q1q2_pjr
end module zm_conv