#define WRF_PORT
#define MODAL_AERO
! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
module module_cu_camzm 2
!---------------------------------------------------------------------------------
! Purpose:
!
! Interface from Zhang-McFarlane convection scheme, includes evaporation of convective
! precip from the ZM scheme
!
! Apr 2006: RBN: Code added to perform a dilute ascent for closure of the CM mass flux
! based on an entraining plume a la Raymond and Blythe (1992)
!
! Author: Byron Boville, from code in tphysbc
!
!---------------------------------------------------------------------------------
use shr_kind_mod
, only: r8 => shr_kind_r8
#ifndef WRF_PORT
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver, pverp
#else
use module_cam_support, only: masterproc, pcols, pver, pverp
#endif
use cldwat, only: cldwat_fice
use physconst, only: cpair, epsilo, gravit, latice, latvap, tmelt, rair, &
cpwv, cpliq, rh2o
#ifndef WRF_PORT
use abortutils, only: endrun
use cam_logfile, only: iulog
#else
use module_cam_support,only: endrun, iulog
#endif
implicit none
save
private ! Make default type private to the module
!
! PUBLIC: interfaces
!
public zmconv_readnl ! read zmconv_nl namelist
public zm_convi ! ZM schemea
public zm_convr ! ZM schemea
public zm_conv_evap ! evaporation of precip from ZM schemea
public convtran ! convective transport
public momtran ! convective momentum transport
!
! Private data
!
real(r8), parameter :: unset_r8 = huge(1.0_r8)
real(r8) :: zmconv_c0_lnd = unset_r8
real(r8) :: zmconv_c0_ocn = unset_r8
real(r8) :: zmconv_ke = unset_r8
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._r8 ! threshold value for cape for deep convection.
real(r8) :: ke ! Tunable evaporation efficiency set from namelist input zmconv_ke
real(r8) :: c0_lnd ! set from namelist input zmconv_c0_lnd
real(r8) :: c0_ocn ! set from namelist input zmconv_c0_ocn
real(r8) tau ! convective time scale
real(r8),parameter :: a = 21.656_r8
real(r8),parameter :: b = 5418._r8
real(r8),parameter :: c1 = 6.112_r8
real(r8),parameter :: c2 = 17.67_r8
real(r8),parameter :: c3 = 243.5_r8
real(r8) :: tfreez
real(r8) :: eps1
logical :: no_deep_pbl ! default = .false.
! no_deep_pbl = .true. eliminates deep convection entirely within PBL
!moved from moistconvection.F90
real(r8) :: rgrav ! reciprocal of grav
real(r8) :: rgas ! gas constant for dry air
real(r8) :: grav ! = gravit
real(r8) :: cp ! = cpres = cpair
integer limcnv ! top interface level limit for convection
real(r8),parameter :: tiedke_add = 0.5_r8
contains
subroutine zmconv_readnl(nlfile) 1,1
#ifndef WRF_PORT
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
#endif
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
#ifndef WRF_PORT
! Local variables
integer :: unitn, ierr
character(len=*), parameter :: subname = 'zmconv_readnl'
namelist /zmconv_nl/ zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'zmconv_nl', status=ierr)
if (ierr == 0) then
read(unitn, zmconv_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
c0_lnd = zmconv_c0_lnd
c0_ocn = zmconv_c0_ocn
ke = zmconv_ke
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast(c0_lnd, 1, mpir8, 0, mpicom)
call mpibcast(c0_ocn, 1, mpir8, 0, mpicom)
call mpibcast(ke, 1, mpir8, 0, mpicom)
#endif
#else
! WRF_PORT currently uses hard-wired values for the namelist input. The
! values could easily be setup to come from the Registry in the future.
! The hard-wired values are the defaults for the fv core. They should be
! verified by somebody knowledgable on the matter.
c0_lnd = 0.0059D0
c0_ocn = 0.0450D0
ke = 1.0E-6
#endif
end subroutine zmconv_readnl
#ifndef WRF_PORT
subroutine zm_convi(limcnv_in, no_deep_pbl_in) 1,7
use dycore, only: dycore_is, get_resolution
#else
subroutine zm_convi(DT,DX,limcnv_in, no_deep_pbl_in) 1,7
#endif
integer, intent(in) :: limcnv_in ! top interface level limit for convection
logical, intent(in), optional :: no_deep_pbl_in ! no_deep_pbl = .true. eliminates ZM convection entirely within PBL
#ifdef WRF_PORT
REAL , INTENT(IN) :: DT, DX
#endif
real(r8) :: DTT
! local variables
character(len=32) :: hgrid ! horizontal grid specifier
#ifdef WRF_PORT
real(r8) :: deltax ! model delta x
real(r8) :: ref_dx ! reference delta x (CAM5 1.9 x 2.5 deg)
real(r8) :: taumin ! minimum allowable tau
real(r8) :: taumax ! maximum allowable tau
#endif
! Initialization of ZM constants
limcnv = limcnv_in
tfreez = tmelt
eps1 = epsilo
rl = latvap
cpres = cpair
rgrav = 1.0_r8/gravit
rgas = rair
grav = gravit
cp = cpres
#ifdef WRF_PORT
!PMA sets parameters
ref_dx = 275000._r8
taumin = 600._r8
taumax = 3600._r8
deltax=dx
#endif
if ( present(no_deep_pbl_in) ) then
no_deep_pbl = no_deep_pbl_in
else
no_deep_pbl = .true.
endif
DTT=DT
! tau=4800. were used in canadian climate center. however, in echam3 t42,
! convection is too weak, thus adjusted to 2400.
#ifndef WRF_PORT
hgrid = get_resolution()
#endif
!PMA: The 1.9x2.5 deg CAM5 uses tau = 3600s. With higher resolution and shorter timestep in CAM,
! tau is set to a smaller value. Williamson (2012) found that without reducing tau in a high-res
! simulation, 'grid-scale storms' can form because the ZM scheme fails to remove CAPE quickly enough
! so that the resolved scheme (micro+macro) has to work very hard. WRF usually is run at much higher
! resolution than CAM and this grid-scale storm problem is apparent, which can cause the model to
! crash. Hence, we make tau to be a function of spatial resolution.
! For WRF, one needs to make dt very small if tau=3600s for the model to run stably.
! tau = 3600._r8
!PMA formula: tau = 3600s when dx = 275km ~ 2.5deg
! tau = 600s when dx = 10km
! tau = 11.32_r8 * deltax + 489.68_r8
!PMA adopts PJR formula
tau = max(taumin, taumax*min(1._r8,deltax/ref_dx))
if ( masterproc ) then
write(iulog,*) 'delta X =',deltax
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,*) 'Convective relaxation time scale (tau) is a tunable parameter in CAM and is a function of spatial resolution.'
write(iulog,*) 'Users are encouraged to consult with the PNNL WRF-CAM5 development team if they want to change tau.'
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,*) 'tuning parameters zm_convi: tau',tau
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,*) 'tuning parameters zm_convi: c0_lnd',c0_lnd, ', c0_ocn', c0_ocn
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,*) 'tuning parameters zm_convi: ke',ke
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,*) 'tuning parameters zm_convi: no_deep_pbl',no_deep_pbl
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
endif
if (masterproc) write(iulog,*)'**** ZM: DILUTE Buoyancy Calculation ****'
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
end subroutine zm_convi
subroutine zm_convr(lchnk ,ncol , & 1,7
t ,qh ,prec ,jctop ,jcbot , &
pblh ,zm ,geos ,zi ,qtnd , &
heat ,pap ,paph ,dpp , &
delt ,mcon ,cme ,cape , &
tpert ,dlf ,pflx ,zdu ,rprd , &
mu ,md ,du ,eu ,ed , &
dp ,dsubcld ,jt ,maxg ,ideep , &
lengath ,ql ,rliq ,landfrac)
!-----------------------------------------------------------------------
!
! 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
!
!-----------------------------------------------------------------------
#ifdef WRF_PORT
use module_cam_support, only: pcnst =>pcnst_runtime
#else
use constituents, only: pcnst
use phys_control, only: cam_physpkg_is
#endif
!
! ************************ 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) ! 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)
real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac
!
! 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) :: cape(pcols) ! w convective available potential energy.
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) 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.
real(r8) landfracg(pcols) ! wg grid slice of landfrac
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._r8
heat(:,:) = 0._r8
mcon(:,:) = 0._r8
rliq(:ncol) = 0._r8
!
! initialize convective tendencies
!
prec(:ncol) = 0._r8
do k = 1,pver
do i = 1,ncol
dqdt(i,k) = 0._r8
dsdt(i,k) = 0._r8
dudt(i,k) = 0._r8
dvdt(i,k) = 0._r8
pflx(i,k) = 0._r8
pflxg(i,k) = 0._r8
cme(i,k) = 0._r8
rprd(i,k) = 0._r8
zdu(i,k) = 0._r8
ql(i,k) = 0._r8
qlg(i,k) = 0._r8
dlf(i,k) = 0._r8
dlg(i,k) = 0._r8
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._r8
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_r8
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_r8
pf(i,k) = paph(i,k)*0.01_r8
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_r8) 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_r8
shat(i,k) = s(i,k)
qhat(i,k) = q(i,k)
end do
end do
do i = 1,ncol
capeg(i) = 0._r8
lclg(i) = 1
lelg(i) = pver
maxg(i) = 1
tlg(i) = 400._r8
dsubcld(i) = 0._r8
end do
#ifndef WRF_PORT
if( cam_physpkg_is('cam3')) then
! For cam3 physics package, call non-dilute
call buoyan(lchnk ,ncol , &
q ,t ,p ,z ,pf , &
tp ,qstp ,tl ,rl ,cape , &
pblt ,lcl ,lel ,lon ,maxi , &
rgas ,grav ,cpres ,msg , &
tpert )
else
#else
! Evaluate Tparcel, qsat(Tparcel), buoyancy and CAPE,
! lcl, lel, parcel launch level at index maxi()=hmax
call buoyan_dilute(lchnk ,ncol , &
q ,t ,p ,z ,pf , &
tp ,qstp ,tl ,rl ,cape , &
pblt ,lcl ,lel ,lon ,maxi , &
rgas ,grav ,cpres ,msg , &
tpert )
#endif
#ifndef WRF_PORT
end if
#endif
!
! 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_r8*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._r8
vg(i,k) = 0._r8
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))
landfracg(i) = landfrac(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._r8
qdifr = 0._r8
if (sg(i,k) > 0._r8 .or. sg(i,k-1) > 0._r8) &
sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
if (qg(i,k) > 0._r8 .or. qg(i,k-1) > 0._r8) &
qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
if (sdifr > 1.E-6_r8) 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_r8* (sg(i,k)+sg(i,k-1))
end if
if (qdifr > 1.E-6_r8) 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_r8* (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 ,landfracg)
!
! 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._r8) then
mb(i) = min(mb(i),0.5_r8/(delt*mumax(i)))
else
mb(i) = 0._r8
endif
end do
! If no_deep_pbl = .true., don't allow convection entirely
! within PBL (suggestion of Bjorn Stevens, 8-2000)
if (no_deep_pbl) then
do i=1,lengath
if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0
end do
end if
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._r8/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
!DIR$ CONCURRENT
do i = 1,lengath
!
! q is updated to compute net precip.
!
q(ideep(i),k) = qh(ideep(i),k,1) + 2._r8*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
!
!DIR$ CONCURRENT
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._r8*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._r8*delt)/1000._r8
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._r8
return
end subroutine zm_convr
!===============================================================================
subroutine zm_conv_evap(ncol,lchnk, & 1,3
t,pmid,pdel,q, &
tend_s, tend_s_snwprd, tend_s_snwevmlt, tend_q, &
prdprec, cldfrc, deltat, &
prec, snow, ntprprd, ntsnprd, flxprec, flxsnow )
!-----------------------------------------------------------------------
! 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 wv_saturation, only: aqsat
#ifndef WRF_PORT
use phys_grid, only: get_rlat_all_p
#endif
!------------------------------Arguments--------------------------------
integer,intent(in) :: ncol, lchnk ! number of columns and chunk index
real(r8),intent(in), dimension(pcols,pver) :: t ! temperature (K)
real(r8),intent(in), dimension(pcols,pver) :: pmid ! midpoint pressure (Pa)
real(r8),intent(in), dimension(pcols,pver) :: pdel ! layer thickness (Pa)
real(r8),intent(in), dimension(pcols,pver) :: q ! water vapor (kg/kg)
real(r8),intent(inout), dimension(pcols,pver) :: tend_s ! heating rate (J/kg/s)
real(r8),intent(inout), dimension(pcols,pver) :: tend_q ! water vapor tendency (kg/kg/s)
real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwprd ! Heating rate of snow production
real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwevmlt ! Heating rate of evap/melting of snow
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
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),intent(out) :: flxprec(pcols,pverp) ! Convective-scale flux of precip at interfaces (kg/m2/s)
real(r8),intent(out) :: flxsnow(pcols,pverp) ! Convective-scale flux of snow at interfaces (kg/m2/s)
real(r8),intent(out) :: ntprprd(pcols,pver) ! net precip production in layer
real(r8),intent(out) :: 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
!-----------------------------------------------------------------------
! convert input precip to kg/m2/s
prec(:ncol) = prec(:ncol)*1000._r8
! determine saturation vapor pressure
call aqsat (t ,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, t, fice, fsnow_conv)
! zero the flux integrals on the top boundary
flxprec(:ncol,1) = 0._r8
flxsnow(:ncol,1) = 0._r8
evpvint(:ncol) = 0._r8
do k = 1, pver
do i = 1, ncol
! Melt snow falling into layer, if necessary.
if (t(i,k) > tmelt) then
flxsntm(i) = 0._r8
snowmlt(i) = flxsnow(i,k) * gravit/ pdel(i,k)
else
flxsntm(i) = flxsnow(i,k)
snowmlt(i) = 0._r8
end if
! relative humidity depression must be > 0 for evaporation
evplimit = max(1._r8 - q(i,k)/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._r8 - 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)-q(i,k)) / 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 / pdel(i,k) + min(prdprec(i,k), 0.)
evplimit = min(evplimit, flxprec(i,k) * gravit / pdel(i,k))
! Total evaporation cannot exceed input precipitation
evplimit = min(evplimit, (prec(i) - evpvint(i)) * gravit / 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._r8) 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._r8
end if
! vertically integrated evaporation
evpvint(i) = evpvint(i) + evpprec(i) * 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._r8) then
work1 = min(max(0._r8,flxsnow(i,k)/flxprec(i,k)),1._r8)
else
work1 = 0._r8
endif
#endif
work2 = max(fsnow_conv(i,k), work1)
if (snowmlt(i).gt.0._r8) work2 = 0._r8
! work2 = fsnow_conv(i,k)
ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i)
tend_s_snwprd (i,k) = prdprec(i,k)*work2*latice
tend_s_snwevmlt(i,k) = - ( evpsnow(i) + snowmlt(i) )*latice
! precipitation fluxes
flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * pdel(i,k)/gravit
flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * 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
tend_s(i,k) =-evpprec(i)*latvap + ntsnprd(i,k)*latice
tend_q(i,k) = evpprec(i)
end do
end do
! set output precipitation rates (m/s)
prec(:ncol) = flxprec(:ncol,pver+1) / 1000._r8
snow(:ncol) = flxsnow(:ncol,pver+1) / 1000._r8
!**********************************************************
!!$ tend_s(:ncol,:) = 0. ! turn heating off
!**********************************************************
end subroutine zm_conv_evap
subroutine convtran(lchnk , & 2,2
doconvtran,q ,ncnst ,mu ,md , &
du ,eu ,ed ,dp ,dsubcld , &
jt ,mx ,ideep ,il1g ,il2g , &
nstep ,fracis ,dqdt ,dpdry )
!-----------------------------------------------------------------------
!
! Purpose:
! Convective transport of trace species
!
! Mixing ratios may be with respect to either dry or moist air
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: P. Rasch
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use constituents, only: cnst_get_type_byind
#ifndef WRF_PORT
use ppgrid
use abortutils, only: endrun
#endif
implicit none
!-----------------------------------------------------------------------
!
! Input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncnst ! number of tracers to transport
logical, intent(in) :: doconvtran(ncnst) ! flag for doing convective transport
real(r8), intent(in) :: q(pcols,pver,ncnst) ! Tracer array including moisture
real(r8), intent(in) :: mu(pcols,pver) ! Mass flux up
real(r8), intent(in) :: md(pcols,pver) ! Mass flux down
real(r8), intent(in) :: du(pcols,pver) ! Mass detraining from updraft
real(r8), intent(in) :: eu(pcols,pver) ! Mass entraining from updraft
real(r8), intent(in) :: ed(pcols,pver) ! Mass entraining from downdraft
real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces
real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc
real(r8), intent(in) :: fracis(pcols,pver,ncnst) ! fraction of tracer that is insoluble
integer, intent(in) :: jt(pcols) ! Index of cloud top for each column
integer, intent(in) :: mx(pcols) ! Index of cloud top for each column
integer, intent(in) :: ideep(pcols) ! Gathering array
integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
integer, intent(in) :: nstep ! Time step index
real(r8), intent(in) :: dpdry(pcols,pver) ! Delta pressure between interfaces
! input/output
real(r8), intent(out) :: dqdt(pcols,pver,ncnst) ! Tracer tendency array
!--------------------------Local Variables------------------------------
integer i ! Work index
integer k ! Work index
integer kbm ! Highest altitude index of cloud base
integer kk ! Work index
integer kkp1 ! Work index
integer km1 ! Work index
integer kp1 ! Work index
integer ktm ! Highest altitude index of cloud top
integer m ! Work index
real(r8) cabv ! Mix ratio of constituent above
real(r8) cbel ! Mix ratio of constituent below
real(r8) cdifr ! Normalized diff between cabv and cbel
real(r8) chat(pcols,pver) ! Mix ratio in env at interfaces
real(r8) cond(pcols,pver) ! Mix ratio in downdraft at interfaces
real(r8) const(pcols,pver) ! Gathered tracer array
real(r8) fisg(pcols,pver) ! gathered insoluble fraction of tracer
real(r8) conu(pcols,pver) ! Mix ratio in updraft at interfaces
real(r8) dcondt(pcols,pver) ! Gathered tend array
real(r8) small ! A small number
real(r8) mbsth ! Threshold for mass fluxes
real(r8) mupdudp ! A work variable
real(r8) minc ! A work variable
real(r8) maxc ! A work variable
real(r8) fluxin ! A work variable
real(r8) fluxout ! A work variable
real(r8) netflux ! A work variable
real(r8) dutmp(pcols,pver) ! Mass detraining from updraft
real(r8) eutmp(pcols,pver) ! Mass entraining from updraft
real(r8) edtmp(pcols,pver) ! Mass entraining from downdraft
real(r8) dptmp(pcols,pver) ! Delta pressure between interfaces
!-----------------------------------------------------------------------
!
small = 1.e-36_r8
! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
mbsth = 1.e-15_r8
! 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
! Loop ever each constituent
do m = 2, ncnst
if (doconvtran(m)) then
if (cnst_get_type_byind(m).eq.'dry') then
do k = 1,pver
do i =il1g,il2g
dptmp(i,k) = dpdry(i,k)
dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
end do
end do
else
do k = 1,pver
do i =il1g,il2g
dptmp(i,k) = dp(i,k)
dutmp(i,k) = du(i,k)
eutmp(i,k) = eu(i,k)
edtmp(i,k) = ed(i,k)
end do
end do
endif
! dptmp = dp
! Gather up the constituent and set tend to zero
do k = 1,pver
do i =il1g,il2g
const(i,k) = q(ideep(i),k,m)
fisg(i,k) = fracis(ideep(i),k,m)
end do
end do
! From now on work only with gathered data
! Interpolate environment tracer values to interfaces
do k = 1,pver
km1 = max(1,k-1)
do i = il1g, il2g
minc = min(const(i,km1),const(i,k))
maxc = max(const(i,km1),const(i,k))
if (minc < 0) then
cdifr = 0._r8
else
cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small)
endif
! If the two layers differ significantly use a geometric averaging
! procedure
if (cdifr > 1.E-6_r8) then
cabv = max(const(i,km1),maxc*1.e-12_r8)
cbel = max(const(i,k),maxc*1.e-12_r8)
chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
else ! Small diff, so just arithmetic mean
chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1))
end if
! Provisional up and down draft values
conu(i,k) = chat(i,k)
cond(i,k) = chat(i,k)
! provisional tends
dcondt(i,k) = 0._r8
end do
end do
! Do levels adjacent to top and bottom
k = 2
km1 = 1
kk = pver
do i = il1g,il2g
mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
if (mupdudp > mbsth) then
conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp
endif
if (md(i,k) < -mbsth) then
cond(i,k) = (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k)
endif
end do
! Updraft from bottom to top
do kk = pver-1,1,-1
kkp1 = min(pver,kk+1)
do i = il1g,il2g
mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
if (mupdudp > mbsth) then
conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* &
const(i,kk)*dptmp(i,kk) )/mupdudp
endif
end do
end do
! Downdraft from top to bottom
do k = 3,pver
km1 = max(1,k-1)
do i = il1g,il2g
if (md(i,k) < -mbsth) then
cond(i,k) = ( md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) &
*dptmp(i,km1) )/md(i,k)
endif
end do
end do
do k = ktm,pver
km1 = max(1,k-1)
kp1 = min(pver,k+1)
do i = il1g,il2g
! version 1 hard to check for roundoff errors
! dcondt(i,k) =
! $ +(+mu(i,kp1)* (conu(i,kp1)-chat(i,kp1))
! $ -mu(i,k)* (conu(i,k)-chat(i,k))
! $ +md(i,kp1)* (cond(i,kp1)-chat(i,kp1))
! $ -md(i,k)* (cond(i,k)-chat(i,k))
! $ )/dp(i,k)
! version 2 hard to limit fluxes
! fluxin = mu(i,kp1)*conu(i,kp1) + mu(i,k)*chat(i,k)
! $ -(md(i,k) *cond(i,k) + md(i,kp1)*chat(i,kp1))
! fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*chat(i,kp1)
! $ -(md(i,kp1)*cond(i,kp1) + md(i,k)*chat(i,k))
! version 3 limit fluxes outside convection to mass in appropriate layer
! these limiters are probably only safe for positive definite quantitities
! it assumes that mu and md already satify a courant number limit of 1
fluxin = mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) &
-(md(i,k) *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1)))
fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) &
-(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k)))
netflux = fluxin - fluxout
if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then
netflux = 0._r8
endif
dcondt(i,k) = netflux/dptmp(i,k)
end do
end do
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!DIR$ NOINTERCHANGE
do k = kbm,pver
km1 = max(1,k-1)
do i = il1g,il2g
if (k == mx(i)) then
! version 1
! dcondt(i,k) = (1./dsubcld(i))*
! $ (-mu(i,k)*(conu(i,k)-chat(i,k))
! $ -md(i,k)*(cond(i,k)-chat(i,k))
! $ )
! version 2
! fluxin = mu(i,k)*chat(i,k) - md(i,k)*cond(i,k)
! fluxout = mu(i,k)*conu(i,k) - md(i,k)*chat(i,k)
! version 3
fluxin = mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k)
fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k))
netflux = fluxin - fluxout
if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then
netflux = 0._r8
endif
! dcondt(i,k) = netflux/dsubcld(i)
dcondt(i,k) = netflux/dptmp(i,k)
else if (k > mx(i)) then
! dcondt(i,k) = dcondt(i,k-1)
dcondt(i,k) = 0._r8
end if
end do
end do
! Initialize to zero everywhere, then scatter tendency back to full array
dqdt(:,:,m) = 0._r8
do k = 1,pver
kp1 = min(pver,k+1)
!DIR$ CONCURRENT
do i = il1g,il2g
dqdt(ideep(i),k,m) = dcondt(i,k)
end do
end do
end if ! for doconvtran
end do
return
end subroutine convtran
!=========================================================================================
subroutine momtran(lchnk, ncol, & 1,2
domomtran,q ,ncnst ,mu ,md , &
du ,eu ,ed ,dp ,dsubcld , &
jt ,mx ,ideep ,il1g ,il2g , &
nstep ,dqdt ,pguall ,pgdall, icwu, icwd, dt, seten )
!-----------------------------------------------------------------------
!
! Purpose:
! Convective transport of momentum
!
! Mixing ratios may be with respect to either dry or moist air
!
! Method:
! Based on the convtran subroutine by P. Rasch
! <Also include any applicable external references.>
!
! Author: J. Richter and P. Rasch
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
#ifndef WRF_PORT
use constituents, only: cnst_get_type_byind
use ppgrid
use abortutils, only: endrun
#endif
implicit none
!-----------------------------------------------------------------------
!
! Input arguments
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: ncnst ! number of tracers to transport
logical, intent(in) :: domomtran(ncnst) ! flag for doing convective transport
real(r8), intent(in) :: q(pcols,pver,ncnst) ! Wind array
real(r8), intent(in) :: mu(pcols,pver) ! Mass flux up
real(r8), intent(in) :: md(pcols,pver) ! Mass flux down
real(r8), intent(in) :: du(pcols,pver) ! Mass detraining from updraft
real(r8), intent(in) :: eu(pcols,pver) ! Mass entraining from updraft
real(r8), intent(in) :: ed(pcols,pver) ! Mass entraining from downdraft
real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces
real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc
real(r8), intent(in) :: dt ! time step in seconds : 2*delta_t
integer, intent(in) :: jt(pcols) ! Index of cloud top for each column
integer, intent(in) :: mx(pcols) ! Index of cloud top for each column
integer, intent(in) :: ideep(pcols) ! Gathering array
integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
integer, intent(in) :: nstep ! Time step index
! input/output
real(r8), intent(out) :: dqdt(pcols,pver,ncnst) ! Tracer tendency array
!--------------------------Local Variables------------------------------
integer i ! Work index
integer k ! Work index
integer kbm ! Highest altitude index of cloud base
integer kk ! Work index
integer kkp1 ! Work index
integer kkm1 ! Work index
integer km1 ! Work index
integer kp1 ! Work index
integer ktm ! Highest altitude index of cloud top
integer m ! Work index
integer ii ! Work index
real(r8) cabv ! Mix ratio of constituent above
real(r8) cbel ! Mix ratio of constituent below
real(r8) cdifr ! Normalized diff between cabv and cbel
real(r8) chat(pcols,pver) ! Mix ratio in env at interfaces
real(r8) cond(pcols,pver) ! Mix ratio in downdraft at interfaces
real(r8) const(pcols,pver) ! Gathered wind array
real(r8) conu(pcols,pver) ! Mix ratio in updraft at interfaces
real(r8) dcondt(pcols,pver) ! Gathered tend array
real(r8) small ! A small number
real(r8) mbsth ! Threshold for mass fluxes
real(r8) mupdudp ! A work variable
real(r8) minc ! A work variable
real(r8) maxc ! A work variable
real(r8) fluxin ! A work variable
real(r8) fluxout ! A work variable
real(r8) netflux ! A work variable
real(r8) momcu ! constant for updraft pressure gradient term
real(r8) momcd ! constant for downdraft pressure gradient term
real(r8) sum ! sum
real(r8) sum2 ! sum2
real(r8) mududp(pcols,pver) ! working variable
real(r8) mddudp(pcols,pver) ! working variable
real(r8) pgu(pcols,pver) ! Pressure gradient term for updraft
real(r8) pgd(pcols,pver) ! Pressure gradient term for downdraft
real(r8),intent(out) :: pguall(pcols,pver,ncnst) ! Apparent force from updraft PG
real(r8),intent(out) :: pgdall(pcols,pver,ncnst) ! Apparent force from downdraft PG
real(r8),intent(out) :: icwu(pcols,pver,ncnst) ! In-cloud winds in updraft
real(r8),intent(out) :: icwd(pcols,pver,ncnst) ! In-cloud winds in downdraft
real(r8),intent(out) :: seten(pcols,pver) ! Dry static energy tendency
real(r8) gseten(pcols,pver) ! Gathered dry static energy tendency
real(r8) mflux(pcols,pverp,ncnst) ! Gathered momentum flux
real(r8) wind0(pcols,pver,ncnst) ! gathered wind before time step
real(r8) windf(pcols,pver,ncnst) ! gathered wind after time step
real(r8) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2
!-----------------------------------------------------------------------
!
! Initialize outgoing fields
pguall(:,:,:) = 0.0_r8
pgdall(:,:,:) = 0.0_r8
! Initialize in-cloud winds to environmental wind
icwu(:ncol,:,:) = q(:ncol,:,:)
icwd(:ncol,:,:) = q(:ncol,:,:)
! Initialize momentum flux and final winds
mflux(:,:,:) = 0.0_r8
wind0(:,:,:) = 0.0_r8
windf(:,:,:) = 0.0_r8
! Initialize dry static energy
seten(:,:) = 0.0_r8
gseten(:,:) = 0.0_r8
! Define constants for parameterization
momcu = 0.4_r8
momcd = 0.4_r8
small = 1.e-36_r8
! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
mbsth = 1.e-15_r8
! 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
! Loop ever each wind component
do m = 1, ncnst !start at m = 1 to transport momentum
if (domomtran(m)) then
! Gather up the winds and set tend to zero
do k = 1,pver
do i =il1g,il2g
const(i,k) = q(ideep(i),k,m)
wind0(i,k,m) = const(i,k)
end do
end do
! From now on work only with gathered data
! Interpolate winds to interfaces
do k = 1,pver
km1 = max(1,k-1)
do i = il1g, il2g
! use arithmetic mean
chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1))
! Provisional up and down draft values
conu(i,k) = chat(i,k)
cond(i,k) = chat(i,k)
! provisional tends
dcondt(i,k) = 0._r8
end do
end do
!
! Pressure Perturbation Term
!
!Top boundary: assume mu is zero
k=1
pgu(:il2g,k) = 0.0_r8
pgd(:il2g,k) = 0.0_r8
do k=2,pver-1
km1 = max(1,k-1)
kp1 = min(pver,k+1)
do i = il1g,il2g
!interior points
mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
+ mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
pgu(i,k) = - momcu * 0.5_r8 * mududp(i,k)
mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
+ md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
pgd(i,k) = - momcd * 0.5_r8 * mddudp(i,k)
end do
end do
! bottom boundary
k = pver
km1 = max(1,k-1)
do i=il1g,il2g
mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
pgu(i,k) = - momcu * mududp(i,k)
mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
pgd(i,k) = - momcd * mddudp(i,k)
end do
!
! In-cloud velocity calculations
!
! Do levels adjacent to top and bottom
k = 2
km1 = 1
kk = pver
kkm1 = max(1,kk-1)
do i = il1g,il2g
mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
if (mupdudp > mbsth) then
conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
endif
if (md(i,k) < -mbsth) then
cond(i,k) = (-ed(i,km1)*const(i,km1)*dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k)
endif
end do
! Updraft from bottom to top
do kk = pver-1,1,-1
kkm1 = max(1,kk-1)
kkp1 = min(pver,kk+1)
do i = il1g,il2g
mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
if (mupdudp > mbsth) then
conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* &
const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
endif
end do
end do
! Downdraft from top to bottom
do k = 3,pver
km1 = max(1,k-1)
do i = il1g,il2g
if (md(i,k) < -mbsth) then
cond(i,k) = ( md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) &
*dp(i,km1)-pgd(i,km1)*dp(i,km1) )/md(i,k)
endif
end do
end do
sum = 0._r8
sum2 = 0._r8
do k = ktm,pver
km1 = max(1,k-1)
kp1 = min(pver,k+1)
do i = il1g,il2g
ii = ideep(i)
! version 1 hard to check for roundoff errors
dcondt(i,k) = &
+(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
-mu(i,k)* (conu(i,k)-chat(i,k)) &
+md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
-md(i,k)* (cond(i,k)-chat(i,k)) &
)/dp(i,k)
end do
end do
! dcont for bottom layer
!
!DIR$ NOINTERCHANGE
do k = kbm,pver
km1 = max(1,k-1)
do i = il1g,il2g
if (k == mx(i)) then
! version 1
dcondt(i,k) = (1./dp(i,k))* &
(-mu(i,k)*(conu(i,k)-chat(i,k)) &
-md(i,k)*(cond(i,k)-chat(i,k)) &
)
end if
end do
end do
! Initialize to zero everywhere, then scatter tendency back to full array
dqdt(:,:,m) = 0._r8
do k = 1,pver
do i = il1g,il2g
ii = ideep(i)
dqdt(ii,k,m) = dcondt(i,k)
! Output apparent force on the mean flow from pressure gradient
pguall(ii,k,m) = -pgu(i,k)
pgdall(ii,k,m) = -pgd(i,k)
icwu(ii,k,m) = conu(i,k)
icwd(ii,k,m) = cond(i,k)
end do
end do
! Calculate momentum flux in units of mb*m/s2
do k = ktm,pver
do i = il1g,il2g
ii = ideep(i)
mflux(i,k,m) = &
-mu(i,k)* (conu(i,k)-chat(i,k)) &
-md(i,k)* (cond(i,k)-chat(i,k))
end do
end do
! Calculate winds at the end of the time step
do k = ktm,pver
do i = il1g,il2g
ii = ideep(i)
km1 = max(1,k-1)
kp1 = k+1
windf(i,k,m) = const(i,k) - (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k)
end do
end do
end if ! for domomtran
end do
! Need to add an energy fix to account for the dissipation of kinetic energy
! Formulation follows from Boville and Bretherton (2003)
! formulation by PJR
do k = ktm,pver
km1 = max(1,k-1)
kp1 = min(pver,k+1)
do i = il1g,il2g
ii = ideep(i)
! calculate the KE fluxes at top and bot of layer
! based on a discrete approximation to b&b eq(35) F_KE = u*F_u + v*F_v at interface
utop = (wind0(i,k,1)+wind0(i,km1,1))/2.
vtop = (wind0(i,k,2)+wind0(i,km1,2))/2.
ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2.
vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2.
fket = utop*mflux(i,k,1) + vtop*mflux(i,k,2) ! top of layer
fkeb = ubot*mflux(i,k+1,1) + vbot*mflux(i,k+1,2) ! bot of layer
! divergence of these fluxes should give a conservative redistribution of KE
ketend_cons = (fket-fkeb)/dp(i,k)
! tendency in kinetic energy resulting from the momentum transport
ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))*0.5/dt
! the difference should be the dissipation
gset2 = ketend_cons - ketend
gseten(i,k) = gset2
end do
end do
! Scatter dry static energy to full array
do k = 1,pver
do i = il1g,il2g
ii = ideep(i)
seten(ii,k) = gseten(i,k)
end do
end do
return
end subroutine momtran
!=========================================================================================
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._r8
end do
end do
!
do i = 1,ncol
lon(i) = pver
knt(i) = 0
lel(i) = pver
mx(i) = lon(i)
cape(i) = 0._r8
hmax(i) = 0._r8
end do
tp(:ncol,:) = t(:ncol,:)
qstp(:ncol,:) = q(:ncol,:)
!!! RBN - Initialize tv and buoy for output.
!!! tv=tv : tpv=tpv : qstp=q : buoy=0.
tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608*q(:ncol,:))/ (1._r8+q(:ncol,:))
tpv(:ncol,:) = tv(:ncol,:)
buoy(:ncol,:) = 0._r8
!
! 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_r8) 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._r8/ (3.5_r8*log(t(i,mx(i)))-log(e)-4.805_r8) + 55._r8
if (tl(i) < t(i,mx(i))) then
plexp(i) = (1._r8/ (0.2854_r8* (1._r8-0.28_r8*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._r8
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._r8+1.608_r8*q(i,k))/ (1._r8+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_r8* (1._r8-0.28_r8*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._r8+1.608_r8*q(i,mx(i)))/ (1._r8+q(i,mx(i)))
buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
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._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
qstp(i,k) = q(i,mx(i))
tp(i,k) = tl(i)* (p(i,k)/pl(i))**(0.2854_r8* (1._r8-0.28_r8*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._r8+ qstp(i,k) / eps1) * rl * eps1 / &
(rd * tp(i,k) ** 2)
a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* &
(1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ &
(rd**2*tp(i,k)**4)-qstp(i,k)* &
(1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ &
(rd*tp(i,k)**3))
a1(i) = 1._r8/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._r8+1.608_r8*qstp(i,k)) / (1._r8+q(i,mx(i)))
buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
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._r8+1.608_r8*q(i,k))/ (1._r8+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_r8* (1._r8-0.28_r8*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._r8+qstp(i,k)/eps1)*rl*eps1/ (rd*tp(i,k)**2)
a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* &
(1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ &
(rd**2*tp(i,k)**4)-qstp(i,k)* &
(1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ &
(rd*tp(i,k)**3))
a1(i) = 1._r8/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._r8+1.608_r8*qstp(i,k))/(1._r8+q(i,mx(i)))
buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
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._r8 .and. buoy(i,k) <= 0._r8) 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,2
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 ,landfrac)
!-----------------------------------------------------------------------
!
! 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
real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac
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)
! RBN For c0mask
real(r8) c0mask(pcols)
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._r8
expnum(i) = 0._r8
expdif(i) = 0._r8
c0mask(i) = c0_ocn * (1._r8-landfrac(i)) + c0_lnd * landfrac(i)
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._r8
i2(i,k) = 0._r8
i3(i,k) = 0._r8
i4(i,k) = 0._r8
mu(i,k) = 0._r8
f(i,k) = 0._r8
eps(i,k) = 0._r8
eu(i,k) = 0._r8
du(i,k) = 0._r8
ql(i,k) = 0._r8
cu(i,k) = 0._r8
evp(i,k) = 0._r8
cmeg(i,k) = 0._r8
qds(i,k) = q(i,k)
md(i,k) = 0._r8
ed(i,k) = 0._r8
sd(i,k) = s(i,k)
qd(i,k) = q(i,k)
mc(i,k) = 0._r8
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._r8 ) then
qst(i,k) = eps1*est(i)/ (p(i,k)-est(i))
else
qst(i,k) = 1.0_r8
end if
!--bee
gamma(i,k) = qst(i,k)*(1._r8 + 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._r8
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._r8
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._r8
totevp(i) = 0._r8
end do
do k = msg + 2,pver
do i = 1,il2g
if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_r8) 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_r8) 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)
!
jt(:) = 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_r8
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*tiedke_add
su(i,k) = s(i,mx(i)) + tiedke_add
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_r8* (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_r8* (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_r8* (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_r8
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._r8
ftemp(i) = 0._r8
if (k < jt(i) .or. k >= jb(i)) then
k1(i,k) = 0._r8
expnum(i) = 0._r8
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._r8 .and. expnum(i) > 0._r8) .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._r8*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* &
ftemp(i)**3 + (-5._r8*k1(i,k)*i2(i,k)*i3(i,k)+ &
5._r8*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_r8 .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._r8) then
mu(i,jb(i)) = 1._r8
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._r8 .and. (k >= jt(i) .and. k < jb(i))) then
zuef(i) = zf(i,k) - zf(i,jb(i))
rmue(i) = (1._r8/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._r8)/zuef(i)
mu(i,k) = (1._r8/eps0(i))* (exp(eps(i,k )*zuef(i))-1._r8)/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._r8) then
if (mu(i,k) < 0.01_r8) then
hu(i,k) = hu(i,jb(i))
mu(i,k) = 0._r8
eu(i,k) = 0._r8
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_r8) then
if (hu(i,k)-hsthat(i,k) < -2000._r8) then
jt(i) = k + 1
doit(i) = .false.
else
jt(i) = k
if (eps0(i) <= 0._r8) doit(i) = .false.
end if
else if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.01_r8) 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._r8) then
mu(i,k) = 0._r8
eu(i,k) = 0._r8
du(i,k) = 0._r8
hu(i,k) = hu(i,jb(i))
end if
if (k == jt(i) .and. eps0(i) > 0._r8) then
du(i,k) = mu(i,k+1)/dz(i,k)
eu(i,k) = 0._r8
mu(i,k) = 0._r8
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_r8
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._r8) 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._r8) then
zdef(i) = zf(i,jd(i)) - zf(i,k)
md(i,k) = -alfa(i)/ (2._r8*eps0(i))*(exp(2._r8*epsm(i)*zdef(i))-1._r8)/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._r8 .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_r8
do k = msg + 1,pver
do i = 1,il2g
if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._r8) 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._r8 .and. jd(i) < jb(i)) then
qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ &
(rl*(1._r8 + 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._r8) 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._r8-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._r8) 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._r8) then
su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._r8+gamhat(i,k)))
qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ &
(rl* (1._r8+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._r8) 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._r8
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
do k = pver,msg + 2,-1
do i = 1,il2g
rprd(i,k) = 0._r8
if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8 .and. mu(i,k) >= 0.0_r8) then
if (mu(i,k) > 0._r8) then
ql1 = 1._r8/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._r8+dz(i,k)*c0mask(i))
else
ql(i,k) = 0._r8
end if
totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1))
rprd(i,k) = c0mask(i)*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._r8) 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._r8) 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._r8 .and. totpcp(i) > 0._r8) 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._r8
ed(i,k) = 0._r8
evp(i,k) = 0._r8
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._r8
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
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
!
!-----------------------------------------------------------------------
#ifndef WRF_PORT
use dycore, only: dycore_is, get_resolution
#endif
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._r8
eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
dtbdt(i) = (1._r8/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._r8/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._r8* (3.5_r8/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
(3.5_r8*log(t(i,mx(i)))-log(eb)-4.805_r8)**2
end do
!
! dtmdt and dqmdt are cumulus heating and drying.
!
do k = msg + 1,pver
do i = il1g,il2g
dtmdt(i,k) = 0._r8
dqmdt(i,k) = 0._r8
end do
end do
!
do k = msg + 1,pver - 1
do i = il1g,il2g
if (k == jt(i)) then
dtmdt(i,k) = (1._r8/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._r8/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._r8
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._r8/p(i,k))** (rd/cp)*(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))
thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k))
dqsdtp(i,k) = qstp(i,k)* (1._r8+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._r8+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._r8/(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))* &
(1.608_r8 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_r8/ &
(1._r8+0.608_r8*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._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,mx(i)))
thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k))
!
! dboydt is the integrand of cape change.
!
dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_r8/ (1._r8+0.608_r8*q(i,mx(i)))*dqbdt(i)- &
dtmdt(i,k)/t(i,k)-0.608_r8/ (1._r8+0.608_r8*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._r8
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._r8* (cape(i)-capelmt)
if (dadt(i) /= 0._r8) 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._r8
dqdt(i,k) = 0._r8
dl(i,k) = 0._r8
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
!
!DIR$ NOINTERCHANGE!
do k = kbm,pver
do i = il1g,il2g
if (k == mx(i)) then
dsdt(i,k) = (1._r8/dsubcld(i))* &
(-mu(i,k)* (su(i,k)-shat(i,k)) &
-md(i,k)* (sd(i,k)-shat(i,k)) &
)
dqdt(i,k) = (1._r8/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
subroutine buoyan_dilute(lchnk ,ncol , & 1,1
q ,t ,p ,z ,pf , &
tp ,qstp ,tl ,rl ,cape , &
pblt ,lcl ,lel ,lon ,mx , &
rd ,grav ,cp ,msg , &
tpert )
!-----------------------------------------------------------------------
!
! Purpose:
! Calculates CAPE the lifting condensation level and the convective top
! where buoyancy is first -ve.
!
! Method: Calculates the parcel temperature based on a simple constant
! entraining plume model. CAPE is integrated from buoyancy.
! 09/09/04 - Simplest approach using an assumed entrainment rate for
! testing (dmpdp).
! 08/04/05 - Swap to convert dmpdz to dmpdp
!
! SCAM Logical Switches - DILUTE:RBN - Now Disabled
! ---------------------
! switch(1) = .T. - Uses the dilute parcel calculation to obtain tendencies.
! switch(2) = .T. - Includes entropy/q changes due to condensate loss and freezing.
! switch(3) = .T. - Adds the PBL Tpert for the parcel temperature at all levels.
!
! References:
! Raymond and Blythe (1992) JAS
!
! Author:
! Richard Neale - September 2004
!
!-----------------------------------------------------------------------
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 (only above lcl, just q below).
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._r8
end do
end do
!
do i = 1,ncol
lon(i) = pver
knt(i) = 0
lel(i) = pver
mx(i) = lon(i)
cape(i) = 0._r8
hmax(i) = 0._r8
end do
tp(:ncol,:) = t(:ncol,:)
qstp(:ncol,:) = q(:ncol,:)
!!! RBN - Initialize tv and buoy for output.
!!! tv=tv : tpv=tpv : qstp=q : buoy=0.
tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608_r8*q(:ncol,:))/ (1._r8+q(:ncol,:))
tpv(:ncol,:) = tv(:ncol,:)
buoy(:ncol,:) = 0._r8
!
! 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_r8) 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
! LCL dilute calculation - initialize to mx(i)
! Determine lcl in parcel_dilute and get pl,tl after parcel_dilute
! Original code actually sets LCL as level above wher condensate forms.
! Therefore in parcel_dilute lcl(i) will be at first level where qsmix < qtmix.
do i = 1,ncol ! Initialise LCL variables.
lcl(i) = mx(i)
tl(i) = t(i,mx(i))
pl(i) = p(i,mx(i))
end do
!
! main buoyancy calculation.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! DILUTE PLUME CALCULATION USING ENTRAINING PLUME !!!
!!! RBN 9/9/04 !!!
call parcel_dilute(lchnk, ncol, msg, mx, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl)
! 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._r8 ! Just change to always allow buoy calculation.
end do
!
! Main buoyancy calculation.
!
do k = pver,msg + 1,-1
do i=1,ncol
if (k <= mx(i) .and. plge600(i)) then ! Define buoy from launch level to cloud top.
tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add ! +0.5K or not?
else
qstp(i,k) = q(i,k)
tp(i,k) = t(i,k)
tpv(i,k) = tv(i,k)
endif
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._r8) 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_dilute
subroutine parcel_dilute (lchnk, ncol, msg, klaunch, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl) 1,6
! Routine to determine
! 1. Tp - Parcel temperature
! 2. qstp - Saturated mixing ratio at the parcel temperature.
!--------------------
implicit none
!--------------------
integer, intent(in) :: lchnk
integer, intent(in) :: ncol
integer, intent(in) :: msg
integer, intent(in), dimension(pcols) :: klaunch(pcols)
real(r8), intent(in), dimension(pcols,pver) :: p
real(r8), intent(in), dimension(pcols,pver) :: t
real(r8), intent(in), dimension(pcols,pver) :: q
real(r8), intent(in), dimension(pcols) :: tpert ! PBL temperature perturbation.
real(r8), intent(inout), dimension(pcols,pver) :: tp ! Parcel temp.
real(r8), intent(inout), dimension(pcols,pver) :: qstp ! Parcel water vapour (sat value above lcl).
real(r8), intent(inout), dimension(pcols) :: tl ! Actual temp of LCL.
real(r8), intent(inout), dimension(pcols) :: pl ! Actual pressure of LCL.
integer, intent(inout), dimension(pcols) :: lcl ! Lifting condesation level (first model level with saturation).
real(r8), intent(out), dimension(pcols,pver) :: tpv ! Define tpv within this routine.
!--------------------
! Have to be careful as s is also dry static energy.
! If we are to retain the fact that CAM loops over grid-points in the internal
! loop then we need to dimension sp,atp,mp,xsh2o with ncol.
real(r8) tmix(pcols,pver) ! Tempertaure of the entraining parcel.
real(r8) qtmix(pcols,pver) ! Total water of the entraining parcel.
real(r8) qsmix(pcols,pver) ! Saturated mixing ratio at the tmix.
real(r8) smix(pcols,pver) ! Entropy of the entraining parcel.
real(r8) xsh2o(pcols,pver) ! Precipitate lost from parcel.
real(r8) ds_xsh2o(pcols,pver) ! Entropy change due to loss of condensate.
real(r8) ds_freeze(pcols,pver) ! Entropy change sue to freezing of precip.
real(r8) mp(pcols) ! Parcel mass flux.
real(r8) qtp(pcols) ! Parcel total water.
real(r8) sp(pcols) ! Parcel entropy.
real(r8) sp0(pcols) ! Parcel launch entropy.
real(r8) qtp0(pcols) ! Parcel launch total water.
real(r8) mp0(pcols) ! Parcel launch relative mass flux.
real(r8) lwmax ! Maximum condesate that can be held in cloud before rainout.
real(r8) dmpdp ! Parcel fractional mass entrainment rate (/mb).
!real(r8) dmpdpc ! In cloud parcel mass entrainment rate (/mb).
real(r8) dmpdz ! Parcel fractional mass entrainment rate (/m)
real(r8) dpdz,dzdp ! Hydrstatic relation and inverse of.
real(r8) senv ! Environmental entropy at each grid point.
real(r8) qtenv ! Environmental total water " " ".
real(r8) penv ! Environmental total pressure " " ".
real(r8) tenv ! Environmental total temperature " " ".
real(r8) new_s ! Hold value for entropy after condensation/freezing adjustments.
real(r8) new_q ! Hold value for total water after condensation/freezing adjustments.
real(r8) dp ! Layer thickness (center to center)
real(r8) tfguess ! First guess for entropy inversion - crucial for efficiency!
real(r8) tscool ! Super cooled temperature offset (in degC) (eg -35).
real(r8) qxsk, qxskp1 ! LCL excess water (k, k+1)
real(r8) dsdp, dqtdp, dqxsdp ! LCL s, qt, p gradients (k, k+1)
real(r8) slcl,qtlcl,qslcl ! LCL s, qt, qs values.
integer rcall ! Number of ientropy call for errors recording
integer nit_lheat ! Number of iterations for condensation/freezing loop.
integer i,k,ii ! Loop counters.
!======================================================================
! SUMMARY
!
! 9/9/04 - Assumes parcel is initiated from level of maxh (klaunch)
! and entrains at each level with a specified entrainment rate.
!
! 15/9/04 - Calculates lcl(i) based on k where qsmix is first < qtmix.
!
!======================================================================
!
! Set some values that may be changed frequently.
!
nit_lheat = 2 ! iterations for ds,dq changes from condensation freezing.
dmpdz=-1.e-3_r8 ! Entrainment rate. (-ve for /m)
!dmpdpc = 3.e-2_r8 ! In cloud entrainment rate (/mb).
lwmax = 1.e-3_r8 ! Need to put formula in for this.
tscool = 0.0_r8 ! Temp at which water loading freezes in the cloud.
qtmix=0._r8
smix=0._r8
qtenv = 0._r8
senv = 0._r8
tenv = 0._r8
penv = 0._r8
qtp0 = 0._r8
sp0 = 0._r8
mp0 = 0._r8
qtp = 0._r8
sp = 0._r8
mp = 0._r8
new_q = 0._r8
new_s = 0._r8
! **** Begin loops ****
do k = pver, msg+1, -1
do i=1,ncol
! Initialize parcel values at launch level.
if (k == klaunch(i)) then
qtp0(i) = q(i,k) ! Parcel launch total water (assuming subsaturated) - OK????.
sp0(i) = entropy(t(i,k),p(i,k),qtp0(i)) ! Parcel launch entropy.
mp0(i) = 1._r8 ! Parcel launch relative mass (i.e. 1 parcel stays 1 parcel for dmpdp=0, undilute).
smix(i,k) = sp0(i)
qtmix(i,k) = qtp0(i)
tfguess = t(i,k)
rcall = 1
call ientropy (rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess)
end if
! Entraining levels
if (k < klaunch(i)) then
! Set environmental values for this level.
dp = (p(i,k)-p(i,k+1)) ! In -ve mb as p decreasing with height - difference between center of layers.
qtenv = 0.5_r8*(q(i,k)+q(i,k+1)) ! Total water of environment.
tenv = 0.5_r8*(t(i,k)+t(i,k+1))
penv = 0.5_r8*(p(i,k)+p(i,k+1))
senv = entropy(tenv,penv,qtenv) ! Entropy of environment.
! Determine fractional entrainment rate /pa given value /m.
dpdz = -(penv*grav)/(rgas*tenv) ! in mb/m since p in mb.
dzdp = 1._r8/dpdz ! in m/mb
dmpdp = dmpdz*dzdp ! /mb Fractional entrainment
! Sum entrainment to current level
! entrains q,s out of intervening dp layers, in which linear variation is assumed
! so really it entrains the mean of the 2 stored values.
sp(i) = sp(i) - dmpdp*dp*senv
qtp(i) = qtp(i) - dmpdp*dp*qtenv
mp(i) = mp(i) - dmpdp*dp
! Entrain s and qt to next level.
smix(i,k) = (sp0(i) + sp(i)) / (mp0(i) + mp(i))
qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i))
! Invert entropy from s and q to determine T and saturation-capped q of mixture.
! t(i,k) used as a first guess so that it converges faster.
tfguess = tmix(i,k+1)
rcall = 2
call ientropy(rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess)
!
! Determine if this is lcl of this column if qsmix <= qtmix.
! FIRST LEVEL where this happens on ascending.
if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then
lcl(i) = k
qxsk = qtmix(i,k) - qsmix(i,k)
qxskp1 = qtmix(i,k+1) - qsmix(i,k+1)
dqxsdp = (qxsk - qxskp1)/dp
pl(i) = p(i,k+1) - qxskp1/dqxsdp ! pressure level of actual lcl.
dsdp = (smix(i,k) - smix(i,k+1))/dp
dqtdp = (qtmix(i,k) - qtmix(i,k+1))/dp
slcl = smix(i,k+1) + dsdp* (pl(i)-p(i,k+1))
qtlcl = qtmix(i,k+1) + dqtdp*(pl(i)-p(i,k+1))
tfguess = tmix(i,k)
rcall = 3
call ientropy (rcall,i,lchnk,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess)
! write(iulog,*)' '
! write(iulog,*)' p',p(i,k+1),pl(i),p(i,lcl(i))
! write(iulog,*)' t',tmix(i,k+1),tl(i),tmix(i,lcl(i))
! write(iulog,*)' s',smix(i,k+1),slcl,smix(i,lcl(i))
! write(iulog,*)'qt',qtmix(i,k+1),qtlcl,qtmix(i,lcl(i))
! write(iulog,*)'qs',qsmix(i,k+1),qslcl,qsmix(i,lcl(i))
endif
!
end if ! k < klaunch
end do ! Levels loop
end do ! Columns loop
!!!!!!!!!!!!!!!!!!!!!!!!!!END ENTRAINMENT LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Could stop now and test with this as it will provide some estimate of buoyancy
!! without the effects of freezing/condensation taken into account for tmix.
!! So we now have a profile of entropy and total water of the entraining parcel
!! Varying with height from the launch level klaunch parcel=environment. To the
!! top allowed level for the existence of convection.
!! Now we have to adjust these values such that the water held in vaopor is < or
!! = to qsmix. Therefore, we assume that the cloud holds a certain amount of
!! condensate (lwmax) and the rest is rained out (xsh2o). This, obviously
!! provides latent heating to the mixed parcel and so this has to be added back
!! to it. But does this also increase qsmix as well? Also freezing processes
xsh2o = 0._r8
ds_xsh2o = 0._r8
ds_freeze = 0._r8
!!!!!!!!!!!!!!!!!!!!!!!!!PRECIPITATION/FREEZING LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Iterate solution twice for accuracy
do k = pver, msg+1, -1
do i=1,ncol
! Initialize variables at k=klaunch
if (k == klaunch(i)) then
! Set parcel values at launch level assume no liquid water.
tp(i,k) = tmix(i,k)
qstp(i,k) = q(i,k)
tpv(i,k) = (tp(i,k) + tpert(i)) * (1._r8+1.608_r8*qstp(i,k)) / (1._r8+qstp(i,k))
end if
if (k < klaunch(i)) then
! Initiaite loop if switch(2) = .T. - RBN:DILUTE - TAKEN OUT BUT COULD BE RETURNED LATER.
! Iterate nit_lheat times for s,qt changes.
do ii=0,nit_lheat-1
! Rain (xsh2o) is excess condensate, bar LWMAX (Accumulated loss from qtmix).
xsh2o(i,k) = max (0._r8, qtmix(i,k) - qsmix(i,k) - lwmax)
! Contribution to ds from precip loss of condensate (Accumulated change from smix).(-ve)
ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._r8,(xsh2o(i,k)-xsh2o(i,k+1)))
!
! Entropy of freezing: latice times amount of water involved divided by T.
!
if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._r8) then ! One off freezing of condensate.
ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._r8,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) ! Gain of LH
end if
if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._r8) then ! Continual freezing of additional condensate.
ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._r8,(qsmix(i,k+1)-qsmix(i,k)))
end if
! Adjust entropy and accordingly to sum of ds (be careful of signs).
new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k)
! Adjust liquid water and accordingly to xsh2o.
new_q = qtmix(i,k) - xsh2o(i,k)
! Invert entropy to get updated Tmix and qsmix of parcel.
tfguess = tmix(i,k)
rcall =4
call ientropy (rcall,i,lchnk,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess)
end do ! Iteration loop for freezing processes.
! tp - Parcel temp is temp of mixture.
! tpv - Parcel v. temp should be density temp with new_q total water.
tp(i,k) = tmix(i,k)
! tpv = tprho in the presence of condensate (i.e. when new_q > qsmix)
if (new_q > qsmix(i,k)) then ! Super-saturated so condensate present - reduces buoyancy.
qstp(i,k) = qsmix(i,k)
else ! Just saturated/sub-saturated - no condensate virtual effects.
qstp(i,k) = new_q
end if
tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+ new_q)
end if ! k < klaunch
end do ! Loop for columns
end do ! Loop for vertical levels.
return
end subroutine parcel_dilute
!-----------------------------------------------------------------------------------------
real(r8) function entropy(TK,p,qtot) 2
!-----------------------------------------------------------------------------------------
!
! TK(K),p(mb),qtot(kg/kg)
! from Raymond and Blyth 1992
!
real(r8), intent(in) :: p,qtot,TK
real(r8) :: qv,qsat,e,esat,L,eref,pref
pref = 1000.0_r8 ! mb
eref = 6.106_r8 ! sat p at tfreez (mb)
L = rl - (cpliq - cpwv)*(TK-tfreez) ! T IN CENTIGRADE
! Replace call to satmixutils.
esat = c1*exp(c2*(TK-tfreez)/(c3+TK-tfreez)) ! esat(T) in mb
qsat=eps1*esat/(p-esat) ! Sat. mixing ratio (in kg/kg).
qv = min(qtot,qsat) ! Partition qtot into vapor part only.
e = qv*p / (eps1 +qv)
entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + &
L*qv/TK - qv*rh2o*log(qv/qsat)
!
return
end FUNCTION entropy
!
!-----------------------------------------------------------------------------------------
SUBROUTINE ientropy (rcall,icol,lchnk,s,p,qt,T,qsat,Tfg) 4,4
!-----------------------------------------------------------------------------------------
!
! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg).
! Inverts entropy, pressure and total water qt
! for T and saturated vapor mixing ratio
!
#ifndef WRF_PORT
use phys_grid, only: get_rlon_p, get_rlat_p
#endif
integer, intent(in) :: icol, lchnk, rcall
real(r8), intent(in) :: s, p, Tfg, qt
real(r8), intent(out) :: qsat, T
real(r8) :: qv,Ts,dTs,fs1,fs2,esat
real(r8) :: pref,eref,L,e
real(r8) :: this_lat,this_lon
integer :: LOOPMAX,i
LOOPMAX = 100 !* max number of iteration loops
! Values for entropy
pref = 1000.0_r8 ! mb ref pressure.
eref = 6.106_r8 ! sat p at tfreez (mb)
! Invert the entropy equation -- use Newton's method
Ts = Tfg ! Better first guess based on Tprofile from conv.
converge: do i=0, LOOPMAX
L = rl - (cpliq - cpwv)*(Ts-tfreez)
esat = c1*exp(c2*(Ts-tfreez)/(c3+Ts-tfreez)) ! Bolton (eq. 10)
qsat = eps1*esat/(p-esat)
qv = min(qt,qsat)
e = qv*p / (eps1 +qv) ! Bolton (eq. 16)
fs1 = (cpres + qt*cpliq)*log( Ts/tfreez ) - rgas*log( (p-e)/pref ) + &
L*qv/Ts - qv*rh2o*log(qv/qsat) - s
L = rl - (cpliq - cpwv)*(Ts-1._r8-tfreez)
esat = c1*exp(c2*(Ts-1._r8-tfreez)/(c3+Ts-1._r8-tfreez))
qsat = eps1*esat/(p-esat)
qv = min(qt,qsat)
e = qv*p / (eps1 +qv)
fs2 = (cpres + qt*cpliq)*log( (Ts-1._r8)/tfreez ) - rgas*log( (p-e)/pref ) + &
L*qv/(Ts-1._r8) - qv*rh2o*log(qv/qsat) - s
dTs = fs1/(fs2 - fs1)
Ts = Ts+dTs
if (abs(dTs).lt.0.001_r8) exit converge
if (i .eq. LOOPMAX - 1) then
#ifndef WRF_PORT
this_lat = get_rlat_p(lchnk, icol)*57.296_r8
this_lon = get_rlon_p(lchnk, icol)*57.296_r8
#else
!Do not worry about the specific lat/lon in WRF
this_lat = 0.
this_lon = 0.
#endif
write(iulog,*) '*** ZM_CONV: IENTROPY: Failed and about to exit, info follows ****'
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,100) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, &
' lat: ',this_lat,' lon: ',this_lon, &
' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._r8*qt, &
' qsat(g/kg) = ', 1000._r8*qsat,', s(J/kg) = ',s
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
write(iulog,*) '*** Please report this crash to Po-Lun.Ma@pnnl.gov ***'
#ifdef WRF_PORT
call wrf_debug(1,iulog)
#endif
call endrun('**** ZM_CONV IENTROPY: Tmix did not converge ****')
end if
enddo converge
! Replace call to satmixutils.
esat = c1*exp(c2*(Ts-tfreez)/(c3+Ts-tfreez))
qsat=eps1*esat/(p-esat)
qv = min(qt,qsat) ! /* check for saturation */
T = Ts
100 format (A,I1,I4,I4,7(A,F6.2))
return
end SUBROUTINE ientropy
end module module_cu_camzm