module hk_conv 3,4
!
! Moist convection. Primarily data used by both Zhang-McFarlane convection
! and Hack shallow convective schemes.
!
! $Id$
!
use shr_kind_mod
, only: r8 => shr_kind_r8
use cam_logfile, only: iulog
use spmd_utils, only: masterproc
use abortutils, only: endrun
implicit none
private
save
!
! Public interfaces
!
public mfinti ! Initialization of data for moist convection
public cmfmca ! Hack shallow convection
public hkconv_readnl ! read hkconv_nl namelist
!
! Private data used for Hack shallow convection
!
real(r8), parameter :: unset_r8 = huge(1.0_r8)
! Namelist variables
real(r8) :: hkconv_c0 = unset_r8
real(r8) :: hkconv_cmftau = unset_r8
real(r8) :: hlat ! latent heat of vaporization
real(r8) :: c0 ! rain water autoconversion coefficient set from namelist input hkconv_c0
real(r8) :: betamn ! minimum overshoot parameter
real(r8) :: rhlat ! reciprocal of hlat
real(r8) :: rcp ! reciprocal of cp
real(r8) :: cmftau ! characteristic adjustment time scale set from namelist input hkconv_cmftau
real(r8) :: rhoh2o ! density of liquid water (STP)
real(r8) :: dzmin ! minimum convective depth for precipitation
real(r8) :: tiny ! arbitrary small num used in transport estimates
real(r8) :: eps ! convergence criteria (machine dependent)
real(r8) :: tpmax ! maximum acceptable t perturbation (degrees C)
real(r8) :: shpmax ! maximum acceptable q perturbation (g/g)
integer :: iloc ! longitude location for diagnostics
integer :: jloc ! latitude location for diagnostics
integer :: nsloc ! nstep for which to produce diagnostics
!
logical :: rlxclm ! logical to relax column versus cloud triplet
real(r8) cp ! specific heat of dry air
real(r8) grav ! gravitational constant
real(r8) rgrav ! reciprocal of grav
real(r8) rgas ! gas constant for dry air
integer limcnv ! top interface level limit for convection
contains
subroutine hkconv_readnl(nlfile) 1,9
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
! Local variables
integer :: unitn, ierr
character(len=*), parameter :: subname = 'hkconv_readnl'
namelist /hkconv_nl/ hkconv_cmftau, hkconv_c0
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'hkconv_nl', status=ierr)
if (ierr == 0) then
read(unitn, hkconv_nl, iostat=ierr)
if (ierr /= 0) then
call endrun(subname // ':: ERROR reading namelist')
end if
end if
close(unitn)
call freeunit(unitn)
! set local variables
cmftau = hkconv_cmftau
c0 = hkconv_c0
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast(cmftau, 1, mpir8, 0, mpicom)
call mpibcast(c0, 1, mpir8, 0, mpicom)
#endif
end subroutine hkconv_readnl
!================================================================================================
subroutine mfinti (rair ,cpair ,gravit ,latvap ,rhowtr,limcnv_in ) 1,2
!-----------------------------------------------------------------------
!
! Purpose:
! Initialize moist convective mass flux procedure common block, cmfmca
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: J. Hack
!
!-----------------------------------------------------------------------
use dycore, only: dycore_is, get_resolution
use spmd_utils, only: masterproc
!------------------------------Arguments--------------------------------
!
! Input arguments
!
real(r8), intent(in) :: rair ! gas constant for dry air
real(r8), intent(in) :: cpair ! specific heat of dry air
real(r8), intent(in) :: gravit ! acceleration due to gravity
real(r8), intent(in) :: latvap ! latent heat of vaporization
real(r8), intent(in) :: rhowtr ! density of liquid water (STP)
integer, intent(in) :: limcnv_in ! top interface level limit for convection
! local variables
character(len=32) :: hgrid ! horizontal grid specifier
!
!-----------------------------------------------------------------------
!
! Initialize physical constants for moist convective mass flux procedure
!
cp = cpair ! specific heat of dry air
hlat = latvap ! latent heat of vaporization
grav = gravit ! gravitational constant
rgas = rair ! gas constant for dry air
rhoh2o = rhowtr ! density of liquid water (STP)
limcnv = limcnv_in
! Initialize free parameters for moist convective mass flux procedure
! cmftau - characteristic adjustment time scale
! c0 - rain water autoconversion coeff (1/m)
if (masterproc) then
write(iulog,*) 'tuning parameters hk_conv: cmftau',cmftau
write(iulog,*) 'tuning parameters hk_conv: c0',c0
endif
dzmin = 0.0_r8 ! minimum cloud depth to precipitate (m)
betamn = 0.10_r8 ! minimum overshoot parameter
tpmax = 1.50_r8 ! maximum acceptable t perturbation (deg C)
shpmax = 1.50e-3_r8 ! maximum acceptable q perturbation (g/g)
rlxclm = .true. ! logical variable to specify that relaxation
! time scale should applied to column as
! opposed to triplets individually
!
! Initialize miscellaneous (frequently used) constants
!
rhlat = 1.0_r8/hlat ! reciprocal latent heat of vaporization
rcp = 1.0_r8/cp ! reciprocal specific heat of dry air
rgrav = 1.0_r8/grav ! reciprocal gravitational constant
!
! Initialize diagnostic location information for moist convection scheme
!
iloc = 1 ! longitude point for diagnostic info
jloc = 1 ! latitude point for diagnostic info
nsloc = 1 ! nstep value at which to begin diagnostics
!
! Initialize other miscellaneous parameters
!
tiny = 1.0e-36_r8 ! arbitrary small number (scalar transport)
eps = 1.0e-13_r8 ! convergence criteria (machine dependent)
!
return
end subroutine mfinti
subroutine cmfmca(lchnk ,ncol , & 1,11
nstep ,ztodt ,pmid ,pdel , &
rpdel ,zm ,tpert ,qpert ,phis , &
pblht ,t ,q ,cmfdt ,dq , &
cmfmc ,cmfdqr ,cmfsl ,cmflq ,precc , &
qc ,cnt ,cnb ,icwmr ,rliq , &
pmiddry ,pdeldry ,rpdeldry)
!-----------------------------------------------------------------------
!
! Purpose:
! Moist convective mass flux procedure:
!
! Method:
! If stratification is unstable to nonentraining parcel ascent,
! complete an adjustment making successive use of a simple cloud model
! consisting of three layers (sometimes referred to as a triplet)
!
! Code generalized to allow specification of parcel ("updraft")
! properties, as well as convective transport of an arbitrary
! number of passive constituents (see q array). The code
! is written so the water vapor field is passed independently
! in the calling list from the block of other transported
! constituents, even though as currently designed, it is the
! first component in the constituents field.
!
! Author: J. Hack
!
! BAB: changed code to report tendencies in cmfdt and dq, instead of
! updating profiles. Cmfdq contains water only, made it a local variable
! made dq (all constituents) the argument.
!
!-----------------------------------------------------------------------
!#######################################################################
!# #
!# Debugging blocks are marked this way for easy identification #
!# #
!#######################################################################
use constituents, only: pcnst
use constituents, only: cnst_get_type_byind
use ppgrid, only: pcols, pver, pverp
use phys_grid, only: get_lat_all_p, get_lon_all_p
use wv_saturation, only: aqsatd, vqsatd
real(r8) ssfac ! supersaturation bound (detrained air)
parameter (ssfac = 1.001_r8)
!------------------------------Arguments--------------------------------
!
! Input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: nstep ! current time step index
real(r8), intent(in) :: ztodt ! 2 delta-t (seconds)
real(r8), intent(in) :: pmid(pcols,pver) ! pressure
real(r8), intent(in) :: pdel(pcols,pver) ! delta-p
real(r8), intent(in) :: pmiddry(pcols,pver) ! pressure
real(r8), intent(in) :: pdeldry(pcols,pver) ! delta-p
real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel
real(r8), intent(in) :: rpdeldry(pcols,pver) ! 1./pdel
real(r8), intent(in) :: zm(pcols,pver) ! height abv sfc at midpoints
real(r8), intent(in) :: tpert(pcols) ! PBL perturbation theta
real(r8), intent(in) :: qpert(pcols,pcnst) ! PBL perturbation specific humidity
real(r8), intent(in) :: phis(pcols) ! surface geopotential
real(r8), intent(in) :: pblht(pcols) ! PBL height (provided by PBL routine)
real(r8), intent(in) :: t(pcols,pver) ! temperature (t bar)
real(r8), intent(in) :: q(pcols,pver,pcnst) ! specific humidity (sh bar)
!
! Output arguments
!
real(r8), intent(out) :: cmfdt(pcols,pver) ! dt/dt due to moist convection
real(r8), intent(out) :: cmfmc(pcols,pverp) ! moist convection cloud mass flux
real(r8), intent(out) :: cmfdqr(pcols,pver) ! dq/dt due to convective rainout
real(r8), intent(out) :: cmfsl(pcols,pver ) ! convective lw static energy flux
real(r8), intent(out) :: cmflq(pcols,pver ) ! convective total water flux
real(r8), intent(out) :: precc(pcols) ! convective precipitation rate
! JJH mod to explicitly export cloud water
real(r8), intent(out) :: qc(pcols,pver) ! dq/dt due to export of cloud water
real(r8), intent(out) :: cnt(pcols) ! top level of convective activity
real(r8), intent(out) :: cnb(pcols) ! bottom level of convective activity
real(r8), intent(out) :: dq(pcols,pver,pcnst) ! constituent tendencies
real(r8), intent(out) :: icwmr(pcols,pver)
real(r8), intent(out) :: rliq(pcols)
!
!---------------------------Local workspace-----------------------------
!
real(r8) pm(pcols,pver) ! pressure
real(r8) pd(pcols,pver) ! delta-p
real(r8) rpd(pcols,pver) ! 1./pdel
real(r8) cmfdq(pcols,pver) ! dq/dt due to moist convection
real(r8) gam(pcols,pver) ! 1/cp (d(qsat)/dT)
real(r8) sb(pcols,pver) ! dry static energy (s bar)
real(r8) hb(pcols,pver) ! moist static energy (h bar)
real(r8) shbs(pcols,pver) ! sat. specific humidity (sh bar star)
real(r8) hbs(pcols,pver) ! sat. moist static energy (h bar star)
real(r8) shbh(pcols,pverp) ! specific humidity on interfaces
real(r8) sbh(pcols,pverp) ! s bar on interfaces
real(r8) hbh(pcols,pverp) ! h bar on interfaces
real(r8) cmrh(pcols,pverp) ! interface constituent mixing ratio
real(r8) prec(pcols) ! instantaneous total precipitation
real(r8) dzcld(pcols) ! depth of convective layer (m)
real(r8) beta(pcols) ! overshoot parameter (fraction)
real(r8) betamx(pcols) ! local maximum on overshoot
real(r8) eta(pcols) ! convective mass flux (kg/m^2 s)
real(r8) etagdt(pcols) ! eta*grav*dt
real(r8) cldwtr(pcols) ! cloud water (mass)
real(r8) rnwtr(pcols) ! rain water (mass)
! JJH extension to facilitate export of cloud liquid water
real(r8) totcond(pcols) ! total condensate; mix of precip and cloud water (mass)
real(r8) sc (pcols) ! dry static energy ("in-cloud")
real(r8) shc (pcols) ! specific humidity ("in-cloud")
real(r8) hc (pcols) ! moist static energy ("in-cloud")
real(r8) cmrc(pcols) ! constituent mix rat ("in-cloud")
real(r8) dq1(pcols) ! shb convective change (lower lvl)
real(r8) dq2(pcols) ! shb convective change (mid level)
real(r8) dq3(pcols) ! shb convective change (upper lvl)
real(r8) ds1(pcols) ! sb convective change (lower lvl)
real(r8) ds2(pcols) ! sb convective change (mid level)
real(r8) ds3(pcols) ! sb convective change (upper lvl)
real(r8) dcmr1(pcols) ! q convective change (lower lvl)
real(r8) dcmr2(pcols) ! q convective change (mid level)
real(r8) dcmr3(pcols) ! q convective change (upper lvl)
real(r8) estemp(pcols,pver) ! saturation vapor pressure (scratch)
real(r8) vtemp1(2*pcols) ! intermediate scratch vector
real(r8) vtemp2(2*pcols) ! intermediate scratch vector
real(r8) vtemp3(2*pcols) ! intermediate scratch vector
real(r8) vtemp4(2*pcols) ! intermediate scratch vector
integer indx1(pcols) ! longitude indices for condition true
logical etagt0 ! true if eta > 0.0
real(r8) sh1 ! dummy arg in qhalf statement func.
real(r8) sh2 ! dummy arg in qhalf statement func.
real(r8) shbs1 ! dummy arg in qhalf statement func.
real(r8) shbs2 ! dummy arg in qhalf statement func.
real(r8) cats ! modified characteristic adj. time
real(r8) rtdt ! 1./ztodt
real(r8) qprime ! modified specific humidity pert.
real(r8) tprime ! modified thermal perturbation
real(r8) pblhgt ! bounded pbl height (max[pblh,1m])
real(r8) fac1 ! intermediate scratch variable
real(r8) shprme ! intermediate specific humidity pert.
real(r8) qsattp ! sat mix rat for thermally pert PBL parcels
real(r8) dz ! local layer depth
real(r8) temp1 ! intermediate scratch variable
real(r8) b1 ! bouyancy measure in detrainment lvl
real(r8) b2 ! bouyancy measure in condensation lvl
real(r8) temp2 ! intermediate scratch variable
real(r8) temp3 ! intermediate scratch variable
real(r8) g ! bounded vertical gradient of hb
real(r8) tmass ! total mass available for convective exch
real(r8) denom ! intermediate scratch variable
real(r8) qtest1 ! used in negative q test (middle lvl)
real(r8) qtest2 ! used in negative q test (lower lvl)
real(r8) fslkp ! flux lw static energy (bot interface)
real(r8) fslkm ! flux lw static energy (top interface)
real(r8) fqlkp ! flux total water (bottom interface)
real(r8) fqlkm ! flux total water (top interface)
real(r8) botflx ! bottom constituent mixing ratio flux
real(r8) topflx ! top constituent mixing ratio flux
real(r8) efac1 ! ratio q to convectively induced chg (btm lvl)
real(r8) efac2 ! ratio q to convectively induced chg (mid lvl)
real(r8) efac3 ! ratio q to convectively induced chg (top lvl)
real(r8) tb(pcols,pver) ! working storage for temp (t bar)
real(r8) shb(pcols,pver) ! working storage for spec hum (sh bar)
real(r8) adjfac ! adjustment factor (relaxation related)
real(r8) rktp
real(r8) rk
#if ( defined DIAGNS )
!
! Following 7 real variables are used in diagnostics calculations
!
real(r8) rh ! relative humidity
real(r8) es ! sat vapor pressure
real(r8) hsum1 ! moist static energy integral
real(r8) qsum1 ! total water integral
real(r8) hsum2 ! final moist static energy integral
real(r8) qsum2 ! final total water integral
real(r8) fac ! intermediate scratch variable
#endif
integer i,k ! longitude, level indices
integer ii ! index on "gathered" vectors
integer len1 ! vector length of "gathered" vectors
integer m ! constituent index
integer ktp ! tmp indx used to track top of convective layer
#if ( defined DIAGNS )
integer n ! vertical index (diagnostics)
integer kp ! vertical index (diagnostics)
integer kpp ! index offset, kp+1 (diagnostics)
integer kpm1 ! index offset, kp-1 (diagnostics)
integer lat(pcols) ! latitude indices
integer lon(pcols) ! longitude indices
#endif
!
!---------------------------Statement functions-------------------------
!
real(r8) qhalf
qhalf(sh1,sh2,shbs1,shbs2) = min(max(sh1,sh2),(shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2))
!
!-----------------------------------------------------------------------
!** BAB initialize output tendencies here
! copy q to dq; use dq below for passive tracer transport
cmfdt(:ncol,:) = 0._r8
cmfdq(:ncol,:) = 0._r8
dq(:ncol,:,2:) = q(:ncol,:,2:)
cmfmc(:ncol,:) = 0._r8
cmfdqr(:ncol,:) = 0._r8
cmfsl(:ncol,:) = 0._r8
cmflq(:ncol,:) = 0._r8
qc(:ncol,:) = 0._r8
rliq(:ncol) = 0._r8
!
#if ( defined DIAGNS )
! Determine chunk latitudes and longitudes
call get_lat_all_p(lchnk, ncol, lat)
call get_lon_all_p(lchnk, ncol, lon)
#endif
!
! Ensure that characteristic adjustment time scale (cmftau) assumed
! in estimate of eta isn't smaller than model time scale (ztodt)
! The time over which the convection is assumed to act (the adjustment
! time scale) can be applied with each application of the three-level
! cloud model, or applied to the column tendencies after a "hard"
! adjustment (i.e., on a 2-delta t time scale) is evaluated
!
if (rlxclm) then
cats = ztodt ! relaxation applied to column
adjfac = ztodt/(max(ztodt,cmftau))
else
cats = max(ztodt,cmftau) ! relaxation applied to triplet
adjfac = 1.0_r8
endif
rtdt = 1.0_r8/ztodt
!
! Move temperature and moisture into working storage
!
do k=limcnv,pver
do i=1,ncol
tb (i,k) = t(i,k)
shb(i,k) = q(i,k,1)
end do
end do
do k=1,pver
do i=1,ncol
icwmr(i,k) = 0._r8
end do
end do
!
! Compute sb,hb,shbs,hbs
!
call aqsatd(tb ,pmid ,estemp ,shbs ,gam , &
pcols ,ncol ,pver ,limcnv ,pver )
!
do k=limcnv,pver
do i=1,ncol
sb (i,k) = cp*tb(i,k) + zm(i,k)*grav + phis(i)
hb (i,k) = sb(i,k) + hlat*shb(i,k)
hbs(i,k) = sb(i,k) + hlat*shbs(i,k)
end do
end do
!
! Compute sbh, shbh
!
do k=limcnv+1,pver
do i=1,ncol
sbh (i,k) = 0.5_r8*(sb(i,k-1) + sb(i,k))
shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
end do
end do
!
! Specify properties at top of model (not used, but filling anyway)
!
do i=1,ncol
sbh (i,limcnv) = sb(i,limcnv)
shbh(i,limcnv) = shb(i,limcnv)
hbh (i,limcnv) = hb(i,limcnv)
end do
!
! Zero vertically independent control, tendency & diagnostic arrays
!
do i=1,ncol
prec(i) = 0.0_r8
dzcld(i) = 0.0_r8
cnb(i) = 0.0_r8
cnt(i) = real(pver+1,r8)
end do
#if ( defined DIAGNS )
!#######################################################################
!# #
!# output initial thermodynamic profile if debug diagnostics #
!# #
do i=1,ncol
if ((lat(i).eq.jloc) .and. (lon(i).eq.iloc) &
.and. (nstep.ge.nsloc)) then
!# #
!# approximate vertical integral of moist static energy #
!# and total preciptable water #
!# #
hsum1 = 0.0_r8
qsum1 = 0.0_r8
do k=limcnv,pver
hsum1 = hsum1 + pdel(i,k)*rgrav*hb(i,k)
qsum1 = qsum1 + pdel(i,k)*rgrav*shb(i,k)
end do
!# #
write(iulog,8010)
fac = grav*864._r8
do k=limcnv,pver
rh = shb(i,k)/shbs(i,k)
write(iulog,8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k),cmfsl(i,k), cmflq(i,k)
write(iulog,8040) tb(i,k),shb(i,k),rh,sb(i,k),hb(i,k),hbs(i,k),ztodt*cmfdt(i,k), &
ztodt*cmfdq(i,k),ztodt*cmfdqr(i,k)
end do
write(iulog, 8000) prec(i)
end if
enddo
#endif
!# #
!# #
!#######################################################################
!
! Begin moist convective mass flux adjustment procedure.
! Formalism ensures that negative cloud liquid water can never occur
!
do 70 k=pver-1,limcnv+1,-1
do 10 i=1,ncol
etagdt(i) = 0.0_r8
eta (i) = 0.0_r8
beta (i) = 0.0_r8
ds1 (i) = 0.0_r8
ds2 (i) = 0.0_r8
ds3 (i) = 0.0_r8
dq1 (i) = 0.0_r8
dq2 (i) = 0.0_r8
dq3 (i) = 0.0_r8
!
! Specification of "cloud base" conditions
!
qprime = 0.0_r8
tprime = 0.0_r8
!
! Assign tprime within the PBL to be proportional to the quantity
! tpert (which will be bounded by tpmax), passed to this routine by
! the PBL routine. Don't allow perturbation to produce a dry
! adiabatically unstable parcel. Assign qprime within the PBL to be
! an appropriately modified value of the quantity qpert (which will be
! bounded by shpmax) passed to this routine by the PBL routine. The
! quantity qprime should be less than the local saturation value
! (qsattp=qsat[t+tprime,p]). In both cases, tpert and qpert are
! linearly reduced toward zero as the PBL top is approached.
!
pblhgt = max(pblht(i),1.0_r8)
if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_r8 ) then
fac1 = max(0.0_r8,1.0_r8-zm(i,k+1)/pblhgt)
tprime = min(tpert(i),tpmax)*fac1
qsattp = shbs(i,k+1) + cp*rhlat*gam(i,k+1)*tprime
shprme = min(min(qpert(i,1),shpmax)*fac1,max(qsattp-shb(i,k+1),0.0_r8))
qprime = max(qprime,shprme)
else
tprime = 0.0_r8
qprime = 0.0_r8
end if
!
! Specify "updraft" (in-cloud) thermodynamic properties
!
sc (i) = sb (i,k+1) + cp*tprime
shc(i) = shb(i,k+1) + qprime
hc (i) = sc (i ) + hlat*shc(i)
vtemp4(i) = hc(i) - hbs(i,k)
dz = pdel(i,k)*rgas*tb(i,k)*rgrav/pmid(i,k)
if (vtemp4(i) > 0.0_r8) then
dzcld(i) = dzcld(i) + dz
else
dzcld(i) = 0.0_r8
end if
10 continue
#if ( defined DIAGNS )
!#######################################################################
!# #
!# output thermodynamic perturbation information #
!# #
do i=1,ncol
if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
write(iulog,8090) k+1,sc(iloc),shc(iloc),hc(iloc)
end if
enddo
!# #
!#######################################################################
#endif
!
! Check on moist convective instability
! Build index vector of points where instability exists
!
len1 = 0
do i=1,ncol
if (vtemp4(i) > 0.0_r8) then
len1 = len1 + 1
indx1(len1) = i
end if
end do
if (len1 <= 0) go to 70
!
! Current level just below top level => no overshoot
!
if (k <= limcnv+1) then
do ii=1,len1
i = indx1(ii)
temp1 = vtemp4(i)/(1.0_r8 + gam(i,k))
cldwtr(i) = max(0.0_r8,(sb(i,k) - sc(i) + temp1))
beta(i) = 0.0_r8
vtemp3(i) = (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k))
end do
else
!
! First guess at overshoot parameter using crude buoyancy closure
! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
! If pre-existing supersaturation in detrainment layer, beta=0
! cldwtr is temporarily equal to hlat*l (l=> liquid water)
!
!cdir nodep
!DIR$ CONCURRENT
do ii=1,len1
i = indx1(ii)
temp1 = vtemp4(i)/(1.0_r8 + gam(i,k))
cldwtr(i) = max(0.0_r8,(sb(i,k)-sc(i)+temp1))
betamx(i) = 1.0_r8 - c0*max(0.0_r8,(dzcld(i)-dzmin))
b1 = (hc(i) - hbs(i,k-1))*pdel(i,k-1)
b2 = (hc(i) - hbs(i,k ))*pdel(i,k )
beta(i) = max(betamn,min(betamx(i), 1.0_r8 + b1/b2))
if (hbs(i,k-1) <= hb(i,k-1)) beta(i) = 0.0_r8
!
! Bound maximum beta to ensure physically realistic solutions
!
! First check constrains beta so that eta remains positive
! (assuming that eta is already positive for beta equal zero)
!
vtemp1(i) = -(hbh(i,k+1) - hc(i))*pdel(i,k)*rpdel(i,k+1)+ &
(1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i))
vtemp2(i) = (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k))
vtemp3(i) = vtemp2(i)
if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
betamx(i) = 0.99_r8*(vtemp1(i)/vtemp2(i))
beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
end if
end do
!
! Second check involves supersaturation of "detrainment layer"
! small amount of supersaturation acceptable (by ssfac factor)
!
!cdir nodep
!DIR$ CONCURRENT
do ii=1,len1
i = indx1(ii)
if (hb(i,k-1) < hbs(i,k-1)) then
vtemp1(i) = vtemp1(i)*rpdel(i,k)
temp2 = gam(i,k-1)*(sbh(i,k) - sc(i) + cldwtr(i)) - &
hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
temp3 = vtemp3(i)*rpdel(i,k)
vtemp2(i) = (ztodt/cats)*(hc(i) - hbs(i,k))*temp2/ &
(pdel(i,k-1)*(hbs(i,k-1) - hb(i,k-1))) + temp3
if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
betamx(i) = ssfac*(vtemp1(i)/vtemp2(i))
beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
end if
else
beta(i) = 0.0_r8
end if
end do
!
! Third check to avoid introducing 2 delta x thermodynamic
! noise in the vertical ... constrain adjusted h (or theta e)
! so that the adjustment doesn't contribute to "kinks" in h
!
!cdir nodep
!DIR$ CONCURRENT
do ii=1,len1
i = indx1(ii)
g = min(0.0_r8,hb(i,k) - hb(i,k-1))
temp1 = (hb(i,k) - hb(i,k-1) - g)*(cats/ztodt)/(hc(i) - hbs(i,k))
vtemp1(i) = temp1*vtemp1(i) + (hc(i) - hbh(i,k+1))*rpdel(i,k)
vtemp2(i) = temp1*vtemp3(i)*rpdel(i,k) + (hc(i) - hbh(i,k) - cldwtr(i))* &
(rpdel(i,k) + rpdel(i,k+1))
if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
if (vtemp2(i) /= 0.0_r8) then
betamx(i) = vtemp1(i)/vtemp2(i)
else
betamx(i) = 0.0_r8
end if
beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
end if
end do
end if
!
! Calculate mass flux required for stabilization.
!
! Ensure that the convective mass flux, eta, is positive by
! setting negative values of eta to zero..
! Ensure that estimated mass flux cannot move more than the
! minimum of total mass contained in either layer k or layer k+1.
! Also test for other pathological cases that result in non-
! physical states and adjust eta accordingly.
!
!cdir nodep
!DIR$ CONCURRENT
do ii=1,len1
i = indx1(ii)
beta(i) = max(0.0_r8,beta(i))
temp1 = hc(i) - hbs(i,k)
temp2 = ((1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i)) - &
beta(i)*vtemp3(i))*rpdel(i,k) - (hbh(i,k+1) - hc(i))*rpdel(i,k+1)
eta(i) = temp1/(temp2*grav*cats)
tmass = min(pdel(i,k),pdel(i,k+1))*rgrav
if (eta(i) > tmass*rtdt .or. eta(i) <= 0.0_r8) eta(i) = 0.0_r8
!
! Check on negative q in top layer (bound beta)
!
if (shc(i)-shbh(i,k) < 0.0_r8 .and. beta(i)*eta(i) /= 0.0_r8) then
denom = eta(i)*grav*ztodt*(shc(i) - shbh(i,k))*rpdel(i,k-1)
beta(i) = max(0.0_r8,min(-0.999_r8*shb(i,k-1)/denom,beta(i)))
end if
!
! Check on negative q in middle layer (zero eta)
!
qtest1 = shb(i,k) + eta(i)*grav*ztodt*((shc(i) - shbh(i,k+1)) - &
(1.0_r8 - beta(i))*cldwtr(i)*rhlat - beta(i)*(shc(i) - shbh(i,k)))* &
rpdel(i,k)
if (qtest1 <= 0.0_r8) eta(i) = 0.0_r8
!
! Check on negative q in lower layer (bound eta)
!
fac1 = -(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
qtest2 = shb(i,k+1) - eta(i)*grav*ztodt*fac1
if (qtest2 < 0.0_r8) then
eta(i) = 0.99_r8*shb(i,k+1)/(grav*ztodt*fac1)
end if
etagdt(i) = eta(i)*grav*ztodt
end do
!
#if ( defined DIAGNS )
!#######################################################################
!# #
do i=1,ncol
if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep >= nsloc)) then
write(iulog,8080) beta(iloc), eta(iloc)
end if
enddo
!# #
!#######################################################################
#endif
!
! Calculate cloud water, rain water, and thermodynamic changes
!
!cdir nodep
!DIR$ CONCURRENT
do 30 ii=1,len1
i = indx1(ii)
icwmr(i,k) = cldwtr(i)*rhlat
cldwtr(i) = etagdt(i)*cldwtr(i)*rhlat*rgrav
! JJH changes to facilitate export of cloud liquid water --------------------------------
totcond(i) = (1.0_r8 - beta(i))*cldwtr(i)
rnwtr(i) = min(totcond(i),c0*(dzcld(i)-dzmin)*cldwtr(i))
ds1(i) = etagdt(i)*(sbh(i,k+1) - sc(i))*rpdel(i,k+1)
dq1(i) = etagdt(i)*(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
ds2(i) = (etagdt(i)*(sc(i) - sbh(i,k+1)) + &
hlat*grav*cldwtr(i) - beta(i)*etagdt(i)*(sc(i) - sbh(i,k)))*rpdel(i,k)
! JJH change for export of cloud liquid water; must use total condensate
! since rainwater no longer represents total condensate
dq2(i) = (etagdt(i)*(shc(i) - shbh(i,k+1)) - grav*totcond(i) - beta(i)* &
etagdt(i)*(shc(i) - shbh(i,k)))*rpdel(i,k)
ds3(i) = beta(i)*(etagdt(i)*(sc(i) - sbh(i,k)) - hlat*grav*cldwtr(i))* &
rpdel(i,k-1)
dq3(i) = beta(i)*etagdt(i)*(shc(i) - shbh(i,k))*rpdel(i,k-1)
!
! Isolate convective fluxes for later diagnostics
!
fslkp = eta(i)*(sc(i) - sbh(i,k+1))
fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) - hlat*cldwtr(i)*rtdt)
fqlkp = eta(i)*(shc(i) - shbh(i,k+1))
fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
!
! Update thermodynamic profile (update sb, hb, & hbs later)
!
tb (i,k+1) = tb(i,k+1) + ds1(i)*rcp
tb (i,k ) = tb(i,k ) + ds2(i)*rcp
tb (i,k-1) = tb(i,k-1) + ds3(i)*rcp
shb(i,k+1) = shb(i,k+1) + dq1(i)
shb(i,k ) = shb(i,k ) + dq2(i)
shb(i,k-1) = shb(i,k-1) + dq3(i)
!
! ** Update diagnostic information for final budget **
! Tracking precipitation, temperature & specific humidity tendencies,
! rainout term, convective mass flux, convective liquid
! water static energy flux, and convective total water flux
! The variable afac makes the necessary adjustment to the
! diagnostic fluxes to account for adjustment time scale based on
! how relaxation time scale is to be applied (column vs. triplet)
!
prec(i) = prec(i) + (rnwtr(i)/rhoh2o)*adjfac
!
! The following variables have units of "units"/second
!
cmfdt (i,k+1) = cmfdt (i,k+1) + ds1(i)*rtdt*adjfac
cmfdt (i,k ) = cmfdt (i,k ) + ds2(i)*rtdt*adjfac
cmfdt (i,k-1) = cmfdt (i,k-1) + ds3(i)*rtdt*adjfac
cmfdq (i,k+1) = cmfdq (i,k+1) + dq1(i)*rtdt*adjfac
cmfdq (i,k ) = cmfdq (i,k ) + dq2(i)*rtdt*adjfac
cmfdq (i,k-1) = cmfdq (i,k-1) + dq3(i)*rtdt*adjfac
! JJH changes to export cloud liquid water --------------------------------
qc (i,k ) = (grav*(totcond(i)-rnwtr(i))*rpdel(i,k))*rtdt*adjfac
cmfdqr(i,k ) = cmfdqr(i,k ) + (grav*rnwtr(i)*rpdel(i,k))*rtdt*adjfac
cmfmc (i,k+1) = cmfmc (i,k+1) + eta(i)*adjfac
cmfmc (i,k ) = cmfmc (i,k ) + beta(i)*eta(i)*adjfac
!
! The following variables have units of w/m**2
!
cmfsl (i,k+1) = cmfsl (i,k+1) + fslkp*adjfac
cmfsl (i,k ) = cmfsl (i,k ) + fslkm*adjfac
cmflq (i,k+1) = cmflq (i,k+1) + hlat*fqlkp*adjfac
cmflq (i,k ) = cmflq (i,k ) + hlat*fqlkm*adjfac
30 continue
!
! Next, convectively modify passive constituents
! For now, when applying relaxation time scale to thermal fields after
! entire column has undergone convective overturning, constituents will
! be mixed using a "relaxed" value of the mass flux determined above
! Although this will be inconsistant with the treatment of the thermal
! fields, it's computationally much cheaper, no more-or-less justifiable,
! and consistent with how the history tape mass fluxes would be used in
! an off-line mode (i.e., using an off-line transport model)
!
do 50 m=2,pcnst ! note: indexing assumes water is first field
if (cnst_get_type_byind(m).eq.'dry') then
pd(:ncol,:) = pdeldry(:ncol,:)
rpd(:ncol,:) = rpdeldry(:ncol,:)
pm(:ncol,:) = pmiddry(:ncol,:)
else
pd(:ncol,:) = pdel(:ncol,:)
rpd(:ncol,:) = rpdel(:ncol,:)
pm(:ncol,:) = pmid(:ncol,:)
endif
!cdir nodep
!DIR$ CONCURRENT
do 40 ii=1,len1
i = indx1(ii)
!
! If any of the reported values of the constituent is negative in
! the three adjacent levels, nothing will be done to the profile
!
if ((dq(i,k+1,m) < 0.0_r8) .or. (dq(i,k,m) < 0.0_r8) .or. (dq(i,k-1,m) < 0.0_r8)) go to 40
!
! Specify constituent interface values (linear interpolation)
!
cmrh(i,k ) = 0.5_r8*(dq(i,k-1,m) + dq(i,k ,m))
cmrh(i,k+1) = 0.5_r8*(dq(i,k ,m) + dq(i,k+1,m))
!
! Specify perturbation properties of constituents in PBL
!
pblhgt = max(pblht(i),1.0_r8)
if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_r8 ) then
fac1 = max(0.0_r8,1.0_r8-zm(i,k+1)/pblhgt)
cmrc(i) = dq(i,k+1,m) + qpert(i,m)*fac1
else
cmrc(i) = dq(i,k+1,m)
end if
!
! Determine fluxes, flux divergence => changes due to convection
! Logic must be included to avoid producing negative values. A bit
! messy since there are no a priori assumptions about profiles.
! Tendency is modified (reduced) when pending disaster detected.
!
botflx = etagdt(i)*(cmrc(i) - cmrh(i,k+1))*adjfac
topflx = beta(i)*etagdt(i)*(cmrc(i)-cmrh(i,k))*adjfac
dcmr1(i) = -botflx*rpd(i,k+1)
efac1 = 1.0_r8
efac2 = 1.0_r8
efac3 = 1.0_r8
!
if (dq(i,k+1,m)+dcmr1(i) < 0.0_r8) then
if ( abs(dcmr1(i)) > 1.e-300_r8 ) then
efac1 = max(tiny,abs(dq(i,k+1,m)/dcmr1(i)) - eps)
else
efac1 = tiny
endif
end if
!
if (efac1 == tiny .or. efac1 > 1.0_r8) efac1 = 0.0_r8
dcmr1(i) = -efac1*botflx*rpd(i,k+1)
dcmr2(i) = (efac1*botflx - topflx)*rpd(i,k)
!
if (dq(i,k,m)+dcmr2(i) < 0.0_r8) then
if ( abs(dcmr2(i)) > 1.e-300_r8 ) then
efac2 = max(tiny,abs(dq(i,k ,m)/dcmr2(i)) - eps)
else
efac2 = tiny
endif
end if
!
if (efac2 == tiny .or. efac2 > 1.0_r8) efac2 = 0.0_r8
dcmr2(i) = (efac1*botflx - efac2*topflx)*rpd(i,k)
dcmr3(i) = efac2*topflx*rpd(i,k-1)
!
if ( (dq(i,k-1,m)+dcmr3(i) < 0.0_r8 ) ) then
if ( abs(dcmr3(i)) > 1.e-300_r8 ) then
efac3 = max(tiny,abs(dq(i,k-1,m)/dcmr3(i)) - eps)
else
efac3 = tiny
endif
end if
!
if (efac3 == tiny .or. efac3 > 1.0_r8) efac3 = 0.0_r8
efac3 = min(efac2,efac3)
dcmr2(i) = (efac1*botflx - efac3*topflx)*rpd(i,k)
dcmr3(i) = efac3*topflx*rpd(i,k-1)
!
dq(i,k+1,m) = dq(i,k+1,m) + dcmr1(i)
dq(i,k ,m) = dq(i,k ,m) + dcmr2(i)
dq(i,k-1,m) = dq(i,k-1,m) + dcmr3(i)
40 continue
50 continue ! end of m=2,pcnst loop
!
! Constituent modifications complete
!
if (k == limcnv+1) go to 60
!
! Complete update of thermodynamic structure at integer levels
! gather/scatter points that need new values of shbs and gamma
!
do ii=1,len1
i = indx1(ii)
vtemp1(ii ) = tb(i,k)
vtemp1(ii+len1) = tb(i,k-1)
vtemp2(ii ) = pmid(i,k)
vtemp2(ii+len1) = pmid(i,k-1)
end do
call vqsatd (vtemp1 ,vtemp2 ,estemp ,vtemp3 , vtemp4 , &
2*len1 ) ! using estemp as extra long vector
!cdir nodep
!DIR$ CONCURRENT
do ii=1,len1
i = indx1(ii)
shbs(i,k ) = vtemp3(ii )
shbs(i,k-1) = vtemp3(ii+len1)
gam(i,k ) = vtemp4(ii )
gam(i,k-1) = vtemp4(ii+len1)
sb (i,k ) = sb(i,k ) + ds2(i)
sb (i,k-1) = sb(i,k-1) + ds3(i)
hb (i,k ) = sb(i,k ) + hlat*shb(i,k )
hb (i,k-1) = sb(i,k-1) + hlat*shb(i,k-1)
hbs(i,k ) = sb(i,k ) + hlat*shbs(i,k )
hbs(i,k-1) = sb(i,k-1) + hlat*shbs(i,k-1)
end do
!
! Update thermodynamic information at half (i.e., interface) levels
!
!DIR$ CONCURRENT
do ii=1,len1
i = indx1(ii)
sbh (i,k) = 0.5_r8*(sb(i,k) + sb(i,k-1))
shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
sbh (i,k-1) = 0.5_r8*(sb(i,k-1) + sb(i,k-2))
shbh(i,k-1) = qhalf(shb(i,k-2),shb(i,k-1),shbs(i,k-2),shbs(i,k-1))
hbh (i,k-1) = sbh(i,k-1) + hlat*shbh(i,k-1)
end do
!
#if ( defined DIAGNS )
!#######################################################################
!# #
!# this update necessary, only if debugging diagnostics requested #
!# #
do i=1,ncol
if (lat(i) == jloc .and. nstep >= nsloc) then
call vqsatd(tb(i,k+1),pmid(i,k+1),es,shbs(i,k+1),gam(i,k+1), &
1)
sb (i,k+1) = sb(i,k+1) + ds1(i)
hb (i,k+1) = sb(i,k+1) + hlat*shb(i,k+1)
hbs(i,k+1) = sb(i,k+1) + hlat*shbs(i,k+1)
kpp = k + 2
if (k+1 == pver) kpp = k + 1
do kp=k+1,kpp
kpm1 = kp-1
sbh(i,kp) = 0.5_r8*(sb(i,kpm1) + sb(i,kp))
shbh(i,kp) = qhalf(shb(i,kpm1),shb(i,kp),shbs(i,kpm1),shbs(i,kp))
hbh(i,kp) = sbh(i,kp) + hlat*shbh(i,kp)
end do
end if
end do
!# #
!# diagnostic output #
!# #
do i=1,ncol
if ((lat(i)== jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
write(iulog, 8060) k
fac = grav*864._r8
do n=limcnv,pver
rh = shb(i,n)/shbs(i,n)
write(iulog,8020)shbh(i,n),sbh(i,n),hbh(i,n),fac*cmfmc(i,n), &
cmfsl(i,n), cmflq(i,n)
!--------------write(iulog, 8050)
!--------------write(iulog, 8030) fac*cmfmc(i,n),cmfsl(i,n), cmflq(i,n)
write(iulog, 8040) tb(i,n),shb(i,n),rh,sb(i,n),hb(i,n), &
hbs(i,n), ztodt*cmfdt(i,n),ztodt*cmfdq(i,n), &
ztodt*cmfdqr(i,n)
end do
write(iulog, 8000) prec(i)
end if
end do
!# #
!# #
!#######################################################################
#endif
!
! Ensure that dzcld is reset if convective mass flux zero
! specify the current vertical extent of the convective activity
! top of convective layer determined by size of overshoot param.
!
60 do i=1,ncol
etagt0 = eta(i).gt.0.0_r8
if ( .not. etagt0) dzcld(i) = 0.0_r8
if (etagt0 .and. beta(i) > betamn) then
ktp = k-1
else
ktp = k
end if
if (etagt0) then
rk=k
rktp=ktp
cnt(i) = min(cnt(i),rktp)
cnb(i) = max(cnb(i),rk)
end if
end do
70 continue ! end of k loop
!
! ** apply final thermodynamic tendencies **
!
!**BAB don't update input profiles
!!$ do k=limcnv,pver
!!$ do i=1,ncol
!!$ t (i,k) = t (i,k) + cmfdt(i,k)*ztodt
!!$ q(i,k,1) = q(i,k,1) + cmfdq(i,k)*ztodt
!!$ end do
!!$ end do
! Set output q tendencies
dq(:ncol,:,1 ) = cmfdq(:ncol,:)
dq(:ncol,:,2:) = (dq(:ncol,:,2:) - q(:ncol,:,2:))/ztodt
!
! Kludge to prevent cnb-cnt from being zero (in the event
! someone decides that they want to divide by this quantity)
!
do i=1,ncol
if (cnb(i) /= 0.0_r8 .and. cnb(i) == cnt(i)) then
cnt(i) = cnt(i) - 1.0_r8
end if
end do
!
do i=1,ncol
precc(i) = prec(i)*rtdt
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) + qc(i,k)*pdel(i,k)/grav
end do
end do
rliq(:ncol) = rliq(:ncol) /1000._r8
#if ( defined DIAGNS )
!#######################################################################
!# #
!# we're done ... show final result if debug diagnostics requested #
!# #
do i=1,ncol
if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
fac = grav*864._r8
write(iulog, 8010)
do k=limcnv,pver
rh = shb(i,k)/shbs(i,k)
write(iulog, 8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k), &
cmfsl(i,k), cmflq(i,k)
write(iulog, 8040) tb(i,k),shb(i,k),rh ,sb(i,k),hb(i,k), &
hbs(i,k), ztodt*cmfdt(i,k),ztodt*cmfdq(i,k), &
ztodt*cmfdqr(i,k)
end do
write(iulog, 8000) prec(i)
!# #
!# approximate vertical integral of moist static energy and #
!# total preciptable water after adjustment and output changes #
!# #
hsum2 = 0.0_r8
qsum2 = 0.0_r8
do k=limcnv,pver
hsum2 = hsum2 + pdel(i,k)*rgrav*hb(i,k)
qsum2 = qsum2 + pdel(i,k)*rgrav*shb(i,k)
end do
!# #
write(iulog,8070) hsum1, hsum2, abs(hsum2-hsum1)/hsum2, &
qsum1, qsum2, abs(qsum2-qsum1)/qsum2
end if
enddo
!# #
!# #
!#######################################################################
#endif
return ! we're all done ... return to calling procedure
#if ( defined DIAGNS )
!
! Formats
!
8000 format(///,10x,'PREC = ',3pf12.6,/)
8010 format('1** TB SHB RH SB', &
' HB HBS CAH CAM PRECC ', &
' ETA FSL FLQ **', /)
8020 format(' ----- ', 9x,3p,f7.3,2x,2p, 9x,-3p,f7.3,2x, &
f7.3, 37x, 0p,2x,f8.2, 0p,2x,f8.2,2x,f8.2, ' ----- ')
8030 format(' ----- ', 0p,82x,f8.2, 0p,2x,f8.2,2x,f8.2, &
' ----- ')
8040 format(' - - - ',f7.3,2x,3p,f7.3,2x,2p,f7.3,2x,-3p,f7.3,2x, &
f7.3, 2x,f8.3,2x,0p,f7.3,3p,2x,f7.3,2x,f7.3,30x, &
' - - - ')
8050 format(' ----- ',110x,' ----- ')
8060 format('1 K =>', i4,/, &
' TB SHB RH SB', &
' HB HBS CAH CAM PREC ', &
' ETA FSL FLQ', /)
8070 format(' VERTICALLY INTEGRATED MOIST STATIC ENERGY BEFORE, AFTER', &
' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/, &
' VERTICALLY INTEGRATED MOISTURE BEFORE, AFTER' &,
' AND PERCENTAGE DIFFERENCE => ',1p,2e15_r8.7_r8,2x,2p,f7.3,/)
8080 format(' BETA, ETA => ', 1p,2e12.3)
8090 format (' k+1, sc, shc, hc => ', 1x, i2, 1p, 3e12.4)
#endif
!
end subroutine cmfmca
end module hk_conv