#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