MODULE CLOUDS 2,3
!@sum CLOUDS column physics of moist conv. and large-scale condensation
!@auth M.S.Yao/A. Del Genio (modifications by Gavin Schmidt)
!@cont MSTCNV,LSCOND
USE CONSTANT
, only : rgas,grav,lhe,lhs,lhm,sha,bysha,pi,by6
* ,by3,tf,bytf,rvap,bygrav,deltx,bymrat,teeny,gamd,rhow
* ,twopi
USE MODEL_COM, only : im,lm,dtsrc,itime,coupled_chem
USE QUSDEF, only : nmom,xymoms,zmoms,zdir
CCC USE RANDOM
IMPLICIT NONE
SAVE
C**** parameters and constants
REAL*8, PARAMETER :: TI=233.16d0 !@param TI pure ice limit
REAL*8, PARAMETER :: CLDMIN=.25d0 !@param CLDMIN min MC/LSC region
!@param WMU critical cloud water content for rapid conversion (g m**-3)
REAL*8, PARAMETER :: WMU=.25
REAL*8, PARAMETER :: WMUL=.5 !@param WMUL WMU over land
REAL*8, PARAMETER :: WMUI=.1d0 !@param WMUI WMU for ice clouds
REAL*8, PARAMETER :: BRCLD=.2d0 !@param BRCLD for cal. BYBR
REAL*8, PARAMETER :: SLHE=LHE*BYSHA
REAL*8, PARAMETER :: SLHS=LHS*BYSHA
!@param CCMUL multiplier for convective cloud cover
!@param CCMUL1 multiplier for deep anvil cloud cover
!@param CCMUL2 multiplier for shallow anvil cloud cover
!@param COETAU multiplier for convective cloud optical thickness
REAL*8, PARAMETER :: CCMUL=1.,CCMUL1=5.,CCMUL2=3.,COETAU=.08d0
REAL*8 :: BYBR,BYDTsrc,XMASS,PLAND
!@var BYBR factor for converting cloud particle radius to effect. radius
!@var XMASS dummy variable
!@var PLAND land fraction
C**** Set-able variables
!@dbparam LMCM max level for originating MC plumes
INTEGER :: LMCM = -1 ! defaults to LS1-1 if not set in rundeck
!@dbparam ISC integer to turn on computation of stratocumulus clouds
INTEGER :: ISC = 0 ! set ISC=1 to compute stratocumulus clouds
!@dbparam U00wtrX multiplies U00ice for critical humidity for water clds
REAL*8 :: U00wtrX = 1.0d0 ! default
!@dbparam U00ice critical humidity for ice cloud condensation
REAL*8 :: U00ice = .7d0 ! default
!@dbparam HRMAX maximum distance an air parcel rises from surface
REAL*8 :: HRMAX = 1000.d0 ! default (m)
!@dbparam RWCldOX multiplies part.size of water clouds over ocean
REAL*8 :: RWCldOX=1.d0
!@dbparam RICldX multiplies part.size of ice clouds at 1000mb
!@+ RICldX changes linearly to 1 as p->0mb
REAL*8 :: RICldX=1.d0 , xRICld
!@dbparam do_blU00 =1 if boundary layer U00 is treated differently
INTEGER :: do_blU00=0 ! default is to disable this
C**** input variables
LOGICAL DEBUG
!@var RA ratio of primary grid box to secondary gridbox
REAL*8, DIMENSION(IM) :: RA
!@var UM,VM,UM1,VM1,U_0,V_0 velocity related variables(UM,VM)=(U,V)*AIRM
REAL*8, DIMENSION(IM,LM) :: UM,VM,UM1,VM1
REAL*8, DIMENSION(IM,LM) :: U_0,V_0
!@var Miscellaneous vertical arrays set in driver
!@var PLE pressure at layer edge
!@var LHP array of precip phase ! may differ from LHX
REAL*8, DIMENSION(LM+1) :: PLE,LHP
REAL*8, DIMENSION(LM) :: PL,PLK,AIRM,BYAM,ETAL,TL,QL,TH,RH,WMX
* ,VSUBL,MCFLX,DGDSM,DPHASE,DTOTW,DQCOND,DGDQM,AQ,DPDT,RH1
* ,FSSL,FSUB,FCONV,FMCL,VLAT,DDMFLX
!@var PL layer pressure (mb)
!@var PLK PL**KAPA
!@var AIRM the layer's pressure depth (mb)
!@var BYAM 1./AIRM
!@var ETAL fractional entrainment rate
!@var TL, QL temperature, specific humidity of the layer
!@var TH potential temperature (K)
!@var RH relative humidity
!@var RH1 relative humidity to compare with the threshold humidity
!@var WMX cloud water mixing ratio (kg/kg)
!@var VSUBL downward vertical velocity due to cumulus subsidence (cm/s)
!@var MCFLX, DGDSM, DPHASE, DQCOND, DGDQM dummy variables
!@var DDMFLX accumulated downdraft mass flux (mb)
!@var AQ time change rate of specific humidity (s**-1)
!@var DPDT time change rate of pressure (mb/s)
!@var FSSL grid fraction for large-scale clouds
!@var FCONV convective fraction
!@var FSUB subsiding fraction
!@var FMCL grid fraction for moist convection
!@var VLAT dummy variable
REAL*8, DIMENSION(LM+1) :: PRECNVL
!@var PRECNVL convective precip entering the layer top
C**** new arrays must be set to model arrays in driver (before MSTCNV)
REAL*8, DIMENSION(LM) :: SDL,WML
!@var SDL vertical velocity in sigma coordinate
!@var WML cloud water mixing ratio (kg/kg)
C**** new arrays must be set to model arrays in driver (after MSTCNV)
REAL*8, DIMENSION(LM) :: TAUMCL,SVLATL,CLDMCL,SVLHXL,SVWMXL
!@var TAUMCL convective cloud optical thickness
!@var SVLATL saved LHX for convective cloud
!@var CLDMCL convective cloud cover
!@var SVLHXL saved LHX for large-scale cloud
!@var SVWMXL saved detrained convective cloud water
REAL*8, DIMENSION(LM) :: CSIZEL
!@var CSIZEL cloud particle radius (micron)
C**** new arrays must be set to model arrays in driver (before LSCOND)
REAL*8, DIMENSION(LM) :: TTOLDL,CLDSAVL,CLDSV1
!@var TTOLDL previous potential temperature
!@var CLDSAVL saved large-scale cloud cover
C**** new arrays must be set to model arrays in driver (after LSCOND)
REAL*8, DIMENSION(LM) :: SSHR,DCTEI,TAUSSL,CLDSSL
!@var SSHR,DCTEI height diagnostics of dry and latent heating by MC
!@var TAUSSL large-scale cloud optical thickness
!@var CLDSSL large-scale cloud cover
!@var SM,QM Vertical profiles of T/Q
REAL*8, DIMENSION(LM) :: SM,QM
REAL*8, DIMENSION(NMOM,LM) :: SMOM,QMOM,SMOMMC,QMOMMC,
* SMOMLS,QMOMLS
!@var KMAX index for surrounding velocity
!@var LP50 50mb level
INTEGER :: KMAX,LP50
!@var PEARTH fraction of land in grid box
!@var TS average surface temperture (C)
!@var RIS, RI1, RI2 Richardson numbers
REAL*8 :: PEARTH,TS,QS,US,VS,RIS,RI1,RI2,DXYPJ
!@var DCL max level of planetary boundary layer
INTEGER :: DCL
C**** output variables
REAL*8 :: PRCPMC,PRCPSS,HCNDSS,WMSUM
!@var PRCPMC precip due to moist convection
!@var PRCPSS precip due to large-scale condensation
!@var HCNDSS heating due to large-scale condensation
!@var WMSUM cloud liquid water path
REAL*8 :: CLDSLWIJ,CLDDEPIJ
!@var CLDSLWIJ shallow convective cloud cover
!@var CLDDEPIJ deep convective cloud cover
INTEGER :: LMCMAX,LMCMIN
!@var LMCMAX upper-most convective layer
!@var LMCMIN lowerest convective layer
!@var AIRXL is convective mass flux (mb)
REAL*8 AIRXL,PRHEAT
!@var RNDSSL stored random number sequences
REAL*8 RNDSSL(3,LM)
CCOMP does not work yet:
CCOMP THREADPRIVATE (RA,UM,VM,U_0,V_0,PLE,PL,PLK,AIRM,BYAM,ETAL
CCOMP* ,TL,QL,TH,RH,WMX,VSUBL,MCFLX,SSHR,DGDSM,DPHASE
CCOMP* ,DTOTW,DQCOND,DCTEI,DGDQM,dxypj,DDMFLX
CCOMP* ,AQ,DPDT,PRECNVL,SDL,WML,SVLATL,SVLHXL,SVWMXL,CSIZEL,RH1
CCOMP* ,TTOLDL,CLDSAVL,TAUMCL,CLDMCL,TAUSSL,CLDSSL,RNDSSL
CCOMP* ,SM,QM,SMOM,QMOM,PEARTH,TS,QS,US,VS,DCL,RIS,RI1,RI2, AIRXL
CCOMP* ,PRCPMC,PRCPSS,HCNDSS,WMSUM,CLDSLWIJ,CLDDEPIJ,LMCMAX
CCOMP* ,LMCMIN,KMAX,DEBUG)
COMMON/CLDPRV/RA,UM,VM,UM1,VM1,U_0,V_0,PLE,PL,PLK,AIRM,BYAM,ETAL
* ,TL,QL,TH,RH,WMX,VSUBL,MCFLX,SSHR,DGDSM,DPHASE,LHP
* ,DTOTW,DQCOND,DCTEI,DGDQM,DXYPJ,DDMFLX,PLAND
* ,AQ,DPDT,PRECNVL,SDL,WML,SVLATL,SVLHXL,SVWMXL,CSIZEL,RH1
* ,TTOLDL,CLDSAVL,TAUMCL,CLDMCL,TAUSSL,CLDSSL,RNDSSL
* ,SM,QM,SMOM,QMOM,PEARTH,TS,QS,US,VS,RIS,RI1,RI2, AIRXL
* ,SMOMMC,QMOMMC,SMOMLS,QMOMLS,CLDSV1,PRHEAT
* ,PRCPMC,PRCPSS,HCNDSS,WMSUM,CLDSLWIJ,CLDDEPIJ,VLAT
* ,FSUB,FCONV,FSSL,FMCL
* ,LMCMAX,LMCMIN,KMAX,DCL,DEBUG ! int/logic last (alignment)
!$OMP THREADPRIVATE (/CLDPRV/)
CONTAINS
SUBROUTINE MSTCNV(IERR,LERR,i_debug,j_debug) 1,22
!@sum MSTCNV moist convective processes (precip, convective clouds,...)
!@auth M.S.Yao/A. Del Genio (modularisation by Gavin Schmidt)
!@ver 1.0 (taken from CB265)
!@calls adv1d,QSAT,DQSATDT,THBAR
IMPLICIT NONE
REAL*8 LHX,MPLUME,MPLUM1,MCLOUD,MPMAX,SENV,QENV
!@var LHX latent heat of evaporation (J/Kg)
!@var MPLUME,MPLUM1 mass of convective plume (mb)
!@var MCLOUD air mass available for re-evaporation of precip
!@var MPMAX convective plume at the detrainment level
!@var SENV,QENV dummy variables
!@var FMC1 grid fraction for moist convection
REAL*8 FMC1
C**** functions
REAL*8 :: QSAT, DQSATDT
!@var QSAT saturation humidity
!@var DQSATDT dQSAT/dT
REAL*8, DIMENSION(0:LM) :: CM !@var CM air mass of subsidence
REAL*8, DIMENSION(IM) :: UMP,VMP,UMDN,VMDN
!@var UMP, VMP momentum carried by convective plumes
!@var UMDN,VMDN dummy variables
!@var DQM,DSM,DQMR,DSMR Vertical profiles of T/Q and changes
REAL*8, DIMENSION(LM) ::
* SMOLD,QMOLD, DQM,DSM,DQMR,DSMR
!@var SMOLD,QMOLD profiles prior to any moist convection
REAL*8, DIMENSION(LM) :: F,CMNEG
REAL*8, DIMENSION(NMOM,LM) :: FMOM
REAL*8, DIMENSION(NMOM,LM) :: SMOMOLD,QMOMOLD
REAL*8, DIMENSION(NMOM,LM) :: DSMOM,DQMOM,DSMOMR,DQMOMR
REAL*8, DIMENSION(NMOM) ::
& SMOMP,QMOMP, SMOMPMAX,QMOMPMAX, SMOMDN,QMOMDN
REAL*8, DIMENSION(LM) ::
* DM,COND,CDHEAT,CCM,SM1,QM1,DMR,ML,SMT,QMT,TPSAV,DDM,CONDP
REAL*8 :: CONDGP,CONDIP
!@var DM change in air mass
!@var COND,CONDGP,CONDIP condensate
!@var CDHEAT heating due to condensation
!@var CCM convective plume mass
!@var SM1, QM1 dummy variables
!@var DMR change of air mass
!@var ML layer air mass
!@var SMT, QMT dummy variables
!@var TPSAV array to save plume temperature
!@var DDM downdraft mass
!@param FCLW fraction of condensate in plume that remains as CLW
REAL*8 :: FCLW
!@var IERRT,LERRT error reports from advection
INTEGER :: IERRT,LERRT
INTEGER, INTENT(IN) :: i_debug,j_debug
INTEGER LDRAFT,LMAX,LMIN,MCCONT,MAXLVL
* ,MINLVL,ITER,IC,LFRZ,NSUB,LDMIN
!@var LDRAFT the layer the downdraft orginates
!@var LEVAP the layer evaporation of precip starts
!@var LMAX, LMIN the lowest, the highest layer of a convective event
!@var MCCONT integer to count convective events
!@var MAXLVL, MINLVL the lowest, the highest layer of convective events
!@var ITER number for iteration
!@var IC integer for cloud types
!@var LFRZ freezing level
!@var nsub = LMAX - LMIN + 1
!@var LDMIN the lowest layer the downdraft drops
REAL*8 TERM1,FMP0,SMO1
* ,QMO1,SMO2,QMO2,SDN,QDN,SUP,QUP,SEDGE,QEDGE,WMDN,WMUP,SVDN
* ,SVUP,WMEDG,SVEDG,DMSE,FPLUME,DFP,FMP2,FRAT1,FRAT2,SMN1
* ,QMN1,SMN2,QMN2,SMP,QMP,TP,GAMA,DQSUM,TNX,TNX1
* ,DQ,DMSE1,FCTYPE,BETAU,ALPHAU
* ,CDHDRT,DDRAFT,DELTA
* ,ALPHA,BETA,CDHM,CDHSUM,CLDM,CLDREF,CONSUM,DQEVP
* ,DQRAT,EPLUME,ETADN,ETAL1,EVPSUM,FCDH
* ,FCDH1,FCLD,FCLOUD,FDDL,FDDP,FENTR,FENTRA,FEVAP,FLEFT
* ,FQCOND,FQEVP,FPRCP,FSEVP,FSSUM,HEAT1
* ,PRCP
* ,QMDN,QMIX,QMPMAX,QMPT,QNX,QSATC,QSATMP
* ,RCLD,RCLDE,SLH,SMDN,SMIX,SMPMAX,SMPT,SUMAJ
* ,SUMDP,DDRUP,EDRAFT
* ,TOLD,TOLD1,TEMWM,TEM,WTEM,WCONST,WORK
* ,FCONV_tmp,FSUB_tmp,FSSL_tmp
* , MNdO,MNdL,MNdI,MCDNCW,MCDNCI !Menon
!@var TERM1 contribution to non-entraining convective cloud
!@var FMP0 non-entraining convective mass
!@var SMO1,QMO1,SMO2,QMO2,SDN,QDN,SUP,QUP,SEDGE,QEDGE dummy variables
!@var WMDN,WMUP,SVDN,SVUP,WMEDG,SVEDG,DDRUP dummy variables
!@var DMSE difference in moist static energy
!@var FPLUME fraction of convective plume
!@var DFP an iterative increment
!@var FMP2,FRAT1,FRAT2,SMN1,QMN1,SMN2,QMN2 dummy variables
!@var SMP, QMP plume's SM, QM
!@var TP plume's temperature (K)
!@var GAMA,DQSUM,TNX,DQ,CONSUM,BETAU,ALPHAU dummy variables
!@var DMSE1 difference in moist static energy
!@var FCTYPE fraction for convective cloud types
!@var CDHDRT,ALPHA,BETA,CDHM,CDHSUM,CLDREF dummy variables
!@var DDRAFT downdraft mass
!@var DELTA fraction of plume stays in the layer
!@var CLDM subsidence due to convection
!@var DQEVP amount of condensate that evaporates in downdrafts
!@var DQRAT fraction for the condensate to evaporate
!@var EPLUME air mass of entrainment
!@var ETADN fraction of the downdraft
!@var ETAL1 fractional entrainment rate
!@var EVPSUM,FCDH,FCDH1,FCLD,FCLOUD,FDDL,FDDP dummy variables
!@var FENTR fraction of entrainment
!@var FENTRA fraction of entrainment
!@var FEVAP fraction of air available for precip evaporation
!@var FLEFT fraction of plume after removing downdraft mass
!@var FQCOND fraction of water vapour that condenses in plume
!@var FQEVP fraction of water vapour that evaporates in downdraft
!@var FPRCP fraction of evaporated precipitation
!@var FSEVP fraction of energy lost to evaporate
!@var FSSUM fraction of energy lost to evaporate
!@var HEAT1 heating due to phase change
!@var PRHEAT energy of condensate
!@var PRCP precipipation
!@var QMDN,QMIX,SMDN,SMIX dummy variables
!@var QMPMAX,SMPMAX detrainment of QMP, SMP
!@var QMPT,SMPT dummy variables
!@var QNX,SUMAJ,SUMDP dummy variables
!@var QSATC saturation vapor mixing ratio
!@var QSATMP plume's saturation vapor mixing ratio
!@var RCLD,RCLDE cloud particle's radius, effective radius
!@var SLH LHX/SHA
!@var EDRAFT entrainment into downdrafts
!@var TOLD,TOLD1 old temperatures
!@var TEMWM,TEM,WTEM,WCONST dummy variables
!@var WORK work done on convective plume
LOGICAL MC1 !@var MC1 true for the first convective event
REAL*8, PARAMETER :: CK1 = 1. !@param CK1 a tunning const.
!@param RHOG,RHOIP density of graupel and ice particles
!@param ITMAX max iteration index
!@param CN0 constant use in computing FLAMW, etc
!@param PN tuning exponential for computing WV
REAL*8, PARAMETER :: CN0=8.d6,PN=1.d0
REAL*8, PARAMETER :: RHOG=400.,RHOIP=100.
INTEGER, PARAMETER :: ITMAX=50
!@var FLAMW,FLAMG,FLAMI lamda for water, graupel and ice, respectively
!@var WMAX specified maximum convective vertical velocity
!@var WV convetive vertical velocity
!@var VT precip terminal velocity
!@var DCW,DCG,DCI critical cloud particle sizes
!@var FG, FI fraction for graupel and ice
!@var FITMAX set to ITMAX
!@var CONDMU convective condensate in Kg/m^3
REAL*8 :: FLAMW,FLAMG,FLAMI,WMAX,WV,DCW,DCG,DCI,FG,FI,FITMAX,DDCW
* ,VT,CONDMU
INTEGER K,L,N !@var K,L,N loop variables
INTEGER ITYPE !@var convective cloud types
!@var IERR,LERR error reports from advection
INTEGER, INTENT(OUT) :: IERR,LERR
!@var DUM, DVM changes of UM,VM
REAL*8, DIMENSION(IM,LM) :: DUM,DVM
REAL*8 THBAR !@var THBAR virtual temperature at layer edge
!@var BELOW_CLOUD logical- is the current level below cloud?
LOGICAL BELOW_CLOUD
C****
C**** MOIST CONVECTION
C****
C**** CONVECTION OCCURS AT THE LOWEST MOIST CONVECTIVELY UNSTABLE
C**** LEVEL AND CONTINUES UNTIL A STABLE LAYER PAIR IS REACHED. RE-
C**** EVAPORATION AND PRECIPITATION ARE COMPUTED AT THE END OF THIS
C**** CYCLE. THE WHOLE PROCESS MAY BE REPEATED FROM A NEW LOWEST
C**** UNSTABLE LEVEL.
C****
ierr=0 ; lerr=0
LMCMIN=0
LMCMAX=0
MCCONT=0
FMC1=0.
FSSL=1.
FITMAX=ITMAX
C**** initiallise arrays of computed ouput
TAUMCL=0
SVWMXL=0
SVLATL=0
VSUBL=0
PRECNVL=0
CLDMCL=0
CLDSLWIJ=0
CLDDEPIJ=0
PRCPMC=0.
TPSAV=0
CSIZEL=RWCLDOX*10.*(1.-PEARTH)+10.*PEARTH ! droplet rad in stem
C**** zero out diagnostics
MCFLX =0.
DGDSM=0.
DPHASE=0.
DTOTW=0.
DQCOND=0.
DGDQM=0.
DDMFLX=0.
C**** save initial values (which will be updated after subsid)
SM1=SM
QM1=QM
FSUB=0.
FCONV=0.
FMCL=0.
VLAT=LHE
C**** SAVE ORIG PROFILES
SMOLD(:) = SM(:)
SMOMOLD(:,:) = SMOM(:,:)
QMOLD(:) = QM(:)
QMOMOLD(:,:) = QMOM(:,:)
C**** OUTER MC LOOP OVER BASE LAYER
DO 600 LMIN=1,LMCM-1
MAXLVL=0
MINLVL=LM
C****
C**** COMPUTE THE CONVECTIVE MASS OF THE NON-ENTRAINING PART
C****
TERM1=-10.*CK1*SDL(LMIN+1)*BYGRAV
FMP0=TERM1*XMASS
IF(FMP0.LE.0.) FMP0=0.
C**** CREATE A PLUME IN THE BOTTOM LAYER
C****
C**** ITERATION TO FIND FPLUME WHICH RESTORES THE ATM TO NEUTRAL STATE
C****
SMO1=SM(LMIN)
QMO1=QM(LMIN)
SMO2=SM(LMIN+1)
QMO2=QM(LMIN+1)
SDN=SMO1*BYAM(LMIN)
SUP=SMO2*BYAM(LMIN+1)
SEDGE=THBAR(SUP,SDN)
QDN=QMO1*BYAM(LMIN)
QUP=QMO2*BYAM(LMIN+1)
WMDN=WML(LMIN)
WMUP=WML(LMIN+1)
SVDN=SDN*(1.+DELTX*QDN-WMDN)
SVUP=SUP*(1.+DELTX*QUP-WMUP)
QEDGE=.5*(QUP+QDN)
WMEDG=.5*(WMUP+WMDN)
SVEDG=SEDGE*(1.+DELTX*QEDGE-WMEDG)
LHX=LHE
SLH=LHX*BYSHA
DMSE=(SVUP-SVEDG)*PLK(LMIN+1)+(SVEDG-SVDN)*PLK(LMIN)+
* SLHE*(QSAT(SUP*PLK(LMIN+1),LHX,PL(LMIN+1))-QDN)
IF(DMSE.GT.-1d-10) GO TO 600
C****
FPLUME=.25
DFP = .25
DO ITER=1,8
DFP=DFP*0.5
FMP2=FPLUME*AIRM(LMIN)
FRAT1=FMP2*BYAM(LMIN+1)
FRAT2=FMP2*BYAM(LMIN+2)
SMN1=SMO1*(1.-FPLUME)+FRAT1*SMO2
QMN1=QMO1*(1.-FPLUME)+FRAT1*QMO2
SMN2=SMO2*(1.-FRAT1)+FRAT2*SM(LMIN+2)
QMN2=QMO2*(1.-FRAT1)+FRAT2*QM(LMIN+2)
SMP=SMO1*FPLUME
QMP=QMO1*FPLUME
TP=SMO1*PLK(LMIN+1)*BYAM(LMIN)
QSATMP=FMP2*QSAT(TP,LHX,PL(LMIN+1))
GAMA=SLH*QSATMP*DQSATDT(TP,LHX)/FMP2
DQSUM=(QMP-QSATMP)/(1.+GAMA)
IF(DQSUM.LE.0.) GO TO 205
FEVAP=.5*FPLUME
MCLOUD=FEVAP*AIRM(LMIN+1)
TNX=SMO2*PLK(LMIN+1)*BYAM(LMIN+1)
QNX=QMO2*BYAM(LMIN+1)
QSATC=QSAT(TNX,LHX,PL(LMIN+1))
DQ=MCLOUD*(QSATC-QNX)/(1.+SLH*QSATC*DQSATDT(TNX,LHX))
IF(DQ.GT.DQSUM) DQ=DQSUM
SMN2=SMN2-SLH*DQ/PLK(LMIN+1)
QMN2=QMN2+DQ
DQSUM=DQSUM-DQ
IF(DQSUM.LE.0.) GO TO 205
MCLOUD=FEVAP*AIRM(LMIN)
TNX=SMO1*PLK(LMIN)*BYAM(LMIN)
QNX=QMO1*BYAM(LMIN)
QSATC=QSAT(TNX,LHX,PL(LMIN))
DQ=MCLOUD*(QSATC-QNX)/(1.+SLH*QSATC*DQSATDT(TNX,LHX)) ! new func'n
IF(DQ.GT.DQSUM) DQ=DQSUM
SMN1=SMN1-SLH*DQ/PLK(LMIN)
QMN1=QMN1+DQ
205 SDN=SMN1*BYAM(LMIN)
SUP=SMN2*BYAM(LMIN+1)
SEDGE=THBAR(SUP,SDN)
QDN=QMN1*BYAM(LMIN)
QUP=QMN2*BYAM(LMIN+1)
SVDN=SDN*(1.+DELTX*QDN-WMDN)
SVUP=SUP*(1.+DELTX*QUP-WMUP)
QEDGE=.5*(QUP+QDN)
SVEDG=SEDGE*(1.+DELTX*QEDGE-WMEDG)
DMSE1=(SVUP-SVEDG)*PLK(LMIN+1)+(SVEDG-SVDN)*PLK(LMIN)+
* SLHE*(QSAT(SUP*PLK(LMIN+1),LHX,PL(LMIN+1))-QDN)
IF (ABS(DMSE1).LE.1.d-3) GO TO 411
IF(DMSE1.GT.1.d-3) FPLUME=FPLUME-DFP
IF(DMSE1.LT.-1.d-3) FPLUME=FPLUME+DFP
END DO
411 IF(FPLUME.LE..001) GO TO 600
C****
C**** ITERATION THROUGH CLOUD TYPES
C****
ITYPE=2 ! always 2 types of clouds:
C IF(LMIN.LE.2) ITYPE=2 ! entraining & non-entraining
FCTYPE=1.
DO 570 IC=1,ITYPE
C**** INITIALLISE VARIABLES USED FOR EACH TYPE
DO L=1,LM
COND(L)=0.
CDHEAT(L)=0.
CONDP(L)=0.
DM(L)=0.
DMR(L)=0.
CCM(L)=0.
DDM(L)=0.
END DO
DUM(1:KMAX,:)=0.
DVM(1:KMAX,:)=0.
DSM(:) = 0.
DSMOM(:,:) = 0.
DSMR(:) = 0.
DSMOMR(:,:) = 0.
DQM(:) = 0.
DQMOM(:,:) = 0.
DQMR(:) = 0.
DQMOMR(:,:) = 0.
MC1=.FALSE.
IF (IC.EQ.1) THEN
WMAX=2.
IF(PLAND.GE..5) WMAX=5.
ELSE
WMAX=1.
IF(PLAND.GE..5) WMAX=2.5
ENDIF
LHX=LHE
MPLUME=MIN(AIRM(LMIN),AIRM(LMIN+1))
IF(MPLUME.GT.FMP2) MPLUME=FMP2
IF(ITYPE.EQ.2) THEN
FCTYPE=1.
IF(MPLUME.GT.FMP0) FCTYPE=FMP0/MPLUME
IF(IC.EQ.2) FCTYPE=1.-FCTYPE
IF(FCTYPE.LT.0.001) GO TO 570
END IF
MPLUM1=MPLUME
MPLUME=MPLUME*FCTYPE
C**** calculate subsiding fraction here. Then we can use FMC1 from the
C**** beginning. The analagous arrays are only set if plume is actually
C**** moved up.
IF (MCCONT.le.0) THEN
FCONV_tmp=MPLUM1/AIRM(LMIN+1)
IF(FCONV_tmp.GT.1.d0) FCONV_tmp=1.d0
FSUB_tmp=1.d0+(AIRM(LMIN+1)-100.d0)/200.d0
IF(FSUB_tmp.GT.1.d0/(FCONV_tmp+1.d-20)-1.d0)
* FSUB_tmp=1.d0/(FCONV_tmp+1.d-20)-1.d0
IF(FSUB_tmp.LT.1.d0) FSUB_tmp=1.d0
IF(FSUB_tmp.GT.5.d0) FSUB_tmp=5.d0
FSSL_tmp=1.d0-(1.d0+FSUB_tmp)*FCONV_tmp
IF(FSSL_tmp.LT.CLDMIN) FSSL_tmp=CLDMIN
IF(FSSL_tmp.GT.1.d0-CLDMIN) FSSL_tmp=1.d0-CLDMIN
FMC1=1.d0-FSSL_tmp+teeny
ELSE
C**** guard against possibility of too big a plume
MPLUME=MIN(0.95d0*AIRM(LMIN)*FMC1,MPLUME)
END IF
C**** adjust MPLUME to take account of restricted area of subsidence
C**** (i.e. MPLUME is now a greater fraction of the relevant airmass.
MPLUME=MIN( MPLUME/FMC1,
* AIRM(LMIN)*0.95d0*QM(LMIN)/(QMOLD(LMIN) + teeny) )
IF(MPLUME.LE..001*AIRM(LMIN)) GO TO 570
FPLUME=MPLUME*BYAM(LMIN)
SMP = SMOLD(LMIN)*FPLUME
SMOMP(xymoms)=SMOMOLD(xymoms,LMIN)*FPLUME
QMP = QMOLD(LMIN)*FPLUME
QMOMP(xymoms)=QMOMOLD(xymoms,LMIN)*FPLUME
TPSAV(LMIN)=SMP*PLK(LMIN)/MPLUME
DMR(LMIN)=-MPLUME
DSMR(LMIN)=-SMP
DSMOMR(xymoms,LMIN)=-SMOMP(xymoms)
DSMOMR(zmoms,LMIN)=-SMOMOLD(zmoms,LMIN)*FPLUME
DQMR(LMIN)=-QMP
DQMOMR(xymoms,LMIN)=-QMOMP(xymoms)
DQMOMR(zmoms,LMIN)=-QMOMOLD(zmoms,LMIN)*FPLUME
DO K=1,KMAX
UMP(K)=UM(K,LMIN)*FPLUME
DUM(K,LMIN)=-UMP(K)
VMP(K)=VM(K,LMIN)*FPLUME
DVM(K,LMIN)=-VMP(K)
ENDDO
C****
C**** RAISE THE PLUME TO THE TOP OF CONVECTION AND CALCULATE
C**** ENTRAINMENT, CONDENSATION, AND SECONDARY MIXING
C****
CDHSUM=0.
CDHDRT=0.
ETADN=0.
LDRAFT=LM
EVPSUM=0.
DDRAFT=0.
LFRZ=0
LMAX=LMIN
220 L=LMAX+1
C**** TEST FOR SUFFICIENT AIR, MOIST STATIC STABILITY AND ENERGY
C IF(L.GT.LMIN+1.AND.SDL(L).GT.0.) GO TO 340
IF(MPLUME.LE..001*AIRM(L)) GO TO 340
SDN=SMP/MPLUME
SUP=SM1(L)*BYAM(L)
QDN=QMP/MPLUME
QUP=QM1(L)*BYAM(L)
WMDN=0.
WMUP=WML(L)
SVDN=SDN*(1.+DELTX*QDN-WMDN)
SVUP=SUP*(1.+DELTX*QUP-WMUP)
IF(LMAX.GT.LMIN) THEN
SEDGE=THBAR(SUP,SDN)
QEDGE=.5*(QUP+QDN)
WMEDG=.5*(WMUP+WMDN)
SVEDG=SEDGE*(1.+DELTX*QEDGE-WMEDG)
LHX=LHE
DMSE=(SVUP-SVEDG)*PLK(L)+(SVEDG-SVDN)*PLK(L-1)+
* SLHE*(QSAT(SUP*PLK(L),LHX,PL(L))-QDN)
IF(DMSE.GT.-1d-10) GO TO 340
END IF
IF(PLK(L-1)*(SVUP-SVDN)+SLHE*(QUP-QDN).GE.0.) GO TO 340
C**** TEST FOR CONDENSATION ALSO DETERMINES IF PLUME REACHES UPPER LAYER
TP=SMP*PLK(L)/MPLUME
TPSAV(L)=TP
IF(TPSAV(L-1).GE.TF.AND.TPSAV(L).LT.TF) LFRZ=L-1
IF(TP.LT.TI) LHX=LHS
QSATMP=MPLUME*QSAT(TP,LHX,PL(L))
IF(QMP.LT.QSATMP) GO TO 340
IF(TP.GE.TF.OR.LHX.EQ.LHS) GO TO 290
LHX=LHS
QSATMP=MPLUME*QSAT(TP,LHX,PL(L))
290 CONTINUE
C**** DEFINE VLAT TO AVOID PHASE DISCREPANCY BETWEEN TWO PLUMES
IF (VLAT(L).EQ.LHS) LHX=LHS
VLAT(L)=LHX
SLH=LHX*BYSHA
MCCONT=MCCONT+1
IF(MCCONT.EQ.1) MC1=.TRUE.
C IF(MC1.AND.L.EQ.LMIN+1) THEN
C FCONV(L)=FCONV_tmp ! these are set here but do not make
C FSSL(L)=FSSL_tmp ! much sense at the moment...
C FSUB(L)=FSUB_tmp
C ENDIF
C****
C**** DEPOSIT PART OF THE PLUME IN LOWER LAYER
C****
IF(MPLUME.GT..95*AIRM(L)) THEN
DELTA=(MPLUME-.95*AIRM(L))/MPLUME
DM(L-1)=DM(L-1)+DELTA*MPLUME
MPLUME=.95*AIRM(L)
DSM(L-1)= DSM(L-1)+DELTA*SMP
SMP = SMP *(1.-DELTA)
DSMOM(xymoms,L-1)=DSMOM(xymoms,L-1)+DELTA*SMOMP(xymoms)
SMOMP(xymoms) = SMOMP(xymoms)*(1.-DELTA)
DQM(L-1)= DQM(L-1)+DELTA*QMP
QMP = QMP *(1.-DELTA)
DQMOM(xymoms,L-1)=DQMOM(xymoms,L-1)+DELTA*QMOMP(xymoms)
QMOMP(xymoms) = QMOMP(xymoms)*(1.-DELTA)
DO K=1,KMAX
DUM(K,L-1)=DUM(K,L-1)+UMP(K)*DELTA
DVM(K,L-1)=DVM(K,L-1)+VMP(K)*DELTA
UMP(K)=UMP(K)-UMP(K)*DELTA
VMP(K)=VMP(K)-VMP(K)*DELTA
ENDDO
END IF
C****
C**** CONVECTION IN UPPER LAYER (WORK DONE COOLS THE PLUME)
C****
WORK=MPLUME*(SUP-SDN)*(PLK(L-1)-PLK(L))/PLK(L-1)
C SMP=SMP-WORK
DSM(L-1)=DSM(L-1)-WORK
CCM(L-1)=MPLUME
WV=WMAX
IF(PL(L).GT.700..AND.PLE(1).GT.700.)
* WV=WMAX*(PLE(1)-PL(L))/(PLE(1)-700.)
IF(PL(L).LT.400.) WV=WMAX*((PL(L)-PLE(LM+1))/(400.-PLE(LM+1)))**PN
DCG=((WV/19.3)*(PL(L)/1000.)**.4)**2.7
DCG=MIN(DCG,1.D-2)
DCI=((WV/1.139)*(PL(L)/1000.)**.4)**9.09
DCI=MIN(DCI,1.D-2)
C****
C**** ENTRAINMENT
C****
IF(IC.EQ.2.OR.(IC.EQ.1.AND.PL(L).GE.800.)) THEN
FENTR=ETAL(L)*FPLUME
IF(FENTR+FPLUME.GT.1.) FENTR=1.-FPLUME
IF(FENTR.LT.teeny) GO TO 293
ETAL1=FENTR/(FPLUME+teeny)
FPLUME=FPLUME+FENTR
EPLUME=MPLUME*ETAL1
C**** Reduce EPLUME so that mass flux is less than mass in box
IF (EPLUME.GT.AIRM(L)*0.975d0-MPLUME) THEN
EPLUME=AIRM(L)*0.975d0-MPLUME
END IF
MPLUME=MPLUME+EPLUME
FENTRA = EPLUME*BYAM(L)
DSMR(L)=DSMR(L)-EPLUME*SUP ! = DSM(L)-SM(L)*FENTRA
DSMOMR(:,L)=DSMOMR(:,L)-SMOM(:,L)*FENTRA
DQMR(L)=DQMR(L)-EPLUME*QUP ! = DQM(L)-QM(L)*FENTRA
DQMOMR(:,L)=DQMOMR(:,L)-QMOM(:,L)*FENTRA
DMR(L)=DMR(L)-EPLUME
SMP=SMP+EPLUME*SUP
SMOMP(xymoms)= SMOMP(xymoms)+ SMOM(xymoms,L)*FENTRA
QMP=QMP+EPLUME*QUP
QMOMP(xymoms)= QMOMP(xymoms)+ QMOM(xymoms,L)*FENTRA
DO K=1,KMAX
UMP(K)=UMP(K)+U_0(K,L)*EPLUME
DUM(K,L)=DUM(K,L)-U_0(K,L)*EPLUME
VMP(K)=VMP(K)+V_0(K,L)*EPLUME
DVM(K,L)=DVM(K,L)-V_0(K,L)*EPLUME
ENDDO
293 CONTINUE
END IF
C****
C**** CHECK THE DOWNDRAFT POSSIBILITY
C****
IF(L-LMIN.LE.1) GO TO 291
IF(ETADN.GT.1d-10) GO TO 291
SMIX=.5*(SUP+SMP/MPLUME)
QMIX=.5*(QUP+QMP/MPLUME)
C WMIX=.5*(WMUP+WMDN)
C SVMIX=SMIX*(1.+DELTX*QMIX-WMIX)
C DMMIX=(SVUP-SVMIX)*PLK(L)+
C * SLHE*(QSAT(SUP*PLK(L),LHX,PL(L))-QMIX)
C IF(DMMIX.LT.1d-10) GO TO 291
IF(SMIX.GE.SUP) GO TO 291
IF(PL(L).GT.700.) GO TO 291
LDRAFT=L
ETADN=BY3
FLEFT=1.-.5*ETADN
DDRAFT=ETADN*MPLUME
FDDP = .5*DDRAFT/MPLUME
FDDL = .5*DDRAFT*BYAM(L)
MPLUME=FLEFT*MPLUME
SMDN=DDRAFT*SMIX ! = SM(L)*FDDL + SMP*FDDP
SMOMDN(xymoms)=SMOM(xymoms,L)*FDDL + SMOMP(xymoms)*FDDP
SMP=FLEFT*SMP
SMOMP(xymoms)= SMOMP(xymoms)*FLEFT
QMDN=DDRAFT*QMIX ! = QM(L)*FDDL + QMP*FDDP
QMOMDN(xymoms)=QMOM(xymoms,L)*FDDL + QMOMP(xymoms)*FDDP
QMP=FLEFT*QMP
QMOMP(xymoms)= QMOMP(xymoms)*FLEFT
DMR(L) = DMR(L)-.5*DDRAFT
DSMR(L)=DSMR(L)-.5*DDRAFT*SUP ! = DSM(L)-SM(L)*FDDL
DSMOMR(:,L)=DSMOMR(:,L) - SMOM(:,L)*FDDL
DQMR(L)=DQMR(L)-.5*DDRAFT*QUP ! = DQM(L)-QM(L)*FDDL
DQMOMR(:,L)=DQMOMR(:,L) - QMOM(:,L)*FDDL
DO K=1,KMAX
UMDN(K)=.5*(ETADN*UMP(K)+DDRAFT*U_0(K,L))
UMP(K)=UMP(K)*FLEFT
DUM(K,L)=DUM(K,L)-.5*DDRAFT*U_0(K,L)
VMDN(K)=.5*(ETADN*VMP(K)+DDRAFT*V_0(K,L))
VMP(K)=VMP(K)*FLEFT
DVM(K,L)=DVM(K,L)-.5*DDRAFT*V_0(K,L)
ENDDO
C****
C**** CONDENSE VAPOR IN THE PLUME AND ADD LATENT HEAT
C****
291 DQSUM=0.
SMPT=SMP
QMPT=QMP
DO 292 N=1,3
TP=SMP*PLK(L)/MPLUME
QSATMP=MPLUME*QSAT(TP,LHX,PL(L))
GAMA=SLH*QSATMP*DQSATDT(TP,LHX)/MPLUME
DQ=(QMP-QSATMP)/(1.+GAMA)
SMP=SMP+SLH*DQ/PLK(L)
QMP=QMP-DQ
292 DQSUM=DQSUM+DQ
FQCOND = 0
IF(DQSUM.GE.0.) THEN
IF (QMPT.gt.teeny) FQCOND = DQSUM/QMPT
QMOMP(xymoms) = QMOMP(xymoms)*(1.-FQCOND)
ELSE ! no change
DQSUM=0.
SMP=SMPT
QMP=QMPT
END IF
COND(L)=DQSUM
CONDMU=100.*COND(L)*PL(L)/(CCM(L-1)*TL(L)*RGAS)
FLAMW=(1000.d0*PI*CN0/(CONDMU+teeny))**.25
FLAMG=(400.d0*PI*CN0/(CONDMU+teeny))**.25
FLAMI=(100.d0*PI*CN0/(CONDMU+teeny))**.25
IF (TP.GE.TF) THEN ! water phase
DDCW=6d-3/FITMAX
IF(PLAND.LT..5) DDCW=1.5d-3/FITMAX
DCW=0.
DO ITER=1,ITMAX-1
VT=(-.267d0+DCW*(5.15D3-DCW*(1.0225D6-7.55D7*DCW)))*
* (1000./PL(L))**.4d0
IF(VT.GE.0..AND.ABS(VT-WV).LT..3) EXIT
IF(VT.GT.WMAX) EXIT
DCW=DCW+DDCW
END DO
CONDP(L)=RHOW*(PI*by6)*CN0*EXP(-FLAMW*DCW)*
* (DCW*DCW*DCW/FLAMW+3.*DCW*DCW/(FLAMW*FLAMW)+
* 6.*DCW/(FLAMW*FLAMW*FLAMW)+6./FLAMW**4)
ELSE IF (TP.LE.TI) THEN ! pure ice phase
CONDP(L)=RHOIP*(PI*by6)*CN0*EXP(-FLAMI*DCI)*
* (DCI*DCI*DCI/FLAMI+3.*DCI*DCI/(FLAMI*FLAMI)+
* 6.*DCI/(FLAMI*FLAMI*FLAMI)+6./FLAMI**4)
ELSE ! mixed phase
FG=(TP-TF+40.)*0.025d0
FI=1.-FG
CONDIP=RHOIP*(PI*by6)*CN0*EXP(-FLAMI*DCI)*
* (DCI*DCI*DCI/FLAMI+3.*DCI*DCI/(FLAMI*FLAMI)+
* 6.*DCI/(FLAMI*FLAMI*FLAMI)+6./FLAMI**4)
CONDGP=RHOG*(PI*by6)*CN0*EXP(-FLAMG*DCG)*
* (DCG*DCG*DCG/FLAMG+3.*DCG*DCG/(FLAMG*FLAMG)+
* 6.*DCG/(FLAMG*FLAMG*FLAMG)+6./FLAMG**4)
CONDP(L)=FG*CONDGP+FI*CONDIP
ENDIF
c convert condp to the same units as cond
CONDP(L)=.01d0*CONDP(L)*CCM(L-1)*TL(L)*RGAS/PL(L)
TAUMCL(L)=TAUMCL(L)+DQSUM*FMC1
CDHEAT(L)=SLH*COND(L)
CDHSUM=CDHSUM+CDHEAT(L)
IF(ETADN.GT.1d-10) CDHDRT=CDHDRT+SLH*COND(L)
C****
C**** UPDATE ALL QUANTITIES CARRIED BY THE PLUME
C****
C MCCONT=MCCONT+1
C IF(MCCONT.EQ.1) MC1=.TRUE.
C IF(MC1.AND.PLE(LMIN)-PLE(L+2).GE.450.) SVLATL(L)=LHX
SVLATL(L)=VLAT(L)
SMPMAX=SMP
SMOMPMAX(xymoms) = SMOMP(xymoms)
QMPMAX=QMP
QMOMPMAX(xymoms) = QMOMP(xymoms)
MPMAX=MPLUME
LMAX = LMAX + 1
IF (LMAX.LT.LM) GO TO 220 ! CHECK FOR NEXT POSSIBLE LMAX
C**** UPDATE CHANGES CARRIED BY THE PLUME IN THE TOP CLOUD LAYER
340 IF(LMIN.EQ.LMAX) GO TO 600
IF(TPSAV(LMAX).GE.TF) LFRZ=LMAX
DM(LMAX)=DM(LMAX)+MPMAX
DSM(LMAX)=DSM(LMAX)+SMPMAX
DSMOM(xymoms,LMAX)=DSMOM(xymoms,LMAX) + SMOMPMAX(xymoms)
DQM(LMAX)=DQM(LMAX)+QMPMAX
DQMOM(xymoms,LMAX)=DQMOM(xymoms,LMAX) + QMOMPMAX(xymoms)
CCM(LMAX)=0.
DO K=1,KMAX
DUM(K,LMAX)=DUM(K,LMAX)+UMP(K)
DVM(K,LMAX)=DVM(K,LMAX)+VMP(K)
ENDDO
CDHM=0.
IF(MINLVL.GT.LMIN) MINLVL=LMIN
IF(MAXLVL.LT.LMAX) MAXLVL=LMAX
IF(LMCMIN.EQ.0) LMCMIN=LMIN
IF(LMCMAX.LT.MAXLVL) LMCMAX=MAXLVL
C****
C**** PROCESS OF DOWNDRAFTING
C****
LDMIN=LMIN
EDRAFT=0.
IF(ETADN.GT.1d-10) THEN
CONSUM=0.
DO 347 L=LMIN,LDRAFT
347 CONSUM=CONSUM+COND(L)
TNX=SMDN*PLK(LMIN)/DDRAFT
QNX=QMDN/DDRAFT
LHX=LHE
IF(TPSAV(LMIN).LT.TF) LHX=LHS
SLH=LHX*BYSHA
QSATC=QSAT(TNX,LHX,PL(LMIN))
DQ=(QSATC-QNX)/(1.+SLH*QSATC*DQSATDT(TNX,LHX))
DQRAT=DQ*DDRAFT/(CONSUM+teeny)
DO 346 L=LDRAFT,1,-1
LHX=LHE
IF (L.GE.LMIN) THEN ! avoids reference to uninitiallised value
IF(TPSAV(L).LT.TF) LHX=LHS
END IF
SLH=LHX*BYSHA
DQEVP=DQRAT*COND(L)
IF(DQEVP.GT.COND(L)) DQEVP=COND(L)
IF (L.LT.LMIN) DQEVP=0.
FSEVP = 0
IF (ABS(PLK(L)*SMDN).gt.teeny) FSEVP = SLH*DQEVP/(PLK(L)*SMDN)
SMDN=SMDN-SLH*DQEVP/PLK(L)
SMOMDN(xymoms)=SMOMDN(xymoms)*(1.-FSEVP)
FQEVP = 0
IF (COND(L).gt.0.) FQEVP = DQEVP/COND(L)
QMDN=QMDN+DQEVP
COND(L)=COND(L)-DQEVP
TAUMCL(L)=TAUMCL(L)-DQEVP*FMC1
CDHEAT(L)=CDHEAT(L)-DQEVP*SLH
EVPSUM=EVPSUM+DQEVP*SLH
C**** ENTRAINMENT OF DOWNDRAFTS
IF(L.LT.LDRAFT.AND.L.GT.1) THEN
DDRUP=DDRAFT
DDRAFT=DDRAFT*(1.+ETAL(L))
IF(DDRUP.GT.DDRAFT) DDRAFT=DDRUP
IF(DDRAFT.GT..95d0*(AIRM(L-1)+DMR(L-1)))
* DDRAFT=.95d0*(AIRM(L-1)+DMR(L-1))
EDRAFT=DDRAFT-DDRUP
IF (EDRAFT.gt.0) THEN ! usual case, entrainment into downdraft
FENTRA=EDRAFT*BYAM(L)
SENV=SM(L)/AIRM(L)
QENV=QM(L)/AIRM(L)
SMDN=SMDN+EDRAFT*SENV
QMDN=QMDN+EDRAFT*QENV
SMOMDN(xymoms)= SMOMDN(xymoms)+ SMOM(xymoms,L)*FENTRA
QMOMDN(xymoms)= QMOMDN(xymoms)+ QMOM(xymoms,L)*FENTRA
DSMR(L)=DSMR(L)-EDRAFT*SENV
DSMOMR(:,L)=DSMOMR(:,L)-SMOM(:,L)*FENTRA
DQMR(L)=DQMR(L)-EDRAFT*QENV
DQMOMR(:,L)=DQMOMR(:,L)-QMOM(:,L)*FENTRA
DMR(L)=DMR(L)-EDRAFT
ELSE ! occasionally detrain into environment if ddraft too big
FENTRA=EDRAFT/DDRUP ! < 0
DSM(L)=DSM(L)-FENTRA*SMDN
DSMOM(xymoms,L)=DSMOM(xymoms,L)-SMOMDN(xymoms)*FENTRA
DQM(L)=DQM(L)-FENTRA*QMDN
DQMOM(xymoms,L)=DQMOM(xymoms,L)-QMOMDN(xymoms)*FENTRA
SMDN=SMDN*(1+FENTRA)
QMDN=QMDN*(1+FENTRA)
SMOMDN(xymoms)= SMOMDN(xymoms)*(1+FENTRA)
QMOMDN(xymoms)= QMOMDN(xymoms)*(1+FENTRA)
DM(L)=DM(L)-EDRAFT
END IF
ENDIF
IF(L.GT.1) DDM(L-1)=DDRAFT
LDMIN=L
C**** ALLOW FOR DOWNDRAFT TO DROP BELOW LMIN, IF IT'S NEGATIVE BUOYANT
IF (L.LE.LMIN.AND.L.GT.1) THEN
SMIX=SMDN/DDRAFT
IF (SMIX.GE.(SM1(L-1)/AIRM(L-1))) GO TO 345
ENDIF
346 CONTINUE
345 DSM(LDMIN)=DSM(LDMIN)+SMDN
DSMOM(xymoms,LDMIN)=DSMOM(xymoms,LDMIN) + SMOMDN(xymoms)
DQM(LDMIN)=DQM(LDMIN)+QMDN
DQMOM(xymoms,LDMIN)=DQMOM(xymoms,LDMIN) + QMOMDN(xymoms)
DO K=1,KMAX
DUM(K,LDMIN)=DUM(K,LDMIN)+UMDN(K)
DVM(K,LDMIN)=DVM(K,LDMIN)+VMDN(K)
ENDDO
DM(LDMIN)=DM(LDMIN)+DDRAFT
END IF
C****
C**** SUBSIDENCE AND MIXING
C****
C**** Calculate vertical mass fluxes (Note CM for subsidence is defined
C**** in opposite sense than normal (positive is down))
DO L=0,LDMIN-1
CM(L) = 0.
END DO
DO L=LDMIN,LMAX
CM(L) = CM(L-1) - DM(L) - DMR(L)
SMT(L)=SM(L) ! Save profiles for diagnostics
QMT(L)=QM(L)
END DO
DO L=LMAX+1,LM
CM(L) = 0.
END DO
C**** simple upwind scheme for momentum
ALPHA=0.
DO 380 L=LDMIN,LMAX
CLDM=CCM(L)
IF(L.LT.LDRAFT.AND.ETADN.GT.1d-10) CLDM=CCM(L)-DDM(L)
IF(MC1) VSUBL(L)=100.*CLDM*RGAS*TL(L)/(PL(L)*GRAV*DTsrc)
BETA=CLDM*BYAM(L+1)
IF(CLDM.LT.0.) BETA=CLDM*BYAM(L)
BETAU=BETA
ALPHAU=ALPHA
IF(BETA.LT.0.) BETAU=0.
IF(ALPHA.LT.0.) ALPHAU=0.
DO K=1,KMAX
UM(K,L)=
* UM(K,L)+RA(K)*(-ALPHAU*UM(K,L)+BETAU*UM(K,L+1)+DUM(K,L))
VM(K,L)=
* VM(K,L)+RA(K)*(-ALPHAU*VM(K,L)+BETAU*VM(K,L+1)+DVM(K,L))
ENDDO
380 ALPHA=BETA
C**** Subsidence uses Quadratic Upstream Scheme for QM and SM
ML(LDMIN:LMAX) = AIRM(LDMIN:LMAX) + DMR(LDMIN:LMAX)
SM(LDMIN:LMAX) = SM(LDMIN:LMAX) + DSMR(LDMIN:LMAX)
SMOM(:,LDMIN:LMAX) = SMOM(:,LDMIN:LMAX) + DSMOMR(:,LDMIN:LMAX)
C****
nsub = lmax-ldmin+1
cmneg(ldmin:lmax)=-cm(ldmin:lmax)
cmneg(lmax)=0. ! avoid roundoff error (esp. for qlimit)
call adv1d(sm(ldmin),smom(1,ldmin), f(ldmin),fmom(1,ldmin),
& ml(ldmin),cmneg(ldmin), nsub,.false.,1, zdir,ierrt,lerrt)
SM(LDMIN:LMAX) = SM(LDMIN:LMAX) + DSM(LDMIN:LMAX)
SMOM(:,LDMIN:LMAX) = SMOM(:,LDMIN:LMAX) + DSMOM(:,LDMIN:LMAX)
ierr=max(ierrt,ierr) ; lerr=max(lerrt+ldmin-1,lerr)
C****
ML(LDMIN:LMAX) = AIRM(LDMIN:LMAX) + DMR(LDMIN:LMAX)
QM(LDMIN:LMAX) = QM(LDMIN:LMAX) + DQMR(LDMIN:LMAX)
QMOM(:,LDMIN:LMAX) = QMOM(:,LDMIN:LMAX) + DQMOMR(:,LDMIN:LMAX)
call adv1d(qm(ldmin),qmom(1,ldmin), f(ldmin),fmom(1,ldmin),
& ml(ldmin),cmneg(ldmin), nsub,.true.,1, zdir,ierrt,lerrt)
QM(LDMIN:LMAX) = QM(LDMIN:LMAX) + DQM(LDMIN:LMAX)
QMOM(:,LDMIN:LMAX) = QMOM(:,LDMIN:LMAX) + DQMOM(:,LDMIN:LMAX)
ierr=max(ierrt,ierr) ; lerr=max(lerrt+ldmin-1,lerr)
C**** diagnostics
DO L=LDMIN,LMAX
FCDH=0.
IF(L.EQ.LMAX) FCDH=CDHSUM-(CDHSUM-CDHDRT)*.5*ETADN+CDHM
FCDH1=0.
IF(L.EQ.LDMIN) FCDH1=(CDHSUM-CDHDRT)*.5*ETADN-EVPSUM
MCFLX(L)=MCFLX(L)+CCM(L)*FMC1
DGDSM(L)=DGDSM(L)+(PLK(L)*(SM(L)-SMT(L))-FCDH-FCDH1)*FMC1
DTOTW(L)=DTOTW(L)+SLHE*(QM(L)-QMT(L)+COND(L))*FMC1
DGDQM(L)=DGDQM(L)+SLHE*(QM(L)-QMT(L))*FMC1
DDMFLX(L)=DDMFLX(L)+DDM(L)*FMC1
IF(QM(L).LT.0.d0) WRITE(0,*) ' Q neg: it,i,j,l,q',
* itime,i_debug,j_debug,l,qm(l)
END DO
C**** save new 'environment' profile for static stability calc.
DO L=1,LM
SM1(L)=SM(L)
QM1(L)=QM(L)
END DO
C****
C**** Partition condensate into precipitation and cloud water
C****
IF(PLE(LMIN)-PLE(LMAX+1).GE.450.) THEN
DO L=LMAX,LMIN,-1
IF(COND(L).LT.CONDP(L)) CONDP(L)=COND(L)
FCLW=0.
IF (COND(L).GT.0) FCLW=(COND(L)-CONDP(L))/COND(L)
SVWMXL(L)=SVWMXL(L)+FCLW*COND(L)*BYAM(L)*FMC1
COND(L)=CONDP(L)
END DO
END IF
C****
C**** REEVAPORATION AND PRECIPITATION
C****
PRCP=COND(LMAX)
PRHEAT=CDHEAT(LMAX)
DPHASE(LMAX)=DPHASE(LMAX)+(CDHSUM-(CDHSUM-CDHDRT)*.5*ETADN+
* CDHM)*FMC1
DO 540 L=LMAX-1,1,-1
IF(PRCP.LE.0.) GO TO 530
FCLOUD=CCMUL*CCM(L)*BYAM(L+1)
IF(PLE(LMIN)-PLE(L+2).GE.450.) FCLOUD=CCMUL1*CCM(L)*BYAM(L+1)
IF(L.LT.LMIN) FCLOUD=CCMUL*CCM(LMIN)*BYAM(LMIN+1)
IF(PLE(LMIN)-PLE(LMAX+1).LT.450.) THEN
IF(L.EQ.LMAX-1) FCLOUD=CCMUL2*CCM(L)*BYAM(L+1)
IF(L.LT.LMIN) FCLOUD=0.
ENDIF
IF(FCLOUD.GT.1.) FCLOUD=1.
FEVAP=.5*CCM(L)*BYAM(L+1)
IF(L.LT.LMIN) FEVAP=.5*CCM(LMIN)*BYAM(LMIN+1)
IF(FEVAP.GT..5) FEVAP=.5
CLDMCL(L+1)=CLDMCL(L+1)+FCLOUD*FMC1
CLDREF=CLDMCL(L+1)
IF(PLE(LMAX+1).GT.700..AND.CLDREF.GT.CLDSLWIJ)
* CLDSLWIJ=CLDREF
IF(PLE(LMIN)-PLE(LMAX+1).GE.450..AND.CLDREF.GT.CLDDEPIJ)
* CLDDEPIJ=CLDREF
C**** FORWARD STEP COMPUTES HUMIDITY CHANGE BY RECONDENSATION
C**** Q = Q + F(TOLD,PRHEAT,QOLD+EVAP)
PRECNVL(L+1)=PRECNVL(L+1)+PRCP*BYGRAV
MCLOUD=0.
IF(L.LE.LMIN) MCLOUD=2.*FEVAP*AIRM(L)
TOLD=SMOLD(L)*PLK(L)*BYAM(L)
TOLD1=SMOLD(L+1)*PLK(L+1)*BYAM(L+1)
C**** decide whether to melt snow/ice
HEAT1=0.
IF(L.EQ.LFRZ.OR.(L.LE.LMIN.AND.TOLD.GE.TF.AND.TOLD1.LT.TF.AND.L.GT
* .LFRZ)) HEAT1=LHM*PRCP*BYSHA
C**** and deal with possible inversions and re-freezing of rain
IF (TOLD.LT.TF.AND.TOLD1.GT.TF) HEAT1=-LHM*PRCP*BYSHA
TNX=TOLD
QNX=QMOLD(L)*BYAM(L)
LHX=LHE
IF(TNX.LT.TF) LHX=LHS
IF(L.GT.LMIN) THEN ! needed to avoid unintiallised value
IF (TPSAV(L).GE.TF) LHX=LHE
END IF
SLH=LHX*BYSHA
DQSUM=0.
DO 510 N=1,3
QSATC=QSAT(TNX,LHX,PL(L))
DQ=(QSATC-QNX)/(1.+SLH*QSATC*DQSATDT(TNX,LHX))
TNX=TNX-SLH*DQ
QNX=QNX+DQ
510 DQSUM=DQSUM+DQ*MCLOUD
IF(DQSUM.LT.0.) DQSUM=0.
IF(DQSUM.GT.PRCP) DQSUM=PRCP
FPRCP=DQSUM/PRCP
PRCP=PRCP-DQSUM
C**** UPDATE TEMPERATURE AND HUMIDITY DUE TO NET REEVAPORATION IN CLOUDS
FSSUM = 0
IF (ABS(PLK(L)*SM(L)).gt.teeny) FSSUM = (SLH*DQSUM+HEAT1)/
* (PLK(L)*SM(L))
SM(L)=SM(L)-(SLH*DQSUM+HEAT1)/PLK(L)
SMOM(:,L) = SMOM(:,L)*(1.-FSSUM)
QM(L)=QM(L)+DQSUM
FCDH1=0.
IF(L.EQ.LDMIN) FCDH1=(CDHSUM-CDHDRT)*.5*ETADN-EVPSUM
DPHASE(L)=DPHASE(L)-(SLH*DQSUM-FCDH1+HEAT1)*FMC1
DQCOND(L)=DQCOND(L)-SLH*DQSUM*FMC1
C**** ADD PRECIPITATION AND LATENT HEAT BELOW
530 PRHEAT=CDHEAT(L)+SLH*PRCP
PRCP=PRCP+COND(L)
540 CONTINUE
C****
IF(PRCP.GT.0.) CLDMCL(1)=CLDMCL(1)+CCM(LMIN)*BYAM(LMIN+1)*FMC1
PRCPMC=PRCPMC+PRCP*FMC1
IF(LMCMIN.GT.LDMIN) LMCMIN=LDMIN
C****
C**** END OF LOOP OVER CLOUD TYPES
C****
MC1=.FALSE. !!!
570 CONTINUE
C****
C**** END OF OUTER LOOP OVER CLOUD BASE
C****
600 CONTINUE
IF(LMCMIN.GT.0) THEN
C**** set fssl array
DO L=1,LM
FSSL(L)=1.-FMC1
C FMCL(L)=FMC1 ! may be generalized
END DO
C**** ADJUSTMENT TO CONSERVE CP*T
SUMAJ=0.
SUMDP=0.
DO L=LMCMIN,LMCMAX
SUMDP=SUMDP+AIRM(L)*FMC1
SUMAJ=SUMAJ+DGDSM(L)
END DO
DO L=LMCMIN,LMCMAX
DGDSM(L)=DGDSM(L)-SUMAJ*AIRM(L)*FMC1/SUMDP
SM(L)=SM(L)-SUMAJ*AIRM(L)/(SUMDP*PLK(L))
END DO
C**** LOAD MASS EXCHANGE ARRAY FOR GWDRAG
AIRXL = 0.
DO L=LMCMIN,LMCMAX
AIRXL = AIRXL+MCFLX(L)
END DO
END IF
C**** CALCULATE OPTICAL THICKNESS
WCONST=WMU*(1.-PEARTH)+WMUL*PEARTH
WMSUM=0.
DO L=1,LMCMAX
TL(L)=(SM(L)*BYAM(L))*PLK(L)
TEMWM=(TAUMCL(L)-SVWMXL(L)*AIRM(L))*1.d2*BYGRAV
IF(TL(L).GE.TF) WMSUM=WMSUM+TEMWM
IF(CLDMCL(L).GT.0.) THEN
TAUMCL(L)=AIRM(L)*COETAU
IF(L.EQ.LMCMAX .AND. PLE(LMCMIN)-PLE(LMCMAX+1).LT.450.)
* TAUMCL(L)=AIRM(L)*.02d0
IF(L.LE.LMCMIN .AND. PLE(LMCMIN)-PLE(LMCMAX+1).GE.450.)
* TAUMCL(L)=AIRM(L)*.02d0
END IF
IF(SVLATL(L).EQ.0.) THEN
SVLATL(L)=LHE
IF ( * (TPSAV(L).eq.0. .and. TL(L).lt.TF) ) SVLATL(L)=LHS
ENDIF
IF(SVWMXL(L).GT.0.) THEN
FCLD=CLDMCL(L)+1.E-20
TEM=1.d5*SVWMXL(L)*AIRM(L)*BYGRAV
WTEM=1.d5*SVWMXL(L)*PL(L)/(FCLD*TL(L)*RGAS)
IF(SVLATL(L).EQ.LHE.AND.SVWMXL(L)/FCLD.GE.WCONST*1.d-3)
* WTEM=1d2*WCONST*PL(L)/(TL(L)*RGAS)
IF(SVLATL(L).EQ.LHS.AND.SVWMXL(L)/FCLD.GE.WMUI*1.d-3)
* WTEM=1d2*WMUI*PL(L)/(TL(L)*RGAS)
IF(WTEM.LT.1.d-10) WTEM=1.d-10
C** Set CDNC for moist conv. clds (const at present)
MNdO = 59.68d0/(RWCLDOX**3)
MNdL = 174.d0
MNdI = 0.06417127d0
MCDNCW=MNdO*(1.-PEARTH)+MNdL*PEARTH
MCDNCI=MNdI
IF(SVLATL(L).EQ.LHE) THEN
! RCLD=(RWCLDOX*10.*(1.-PEARTH)+7.0*PEARTH)*(WTEM*4.)**BY3
RCLD=100.d0*(WTEM/(2.d0*BY3*TWOPI*MCDNCW))**BY3
ELSE
! RCLD=25.0*(WTEM/4.2d-3)**BY3 * (1.+pl(l)*xRICld)
RCLD=100.d0*(WTEM/(2.d0*BY3*TWOPI*MCDNCI))**BY3
* *(1.+pl(l)*xRICld)
END IF
RCLDE=RCLD/BYBR ! effective droplet radius in anvil
CSIZEL(L)=RCLDE ! effective droplet radius in anvil
TAUMCL(L)=1.5*TEM/(FCLD*RCLDE+1.E-20)
IF(TAUMCL(L).GT.100.) TAUMCL(L)=100.
END IF
END DO
RETURN
END SUBROUTINE MSTCNV
SUBROUTINE LSCOND(IERR,WMERR,LERR,i_debug,j_debug) 1,17
!@sum LSCOND column physics of large scale condensation
!@auth M.S.Yao/A. Del Genio (modularisation by Gavin Schmidt)
!@ver 1.0 (taken from CB265)
!@calls CTMIX,QSAT,DQSATDT,THBAR
IMPLICIT NONE
!@var IERR,WMERR,LERR error reporting
INTEGER, INTENT(OUT) :: IERR,LERR
REAL*8, INTENT(OUT) :: WMERR
INTEGER, INTENT(IN) :: i_debug,j_debug
logical :: debug_out
REAL*8 LHX
C**** functions
REAL*8 :: QSAT, DQSATDT,ERMAX
!@var QSAT saturation humidity
!@var DQSATDT dQSAT/dT
!@param CM00 upper limit for autoconversion rate
!@param AIRM0 scaling factor for computing rain evaporation
!@param HEFOLD e-folding length for computing HDEP
!@param COEFM coefficient for ratio of cloud water amount and WCONST
!@param COEFT coefficient used in computing PRATM
!@param COESIG coefficient for equ. 23 of Del Genio et al. (1996)
!@param COEEC coefficient for computing cloud evaporation
!@param ERP exponential power for computing ER
C**** Adjust COEFT and COEFM to change proportion of super-cooled rain
C**** to snow. Increasing COEFT reduces temperature range of super
C**** -cooled rain, increasing COEFM enhances probability of snow.
REAL*8, PARAMETER :: CM00=1.d-4, AIRM0=100.d0, GbyAIRM0=GRAV/AIRM0
REAL*8, PARAMETER :: HEFOLD=500.,COEFM=10.,COEFT=2.5
REAL*8, PARAMETER :: COESIG=1d-3,COEEC=1000.
INTEGER, PARAMETER :: ERP=2
REAL*8, DIMENSION(IM) :: UMO1,UMO2,UMN1,UMN2 !@var dummy variables
REAL*8, DIMENSION(IM) :: VMO1,VMO2,VMN1,VMN2 !@var dummy variables
!@var Miscellaneous vertical arrays
REAL*8, DIMENSION(LM) ::
* QSATL,RHF,ATH,SQ,ER,QHEAT,
* CAREA,PREP,RH00,EC,WMXM
!@var QSATL saturation water vapor mixing ratio
!@var RHF environmental relative humidity
!@var ATH change in potential temperature
!@var SQ ERMAX dummy variables
!@var ER evaporation of precip
!@var QHEAT change of latent heat
!@var CAREA fraction of clear region
!@var PREP precip conversion rate
!@var RH00 threshold relative humidity
!@var EC cloud evaporation rate
!@var WMXM cloud water mass (mb)
REAL*8, DIMENSION(LM+1) :: PREBAR,PREICE
!@var PREBAR,PREICE precip entering layer top for total, snow
REAL*8 AIRMR,BETA,BMAX
* ,CBF,CBFC0,CK,CKIJ,CK1,CK2,CKM,CKR,CM,CM0,CM1,DFX,DQ,DQSDT
* ,DQSUM,DQUP,DRHDT,DSE,DSEC,DSEDIF,DWDT,DWDT1,DWM,ECRATE,EXPST
* ,FCLD,FMASS,FMIX,FPLUME,FPMAX,FQTOW,FRAT,FUNI
* ,FUNIL,FUNIO,HCHANG,HDEP,HPHASE,OLDLAT,OLDLHX,PFR,PMI,PML
* ,HPBL,PRATIO,QCONV,QHEATC,QLT1,QLT2,QMN1,QMN2,QMO1,QMO2,QNEW
* ,QNEWU,QOLD,QOLDU,QSATC,QSATE,RANDNO,RCLDE,RHI,RHN,RHO,RHT1
* ,RHW,SEDGE,SIGK,SLH,SMN1,SMN2,SMO1,SMO2,TEM,TEMP,TEVAP,THT1
* ,THT2,TLT1,TNEW,TNEWU,TOLD,TOLDU,TOLDUP,VDEF,WCONST,WMN1,WMN2
* ,WMNEW,WMO1,WMO2,WMT1,WMT2,WMX1,WTEM,VVEL,RCLD,FCOND
* ,PRATM,SMN12,SMO12
real*8 SNdO,SNdL,SNdI,SCDNCW,SCDNCI
!@var BETA,BMAX,CBFC0,CKIJ,CK1,CK2,PRATM dummy variabls
!@var SMN12,SMO12 dummy variables
!@var AIRMR
!@var CBF enhancing factor for precip conversion
!@var CK ratio of cloud top jumps in moist static energy and total water
!@var CKM CTEI threshold in MacVean and Mason theory
!@var CKR CTEI threshold in Randall theory
!@var CM conversion rate for large cloud water content
!@var CM0,CM1 limiting autoconversion rate for large cloud water content
!@var DFX iteration increment
!@var DQ condensed water vapor
!@var DQSDT derivative of saturation vapor pressure w.r.t. temperature
!@var DQSUM,DQUP dummy variables
!@var DRHDT time change of relative humidity
!@var DSE moist static energy jump at cloud top
!@var DSEC critical DSE for CTEI to operate
!@var DSFDIF DSE-DSEC
!@var DWDT,DWDT1 time change rates of cloud water
!@var ECRATE cloud droplet evaporation rate
!@var EXPST exponential term in determining the fraction in CTEI
!@var FCLD cloud fraction
!@var FMIX,FRAT fraction of mixing in CTEI
!@var FMASS mass of mixing in CTEI
!@var FPLUME fraction of mixing in CTEI
!@var FPMAX max fraction of mixing in CTEI
!@var FQTOW fraction of water vapour that goes to CLW
!@var FUNI the probablity for ice cloud to form
!@var FUNIL,FUNIO FUNI over land, ocean
!@var HCHANG,HPHASE latent heats for changing phase
!@var HDEP,HPBL layer depth (m) (note change of unit km-->m)
!@var OLDLAT,OLDLHX previous LHX
!@var PFR PROBABLITY OF GLACIATION OF SUPER-COOLED WATER
!@var PMI icy precip entering the layer top
!@var PML layer's cloud water devided by GRAV
!@var PRATIO PMI/PML
!@var QCONV convergence of latent heat
!@var QHEATC,QLT1,QLT2,QMN1,QMN2,QMO1,QMO2 dummy variables
!@var QNEW,QNEWU updated specific humidity
!@var QOLD,QOLDU previous specific humidity
!@var QSATC saturation vapor mixing ratio
!@var QSATE saturation vapor mixing ratio w.r.t. water
!@var RANDNO random number
!@var RCLD,RCLDE cloud particle's radius, effective radius
!@var RHI relative humidity w.r.t. ice
!@var RHN dummy variable
!@var RHO air density
!@var RHT1,RHW,SIGK dummy variables
!@var SEDGE potential temperature at layer edge
!@var SMN1,SMN2,SMO1,SMO2,TEM,TEMP,TEVAP,THT1,THT2,TLT1 dummy variables
!@var TNEW,TNEWU updated tempertures
!@var TOLD,TOLDU,TOLDUP previous temperatures
!@var VDEF = VVEL - VSUB
!@var WCONST,WMN1,WMN2,WMO1,WMO2,WMT1,WMT2,WMX1 dummy variables
!@var WMNEW updated cloud water mixing ratio
!@var WTEM cloud water density (g m**-3)
!@var VVEL vertical velocity (cm/s)
!@var FCOND dummy variables
INTEGER LN,ITER !@var LN,ITER loop variables
LOGICAL BANDF !@var BANDF true if Bergero-Findeisen proc. occurs
LOGICAL FORM_CLOUDS !@var FORM_CLOUDS true if clouds are formed
INTEGER K,L,N !@var K,L,N loop variables
REAL*8 THBAR !@var THBAR potential temperature at layer edge
C****
C**** LARGE-SCALE CLOUDS AND PRECIPITATION
C**** THE LIQUID WATER CONTENT IS PREDICTED
C****
IERR=0
C****
debug_out=.false.
PRCPSS=0.
HCNDSS=0.
CKIJ=1.
C**** initialise vertical arrays
ER=0.
EC=0.
PREP=0.
PREBAR=0.
LHP=0.
QHEAT=0.
CLDSSL=0
TAUSSL=0
DO L=1,LP50
CAREA(L)=1.-CLDSAVL(L)
C IF(RH(L).LE.1.) CAREA(L)=DSQRT((1.-RH(L))/(1.-RH00(L)+teeny))
C IF(RH(L).LE.1.) CAREA(L)=DSQRT((1.-RH(L))/(1.-U00ice+teeny))
C IF(CAREA(L).GT.1.) CAREA(L)=1.
C IF(RH(L).GT.1.) CAREA(L)=0.
IF(WMX(L).LE.0.) CAREA(L)=1.
END DO
DQUP=0.
TOLDUP=TL(LP50)
PREICE(LP50+1)=0.
WCONST=WMU*(1.-PEARTH)+WMUL*PEARTH
SSHR=0.
DCTEI=0.
C****
C**** MAIN L LOOP FOR LARGE-SCALE CONDENSATION, PRECIPITATION AND CLOUDS
C****
DO L=LP50,1,-1
TOLD=TL(L)
QOLD=QL(L)
OLDLHX=SVLHXL(L)
OLDLAT=SVLATL(L)
C**** COMPUTE VERTICAL VELOCITY IN CM/S
TEMP=100.*RGAS*TL(L)/(PL(L)*GRAV)
IF(L.EQ.1) THEN
VVEL=-SDL(L+1)*TEMP
ELSE IF(L.EQ.LM) THEN
VVEL=-SDL(L)*TEMP
ELSE
VVEL=-.5*(SDL(L)+SDL(L+1))*TEMP
END IF
C**** COMPUTE THE LIMITING AUTOCONVERSION RATE FOR CLOUD WATER CONTENT
CM0=CM00
VDEF=VVEL-VSUBL(L)
IF(VDEF.GT.0.) CM0=CM00*10.**(-VDEF)
c FCLD=1.-CAREA(L)+1.E-20 !!!
FCLD=(1.-CAREA(L))*FSSL(L)+teeny
C**** COMPUTE THE PROBABILITY OF ICE FORMATION, FUNI, AND
C**** THE PROBABLITY OF GLACIATION OF SUPER-COOLED WATER, PFR
C**** DETERMINE THE PHASE MOISTURE CONDENSES TO
C**** DETERMINE THE POSSIBILITY OF B-F PROCESS
BANDF=.FALSE.
LHX=LHE
CBF=1. + EXP(
IF (TL(L).LE.TI) THEN ! below -40: force ice
LHX=LHS
ELSEIF (TL(L).GE.TF) THEN ! above freezing: force water
LHX=LHE
ELSE ! in between: compute probability
IF(TL(L).GT.269.16) THEN ! OC/SI/LI clouds: water above -4
FUNIO=0.
ELSE
FUNIO=1.-EXP( END IF
IF(TL(L).GT.263.16) THEN ! land clouds water: above -10
FUNIL=0.
ELSE
FUNIL=1.-EXP( END IF
FUNI=FUNIO*(1.-PEARTH)+FUNIL*PEARTH
RANDNO=RNDSSL(1,L) ! RANDNO=RANDU(XY)
IF(RANDNO.LT.FUNI) LHX=LHS
C**** special case 1) if ice previously then stay as ice (if T<Tf)
IF((OLDLHX.EQ.LHS.OR.OLDLAT.EQ.LHS).AND.TL(L).LT.TF) THEN
IF(LHX.EQ.LHE) BANDF=.TRUE.
LHX=LHS
ENDIF
IF (L.LT.LP50) THEN
C**** Decide whether precip initiates B-F process
PML=WMX(L)*AIRM(L)*BYGRAV
PMI=PREICE(L+1)*DTsrc
RANDNO=RNDSSL(2,L) ! RANDNO=RANDU(XY)
C**** Calculate probability of ice precip seeding a water cloud
IF (LHX.EQ.LHE.AND.PMI.gt.0) THEN
PRATIO=MIN(PMI/(PML+1.E-20),10d0)
CBFC0=.5*CM0*CBF*DTsrc
PFR=(1.-EXP( IF(PFR.GT.RANDNO) THEN
BANDF=.TRUE.
LHX=LHS
END IF
END IF
C**** If liquid rain falls into an ice cloud, B-F must occur
IF (LHP(L+1).EQ.LHE .AND. LHX.EQ.LHS .AND. PML.GT.0.)
* BANDF=.TRUE.
END IF
END IF
IF(LHX.EQ.LHS .AND. (OLDLHX.EQ.LHE.OR.OLDLAT.EQ.LHE)) BANDF=.TRUE.
C**** COMPUTE RELATIVE HUMIDITY
QSATL(L)=QSAT(TL(L),LHX,PL(L))
RH1(L)=QL(L)/QSATL(L)
IF(LHX.EQ.LHS) THEN
QSATE=QSAT(TL(L),LHE,PL(L))
RHW=.00536d0*TL(L)-.276d0
RH1(L)=QL(L)/QSATE
IF(TL(L).LT.238.16) RH1(L)=QL(L)/(QSATE*RHW)
END IF
C**** PHASE CHANGE OF CLOUD WATER CONTENT
HCHANG=0.
IF(OLDLHX.EQ.LHE.AND.LHX.EQ.LHS) HCHANG= WML(L)*LHM
IF(OLDLHX.EQ.LHS.AND.LHX.EQ.LHE) HCHANG=-WML(L)*LHM
IF(OLDLAT.EQ.LHE.AND.LHX.EQ.LHS) HCHANG=HCHANG+SVWMXL(L)*LHM
IF(OLDLAT.EQ.LHS.AND.LHX.EQ.LHE) HCHANG=HCHANG-SVWMXL(L)*LHM
SVLHXL(L)=LHX
TL(L)=TL(L)+HCHANG/(SHA*FSSL(L)+teeny)
TH(L)=TL(L)/PLK(L)
ATH(L)=(TH(L)-TTOLDL(L))*BYDTsrc
C**** COMPUTE RH IN THE CLOUD-FREE AREA, RHF
RHI=QL(L)/QSAT(TL(L),LHS,PL(L))
! this formulation is used for consistency with current practice
RH00(L)=U00wtrX*U00ice
IF(LHX.EQ.LHS) RH00(L)=U00ice
C**** Option to treat boundary layer differently
IF (do_blU00.eq.1) then
IF (L.LE.DCL) THEN ! boundary layer clouds
C**** calculate total pbl depth
HPBL=0.
DO LN=1,DCL
HPBL=HPBL+AIRM(LN)*TL(LN)*RGAS/(GRAV*PL(LN))
END DO
C**** Scale HPBL by HRMAX to provide tuning control for PBL clouds
HDEP = MIN(HPBL,HRMAX*(1.-EXP(-HPBL/HEFOLD)))
C**** Special conditions for boundary layer contained wholly in layer 1
IF (DCL.LE.1) THEN
IF (RIS.GT.1.) HDEP=10d0
IF (RIS.LE.1..AND.RI1.GT.1.) HDEP=50d0
IF (RIS.LE.1..AND.RI1.LE.1..AND.RI2.GT.1.) HDEP=100d0
END IF
C**** Estimate critical rel. hum. based on parcel lifting argument
RH00(L)=1.-GAMD*LHE*HDEP/(RVAP*TS*TS)
IF(RH00(L).LT.0.) RH00(L)=0.
END IF
END IF
C****
IF(RH00(L).GT.1.) RH00(L)=1.
RHF(L)=RH00(L)+(1.-CAREA(L))*(1.-RH00(L))
C**** Set precip phase to be the same as the cloud, unless precip above
C**** is ice and temperatures after ice melt would still be below TFrez
LHP(L)=LHX
IF (LHP(L+1).eq.LHS .and.
* TL(L).lt.TF+DTsrc*LHM*PREICE(L+1)*GRAV*BYAM(L)*BYSHA/(FSSL(L)
* +teeny)) LHP(L)=LHP(L+1)
C***Setting constant values of CDNC over land and ocean to get RCLD=f(CDNC,LWC)
SNdO = 59.68d0/(RWCLDOX**3)
SNdL = 174.d0
SNdI = 0.06417127d0
SCDNCW=SNdO*(1.-PEARTH)+SNdL*PEARTH
SCDNCI=SNdI
C**** COMPUTE THE AUTOCONVERSION RATE OF CLOUD WATER TO PRECIPITATION
IF(WMX(L).GT.0.) THEN
RHO=1d5*PL(L)/(RGAS*TL(L))
TEM=RHO*WMX(L)/(WCONST*FCLD+teeny)
IF(LHX.EQ.LHS ) TEM=RHO*WMX(L)/(WMUI*FCLD+teeny)
TEM=TEM*TEM
IF(TEM.GT.10.) TEM=10.
CM1=CM0
IF(BANDF) CM1=CM0*CBF
IF(LHX.EQ.LHS) CM1=CM0
CM=CM1*(1.-1./EXP(TEM*TEM))+100.*(PREBAR(L+1)+
* PRECNVL(L+1)*BYDTsrc)
IF(CM.GT.BYDTsrc) CM=BYDTsrc
PREP(L)=WMX(L)*CM
IF(TL(L).LT.TF.AND.LHX.EQ.LHE) THEN ! check snowing pdf
PRATM=1d5*COEFM*WMX(L)*PL(L)/(WCONST*FCLD*TL(L)*RGAS+teeny)
PRATM=MIN(PRATM,1d0)*(1.-EXP(MAX(-1d2,(TL(L)-TF)/COEFT)))
IF(PRATM.GT.RNDSSL(3,L)) LHP(L)=LHS
END IF
ELSE
CM=0.
END IF
C**** DECIDE WHETHER TO FORM CLOUDS
C**** FORM CLOUDS ONLY IF RH GT RH00
IF (RH1(L).LT.RH00(L)) THEN
FORM_CLOUDS=.FALSE.
ELSE ! COMPUTE THE CONVERGENCE OF AVAILABLE LATENT HEAT
SQ(L)=LHX*QSATL(L)*DQSATDT(TL(L),LHX)*BYSHA
TEM=-LHX*DPDT(L)/PL(L)
QCONV=LHX*AQ(L)-RH(L)*SQ(L)*SHA*PLK(L)*ATH(L)
* -TEM*QSATL(L)*RH(L)
FORM_CLOUDS= (QCONV.GT.0. .OR. WMX(L).GT.0.)
END IF
C****
ERMAX=LHX*PREBAR(L+1)*GRAV*BYAM(L)
IF (FORM_CLOUDS) THEN
C**** COMPUTE EVAPORATION OF RAIN WATER, ER
RHN=MIN(RH(L),RHF(L))
IF(WMX(L).GT.0.) THEN
ER(L)=(1.-RHN)**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0
ELSE ! WMX(l).le.0.
IF(PREICE(L+1).GT.0..AND.TL(L).LT.TF) THEN
ER(L)=(1.-RHI)**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0
ELSE
ER(L)=(1.-RH(L))**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0
END IF
END IF
ER(L)=MAX(0d0,MIN(ER(L),ERMAX))
C**** COMPUTATION OF CLOUD WATER EVAPORATION
IF (CAREA(L).GT.0.) THEN
WTEM=1d5*WMX(L)*PL(L)/(FCLD*TL(L)*RGAS+teeny)
IF(LHX.EQ.LHE.AND.WMX(L)/FCLD.GE.WCONST*1d-3)
* WTEM=1d2*WCONST*PL(L)/(TL(L)*RGAS)
IF(LHX.EQ.LHS.AND.WMX(L)/FCLD.GE.WMUI*1d-3)
* WTEM=1d2*WMUI*PL(L)/(TL(L)*RGAS)
IF(WTEM.LT.1d-10) WTEM=1d-10
IF(LHX.EQ.LHE) THEN
! RCLD=1d-6*(RWCLDOX*10.*(1.-PEARTH)+7.*PEARTH)*(WTEM*4.)**BY3
RCLD=1d-6*100.d0*(WTEM/(2.d0*BY3*TWOPI*SCDNCW))**BY3
ELSE
! RCLD=25.d-6*(WTEM/4.2d-3)**BY3 * (1.+pl(l)*xRICld)
RCLD=100.d-6*(WTEM/(2.d0*BY3*TWOPI*SCDNCI))**BY3
* *(1.+pl(l)*xRICld)
END IF
CK1=1000.*LHX*LHX/(2.4d-2*RVAP*TL(L)*TL(L))
CK2=1000.*RGAS*TL(L)/(2.4d-3*QSATL(L)*PL(L))
TEVAP=COEEC*(CK1+CK2)*RCLD*RCLD
WMX1=WMX(L)-PREP(L)*DTsrc
ECRATE=(1.-RHF(L))/(TEVAP*FCLD+teeny)
IF(ECRATE.GT.BYDTsrc) ECRATE=BYDTsrc
EC(L)=WMX1*ECRATE*LHX
END IF
C**** COMPUTE NET LATENT HEATING DUE TO STRATIFORM CLOUD PHASE CHANGE,
C**** QHEAT, AND NEW CLOUD WATER CONTENT, WMNEW
DRHDT=2.*CAREA(L)*CAREA(L)*(1.-RH00(L))*(QCONV+ER(L))/LHX/
* (WMX(L)/(FCLD+teeny)+2.*CAREA(L)*QSATL(L)*(1.-RH00(L))
* +teeny)
IF(ER(L).EQ.0.AND.WMX(L).LE.0.) DRHDT=0.
QHEAT(L)=FSSL(L)*(QCONV-LHX*DRHDT*QSATL(L))/(1.+RH(L)*SQ(L))
DWDT=QHEAT(L)/LHX-PREP(L)+CAREA(L)*FSSL(L)*ER(L)/LHX
WMNEW =WMX(L)+DWDT*DTsrc
IF(WMNEW.LT.0.) THEN
WMNEW=0.
QHEAT(L)=(-WMX(L)*BYDTsrc+PREP(L))*LHX-CAREA(L)*FSSL(L)*ER(L)
END IF
ELSE
C**** UNFAVORABLE CONDITIONS FOR CLOUDS TO EXIT, PRECIP OUT CLOUD WATER
IF(WMX(L).GT.0.) PREP(L)=WMX(L)*BYDTsrc
ER(L)=(1.-RH(L))**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0
IF(PREICE(L+1).GT.0..AND.TL(L).LT.TF)
* ER(L)=(1.-RHI)**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0
ER(L)=MAX(0d0,MIN(ER(L),ERMAX))
QHEAT(L)=-CAREA(L)*FSSL(L)*ER(L)
WMNEW=0.
END IF
C**** PHASE CHANGE OF PRECIPITATION, FROM ICE TO WATER
C**** This occurs if 0 C isotherm is crossed, or ice is falling into
C**** a super-cooled water cloud (that has not had B-F occur).
C**** Note: on rare occasions we have ice clouds even if T>0
C**** In such a case, no energy of phase change is needed.
HPHASE=0.
IF (LHP(L+1).EQ.LHS.AND.LHP(L).EQ.LHE.AND.PREICE(L+1).GT.0) THEN
HPHASE=HPHASE+LHM*PREICE(L+1)*GRAV*BYAM(L)
PREICE(L+1)=0.
ENDIF
C**** PHASE CHANGE OF PRECIP, FROM WATER TO ICE
IF (LHP(L+1).EQ.LHE.AND.LHP(L).EQ.LHS.AND.PREBAR(L+1).GT.0)
* HPHASE=HPHASE-LHM*PREBAR(L+1)*GRAV*BYAM(L)
C**** Make sure energy is conserved for transfers between P and CLW
IF (LHP(L).NE.LHX)
* HPHASE=HPHASE+(ER(L)*CAREA(L)*FSSL(L)/LHX-PREP(L))*LHM
C**** COMPUTE THE PRECIP AMOUNT ENTERING THE LAYER TOP
IF (ER(L).eq.ERMAX) THEN ! to avoid round off problem
PREBAR(L)=PREBAR(L+1)*(1.-CAREA(L)*FSSL(L))+
* AIRM(L)*PREP(L)*BYGRAV
ELSE
PREBAR(L)=MAX(0d0,PREBAR(L+1)+
* AIRM(L)*(PREP(L)-ER(L)*CAREA(L)*FSSL(L)/LHX)*BYGRAV)
END IF
C**** UPDATE NEW TEMPERATURE AND SPECIFIC HUMIDITY
QNEW =QL(L)-DTsrc*QHEAT(L)/(LHX*FSSL(L)+teeny)
IF(QNEW.LT.0.) THEN
QNEW=0.
QHEAT(L)=QL(L)*LHX*BYDTsrc*FSSL(L)
DWDT1=QHEAT(L)/LHX-PREP(L)+CAREA(L)*FSSL(L)*ER(L)/LHX
WMNEW=WMX(L)+DWDT1*DTsrc
C**** IF WMNEW .LT. 0., THE COMPUTATION IS UNSTABLE
IF(WMNEW.LT.0.) THEN
IERR=1
LERR=L
WMERR=WMNEW
WMNEW=0.
END IF
END IF
C**** Only Calculate fractional changes of Q to W
FQTOW=0. ! Q->CLW
IF (FSSL(L).gt.0) THEN
IF (QHEAT(L)+CAREA(L)*FSSL(L)*ER(L).gt.0) THEN
IF (LHX*QL(L)+DTsrc*CAREA(L)*ER(L).gt.0.) FQTOW=(QHEAT(L
* )+CAREA(L)*FSSL(L)*ER(L))*DTsrc/((LHX*QL(L)+DTsrc*CAREA(L)
* *ER(L))*FSSL(L))
END IF
END IF
QL(L)=QNEW
C**** adjust gradients down if Q decreases
QMOM(:,L)= QMOM(:,L)*(1.-FQTOW)
WMX(L)=WMNEW
! if(abs(DTsrc*(QHEAT(L)-HPHASE)).gt.100.*SHA*FSSL(L)) then ! warn
! write(0,*) 'it,i,j,l,tlold,dtl,qht,hph',itime,i_debug,j_debug,
! * L,TL(L),DTsrc*(QHEAT(L)-HPHASE)/(SHA*FSSL(L)+teeny),
! * QHEAT(L),HPHASE
! debug_out=.true.
! end if
TL(L)=TL(L)+DTsrc*(QHEAT(L)-HPHASE)/(SHA*FSSL(L)+teeny)
TH(L)=TL(L)/PLK(L)
TNEW=TL(L)
QSATC=QSAT(TL(L),LHX,PL(L))
RH(L)=QL(L)/QSATC
C**** CONDENSE MORE MOISTURE IF RELATIVE HUMIDITY .GT. 1
IF(RH(L).GT.1.) THEN
SLH=LHX*BYSHA
DQSUM=0.
DO N=1,3
IF(N.NE.1) QSATC=QSAT(TL(L),LHX,PL(L))
DQ=(QL(L)-QSATC)/(1.+SLH*QSATC*DQSATDT(TL(L),LHX))
TL(L)=TL(L)+SLH*DQ
QL(L)=QL(L)-DQ
DQSUM=DQSUM+DQ
END DO
IF(DQSUM.GT.0.) THEN
WMX(L)=WMX(L)+DQSUM*FSSL(L)
FCOND=DQSUM/QNEW
C**** adjust gradients down if Q decreases
QMOM(:,L)= QMOM(:,L)*(1.-FCOND)
ELSE
TL(L)=TNEW
QL(L)=QNEW
END IF
RH(L)=QL(L)/QSAT(TL(L),LHX,PL(L))
TH(L)=TL(L)/PLK(L)
TNEW=TL(L)
! if (debug_out) write(0,*) 'after condensation: l,tlnew,',l,tl(l)
END IF
IF(RH(L).LE.1.) CAREA(L)=DSQRT((1.-RH(L))/(1.-RH00(L)+teeny))
IF(CAREA(L).GT.1.) CAREA(L)=1.
IF(RH(L).GT.1.) CAREA(L)=0.
IF(WMX(L).LE.0.) CAREA(L)=1.
IF(CAREA(L).LT.0.) CAREA(L)=0.
RHF(L)=RH00(L)+(1.-CAREA(L))*(1.-RH00(L))
IF(RH(L).LE.RHF(L).AND.RH(L).LT..999999.AND.WMX(L).gt.0.) THEN
C**** PRECIP OUT CLOUD WATER IF RH LESS THAN THE RH OF THE ENVIRONMENT
PREBAR(L)=PREBAR(L)+WMX(L)*AIRM(L)*BYGRAV*BYDTsrc
IF(LHP(L).EQ.LHS .AND. LHX.EQ.LHE) THEN
HCHANG=WMX(L)*LHM
TL(L)=TL(L)+HCHANG/(SHA*FSSL(L)+teeny)
! if(debug_out) write(0,*) 'after rain out: l,tlnew',l,tl(l)
TH(L)=TL(L)/PLK(L)
END IF
WMX(L)=0.
END IF
C**** set phase of condensation for next box down
PREICE(L)=0.
IF (PREBAR(L).gt.0 .AND. LHP(L).EQ.LHS) PREICE(L)=PREBAR(L)
IF (PREBAR(L).le.0) LHP(L)=0.
C**** COMPUTE THE LARGE-SCALE CLOUD COVER
IF(RH(L).LE.1.) CAREA(L)=DSQRT((1.-RH(L))/(1.-RH00(L)+teeny))
IF(CAREA(L).GT.1.) CAREA(L)=1.
IF(RH(L).GT.1.) CAREA(L)=0.
IF(WMX(L).LE.0.) CAREA(L)=1.
IF(CAREA(L).LT.0.) CAREA(L)=0.
CLDSSL(L)=FSSL(L)*(1.-CAREA(L))
CLDSAVL(L)=1.-CAREA(L)
TOLDUP=TOLD
C**** ACCUMULATE SOME DIAGNOSTICS
HCNDSS=HCNDSS+FSSL(L)*(TNEW-TOLD)*AIRM(L)
SSHR(L)=SSHR(L)+FSSL(L)*(TNEW-TOLD)*AIRM(L)
END DO ! end of loop over L
PRCPSS=MAX(0d0,PREBAR(1)*GRAV*DTsrc) ! fix small round off err
C****
C**** CLOUD-TOP ENTRAINMENT INSTABILITY
C****
DO L=LP50-1,1,-1
SM(L)=TH(L)*AIRM(L)
QM(L)=QL(L)*AIRM(L)
WMXM(L)=WMX(L)*AIRM(L)
SM(L+1)=TH(L+1)*AIRM(L+1)
QM(L+1)=QL(L+1)*AIRM(L+1)
WMXM(L+1)=WMX(L+1)*AIRM(L+1)
TOLD=TL(L)
TOLDU=TL(L+1)
QOLD=QL(L)
QOLDU=QL(L+1)
FCLD=(1.-CAREA(L))*FSSL(L)+teeny
IF(CAREA(L).EQ.1. .OR. (CAREA(L).LT.1..AND.CAREA(L+1).LT.1.))
* CYCLE
SEDGE=THBAR(TH(L+1),TH(L))
DSE=(TH(L+1)-SEDGE)*PLK(L+1)+(SEDGE-TH(L))*PLK(L)+
* SLHE*(QL(L+1)-QL(L))
DWM=QL(L+1)-QL(L)+(WMX(L+1)-WMX(L))/FCLD
DQSDT=DQSATDT(TL(L),LHE)*QL(L)/(RH(L)+1d-30)
BETA=(1.+BYMRAT*TL(L)*DQSDT)/(1.+SLHE*DQSDT)
CKM=(1.+SLHE*DQSDT)*(1.+(1.-DELTX)*TL(L)/SLHE)/
* (2.+(1.+BYMRAT*TL(L)/SLHE)*SLHE*DQSDT)
CKR=TL(L)/(BETA*SLHE)
CK=DSE/(SLHE*DWM)
SIGK=0.
IF(CKR.GT.CKM) CYCLE
IF(CK.GT.CKR) SIGK=COESIG*((CK-CKR)/(CKM-CKR+teeny))**5
EXPST=EXP(-SIGK*DTsrc)
IF(L.LE.1) CKIJ=EXPST
DSEC=DWM*TL(L)/BETA
IF(CK.LT.CKR) CYCLE
FPMAX=MIN(1d0,1.-EXPST)
IF(FPMAX.LE.0.) CYCLE
IF(DSE.GE.DSEC) CYCLE
C**** MIXING TO REMOVE CLOUD-TOP ENTRAINMENT INSTABILITY
AIRMR=(AIRM(L+1)+AIRM(L))*BYAM(L+1)*BYAM(L)
SMO1=SM(L)
QMO1=QM(L)
WMO1=WMXM(L)
SMO2=SM(L+1)
QMO2=QM(L+1)
WMO2=WMXM(L+1)
SMO12=SMO1*PLK(L)+SMO2*PLK(L+1)
DO K=1,KMAX
UMO1(K)=UM(K,L)
VMO1(K)=VM(K,L)
UMO2(K)=UM(K,L+1)
VMO2(K)=VM(K,L+1)
ENDDO
FPLUME=FPMAX*FSSL(L)
DFX=FPLUME ! not FPMAX
DO ITER=1,9
DFX=DFX*0.5
FMIX=FPLUME*FCLD
FMASS=FMIX*AIRM(L)
FMASS=MIN(FMASS,(AIRM(L+1)*AIRM(L))/(AIRM(L+1)+AIRM(L)))
FMIX=FMASS*BYAM(L)
FRAT=FMASS*BYAM(L+1)
SMN1=SMO1*(1.-FMIX)+FRAT*SMO2
QMN1=QMO1*(1.-FMIX)+FRAT*QMO2
WMN1=WMO1*(1.-FMIX)+FRAT*WMO2
SMN2=SMO2*(1.-FRAT)+FMIX*SMO1
QMN2=QMO2*(1.-FRAT)+FMIX*QMO1
WMN2=WMO2*(1.-FRAT)+FMIX*WMO1
THT1=SMN1*BYAM(L)/PLK(L)
QLT1=QMN1*BYAM(L)
TLT1=THT1*PLK(L)
LHX=SVLHXL(L)
RHT1=QLT1/(QSAT(TLT1,LHX,PL(L)))
WMT1=WMN1*BYAM(L)
THT2=SMN2*BYAM(L+1)/PLK(L+1)
QLT2=QMN2*BYAM(L+1)
WMT2=WMN2*BYAM(L+1)
SEDGE=THBAR(THT2,THT1)
DSE=(THT2-SEDGE)*PLK(L+1)+(SEDGE-THT1)*PLK(L)+SLHE*(QLT2-QLT1)
DWM=QLT2-QLT1+(WMT2-WMT1)/FCLD
DQSDT=DQSATDT(TLT1,LHE)*QLT1/(RHT1+1d-30)
BETA=(1.+BYMRAT*TLT1*DQSDT)/(1.+SLHE*DQSDT)
CKM=(1.+SLHE*DQSDT)*(1.+(1.-DELTX)*TLT1/SLHE)/
* (2.+(1.+BYMRAT*TLT1/SLHE)*SLHE*DQSDT)
DSEC=DWM*TLT1/BETA
DSEDIF=DSE-DSEC
IF(DSEDIF.GT.1d-3) FPLUME=FPLUME-DFX
IF(DSEDIF.LT.-1d-3) FPLUME=FPLUME+DFX
IF(ABS(DSEDIF).LE.1d-3.OR.FPLUME.GT.FPMAX) EXIT
END DO
C**** UPDATE TEMPERATURE, SPECIFIC HUMIDITY AND MOMENTUM DUE TO CTEI
SMN12=SMN1*PLK(L)+SMN2*PLK(L+1)
SMN1=SMN1-(SMN12-SMO12)*AIRM(L)/((AIRM(L)+AIRM(L+1))*PLK(L))
SMN2=SMN2-(SMN12-SMO12)*AIRM(L+1)/((AIRM(L)+AIRM(L+1))*PLK(L+1))
TH(L)=SMN1*BYAM(L)
TL(L)=TH(L)*PLK(L)
! if(debug_out) write(0,*) 'after CTEI: l,tlnew',l,tl(l)
QL(L)=QMN1*BYAM(L)
LHX=SVLHXL(L)
RH(L)=QL(L)/QSAT(TL(L),LHX,PL(L))
WMX(L)=WMN1*BYAM(L)
TH(L+1)=SMN2*BYAM(L+1)
QL(L+1)=QMN2*BYAM(L+1)
WMX(L+1)=WMN2*BYAM(L+1)
CALL CTMIX (SM(L),SMOM(1,L),FMASS*AIRMR,FMIX,FRAT)
CALL CTMIX (QM(L),QMOM(1,L),FMASS*AIRMR,FMIX,FRAT)
C****
DO K=1,KMAX
UMN1(K)=(UMO1(K)*(1.-FMIX)+FRAT*UMO2(K))
VMN1(K)=(VMO1(K)*(1.-FMIX)+FRAT*VMO2(K))
UMN2(K)=(UMO2(K)*(1.-FRAT)+FMIX*UMO1(K))
VMN2(K)=(VMO2(K)*(1.-FRAT)+FMIX*VMO1(K))
UM(K,L)=UM(K,L)+(UMN1(K)-UMO1(K))*RA(K)
VM(K,L)=VM(K,L)+(VMN1(K)-VMO1(K))*RA(K)
UM(K,L+1)=UM(K,L+1)+(UMN2(K)-UMO2(K))*RA(K)
VM(K,L+1)=VM(K,L+1)+(VMN2(K)-VMO2(K))*RA(K)
END DO
C**** RE-EVAPORATION OF CLW IN THE UPPER LAYER
QL(L+1)=QL(L+1)+WMX(L+1)/(FSSL(L)+teeny)
TH(L+1)=TH(L+1)-(LHX*BYSHA)*WMX(L+1)/(PLK(L+1)*FSSL(L)+teeny)
TL(L+1)=TH(L+1)*PLK(L+1)
! if(debug_out) write(0,*) 'after re-evap: l,tlnew',l+1,tl(l+1)
RH(L+1)=QL(L+1)/QSAT(TL(L+1),LHX,PL(L+1))
WMX(L+1)=0.
IF(RH(L).LE.1.) CAREA(L)=DSQRT((1.-RH(L))/(1.-RH00(L)+teeny))
IF(CAREA(L).GT.1.) CAREA(L)=1.
IF(RH(L).GT.1.) CAREA(L)=0.
CLDSSL(L)=FSSL(L)*(1.-CAREA(L))
CLDSAVL(L)=1.-CAREA(L)
TNEW=TL(L)
TNEWU=TL(L+1)
QNEW=QL(L)
QNEWU=QL(L+1)
HCNDSS=HCNDSS+FSSL(L)*(TNEW-TOLD)*AIRM(L)+
* FSSL(L+1)*(TNEWU-TOLDU)*AIRM(L+1)
SSHR(L)=SSHR(L)+FSSL(L)*(TNEW-TOLD)*AIRM(L)
SSHR(L+1)=SSHR(L+1)+FSSL(L+1)*(TNEWU-TOLDU)*AIRM(L+1)
DCTEI(L)=DCTEI(L)+FSSL(L)*(QNEW-QOLD)*AIRM(L)*LHX*BYSHA
DCTEI(L+1)=DCTEI(L+1)+FSSL(L+1)*(QNEWU-QOLDU)*AIRM(L+1)*LHX*BYSHA
END DO
C**** COMPUTE CLOUD PARTICLE SIZE AND OPTICAL THICKNESS
WMSUM=0.
DO L=1,LP50
FCLD=CLDSSL(L)+teeny
WTEM=1.d5*WMX(L)*PL(L)/(FCLD*TL(L)*RGAS+teeny)
LHX=SVLHXL(L)
IF(LHX.EQ.LHS.AND.WMX(L)/FCLD.GE.WMUI*1d-3)
* WTEM=1d5*WMUI*1.d-3*PL(L)/(TL(L)*RGAS)
IF(WTEM.LT.1d-10) WTEM=1.d-10
C***Setting constant values of CDNC over land and ocean to get RCLD=f(CDNC,LWC)
SNdO = 59.68d0/(RWCLDOX**3)
SNdL = 174.d0
SNdI = 0.06417127d0
SCDNCW=SNdO*(1.-PEARTH)+SNdL*PEARTH
SCDNCI=SNdI
IF(LHX.EQ.LHE) THEN
! RCLD=(RWCLDOX*10.*(1.-PEARTH)+7.0*PEARTH)*(WTEM*4.)**BY3
RCLD=100.d0*(WTEM/(2.d0*BY3*TWOPI*SCDNCW))**BY3
QHEATC=(QHEAT(L)+FSSL(L)*CAREA(L)*(EC(L)+ER(L)))/LHX
IF(RCLD.GT.20..AND.PREP(L).GT.QHEATC) RCLD=20.
ELSE
! RCLD=25.0*(WTEM/4.2d-3)**BY3 * (1.+pl(l)*xRICld)
RCLD=100.d0*(WTEM/(2.d0*BY3*TWOPI*SCDNCI))**BY3
* *(1.+pl(l)*xRICld)
ENDIF
RCLDE=RCLD/BYBR
CSIZEL(L)=RCLDE
TEM=AIRM(L)*WMX(L)*1.d2*BYGRAV
TAUSSL(L)=1.5d3*TEM/(FCLD*RCLDE+teeny)
IF(TAUSSL(L).GT.100.) TAUSSL(L)=100.
IF(LHX.EQ.LHE) WMSUM=WMSUM+TEM
END DO
C**** CALCULATE OPTICAL THICKNESS
DO L=1,LP50
CLDSV1(L)=CLDSSL(L)
IF(WMX(L).LE.0.) SVLHXL(L)=0.
IF(TAUMCL(L).EQ.0..OR.CKIJ.NE.1.) THEN
BMAX=1.-EXP( IF(CLDSV1(L).GE..95d0) BMAX=CLDSV1(L)
IF(L.EQ.1.OR.L.LE.DCL) THEN
CLDSSL(L)=CLDSSL(L)+(BMAX-CLDSSL(L))*CKIJ
TAUSSL(L)=TAUSSL(L)*CLDSV1(L)/(CLDSSL(L)+teeny)
ENDIF
IF(TAUSSL(L).LE.0.) CLDSSL(L)=0.
IF(L.GT.DCL .AND. TAUMCL(L).LE.0.) THEN
CLDSSL(L)=CLDSSL(L)**(2.*BY3)
TAUSSL(L)=TAUSSL(L)*CLDSV1(L)**BY3
END IF
END IF
END DO
RETURN
END SUBROUTINE LSCOND
END MODULE CLOUDS
C-----------------------------------------------------------------------
SUBROUTINE ISCCP_CLOUD_TYPES(pfull,phalf,qv,cc,conv,dtau_s,dtau_c 1,5
* ,skt,at,dem_s,dem_c,itrop,fq_isccp,ctp,tauopt,nbox,jerr)
!@sum ISCCP_CLOUD_TYPES calculate isccp cloud diagnostics in a column
!@auth G. Tselioudis (modifications by Gavin Schmidt)
!@ver 1.0
USE CONSTANT, only : bygrav, wtmair=>mair, bymrat,avog
USE RANDOM, only : rinit,rfinal,randu
USE MODEL_COM, only : nlev=>lm,qcheck
implicit none
!@var overlap type: 1=max, 2=rand, 3=max/rand
!@var top_height 1 = adjust top height, that is compute infrared
!@+ brightness temperature and adjust cloud top pressure accordingly
!@+ 2 = do not adjust top height, that is cloud top
!@+ pressure is the actual cloud top pressure in the model
INTEGER, PARAMETER :: top_height=1, overlap=3
!@var emsfc_lw longwave emissivity of surface at 10.5 microns
REAL*8, PARAMETER :: emsfc_lw=0.99d0
INTEGER,PARAMETER :: ncol =20 !@var ncol number of subcolumns
REAL*8, PARAMETER :: byncol = 1d0/ncol
!@var pc1bylam Planck constant c1 by wavelength (10.5 microns)
REAL*8, PARAMETER :: pc1bylam = 1.439d0/10.5d-4
!@var t0 ave temp (K)
REAL*8, PARAMETER :: t0 = 296.
!@var t0bypstd ave temp by sea level press
REAL*8, PARAMETER :: t0bypstd = t0/1.013250d6
real*8, parameter :: bywc = 1./2.56d0 , byic= 1./2.13d0
! -----
! Input
! -----
!@var pfull pressure of full model levels (Pascals)
REAL*8 pfull(nlev) ! pfull(1) is top level, pfull(nlev) is bottom
!@var phalf pressure of half model levels (Pascals)
REAL*8 phalf(nlev+1) ! phalf(1) is top, phalf(nlev+1) is surf pres
!@var qv water vapor specific humidity (kg vapor/ kg air)
REAL*8 qv(nlev)
!@var cc input cloud cover in each model level (fraction)
REAL*8 cc(nlev) ! NOTE: This is the HORIZONTAL area of
! each grid box covered by clouds
!@var conv input convective cloud cover in each model level(frac)
REAL*8 conv(nlev) ! NOTE: This is the HORIZONTAL area of
! each box covered by convective clouds
!@var dtau_s mean 0.67 micron optical depth of stratiform cloud level
!@var dtau_c mean 0.67 micron optical depth of convective cloud level
REAL*8 dtau_s(nlev), dtau_c(nlev)
! NOTE: these cloud optical depths are only for the
! cloudy part of the grid box, they are not weighted
! with the 0 cloud optical depth of the clear
! part of the grid box
INTEGER :: itrop !@var itrop tropopause level (WMO definition)
! The following input variables are used only if top_height = 1
REAL*8 skt !@var skt skin Temperature (K)
REAL*8 at(nlev) !@var at temperature in each model level (K)
!@var dem_s 10.5 micron longwave emissivity of stratiform clouds
!@var dem_c 10.5 micron longwave emissivity of convective clouds
REAL*8 dem_s(nlev),dem_c(nlev)
! ------
! Output
! ------
!@var fq_isccp the fraction of the model grid box covered by each
!@+ of the 49 ISCCP D level cloud types
REAL*8 fq_isccp(7,7)
!@var ctp cloud top pressure averaged over grid box (mb)
REAL*8 ctp
!@var tauopt cloud optical thickness averaged over grid box
REAL*8 tauopt
!@var nbox number of boxes with clouds (used as a flag)
INTEGER nbox
! ------
! Working variables
! ------
!@var frac_out boxes gridbox divided up into
REAL*8 frac_out(ncol,nlev) ! Equivalent of BOX in orig. version,
! but indexed by column then row,
! rather than by row then column
!@var tca total cloud cover (fraction) with extra layer of zeros on top
REAL*8 tca(ncol,0:nlev) ! in this version this just contains the
! values input from cc but replicated across ncol
!@var cca convect. cloud cover each model level(frac)
REAL*8 cca(ncol,nlev) ! from conv but replicated across ncol
!@var threshold pointer to position in gridbox
!@var maxocc Flag for max overlapped conv cld
!@var maxpsc Flag for max overlapped strat cld
!@var boxpos ordered pointer to position in gridbox
!@var threshold_min min. allowed value of threshold
REAL*8, dimension(ncol) :: threshold,maxocc,maxosc,boxpos,
* threshold_min
!@var dem,bb,bbs working variables for 10.5 micron longwave emissivity
!@+ in part of gridbox under consideration
real*8 dem(ncol),bb(nlev),bbs
!@var dtautmp temporary variables for dtau,dem
real*8, dimension(ncol) :: dtautmp
real*8 ptrop,attrop,atmax,atmin,btcmin,transmax
integer ilev,ibox,ipres,itau
integer acc(nlev,ncol),match(nlev-1),nmatch,levmatch(ncol)
!variables needed for water vapor continuum absorption
real*8 fluxtop_clrsky,trans_layers_above_clrsky,taumin
real*8 dem_wv(nlev)
real*8 press, dpress, atmden, rvh20, wk, rhoave, rh20s, rfrgn
real*8 tmpexp,tauwv, XX
real*8, dimension(ncol) :: tau,tb,ptop,emcld,fluxtop,
* trans_layers_above
real*8, parameter :: isccp_taumin = 0.1d0
INTEGER JERR
character*1 cchar(6),cchar_realtops(6)
data cchar / ' ','-','1','+','I','+'/
data cchar_realtops / ' ',' ','1','1','I','I'/
JERR=0
! assign 2d tca array using 1d input array cc
do ilev=0,nlev
do ibox=1,ncol
if (ilev.eq.0) then
tca(ibox,ilev)=0
else
tca(ibox,ilev)=cc(ilev)
endif
enddo
enddo
! assign 2d cca array using 1d input array conv
do ilev=1,nlev
do ibox=1,ncol
cca(ibox,ilev)=conv(ilev)
enddo
enddo
C**** In order to ensure that the model does not go down a different
C**** path depending on whether this routine is used, we save current
C**** seed and reset it afterwards
cc CALL RFINAL(SEED)
if (top_height .eq. 1) then
ptrop = pfull(itrop)
attrop = at(itrop)
atmin = 400.
atmax = 0.
do ilev=1,nlev-1
if (at(ilev) .gt. atmax) atmax=at(ilev)
if (at(ilev) .lt. atmin) atmin=at(ilev)
end do
end if
! ---------------------------------------------------!
! find unpermittable data.....
!
do ilev=1,nlev
if (cc(ilev) .lt. 0.) then
print*, ' error = cloud fraction less than zero'
JERR=1 ; return; ! stop
end if
if (cc(ilev) .gt. 1.) then
print*,' error = cloud fraction greater than 1'
JERR=1 ; return; ! stop
end if
if (conv(ilev) .lt. 0.) then
print*,' error = convective cloud fraction less than zero'
JERR=1 ; return; ! stop
end if
if (conv(ilev) .gt. 1.) then
print*,' error = convective cloud fraction greater than 1'
JERR=1 ; return; ! stop
end if
if (dtau_s(ilev) .lt. 0.) then
print*,' error = stratiform cloud opt. depth less than zero'
JERR=1 ; return; ! stop
end if
if (dem_s(ilev) .lt. 0.) then
print*,' error = stratiform cloud emissivity less than zero'
JERR=1 ; return; ! stop
end if
if (dem_s(ilev) .gt. 1.) then
print*,' error = stratiform cloud emissivity greater than 1'
JERR=1 ; return; ! stop
end if
if (dtau_c(ilev) .lt. 0.) then
print*, ' error = convective cloud opt. depth less than zero'
JERR=1 ; return; ! stop
end if
if (dem_c(ilev) .lt. 0.) then
print*,' error = convective cloud emissivity less than zero'
JERR=1 ; return; ! stop
end if
if (dem_c(ilev) .gt. 1.) then
print*,' error = convective cloud emissivity greater than 1'
JERR=1 ; return; ! stop
end if
end do
! ---------------------------------------------------!
! Initialise working variables
! ---------------------------------------------------!
! Initialised frac_out to zero
frac_out(:,:)=0.0
do ibox=1,ncol
boxpos(ibox)=(ibox-.5)/ncol
enddo
! ---------------------------------------------------!
! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
! frac_out is the array that contains the information
! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
! convective cloud
!loop over vertical levels
do ilev = 1,nlev
! Initialise threshold
if (ilev.eq.1) then
do ibox=1,ncol
! if max overlap
if (overlap.eq.1) then
! select pixels spread evenly
! across the gridbox
threshold(ibox)=boxpos(ibox)
else
! select random pixels from the non-convective
! part the gridbox ( some will be converted into
! convective pixels below )
threshold(ibox)=
* cca(ibox,ilev)+(1-cca(ibox,ilev))*randu(xx)
endif
enddo
endif
do ibox=1,ncol
! All versions
if (boxpos(ibox).le.cca(ibox,ilev)) then
maxocc(ibox) = 1
else
maxocc(ibox) = 0
end if
select case (overlap)
case (1) ! Max overlap
threshold_min(ibox)=cca(ibox,ilev)
maxosc(ibox)=1
case (2) ! Random overlap
threshold_min(ibox)=cca(ibox,ilev)
maxosc(ibox)=0
case (3) ! Max/Random overlap
threshold_min(ibox)=max(cca(ibox,ilev),
& min(tca(ibox,ilev-1),tca(ibox,ilev)))
if (threshold(ibox).lt.min(tca(ibox,ilev-1),tca(ibox,ilev))
& .and.(threshold(ibox).gt.cca(ibox,ilev))) then
maxosc(ibox)= 1
else
maxosc(ibox)= 0
end if
end select
! Reset threshold
threshold(ibox)=
!if max overlapped conv cloud
& maxocc(ibox) * boxpos(ibox) +
!else
& (1-maxocc(ibox)) * (
!if max overlapped strat cloud
& (maxosc(ibox)) * (
!threshold=boxpos
& threshold(ibox) ) +
!else
& (1-maxosc(ibox)) * (
!threshold_min=random[thrmin,1]
& threshold_min(ibox)+(1-threshold_min(ibox))*RANDU(XX)))
end do
! Fill frac_out with 1's where tca is greater than the threshold
do ibox=1,ncol
if (tca(ibox,ilev).gt.threshold(ibox)) then
frac_out(ibox,ilev)=1
else
frac_out(ibox,ilev)=0
end if
end do
! Code to partition boxes into stratiform and convective parts
! goes here
do ibox=1,ncol
if (threshold(ibox).le.cca(ibox,ilev)) then
! = 2 IF threshold le cca(ibox)
frac_out(ibox,ilev) = 2
else
! = the same IF NOT threshold le cca(ibox)
frac_out(ibox,ilev) = frac_out(ibox,ilev)
end if
end do
end do !loop over nlev
!
! ---------------------------------------------------!
! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
! put into vector tau
!initialize tau to zero
tau(:)=0.
!compute total cloud optical depth for each column
do ilev=1,nlev
!increment tau for each of the boxes
do ibox=1,ncol
if (frac_out(ibox,ilev).eq.1) then
dtautmp(ibox)= dtau_s(ilev)
else if (frac_out(ibox,ilev).eq.2) then
dtautmp(ibox)= dtau_c(ilev)
else
dtautmp(ibox)= 0.
end if
tau(ibox)=tau(ibox)+dtautmp(ibox)
end do
end do
! ---------------------------------------------------!
! COMPUTE INFRARED BRIGHTNESS TEMPERATURES
! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
!
! again this is only done if top_height = 1
!
! fluxtop is the 10.5 micron radiance at the top of the
! atmosphere
! trans_layers_above is the total transmissivity in the layers
! above the current layer
! fluxtop_clrsky and trans_layers_above_clrsky are the clear
! sky versions of these quantities.
if (top_height .eq. 1) then
!---------------------------------------------------------------
!
! DO CLEAR SKY RADIANCE CALCULATION FIRST
!
!compute water vapor continuum emissivity
!this treatment follows Schwarkzopf and Ramasamy
!JGR 1999,vol 104, pages 9467-9499.
!the emissivity is calculated at a wavenumber of 955 cm-1,
!or 10.47 microns
do ilev=1,nlev
!press and dpress are dyne/cm2 = Pascals *10
press = pfull(ilev)*10.
dpress = (phalf(ilev+1)-phalf(ilev))*10.
!atmden = g/cm2 = kg/m2 / 10
atmden = dpress*bygrav
rvh20 = qv(ilev)*bymrat !wtmair/wtmh20
wk = rvh20*avog*atmden/wtmair
c rhoave = (press/pstd)*(t0/at(ilev))
rhoave = (press/at(ilev))*t0bypstd
rh20s = rvh20*rhoave
rfrgn = rhoave-rh20s
tmpexp = exp(-0.02d0*(at(ilev)-t0))
tauwv = wk*1d-20*( & (3.41817d-7*rfrgn) )*0.98d0
dem_wv(ilev) = 1. - exp(-tauwv)
end do
!initialize variables
fluxtop_clrsky = 0.
trans_layers_above_clrsky=1.
do ilev=1,nlev
! Black body emission at temperature of the layer
bb(ilev)=1 / ( exp(pc1bylam/at(ilev)) - 1. )
! increase TOA flux by flux emitted from layer
! times total transmittance in layers above
fluxtop_clrsky = fluxtop_clrsky
& + dem_wv(ilev) * bb(ilev) * trans_layers_above_clrsky
! update trans_layers_above with transmissivity
! from this layer for next time around loop
trans_layers_above_clrsky=
& trans_layers_above_clrsky*(1.-dem_wv(ilev))
end do !loop over level
!add in surface emission
bbs=1/( exp(pc1bylam/skt) - 1. )
fluxtop_clrsky = fluxtop_clrsky + emsfc_lw * bbs
& * trans_layers_above_clrsky
! END OF CLEAR SKY CALCULATION
!
!----------------------------------------------------------------
fluxtop(:)=0.
trans_layers_above(:)=1.
do ilev=1,nlev
do ibox=1,ncol
! emissivity for point in this layer
if (frac_out(ibox,ilev).eq.1) then
dem(ibox)= 1. -
& ( else if (frac_out(ibox,ilev).eq.2) then
dem(ibox)= 1. -
& ( else
dem(ibox)= dem_wv(ilev)
end if
! increase TOA flux by flux emitted from layer
! times total transmittance in layers above
fluxtop(ibox) = fluxtop(ibox)
& + dem(ibox) * bb(ilev)
& * trans_layers_above(ibox)
! update trans_layers_above with transmissivity
! from this layer for next time around loop
trans_layers_above(ibox)=
& trans_layers_above(ibox)*(1.-dem(ibox))
end do ! ibox
end do