subroutine adnlmsas(im,ix,km,jcap,delt,del,sl,slk,rcs,& 2,25
slimsk,inew,xkt2,ncloud,kcheck,ps, &
t0,q0,cwm0,u0,v0,dot0, &
t1,q1,cwm1,u1,v1,rn1, &
cldwrk,kbot,ktop,jmin,kuo,kb, &
ad_t0,ad_q0,ad_cwm0,ad_u0,ad_v0,ad_dot0, &
ad_t1,ad_q1,ad_cwm1,ad_u1,ad_v1,ad_rn1, &
adjoint,mype)
!C$$$ SUBPROGRAM DOCUMENTATION BLOCK
!C . . . .
!C SUBPROGRAM: SASCNV COMPUTES CONVECTIVE HEATING AND MOISNG
!C PRGMMR: HUA-LU PAN ORG: W/NMC23 DATE: 92-03-01
!C
!C ABSTRACT: COMPUTES CONVECTIVE HEATING AND MOISTENING USING A ONE
!C CLOUD TYPE ARAKAWA-SCHUBERT CONVECTION SCHEME ORIGINALLY DEVELOPED
!C BY GEORG GRELL. THE SCHEME INCLUDES UPDRAFT AND DOWNDRAFT EFFECTS.
!C THE CLOSURE IS THE CLOUD WORK FUNCTION. BOTH UPDRAFT AND DOWNDRAFT
!C ARE ASSUMED TO BE SATURATED AND THE HEATING AND MOISTENING ARE
!C ACCOMPLISHED BY THE COMPENSATING ENVIRONMENT. THE NAME COMES FROM
!C "SIMPLIFIED ARAKAWA-SCHUBERT CONVECTION PARAMETERIZATION".
!C
!C PROGRAM HISTORY LOG:
!C 92-03-01 HUA-LU PAN
!C
!C USAGE: CALL SASCNV(IM,IX,KM,JCAP,DELT,DEL,SL,SLK,PS,QN,TN,
!C & Q1,T0,RN,KBOT,KTOP,KUO,SPD,SLIMSK)
!C
!C INPUT ARGUMENT LIST:
!C IM - INTEGER NUMBER OF POINTS
!C IX - LEADING DIMENSION OF QN,TN,Q1,T0,SPD
!C KM - INTEGER NUMBER OF LEVELS
!C JCAP - INTEGER SPECTRAL TRUNCATION
!C DT - REAL TIME STEP IN SECONDS
!C DEL - REAL (KM) SIGMA LAYER THICKNESS
!C SL - REAL (KM) SIGMA VALUES
!C SLK - REAL (KM) SIGMA VALUES TO THE KAPPA
!C PS - REAL (IM) SURFACE PRESSURE IN KILOPASCALS (CB)
!C QN - REAL (KM,IX) PREVIOUS SPECIFI!C HUMIDITY IN KG/KG
!C TN - REAL (KM,IX) PREVIOUS TEMPERATURE IN KELVIN
!C Q1 - REAL (KM,IX) CURRENT SPECIFI!C HUMIDITY IN KG/KG
!C T0 - REAL (KM,IX) CURRENT TEMPERATURE IN KELVIN
!C SPD - REAL (KM,IX) CURRENT WIND SPEED
!C SLIMSK - REAL (IM) LAND(1),SEA(0), ICE(2) FLAG
!C
!C OUTPUT ARGUMENT LIST:
!C Q1 - REAL (KM,IX) ADJUSTED SPECIFIC HUMIDITY IN KG/KG
!C T0 - REAL (KM,IX) ADJUSTED TEMPERATURE IN KELVIN
!C RN - REAL (IM) CONVECTIVE RAIN IN METERS
!C KBOT - INTEGER (IM) CLOUD BOTTOM LEVEL
!C KTOP - INTEGER (IM) CLOUD TOP LEVEL
!C KUO - INTEGER (IM) BIT FLAG INDICATING DEEP CONVECTION
!C
!C SUBPROGRAMS CALLED:
!C FESB - FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
!C
!C REMARKS: FUNCTION FESB IS INLINED BY FPP.
!C NONSTANDARD AUTOMATIC ARRAYS ARE USED.
!C
!C ATTRIBUTES:
!C LANGUAGE: FORTRAN 77.
!C MACHINE: CRAY.
!C
!C$$$
use type_kinds
, only: fp_kind
implicit none
include 'constant.h'
logical adjoint0
real(fp_kind) g,elocp,el2orc,eps,epsm1,terr,c0,delta,fact1,fact2
real(fp_kind) c1,c2,c3,p0,p1,p2,p3
real(fp_kind) tiny,qsmall
parameter (g=grav)
parameter(elocp=hvap/cp)
parameter(el2orc=hvap*hvap/(rv*cp),eps=rd/rv,epsm1=rd/rv-1.)
parameter(terr=0.,c0=.002,delta=0.6077338)
parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
parameter (qsmall=1.e-10)
parameter (tiny=1.e-15)
parameter (adjoint0=.false.)
! Constants for es(T) calculation
parameter (c1=0.61078,c2=17.67,c3=243.5)
! Constanst for precipitation efficiency calculation (equation below)
parameter (p0=1.591, p1=-0.639, p2=0.0953, p3=-0.00496)
! term = 1.e3*vshear/dz
! e1 = 1.591-.639*term+.0953*(term**2)-.00496*(term**3)
! e1 = p0 + p1 *term+ p2 * term**2 + p3 * term**3
integer,dimension(ix):: kb,kbcon,lmin,jmin,ktcon,kbdtr,kt2,kcheck,&
ktcon0,kbot,ktop,kuo
real(fp_kind),dimension(km):: del,slk,sl
real(fp_kind),dimension(ix):: ps,slimsk,cldwrk,rn1,rcs,ad_rn1
real(fp_kind),dimension(km,ix):: t0,q0,t1,q1,cwm1,u1,v1,dot0,&
u0,v0,cwm0,ad_t0,ad_q0,ad_cwm0,ad_u0,ad_v0,ad_dot0,ad_t1,ad_q1,&
ad_u1,ad_v1,ad_cwm1
!c
!c local variables and arrays
logical adjoint,totflg
logical,dimension(ix):: cnvflg,dwnflg,dwnflg2,flg
logical,dimension(km,ix):: flgk
integer km,jcap,im,ix,mype,i,k,indx
integer onemf,onemfu,inew,ncloud
integer kbmax,kmax,kbm,jmn
real(fp_kind) term1,dp,dqsdp,desdt,pprime,qs,a,ad_rfact
real(fp_kind) delt,alphas,betal,edtmaxs,mbdt,qevap
real(fp_kind) betas,pdpdwn,pdetrn,evfact,evfactl,dtmax,edtmaxl,dt2
real(fp_kind) dtmin,dvq1,dv3v,detau,dvu1,xtemp,xdby,dvv1,dellat
real(fp_kind) dv2q,dv3q,dv1q,dv1u,dv1v,dv2v,dv2u,dv3u,xqrch,delqev
real(fp_kind) alpha,factor,termq,esl,gamma
real(fp_kind) dqsdt,dt,po,dq,ad_term1,ad_term2,ad_evef,ad_dellat
real(fp_kind) ad_rain,ad_term3,ad_qcond,ad_qevap,qevap0,ad_term4
real(fp_kind) qevap1,dz,w3l,w2l,w4l,w2s,w1s,w1l,xlambu
real(fp_kind) alphal,fjcap,fkm,val,ad_es0,adt,es,term2,dlnsig,es0
real(fp_kind) w4s,w3s,ad_es,ad_qs,ad_desdt,term,ad_dvu1,ad_dvv1
real(fp_kind) ad_pprime,ad_etah,ad_xpw,ad_esl,ad_dqsdt,ad_dqsdp
real(fp_kind) ad_dv1u,ad_dv2u,ad_dv3q,ad_dv1q,ad_dv2q,ad_dv3u
real(fp_kind) ad_dv1,ad_dvq1,ad_dv3v,ad_dv1v,ad_dv2v,ad_qlk
real(fp_kind) ad_term5,ad_dhh,ad_term6,ad_xpwd,ad_gamma,ad_dh,ad_dqs
real(fp_kind) ad_ratio,ad_dz,ad_dg,ad_dt,ad_onemf,ad_termq,ad_factor
real(fp_kind) ad_xqc,ad_xqrch,ad_xtemp,ad_detad,ad_dz1,ad_xdby
real(fp_kind) ad_rfactg,x,funsqr,ftanh,dftanh,pw0
real(fp_kind) ad_qrch,ad_temp,ad_termv,ad_termu
real(fp_kind) ad_onemfu,ad_heol2,ad_qc,ad_fuv,ad_shear,rterm
real(fp_kind) ad_dv3,dv3,ad_dq,ad_term,ad_e1,ad_dv2,ad_detau,ad_tem1
real(fp_kind) ad_tem2,ad_vol2,ad_uol2
real(fp_kind) rain,rn0,qcond,ad_rn0,ad_delqev,delq2,evef,xpwd,w1,xqc
real(fp_kind) xpw,w2,ratio,fixed,w3,w4,qrch,dz1,vol2,temp,etah,rfact
real(fp_kind) term3,qlk,qc,tem1,tem2,dif,heol2,uol2,fuv,term4
real(fp_kind) edtmax,dhh,term6,detad,dg,dv1,dv2,aup,adw,termu,termv
real(fp_kind) term5,dqs,shear,dh,e1,beta
real(fp_kind) r1200,r3600,one,zero
real(fp_kind),dimension(15):: pcrit,acritt,acrit
real(fp_kind),dimension(ix):: pdot,acrtfct,acrtfct1,acrtfct0,&
dtconv0,dtconv,deltv,acrt,psfc,hmax,delq,hkbo,qkbo,pbcdif,&
hmin,pwavo,aa1,vshear,edt,edto,pwevo,hcdo,qcdo,ddp,pp2,adet,&
aatmp,xaa0,f,xk,xmb,edtx,xqcd,hsbar,xmbmax,xlamb,xlamd,&
xlamdet,xlamdet0,xlamdet1,delhbar,delqvar,deltbar,rntot,&
xkt2,ad_y1,ad_y2,ad_y3,ad_y4,ad_y5,ad_y6,ad_y7,ad_y8,ad_y9
real(fp_kind),dimension(ix):: ucdo,vcdo,ukbo,vkbo,qlko_ktcon,&
dellal,dellal1,ad_dellal1,ad_dellal,ad_ukbo,ad_vkbo,ad_ucdo,&
ad_vcdo,ad_qlko_ktcon,ad_xlamdet1,ad_xlamdet0,dpscl,tdpscl,&
dpcld,tdpcld,tcwfup,tcwfdn,tdetrn,tdpdwn,aa10,edt0,edto0,&
edtx0,xmb0,rn11,edto1,edto10,edtx1,edtx10,ad_hkbo,ad_xlamb,&
ad_xlamdet,ad_aa1,ad_pwavo,ad_tcwfup,ad_tcwfdn,ad_vshear,&
ad_edt,ad_edto,ad_edtx,ad_xlamd,ad_pwevo,ad_qcdo,ad_hcdo,&
ad_xaa0,ad_xpwav,ad_xpwev,ad_xhcd,ad_pdot,ad_acrtfct1,&
ad_acrtfct,ad_dtconv,ad_f,ad_xk,ad_xmb,ad_rntot,ad_edto1,&
ad_edtx1,xhkb,xqkb,xpwav,xpwev,xhcd
real(fp_kind),dimension(km,ix):: to2,qo2,to3,qo3,p,to,qo,uo,vo,&
qeso,tvo,dbyo,zo,heo,heso,qrck,dellah,dellaq,hco,qcko,eta,&
etad,qrcdo,pwo,pwdo,qrcd,hcko
real(fp_kind),dimension(km,ix):: qol,uol,vol,tol,heol,hesol,qesol,qol0,&
xqo,xto,xtvo,xzo,xqo0,xheo,xheso,xqeso,xqrcd,xhcko,xqcko,xtol,&
xqol,xqesol,xqcko0,xqrch0,eta0,hcko0,hckod,hcko1,qcko0,dqk,qcko00,&
xtemp0,qeso2,rn0k,delqevk,ad_to,ad_qo,ad_qeso,ad_tvo,&
ad_zo,ad_heo,ad_heso,ad_tol,ad_qol,ad_hesol,ad_heol,&
ad_qesol,ad_eta,ad_hcko,ad_qcko,ad_hcko1,&
ad_dbyo,ad_hckod,ad_pwo,ad_uo,ad_vo,ad_uol,&
ad_vol,ad_etad,ad_qrcdo,ad_pwdo,ad_dellah,ad_dellaq,ad_xto,&
ad_xqo,ad_xqeso,ad_xtvo,ad_xzo,ad_xtol,ad_xqol,ad_xqesol,&
ad_xheo,ad_xheso,ad_xhcko,ad_xqcko,ad_xqrcd,ad_to2
real(fp_kind),dimension(km,ix):: ad_qo2,ad_qeso2,ad_to3,ad_qo3,ad_x1,&
ad_x2,ad_x3,ad_x7,ad_x4,ad_x5,ad_x6,ad_x8,ad_x9,ad_x10,ad_x11,&
ad_x12,ad_x13,ad_x14,sumz,sumh,ad_cwmo,ucko,vcko,etau,etau0,&
ucko0,vcko0,uckod,vckod,ad_dellalk,dellalk,ad_cwmo2,dellau,&
dellav,uo2,vo2,cwmo,cwmo2,qo0,ad_qo0,ad_xqol0,ad_uo2,ad_vo2,&
xqol0,ad_xqo0,ad_ucko,ad_vcko,ad_uckod,ad_vckod,ad_dellau,&
ad_dellav,ad_etau,ad_sumz,ad_sumh,ad_qol0
save pcrit, acritt
data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
350.,300.,250.,200.,150./
data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
.3151,.3677,.41,.5255,.7663,1.1686,1.6851/
!C GDAS DERIVED ACRIT
!C DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688,
!C & .743,.813,.886,.947,1.138,1.377,1.896/
!C-----------------------------------------------------------------------
!c
!c
!c Define smoothing functions
funsqr(x,a,k) = 1.0 - 0.5**((x/a)**(2*k))
ftanh(x) = 0.5*(1.0 + tanh(x))
dftanh(x) = 0.5/(cosh(x)**2)
!c
do i = 1,ix
do k = 1,km
ad_cwmo(k,i)=0.0
ad_cwmo2(k,i)=0.0
ad_dellalk(k,i)=0.0
ad_t0(k,i) = 0.0
ad_q0(k,i) = 0.0
ad_u0(k,i) = 0.0
ad_v0(k,i) = 0.0
ad_cwm0(k,i)=0.0
ad_dot0(k,i) = 0.0
ad_tvo(k,i) = 0.0
ad_qeso(k,i) = 0.0
ad_qo(k,i) = 0.0
ad_to(k,i) = 0.0
ad_uo(k,i) = 0.0
ad_vo(k,i) = 0.0
ad_zo(k,i) = 0.0
ad_heso(k,i) = 0.0
ad_heo(k,i) = 0.0
ad_tol(k,i) = 0.0
ad_qol(k,i) = 0.0
ad_uol(k,i) = 0.0
ad_vol(k,i) = 0.0
ad_heol(k,i) = 0.0
ad_hesol(k,i) = 0.0
ad_qesol(k,i) = 0.0
ad_eta(k,i) = 0.0
ad_hcko(k,i) = 0.0
ad_qcko(k,i) = 0.0
ad_hcko1(k,i) = 0.0
ad_dbyo(k,i) = 0.0
ad_hckod(k,i) = 0.0
ad_qcko(k,i) = 0.0
ad_pwo(k,i) = 0.0
ad_etad(k,i) = 0.0
ad_qrcdo(k,i) = 0.0
ad_pwdo(k,i) = 0.0
ad_dellah(k,i) = 0.0
ad_dellaq(k,i) = 0.0
ad_xto(k,i) = 0.0
ad_xqo(k,i) = 0.0
ad_xqeso(k,i) = 0.0
ad_xtvo(k,i) = 0.0
ad_xzo(k,i) = 0.0
ad_xtol(k,i) = 0.0
ad_xqol(k,i) = 0.0
ad_xqesol(k,i) = 0.0
ad_xheo(k,i) = 0.0
ad_xheso(k,i) = 0.0
ad_xhcko(k,i) = 0.0
ad_xqcko(k,i) = 0.0
ad_xqrcd(k,i) = 0.0
ad_to2(k,i) = 0.0
ad_qo2(k,i) = 0.0
ad_qeso2(k,i) = 0.0
ad_to3(k,i) = 0.0
ad_qo3(k,i) = 0.0
ad_qo0(k,i)=0.0
ad_xqol0(k,i)=0.0
ad_uo2(k,i)=0.0
ad_vo2(k,i)=0.0
ad_xqo0(k,i)=0.0
ad_ucko(k,i)=0.0
ad_vcko(k,i)=0.0
ad_uckod(k,i)=0.0
ad_vckod(k,i)=0.0
ad_dellau(k,i)=0.0
ad_dellav(k,i)=0.0
ad_etau(k,i)=0.0
ad_sumz(k,i)=0.0
ad_sumh(k,i)=0.0
ad_qol0(k,i)=0.0
!c
p(k,i) = 0.0
to(k,i) = 0.0
qo(k,i) = 0.0
to2(k,i) = 0.0
qo2(k,i) = 0.0
qeso(k,i) = 0.0
qeso2(k,i) = 0.0
to3(k,i) = 0.0
qo3(k,i) = 0.0
tvo(k,i) = 0.0
dbyo(k,i) = 0.0
zo(k,i) = 0.0
heo(k,i) = 0.0
heso(k,i) = 0.0
xqo(k,i) = 0.0
xqeso(k,i) = 0.0
xto(k,i) = 0.0
xtvo(k,i) = 0.0
xzo(k,i) = 0.0
xheo(k,i) = 0.0
xheso(k,i) = 0.0
tol(k,i) = 0.0
qol(k,i) = 0.0
qesol(k,i) = 0.0
heol(k,i) = 0.0
hesol(k,i) = 0.0
uol(k,i) = 0.0
vol(k,i) = 0.0
xtol(k,i) = 0.0
xqol(k,i) = 0.0
xqesol(k,i) = 0.0
qrcd(k,i) = 0.0
dellah(k,i) = 0.0
dellaq(k,i) = 0.0
hcko(k,i) = 0.0
qcko(k,i) = 0.0
eta(k,i) = 0.0
hcko0(k,i) = 0.0
qcko0(k,i) = 0.0
qcko00(k,i)= 0.0
dqk(k,i) = 0.0
eta0(k,i) = 0.0
etad(k,i) = 0.0
xhcko(k,i) = 0.0
xqcko(k,i) = 0.0
xqcko0(k,i) = 0.0
xtemp0(k,i) = 0.0
xqrch0(k,i) = 0.0
qrcdo(k,i) = 0.0
pwo(k,i) = 0.0
pwdo(k,i) = 0.0
hcko1(k,i) = 0.0
hckod(k,i) = 0.0
xqrcd(k,i) = 0.0
rn0k(k,i) = 0.0
delqevk(k,i) = 0.0
flgk(k,i) = .true.
!!new
sumz(k,i) = 0.0
sumh(k,i) = 0.0
ucko(k,i) = 0.0
vcko(k,i) = 0.0
etau(k,i) = 0.0
uckod(k,i)= 0.0
vckod(k,i)=0.0
hckod(k,i)=0.0
end do
end do
do i = 1,ix
ad_xlamdet1(i)=0.0
ad_xlamdet0(i)=0.0
ad_ukbo(i)=0.0
ad_vkbo(i)=0.0
ad_ucdo(i)=0.0
ad_vcdo(i)=0.0
ad_qlko_ktcon(i)=0.0
ad_dellal(i)=0.0
ad_dellal1(i)=0.0
ad_hkbo(i) = 0.0
ad_xlamb(i) = 0.0
ad_xlamdet(i) = 0.0
ad_aa1(i) = 0.0
ad_pwavo(i) = 0.0
ad_tcwfup(i) = 0.0
ad_tcwfdn(i) = 0.0
ad_vshear(i) = 0.0
ad_edt(i) = 0.0
ad_edto(i) = 0.0
ad_edtx(i) = 0.0
ad_edto1(i) = 0.0
ad_edtx1(i) = 0.0
ad_xlamd(i) = 0.0
ad_pwevo(i) = 0.0
ad_qcdo(i) = 0.0
ad_hcdo(i) = 0.0
ad_xaa0(i) = 0.0
ad_xpwav(i) = 0.0
ad_xpwev(i) = 0.0
ad_xhcd(i) = 0.0
ad_pdot(i) = 0.0
ad_acrtfct1(i) = 0.0
ad_acrtfct(i) = 0.0
ad_dtconv(i) = 0.0
ad_f(i) = 0.0
ad_xk(i) = 0.0
ad_xmb(i) = 0.0
ad_rntot(i) = 0.0
dpscl(i) = 0.0
dpcld(i) = 0.0
tdpscl(i) = 0.0
tdpcld(i) = 0.0
tdetrn(i) = 0.0
tcwfup(i) = 0.0
tcwfdn(i) = 0.0
acrtfct(i) = 0.0
acrtfct0(i) = 0.0
acrtfct1(i) = 0.0
acrt(i) = 0.0
psfc(i) = 0.0
hmax(i) = 0.0
kb(i) = 0
hkbo(i) = 0.0
qkbo(i) = 0.0
kbcon(i) = 0
pbcdif(i) = 0.0
hmin(i) = 0.0
lmin(i) = 0
jmin(i) = 0
pwavo(i) = 0.0
aa1(i) = 0.0
vshear(i) = 0.0
edt(i) = 0.0
edt0(i) = 0.0
edto(i) = 0.0
edto1(i) = 0.0
edto10(i) = 0.0
pwevo(i) = 0.0
hcdo(i) = 0.0
qcdo(i) = 0.0
xhkb(i) = 0.0
xqkb(i) = 0.0
xpwav(i) = 0.0
xpwev(i) = 0.0
xhcd(i) = 0.0
xaa0(i) = 0.0
f(i) = 0.0
xk(i) = 0.0
xmb(i) = 0.0
xmb0(i) = 0.0
ktcon(i) = 0
edtx(i) = 0.0
edtx1(i) = 0.0
edtx10(i) = 0.0
xqcd(i) = 0.0
hsbar(i) = 0.0
xmbmax(i) = 0.0
xlamb(i) = 0.0
xlamdet(i) = 0.0
xlamdet0(i)=0.0
xlamdet1(i)=0.0
xlamd(i) = 0.0
kbdtr(i) = 0
pdot(i) = 0.0
cldwrk(i) = 0.0
kbot(i) = km+1
ktop(i) = 0
kuo(i) = 0
rn1(i) = 0.0
aa10(i) = 0.0
rn11(i) = 0.0
qlko_ktcon(i)=0.0
end do
! Initialize variables
r1200=1200.
r3600=3600.
zero=0.0
one=1.0
!c
!C INITIALIZE ARRAYS
!C
DO I = 1,IM
CNVFLG(I) = .TRUE.
dwnflg(i) = .true.
dwnflg2(i) = .true.
flg(i) = .true.
DTCONV(I) = 3600.
dtconv0(i) = dtconv(i)
xmbmax(i) = 0.1
ENDDO
DO K = 1, 15
ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
ENDDO
DT2 = DELT
dtmin = max(dt2,r1200)
dtmax = max(dt2,r3600)
!C MODEL TUNABLE PARAMETERS ARE ALL HERE
MBDT = 10.
EDTMAXl = 0.3
EDTMAXs = 0.3
ALPHAl = 0.5
ALPHAs = 0.5
BETAl = 0.15
betas = 0.15
BETAl = 0.05
betas = 0.05
!c EVEF = 0.07
!! evfact = 0.3
!! evfactl = 0.3
evfact = 0.7
evfactl = .7
PDPDWN = 0.
PDETRN = 200.
xlambu = 1.e-4
fjcap = (float(jcap) / 126.) ** 2
val = 1.
fjcap = max(fjcap,val)
!c W1l = -2.E-3 * (JCAP / 62.) **2
!c W2l = -1.E-2 * (JCAP / 62.) **2
!c W3l = -1.E-2 * (JCAP / 62.) **2
!c W4l = -1.E-3 * (JCAP / 62.) **2
!c W1s = -1.E-4 * (JCAP / 62.)
!c W2s = -1.E-3 * (JCAP / 62.)
!c W3s = -1.E-3 * (JCAP / 62.)
!c W4s = -1.E-4 * (JCAP / 62.)
fkm = (float(km) / 28.) ** 2
fkm = max(fkm,one)
W1l = -8.E-3
W2l = -4.E-2
W3l = -5.E-3
W4l = -5.E-4
W1s = -2.E-4
W2s = -2.E-3
W3s = -1.E-3
W4s = -2.E-5
!C
!C DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
!C AND THE MAXIMUM THETAE FOR UPDRAFT
!C
KBMAX = KM
KBM = KM
KMAX = KM
DO K = 1, KM
IF (SL(K) .GT. 0.45) KBMAX = K + 1
IF (SL(K) .GT. 0.7) KBM = K + 1
IF (SL(K) .GT. 0.04) KMAX = K + 1
ENDDO
!C kbm=10, kbmax=14, kmax=25 (=23 if sl(k).gt.06)
!C
!C CONVERT SURFACE PRESSURE TO MB FROM CB
!C
DO I = 1, IM
PSFC(I) = PS(I) * 10.
ENDDO
do i = 1,im
do k = 1,km
P(k,i) = PSFC(I) * SL(K)
to(k,i) = t0(k,i)
!! qo(k,i) = q0(k,i)
if (q0(k,i)>0.) then
qo(k,i) = q0(k,i)
else
qo(k,i) = qsmall
endif
uo(k,i) = u0(k,i)
vo(k,i) = v0(k,i)
cwmo(k,i)= cwm0(k,i)
qo0(k,i) = q0(k,i)
ENDDO
ENDDO
!C
!C COLUMN VARIABLES
!C P IS PRESSURE OF THE LAYER (MB)
!C T IS TEMPERATURE AT T-DT (K)..TN
!C Q IS MIXING RATIO AT T-DT (KG/KG)..QN
!C TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
!C QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
!C
do i = 1,im
DO K = 1, KMAX
call adfpvsx(to(k,i),es0,adt,ad_es0,adjoint0)
es = 10*es0
QESO(k,i) = EPS*es / (
!new code bounds qeso to be >= 1.e-8, qo>=1.e-10
TVO(k,i) = TO(k,i) + DELTA * TO(k,i) * QO(k,i)
ENDDO
ENDDO
!C
!C HYDROSTATIC HEIGHT ASSUME ZERO TERR
!C
DLNSIG = LOG(SL(1))
DO I = 1, IM
ZO(1,i) = TERR - DLNSIG * RD / G * TVO(1,i)
ENDDO
!cround
!c Reverse order of this loop leads to slightly different
!c final results. Change order of loop in sascnv2 to
!c restore similar results.
do i = 1,im
DO K = 2, KMAX
DLNSIG = LOG(SL(K) / SL(K-1))
term1 = dlnsig * rd / g
term2 = 0.5 * (tvo(k,i) + tvo(k-1,i))
zo(k,i) = zo(k-1,i) - term1*term2
ENDDO
ENDDO
!c
!C COMPUTE MOIST STATIC ENERGY
do i = 1,im
DO K = 1, KMAX
HEO(k,i) = G * ZO(k,i) + CP * TO(k,i) + HVAP * QO(k,i)
HESO(k,i) = G * ZO(k,i) + CP * TO(k,i) + HVAP * QESO(k,i)
ENDDO
ENDDO
!C
!C DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
!C THIS IS THE LEVEL WHERE UPDRAFT STARTS
!C
DO I = 1, IM
HMAX(I) = HEO(1,i)
KB(I) = 1
ENDDO
do i = 1,im
DO K = 2, KBM
IF (HEO(k,i).GT.HMAX(I)) THEN
KB(I) = K
HMAX(I) = HEO(k,i)
ENDIF
ENDDO
ENDDO
!C
!C SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
!C
!! DO I = 1, IM
!! HMIN(I) = HESO(1,i)
!! LMIN(I) = KBMAX
!! JMIN(I) = KBMAX
!! ENDDO
!! do i = 1,im
!! DO K = 2, KBMAX
!! IF (HESO(k,i).LT.HMIN(I)) THEN
!! LMIN(I) = K + 1
!! HMIN(I) = HESO(k,i)
!! ENDIF
!! ENDDO
!! ENDDO
!c
do i = 1,im
DO K = 1, KMAX - 1
DZ = 0.5 * (ZO(k+1,i) - ZO(k,i))
DP = 0.5 * ( call adfpvsx(to(k+1,i),es0,adt,ad_es0,adjoint0)
es = 10*es0
PPRIME = P(k+1,i) + EPSM1 * ES
QS = EPS * ES / PPRIME
DQSDP = - QS / PPRIME
DESDT = ES * (FACT1 / TO(k+1,i) + FACT2 / (TO(k+1,i)**2))
DQSDT = QS * P(k+1,i) * DESDT / (ES * PPRIME)
GAMMA = EL2ORC * QESO(k+1,i) / (TO(k+1,i)**2)
DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
DQ = DQSDT * DT + DQSDP * DP
TOL(k,i) = TO(k+1,i) + DT
QOL0(k,i) = QO(k+1,i) + DQ
if (qol0(k,i)>0.) then
qol(k,i) = qol0(k,i)
else
qol(k,i) = qsmall
endif
end do
end do
!c
do i = 1,im
do k = 1,kmax-1
PO = 0.5 * ( call adfpvsx(tol(k,i),es0,adt,ad_es0,adjoint0)
esl = 10.*es0
QESOL(k,i) = EPS*esl / (PO + EPSM1*esl)
HEOL(k,i) = 0.5 * G * (ZO(k,i) + ZO(k+1,i)) + &
CP * TOL(k,i) + HVAP * QOL(k,i)
HESOL(k,i) = 0.5 * G * (ZO(k,i) + ZO(k+1,i)) + &
CP * TOL(k,i) + HVAP * QESOL(k,i)
uol(k,i) = 0.5 * (uo(k,i) + uo(k+1,i))
vol(k,i) = 0.5 * (vo(k,i) + vo(k+1,i))
ENDDO
ENDDO
k = kmax
do i = 1, im
heol(k,i) = heo(k,i)
hesol(k,i) = heso(k,i)
enddo
!C
!C LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
!C
DO I = 1, IM
INDX = KB(I)
HKBO(I) = HEOL(INDX,I)
ukbo(i) = uol(indx,i)
vkbo(i) = vol(indx,i)
KBCON(I) = KMAX
flg(i) = .true.
ENDDO
do i = 1,im
DO K = 1, KBMAX
IF (FLG(I) .AND. K.GT.KB(I)) THEN
HSBAR(I) = HESOL(k,i)
IF (HKBO(I).GT.HSBAR(I)) THEN
FLG(I) = .FALSE.
KBCON(I) = K
ENDIF
ENDIF
ENDDO
ENDDO
DO I = 1, IM
dpscl(I) = -P(KBCON(I),I) + P(KB(I),I)
!! if (abs(dpscl(i)).gt.tiny) then
!! tdpscl(i) = funsqr(150.,dpscl(i),1)
!! else
!! tdpscl(i) = 1
!! endif
tdpscl(i) = 1.0
PDOT(I) = 10.* DOT0(KBCON(I),I)
IF (dpscl(I).GT.150.) tdpscl(i) = 0.0
IF (KBCON(I).EQ.KMAX) tdpscl(i) = 0.0
ENDDO
!c
!C FOUND LFC, CAN DEFINE REST OF VARIABLES
!C
!C DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
!C
DO I = 1, IM
alpha = alphas
if (nint(slimsk(i)) .eq. 1) alpha = alphal
IF (KB(I).EQ.1) THEN
DZ = 0.5 * (ZO(KBCON(I),I) + ZO(KBCON(I)-1,I)) - ZO(1,i)
ELSE
DZ = 0.5 * (ZO(KBCON(I),I) + ZO(KBCON(I)-1,I)) &
- 0.5 * (ZO(KB(I),I) + ZO(KB(I)-1,I))
ENDIF
IF (KBCON(I).NE.KB(I)) THEN
XLAMB(I) = -LOG(ALPHA) / DZ
ELSE
XLAMB(I) = 0.
ENDIF
ENDDO
!c
!C DETERMINE UPDRAFT MASS FLUX
do i = 1,im
DO K = 1, KMAX
ETA(k,i) = 1.
etau(k,i)= 1.
ENDDO
ENDDO
do i = 1,im
DO K = KBMAX, 2, -1
if (k.lt.kbcon(i) .and. k.ge.kb(i)) then
DZ = 0.5 * (ZO(k+1,i) - ZO(k-1,i))
ETA(k,i) = ETA(k+1,i) * EXP(-XLAMB(I) * DZ)
etau(k,i)= eta(k,i)
ENDIF
ENDDO
ENDDO
DO I = 1, IM
if (kb(i).eq.1 .and. kbcon(i).gt.1) then
DZ = 0.5 * (ZO(2,i) - ZO(1,i))
ETA(1,i) = ETA(2,i) * EXP(-XLAMB(I) * DZ)
etau(1,i)= eta(1,i)
ENDIF
ENDDO
!C
!C WORK UP UPDRAFT CLOUD PROPERTIES
!C
DO I = 1, IM
INDX = KB(I)
HCKO(INDX,I) = HKBO(I)
ucko(indx,i) = ukbo(i)
vcko(indx,i) = vkbo(i)
ENDDO
!C
!C CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
!C
do i = 1,im
DO K = 2, KMAX - 1
if (k.gt.kb(i) .and. k.le.kbcon(i)) then
FACTOR = ETA(k-1,i) / ETA(k,i)
ONEMF = 1. - FACTOR
HCKO(k,i) = FACTOR * HCKO(k-1,i) + ONEMF * &
0.5 * (HEOL(k,i) + HEOL(k+1,i))
UCKO(k,i) = FACTOR * UCKO(k-1,i) + ONEMF * &
0.5 * (uol(k,i) + uol(k+1,i))
VCKO(k,i) = FACTOR * VCKO(k-1,i) + ONEMF * &
0.5 * (vol(k,i) + vol(k+1,i))
ENDIF
if (k.gt.kbcon(i)) then
hcko(k,i) = hcko(kbcon(i),i)
ucko(k,i) = ucko(kbcon(i),i)
vcko(k,i) = vcko(kbcon(i),i)
ENDIF
ENDDO
ENDDO
!c
!c Load eta into new work array to preserve previous value
do i = 1,im
do k = 1,km
eta0(k,i) = eta(k,i)
hcko0(k,i) = hcko(k,i)
ucko0(k,i) = ucko(k,i)
vcko0(k,i) = vcko(k,i)
etau0(k,i) = etau(k,i)
end do
end do
!c
!C DETERMINE CLOUD TOP
DO I = 1, IM
FLG(I) = .true.
KTCON(I) = 1
ENDDO
do i = 1,im
DO K = 2, KMAX
dif = hcko(k,i) - hesol(k,i)
if (dif.lt.0. .and. flg(i) .and. k.gt.kbcon(i)) then
KTCON(I) = K
FLG(I) = .FALSE.
ENDIF
ENDDO
ENDDO
!C
!C Make sure that jmin is within the cloud
!C
!new
!c Search for downdraft origin level above theta-e minimum
do i=1,im
HMIN(i) = HEOl(kbcon(i),i)
LMIN(i) = KBMAX
JMIN(i) = KBMAX
if(cnvflg(i)) then
DO K = kbcon(i), KBMAX
IF(HEOl(k,i).LT.HMIN(i)) THEN
LMIN(i) = K + 1
HMIN(i) = HEOl(k,i)
ENDIF
ENDDO
endif
end do
!c Make sure jmin is within the cloud
DO I = 1, IM
JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
if (ktcon(i).lt.kbcon(i)) jmin(i) = 1
ENDDO
!c
DO I = 1, IM
dpcld(i) = p(kbcon(i),i) - p(ktcon(i),i)
!! tdpcld(i) = funsqr(dpcld(i),150.,1)
tdpcld(i) = 1.0
IF (dpcld(i).LT.150.) tdpcld(i) = 0.0
ENDDO
!!new
! KTCON is model level of deepest cloud in ensemble.
! Now randomly select cloud from ensemble of clouds
! ranging from one level above the downdraft origin to
! to level of neutral buoyancy (ktcon)
do i=1,im
do k=2,kmax-1
if (k.gt.jmin(i) .and. k.le.ktcon(i)) then
sumz(k,i) = sumz(k-1,i) + .5 * (zo(k+1,i) - zo(k-1,i))
sumh(k,i) = sumh(k-1,i) + .5 * (zo(k+1,i) - zo(k-1,i)) &
* heol(k,i)
endif
end do
end do
! Select cloud from ensemble
do i=1,im
kt2(i) = nint(xkt2(i)*float(ktcon(i)-jmin(i))-.5) + jmin(i) + 1
tem1 = hcko(jmin(i),i) - hesol(kt2(i),i)
tem2 = sumz(kt2(i),i) * hesol(kt2(i),i) - sumh(kt2(i),i)
if (abs(tem2) .gt. 0.000001) then
xlamdet0(i) = tem1 / tem2
else
xlamdet0(i) = 0.0
cnvflg = .false.
endif
if (xlamdet0(i).ge.0.) then
xlamdet1(i) = xlamdet0(i)
else
xlamdet1(i) = 0.0
endif
if (sumz(kt2(i),i) .gt. 0.000001) then
term = 2.3/sumz(kt2(i),i)
if (xlamdet1(i).lt.term) then
xlamdet(i) = xlamdet1(i)
else
xlamdet(i) = term
endif
else
xlamdet(i) = xlamdet1(i)
endif
end do
! Apply logic to determine if convection is possible
do i=1,im
if (kt2(i).ge.ktcon(i)) dwnflg(i) = .false.
if(xlamdet(i).le.1.e-30 .or. &
hcko(jmin(i),i)-hesol(kt2(i),i).le.1.e-30) &
dwnflg(i) = .false.
do k = jmin(i), kt2(i)
if(dwnflg(i) .and. &
heol(k,i).gt.hesol(kt2(i),i)) dwnflg(i)=.false.
enddo
end do
!C
!C DETRAINING CLOUD
!C
DO I = 1, IM
!! if (abs(dpcld(i)).gt.tiny) then
!! tdetrn(i) = funsqr(pdetrn,dpcld(i),1)
!! else
!! tdetrn(i) = 1.0
!! endif
tdetrn(i) = 1.0
IF (dpcld(i).LT.PDPDWN) DWNFLG2(I)=.FALSE.
!! if (abs(pdpdwn).gt.tiny) then
!! tdpdwn(i) = funsqr(dpcld(i),pdpdwn,1)
!! else
!! tdpdwn(i) = 1.0
!! endif
tdpdwn(i) = 1.0
if (dwnflg2(i) .and. jmin(i).le.kbcon(i)) then
dwnflg2(i) = .false.
tdetrn(i) = 0.0
tdpdwn(i) = 0.0
endif
ENDDO
!c
!C CLOUD PROPERTY ABOVE CLOUD TOP IS MODIFIED BY THE DETRAINMENT PROCESS
do i=1,im
if (dwnflg(i)) then
do k=2,kmax-1
if (k.gt.jmin(i) .and. k.le.kt2(i)) then
DZ = .5 * (ZO(k+1,i) - ZO(k-1,i))
!c to simplify matter, we will take the linear approach here
eta(k,i) = eta(k-1,i) * (1. + xlamdet(i) * dz)
etau(k,i) = etau(k-1,i) * (1.+(xlamdet(i)+xlambu)* dz)
endif
end do
else
do k=2,kmax-1
if (k.gt.jmin(i) .and. k.le.ktcon(i)) then
DZ = .5 * (ZO(k+1,i) - ZO(k-1,i))
etau(k,i) = etau(k-1,i) * (1. + xlambu * dz)
endif
end do
endif
ENDDO
do i=1,im
if (dwnflg(i)) then
ktcon0(i) = ktcon(i)
ktcon(i) = kt2(i)
else
ktcon0(i) = ktcon(i)
endif
end do
!c
!c Load hcko at cloud base into hckod
do i = 1,im
if (dwnflg(i)) then
indx = kbcon(i)
hckod(indx,i) = hcko(indx,i)
uckod(indx,i) = ucko(indx,i)
vckod(indx,i) = vcko(indx,i)
endif
end do
!c
!c Adjustment to hckod by detraining cloud
do i = 1,im
do k = 2,kmax-1
if (k.gt.kbcon(i) .and. k.le.ktcon(i)) then
FACTOR = ETA(k-1,i) / ETA(k,i)
ONEMF = 1. - FACTOR
fuv = etau(k-1,i)/etau(k,i)
onemfu = 1. - fuv
heol2 = 0.5*(heol(k,i) + heol(k+1,i))
HCKOd(k,i) = FACTOR*HCKOd(k-1,i) + ONEMF*heol2
uol2 = 0.5*(uol(k,i)+uol(k+1,i))
uckod(k,i) = fuv * uckod(k-1,i) + onemfu*uol2
vol2 = 0.5*(vol(k,i)+vol(k+1,i))
vckod(k,i) = fuv * vckod(k-1,i) + onemfu*vol2
ENDIF
end do
end do
!c
!c Load hckod and hcko into hcko1 and calculate dbyo
do k = 2,kmax-1
do i = 1,im
if (k.gt.kb(i).and.k.le.kbcon(i)) then
hcko1(k,i) = hcko(k,i)
elseif (k.gt.kbcon(i).and.k.le.ktcon(i)) then
if (.not.dwnflg(i)) then
hcko1(k,i) = hcko(k,i)
else
hcko1(k,i) = hckod(k,i)
endif
elseif (k.gt.ktcon(i)) then
hcko1(k,i) = hcko(k,i)
endif
if (k.gt.kb(i)) then
dbyo(k,i) = hcko1(k,i) - hesol(k,i)
endif
end do
end do
!C
!C COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
!C
DO I = 1, IM
indx = kb(i)
qcko(indx,i) = qol(indx,i)
ENDDO
do k = 1,km
do i = 1,im
qcko0(k,i) = qcko(k,i)
end do
end do
do i = 1,im
DO K = 1, KMAX
if (k.gt.kb(i) .and. k.lt.ktcon(i)) then
FACTOR = ETA(k-1,i) / ETA(k,i)
ONEMF = 1. - FACTOR
temp = FACTOR * QCKO(k-1,i) + ONEMF * &
0.5 * (QOL(k,i) + QOL(k+1,i))
qcko0(k,i) = temp
GAMMA = EL2ORC * QESOL(k,i) / (TOL(k,i)**2)
QRCH = QESOL(k,i) &
+ GAMMA * DBYO(k,i) / (HVAP * (1. + GAMMA))
dq = eta(k,i) * (temp - qrch)
dqk(k,i) = dq
!C
!C BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
!C
IF (DQ.GT.0.) THEN
dz = 0.5 * (zo(k+1,i) - zo(k-1,i))
dz1 = (zo(k,i) - zo(k-1,i))
ETAH = 0.5 * (ETA(k,i) + ETA(k-1,i))
QLK = DQ / (ETA(k,i) + ETAH * C0 * DZ)
AA1(I) = AA1(I) - DZ1 * G * QLK
PWO(k,i) = ETAH * C0 * DZ * QLK
PWAVO(I) = PWAVO(I) + PWO(k,i)
qc = qlk + qrch
QCKO(k,i) = QC
else
qcko(k,i) = temp
ENDIF
ENDIF
ENDDO
ENDDO
! This section is for cloud liquid water
if (ncloud.gt.0) then
do i=1,im
k = ktcon(i)
if (cnvflg(i)) then
GAMMA = EL2ORC * QESOl(K,i) / (TOl(K,i)**2)
QRCH = QESOl(K,i) &
+ GAMMA * DBYO(K,i) / (HVAP * (1. + GAMMA))
DQ = QCKO(K-1,i) - QRCH
qcko00(k-1,i) = qcko(k-1,i)
!C
!C CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
!C
IF(DQ.GT.0.) THEN
qlko_ktcon(i) = dq
QCKO(K-1,i) = QRCH
ENDIF
endif
end do
endif
!C
!C CALCULATE CLOUD WORK FUNCTION AT T+DT
!C
do i = 1,im
DO K = 1, KMAX
if (k.gt.kbcon(i) .and. k.le.ktcon(i)) then
DZ1 = ZO(k,i) - ZO(k-1,i)
GAMMA = EL2ORC * QESOL(k-1,i) / (TOL(k-1,i)**2)
RFACT = 1. + DELTA * CP * GAMMA &
* TOL(k-1,i) / HVAP
term1 = dz1*g/cp
term2 = 1./tol(k-1,i)
term3 = dbyo(k-1,i)
term4 = 1./(1.+gamma)
term5 = rfact
aa1(i) = aa1(i) + term1*term2*term3*term4*term5
dqs = qesol(k-1,i) - qol(k-1,i)
if (dqs.gt.0.) then
AA1(I) = AA1(I)+ DZ1 * G * DELTA * dqs
endif
ENDIF
ENDDO
ENDDO
!c
!c Set cloud mass flux threshold function based on updaft cwf.
do i = 1,im
term = 1.e3*aa1(i)
aa10(i) = aa1(i)
if (abs(term) .lt. 100.) then
tcwfup(i) = ftanh(term)
else
if (term .lt. -100.) then
tcwfup(i) = 0.0
else
tcwfup(i) = 1.0
endif
endif
end do
!C
!C------- DOWNDRAFT CALCULATIONS
!C
!C
!C--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
!C
do i = 1,im
DO K = 1, KMAX
if (k.ge.kb(i) .and. k.le.ktcon(i)) then
termu = (uol(k+1,i) - uol(k,i))**2
termv = (vol(k+1,i) - vol(k,i))**2
if (termu+termv.gt.tiny) then
shear = rcs(i) * sqrt(termu+termv)
else
shear = 0.0
endif
vshear(i)= vshear(i) + shear
ENDIF
ENDDO
ENDDO
DO I = 1, IM
dz = zo(ktcon(i),i)-zo(kb(i),i)
if (abs(dz).gt.tiny) then
term = 1.E3 * VSHEAR(I) / dz
e1 = p0 + p1*term + p2*term**2 + p3*term**3
edt(i) = 1.-e1
else
edt(i) = 0.
endif
edt0(i) = edt(i)
end do
!c
!c Impose upper and lower bounds on precipitation efficiency
do i = 1,im
if (edt(i).gt.0.9) then
edt(i) = 0.9
elseif (edt(i).lt.0.0) then
edt(i) = 0.0
endif
EDTO(I) = EDT(I)
EDTX(I) = EDT(I)
ENDDO
!c
!C DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
DO I = 1, IM
KBDTR(I) = KBCON(I)
beta = betas
if (nint(slimsk(i)) .eq. 1) beta = betal
KBDTR(I) = KBCON(I)
KBDTR(I) = MAX(KBDTR(I),1)
XLAMD(I) = 0.
IF (KBDTR(I).GT.1) THEN
DZ = 0.5 * ZO(KBDTR(I),I) + 0.5 * ZO(kbdtr(i)-1,i) &
- ZO(1,i)
XLAMD(I) = LOG(BETA) / DZ
ENDIF
ENDDO
!C DETERMINE DOWNDRAFT MASS FLUX
do i = 1,im
DO K = 1, KMAX
ETAD(k,i) = 1.
ENDDO
ENDDO
do i = 1,im
DO K = KBMAX, 2, -1
if (k.lt.kbdtr(i)) then
DZ = 0.5 * (ZO(k+1,i) - ZO(k-1,i))
ETAD(k,i) = ETAD(k+1,i) * EXP(XLAMD(I) * DZ)
ENDIF
ENDDO
ENDDO
K = 1
DO I = 1, IM
if (kbdtr(i).gt.1) then
DZ = 0.5 * (ZO(2,i) - ZO(1,i))
ETAD(k,i) = ETAD(k+1,i) * EXP(XLAMD(I) * DZ)
ENDIF
ENDDO
!C
!C--- DOWNDRAFT MOISTURE PROPERTIES
!C
do i = 1,im
do k = 1,km
qrcdo(k,i) = 0.0
pwdo(k,i) = 0.0
end do
end do
DO I = 1, IM
pwevo(i) = 0.
JMN = JMIN(I)
HCDO(I) = HEOL(JMN,I)
QRCDO(JMN,I) = QESOL(JMN,I)
ucdo(i) = uol(jmn,i)
vcdo(i) = vol(jmn,i)
ENDDO
do i = 1,im
DO K = KMAX-1, 1, -1
if (k.lt.jmin(i)) then
DQ = QESOL(k,i)
DT = TOL(k,i)
GAMMA = EL2ORC * DQ / DT**2
DH = HCDO(I) - HESOL(k,i)
term5 = (1./hvap) * gamma
term6 = 1./(1.+gamma)
qrcdo(k,i) = dq + term5*term6*dh
DETAD = ETAD(k+1,i) - ETAD(k,i)
if (k.lt.jmin(i)-1) then
term1 = etad(k+1,i)*qrcdo(k+1,i)
term2 = etad(k,i) * qrcdo(k,i)
term3 = detad*0.5*(qrcdo(k,i)+qrcdo(k+1,i))
pwdo(k,i) = term1 - term2 - term3
else
term1 = etad(k+1,i)*qol(k+1,i)
term2 = etad(k,i) * qrcdo(k,i)
term3 = detad*0.5*(qrcdo(k,i)+qesol(k+1,i))
pwdo(k,i) = term1 - term2 - term3
endif
PWEVO(I) = PWEVO(I) + PWDO(k,i)
ENDIF
ENDDO
ENDDO
do i = 1,im
qcdo(i) = qrcdo(1,i)
end do
!C
!C--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
!C--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
!C--- EVAPORATE (PWEV)
!C
DO I = 1, IM
edtmax = edtmaxl
if (nint(slimsk(i)) .eq. 0) edtmax = edtmaxs
IF (DWNFLG2(I)) THEN
IF (PWEVO(I).LT.0.) THEN
EDTO1(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
edto10(i) = edto1(i)
if (edto1(i).gt.edtmax) then
edto1(i) = edtmax
endif
ELSE
EDTO1(I) = 0.
ENDIF
ELSE
EDTO1(I) = 0.
ENDIF
ENDDO
!C
!C
!C--- DOWNDRAFT CLOUDWORK FUNCTIONS
!C
!C
!cround reverse order of loop --> get slightly different result
do i = 1,im
DO K = KMAX-1, 1, -1
IF (DWNFLG2(I) .AND. K.LT.JMIN(I)) THEN
GAMMA = EL2ORC * QESOL(k+1,i) / TOL(k+1,i)**2
DHH = HCDO(I)
DT = TOL(k+1,i)
DG = GAMMA
DH = HESOL(k+1,i)
dz = zo(k,i) - zo(k+1,i)
term1 = edto1(i)
term2 = dz
term3 = g/(cp*dt)
term4 = dhh-dh
term5 = 1./(1.+dg)
term6 = 1. + delta*cp*dg*dt/hvap
aa1(i) = aa1(i) &
+ term1*term2*term3*term4*term5*term6
dqs = qesol(k+1,i)-qol(k+1,i)
if (dqs.gt.0.) then
AA1(I)= AA1(I) + EDTO1(I)*DZ*G*DELTA*dqs
endif
ENDIF
ENDDO
ENDDO
!c
!c Compute threshold function based on final cloud work function
do i = 1,im
term = 1.e3*aa1(i)
if (abs(term) .lt. 100.) then
tcwfdn(i) = ftanh(term)
else
if (term .lt. -100.) then
tcwfdn(i) = 0.0
else
tcwfdn(i) = 1.0
endif
endif
end do
!C
!C
!C--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
!C--- WILL DO TO THE ENVIRONMENT?
!C
do i = 1,im
DO K = 1, KMAX
DELLAH(k,i) = 0.
DELLAQ(k,i) = 0.
dellau(k,i) = 0.
dellav(k,i) = 0.
ENDDO
ENDDO
DO I = 1, IM
if (ktcon(i).ne.1) then
DP = 100. * PSFC(I) * DEL(1)
DELLAH(1,i) = EDTO1(I) * ETAD(1,i) * (HCDO(I) &
- HEOL(1,i)) * G / DP
DELLAQ(1,i) = EDTO1(I) * ETAD(1,i) * (QCDO(I) &
- QOL(1,i)) * G / DP
DELLAU(1,i) = EDTO1(I) * ETAD(1,i) * (UCDO(i) &
- UOl(1,i)) * G / DP
DELLAV(1,i) = EDTO1(I) * ETAD(1,i) * (VCDO(i) &
- VOl(1,i)) * G / DP
endif
ENDDO
!C
!C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
!C
!cround reverse order of loop yields different results
do i = 1,im
DO K = 2, KMAX-1
if ( k.lt.ktcon(i) .and. ktcon(i).ne.1 ) then
AUP = 1.
IF (K.LE.KB(I)) AUP = 0.
ADW = 1.
IF (K.GT.JMIN(I)) ADW = 0.
DV1 = HEOL(k,i)
DV2 = 0.5 * (HEOL(k,i) + HEOL(k+1,i))
DV3 = HEOL(k-1,i)
DV1Q = QOL(k,i)
DV2Q = 0.5 * (QOL(k,i) + QOL(k+1,i))
DV3Q = QOL(k-1,i)
DV1u = uOL(k,i)
DV2u = 0.5 * (uOL(k,i) + uOL(k+1,i))
DV3u = uOL(k-1,i)
DV1v = vOL(k,i)
DV2v = 0.5 * (vOL(k,i) + vOL(k+1,i))
DV3v = vOL(k-1,i)
DP = 100. * PSFC(I) * DEL(K)
DETAu = ETA(k,i) - ETA(k-1,i)
DETAD = ETAD(k,i) - ETAD(k-1,i)
term1 = aup*eta(k,i) - adw*edto1(i)*etad(k,i)
term2 = aup*detau
term3 = aup*eta(k-1,i) - adw*edto1(i)*etad(k-1,i)
term4 = adw*edto1(i)*detad
termq = 0.5*(qrcdo(k,i)+qrcdo(k-1,i))
dellah(k,i) = dellah(k,i) + (g/dp) * &
(term1*dv1 - term3*dv3 - term2*dv2 + term4*hcdo(i))
dellaq(k,i) = dellaq(k,i) + (g/dp) * &
( term1*dv1q - term3*dv3q - term2*dv2q + term4*termq )
dellau(k,i) = dellau(k,i) + (g/dp) * &
(term1*dv1u - term3*dv3u - term2*dv2u + term4*ucdo(i))
dellav(k,i) = dellav(k,i) + (g/dp) * &
(term1*dv1v - term3*dv3v - term2*dv2v + term4*vcdo(i))
ENDIF
ENDDO
ENDDO
!C
!C------- CLOUD TOP
!C
DO I = 1, IM
if (ktcon(i).ne.1) then
INDX = KTCON(I)
DP = 100. * PSFC(I) * DEL(INDX)
DV1 = HEOL(INDX-1,I)
DELLAH(INDX,I) = ETA(INDX-1,I) * &
(HCKO1(INDX-1,I) - DV1) * G / DP
DVQ1 = QOL(INDX-1,I)
DELLAQ(INDX,I) = ETA(INDX-1,I) * &
(QCKO(INDX-1,I) - DVQ1) * G / DP
DVu1 = uOL(INDX-1,I)
DELLAu(INDX,I) = ETA(INDX-1,I) * &
(uCKO(INDX-1,I) - DVu1) * G / DP
DVv1 = vOL(INDX-1,I)
DELLAv(INDX,I) = ETA(INDX-1,I) * &
(vCKO(INDX-1,I) - DVv1) * G / DP
!cloud liquid water
dellal(i) = eta(indx-1,i) * qlko_ktcon(i) *g/dp
endif
ENDDO
!C
!C------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
!C
do i = 1,im
DO K = 1, KMAX
if ( k.le.ktcon(i) .and. ktcon(i).ne.1 ) then
xqo0(k,i) = DELLAQ(k,i) * MBDT + Qo0(k,i)
if (xqo0(k,i)>0.) then
xqo(k,i) = xqo0(k,i)
else
xqo(k,i) = qsmall
endif
DELLAT = (DELLAH(k,i) - HVAP * DELLAQ(k,i)) / CP
xto(k,i) = DELLAT * MBDT + To(k,i)
else
xqo(k,i) = Qo0(k,i)
xto(k,i) = To(k,i)
ENDIF
ENDDO
ENDDO
!c
!C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!C
!C--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
!C--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
!C--- WOULD HAVE ON THE STABILITY,
!C--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
!C--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
!C--- DESTABILIZATION.
!C
!C--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
!C
do i = 1,im
DO K = 1, KMAX
call adfpvsx(xto(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
xqeso(k,i) = EPS*es / ( xtvo(k,i) = xto(k,i) + DELTA * xto(k,i) * xqo(k,i)
ENDDO
ENDDO
!C
!C HYDROSTATIC HEIGHT ASSUME ZERO TERR
!C
DLNSIG = LOG(SL(1))
DO I = 1, IM
xzo(1,i) = TERR - DLNSIG * RD / G * xtvo(1,i)
ENDDO
do i = 1,im
DO K = 2, KMAX
DLNSIG = LOG(SL(K) / SL(K-1))
term1 = dlnsig * rd / g
term2 = 0.5 * (xtvo(k,i) + xtvo(k-1,i))
xzo(k,i) = xzo(k-1,i) - term1*term2
ENDDO
ENDDO
!C
!C--- MOIST STATIC ENERGY
!C
do i = 1,im
DO K = 1, KMAX - 1
DZ = 0.5 * (xzo(k+1,i) - xzo(k,i))
DP = 0.5 * ( call adfpvsx(xto(k+1,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
PPRIME = P(k+1,i) + EPSM1 * ES
QS = EPS * ES / PPRIME
DQSDP = - QS / PPRIME
DESDT = ES * (FACT1 / xto(k+1,i) + FACT2 / (xto(k+1,i)**2))
DQSDT = QS * P(k+1,i) * DESDT / (ES * PPRIME)
GAMMA = EL2ORC * xqeso(k+1,i) / (xto(k+1,i)**2)
DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
DQ = DQSDT * DT + DQSDP * DP
xtol(k,i) = xto(k+1,i) + DT
xqol0(k,i) = xqo(k+1,i) + DQ
if (xqol0(k,i)>0.) then
xqol(k,i) = xqol0(k,i)
else
xqol(k,i) = qsmall
endif
end do
end do
do i = 1,im
do k = 1,kmax-1
PO = 0.5 * ( call adfpvsx(xtol(k,i),es0,adt,ad_es0,adjoint0)
esl = 10.*es0
xqesol(k,i) = EPS*esl / (PO + EPSM1*esl)
xheo(k,i) = 0.5 * G * (xzo(k,i) + xzo(k+1,i)) + &
CP * xtol(k,i) + HVAP * xqol(k,i)
xheso(k,i) = 0.5 * G * (xzo(k,i) + xzo(k+1,i)) + &
CP * xtol(k,i) + HVAP * xqesol(k,i)
ENDDO
ENDDO
k = kmax
do i = 1, im
xheo(k,i) = G*xzo(k,i) + CP*xtol(k,i) + HVAP*xqol(k,i)
xheso(k,i) = G*xzo(k,i) + CP*xtol(k,i) + HVAP*xqesol(k,i)
enddo
DO I = 1, IM
INDX = KB(I)
xhcko(INDX,I) = xheo(INDX,I)
xqcko(INDX,I) = xqol(INDX,I)
ENDDO
!C
!C
!C**************************** STATIC CONTROL
!C
!C
!C------- MOISTURE AND CLOUD WORK FUNCTIONS
!C
do i = 1,im
DO K = 2, KMAX - 1
if (k.gt.kb(i) .and. k.le.ktcon(i)) then
FACTOR = ETA(k-1,i) / ETA(k,i)
ONEMF = 1. - FACTOR
xhcko(k,i) = FACTOR * xhcko(k-1,i) + ONEMF * &
0.5 * (xheo(k,i) + xheo(k+1,i))
ENDIF
ENDDO
ENDDO
!c
do i = 1,im
xaa0(i) = 0.0
xpwav(i) = 0.0
end do
do i = 1,im
do k = 2,kmax-1
if (k.gt.kb(i) .and. k.lt.ktcon(i)) then
!c
!c xtemp calculation
factor = eta(k-1,i)/eta(k,i)
onemf = 1. - factor
termq = 0.5*(xqol(k,i) + xqol(k+1,i))
xtemp = factor*xqcko(k-1,i) + onemf*termq
xtemp0(k,i) = xtemp
!c
!c xqrch calculation
xdby = xhcko(k,i) - xheso(k,i)
if (xdby.gt.0.) then
gamma = el2orc*xqesol(k,i)/xtol(k,i)**2
xqrch = xqesol(k,i) + gamma*xdby/(hvap*(1.+gamma))
else
xqrch = xqesol(k,i)
endif
xqrch0(k,i) = xqrch
!c
!c dq calculation
dq = eta(k,i)*(xtemp-xqrch)
!c
!c Actions dependent upon sign of dq
if (dq.gt.0.) then
dz = 0.5 * (xzo(k+1,i) - xzo(k-1,i))
dz1 = xzo(k,i) - xzo(k-1,i)
etah = 0.5 * (eta(k,i) + eta(k-1,i))
qlk = dq / (eta(k,i) + etah * c0 * dz)
xqc = qlk + xqrch
xqcko(k,i) = xqc
xaa0(i) = xaa0(i) - dz1*g*qlk
xpw = etah * c0 * dz * qlk
xpwav(i) = xpwav(i) + xpw
else
xqcko(k,i) = xtemp
endif
!c
endif
end do
end do
!c
do i = 1,im
do k = 2,kmax-1
if (k.gt.kbcon(i) .and. k.le.ktcon(i)) then
DZ1 = xzo(k,i) - xzo(k-1,i)
GAMMA = EL2ORC * xqesol(k-1,i) / (xtol(k-1,i)**2)
RFACT = 1. + DELTA * CP * GAMMA &
* xtol(k-1,i) / HVAP
XDBY = xhcko(k-1,i) - xheso(k-1,i)
term1 = dz1 * (g/cp)
term2 = 1./xtol(k-1,i)
term3 = xdby
term4 = 1./(1.+gamma)
term5 = rfact
xaa0(i) = xaa0(i) &
+ term1*term2*term3*term4*term5
dqs = xqesol(k-1,i) - xqol(k-1,i)
if (dqs.gt.0.) then
XAA0(I) = XAA0(I)+ DZ1 * G * DELTA * dqs
endif
ENDIF
ENDDO
ENDDO
!C
!C------- DOWNDRAFT CALCULATIONS
!C
!C
!C--- DOWNDRAFT MOISTURE PROPERTIES
!C
DO I = 1, IM
XPWEV(I) = 0.
ENDDO
DO I = 1, IM
IF (DWNFLG2(I)) THEN
JMN = JMIN(I)
XHCD(I) = xheo(JMN,I)
XQRCD(JMN,I) = xqesol(JMN,I)
ENDIF
ENDDO
do i = 1,im
DO K = KMAX-1, 1, -1
IF (DWNFLG2(I) .AND. K.LT.JMIN(I)) THEN
dq = xqesol(k,i)
dt = xtol(k,i)
gamma = el2orc*dq/dt**2
dh = xhcd(i) - xheso(k,i)
term5 = (1./hvap) * gamma
term6 = 1./(1.+gamma)
xqrcd(k,i) = dq + term5*term6*dh
detad = etad(k+1,i) - etad(k,i)
term1 = etad(k+1,i)*xqrcd(k+1,i)
term2 = etad(k,i)*xqrcd(k,i)
term3 = detad * 0.5 * (xqrcd(k,i) + xqrcd(k+1,i))
xpwd = term1 - term2 - term3
XPWEV(I) = XPWEV(I) + XPWD
ENDIF
ENDDO
ENDDO
!c
!C
DO I = 1, IM
edtmax = edtmaxl
if (nint(slimsk(i)) .eq. 0) edtmax = edtmaxs
IF (DWNFLG2(I)) THEN
IF (XPWEV(I).GE.0.) THEN
EDTX1(I) = 0.
ELSE
EDTX1(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
edtx10(i) = edtx1(i)
if (edtx1(i).gt.edtmax) then
edtx1(i) = edtmax
endif
ENDIF
ELSE
EDTX1(I) = 0.
ENDIF
ENDDO
!C
!C
!C
!C--- DOWNDRAFT CLOUDWORK FUNCTIONS
!C
!C
!cround reverse order of loop yields slightly difference results
!c difference only in q and rn, not t
do i = 1,im
DO K = KMAX-1, 1, -1
IF (DWNFLG2(I) .AND. K.LT.JMIN(I)) THEN
GAMMA = EL2ORC * xqesol(k+1,i) / xtol(k+1,i)**2
DHH = XHCD(I)
DT = xtol(k+1,i)
DG = GAMMA
DH = xheso(k+1,i)
dz = xzo(k,i) - xzo(k+1,i)
term1 = edtx1(i)
term2 = dz
term3 = g/(cp*dt)
term4 = dhh-dh
term5 = 1./(1.+dg)
term6 = 1. + delta*cp*dg*dt/hvap
dqs = xqesol(k+1,i)-xqol(k+1,i)
xaa0(i) = xaa0(i) &
+ term1*term2*term3*term4*term5*term6
if (dqs.gt.0.) then
XAA0(I) = XAA0(I) + EDTX1(I)*DZ*G*DELTA*dqs
endif
ENDIF
ENDDO
ENDDO
!C
!C CALCULATE CRITICAL CLOUD WORK FUNCTION
!C
DO I = 1, IM
ACRT(I) = 0.
IF ( ACRT(I) = ACRIT(15)*(975.-P(KTCON(I),I)) &
/(975.-PCRIT(15))
ELSEIF ( ACRT(I) = ACRIT(1)
ELSE
k = int((850. - p(ktcon(i),i))/50.) + 2
K = MIN(K,15)
K = MAX(K,2)
ACRT(I) = ACRIT(K)+(ACRIT(K-1)-ACRIT(K))* &
( ENDIF
ENDDO
DO I = 1, IM
ACRTFCT(I) = 1.
if (nint(slimsk(i)) .eq. 1) then
w1 = w1l
w2 = w2l
w3 = w3l
w4 = w4l
else
w1 = w1s
w2 = w2s
w3 = w3s
w4 = w4s
endif
IF (PDOT(I).LE.W4) THEN
ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
ELSEIF (PDOT(I).GE.-W4) THEN
ACRTFCT(I) = -1*(PDOT(I) + W4) / (W4 - W3)
ELSE
ACRTFCT(I) = 0.
ENDIF
acrtfct0(i) = acrtfct(i)
!c
if (acrtfct(i).lt.-1.) then
acrtfct(i) = -1.
endif
if (acrtfct(i).gt.1.) then
acrtfct(i) = 1.
endif
!c
ACRTFCT1(I) = 1. - ACRTFCT(I)
!c
DTCONV(I) = DT2 + (1800. - DT2) * &
(PDOT(I) - W2) / (W1 - W2)
dtconv0(i) = dtconv(i)
if (dtconv(i).lt.dtmin) then
dtconv(i) = dtmin
endif
if (dtconv(i).gt.dtmax) then
dtconv(i) = dtmax
endif
ENDDO
!C
!C--- LARGE SCALE FORCING
!C
DO I= 1, IM
F(I) = (AA1(I) - ACRT(I) * ACRTFCT1(I)) / DTCONV(I)
XK(I) = (XAA0(I) - AA1(I)) / MBDT
!C
!C--- KERNEL, CLOUD BASE MASS FLUX
!C
if (abs(xk(i)).gt.tiny) then
ratio = -f(i)/xk(i)
else
ratio = 0.0
endif
if ( fixed = tdpscl(i)*tdpcld(i)*tdetrn(i)*tdpdwn(i)
xmb(i) = fixed * ratio*tcwfup(i)*tcwfdn(i)
else
xmb(i) = 0.0
endif
xmb0(i) = xmb(i)
if (xmb(i).gt.xmbmax(i)) then
xmb(i) = xmbmax(i)
endif
ENDDO
!c
!C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!C
!C--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
!C--- MULTIPLIED BY THE MASS FLUX NECESSARY TO KEEP THE
!C--- EQUILIBRIUM WITH THE LARGER-SCALE.
!C
!cround
!c reverse order of loop, different q and rn
do i = 1,im
do k = 1,km
dellal1(i) = 0.0
dellalk(k,i) = 0.0
dellalk(k,i) = dellal(i)
if (k.le.ktcon(i) .and. ktcon(i).ne.1) then
dellat = (dellah(k,i)-hvap*dellaq(k,i))/cp
to2(k,i) = to(k,i) + dellat*xmb(i)*dt2
qo2(k,i) = qo0(k,i) + dellaq(k,i)*xmb(i)*dt2
uo2(k,i) = uo(k,i) + dellau(k,i)*xmb(i)*dt2
vo2(k,i) = vo(k,i) + dellav(k,i)*xmb(i)*dt2
call adfpvsx(to2(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
qeso2(k,i) = eps*es / ( if (k.eq.ktcon(i) .and. ncloud.gt.0 .and. cnvflg(i)) then
dp = 100. * psfc(i) * del(k)
cwmo2(k,i) = cwmo(k,i) + dellalk(k,i) * xmb(i) * dt2
dellal1(i) = dellalk(k,i)*xmb(i)*dp/g
else
cwmo2(k,i) = cwmo(k,i)
dellal1(i) = dellalk(k,i)
endif
else
to2(k,i) = to(k,i)
qo2(k,i) = qo(k,i)
uo2(k,i) = uo(k,i)
vo2(k,i) = vo(k,i)
cwmo2(k,i)= cwmo(k,i)
call adfpvsx(to2(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
qeso2(k,i) = eps*es / ( dellal1(i) = dellalk(k,i)
ENDIF
dellal(i) = dellal1(i)
dellal1(i)= 0.0
ENDDO
ENDDO
!c
DO I = 1, IM
rntot(i) = 0.
ENDDO
!cround different results when reverse loop
do i = 1,im
DO K = KMAX, 1, -1
if (k.le.ktcon(i) .and. ktcon(i).ne.1) then
AUP = 1.
IF (K.Le.KB(I)) AUP = 0.
ADW = 1.
IF (K.GT.JMIN(I)) ADW = 0.
rain = AUP * PWO(k,i) + ADW * EDTO1(I) * PWDO(k,i)
RNtot(I) = RNtot(I) + rain * xmb(i) * .001 * dt2
ENDIF
ENDDO
ENDDO
!c
!cround very slight difference in q and r
do i = 1,im
rn0 = 0.0
delqev = 0.0
flg(i) = .true.
!c
do k = km,1,-1
qevap = 0.
qcond = 0.
delq2 = 0.
!c
!c Accumulate precipitation from cloud level
if (k.le.ktcon(i) .and. ktcon(i).ne.1) then
aup = 1.
if (k.le.kb(i)) aup = 0.
adw = 1.
if (k.gt.jmin(i)) adw = 0.
rain = aup*pwo(k,i) + adw*edto1(i)*pwdo(k,i)
rn0 = rn0 + rain*xmb(i)*.001*dt2
endif
!c
!c Save layer specific quantities for use in adjoint
rn0k(k,i) = rn0
flgk(k,i) = flg(i)
delqevk(k,i) = delqev
!c
!c Check if any falling precipitation will evaporate
if (flg(i) .and. k.le.ktcon(i)) then
!! if (nint(slimsk(i)) .eq. 1) then
!! evef = 0.07
!! else
!! evef = edt(i)*evfact
!! endif
evef = edt(i)*evfact
term1 = evef
term2 = qo2(k,i) - qeso2(k,i)
term3 = el2orc*qeso2(k,i)/to2(k,i)**2
term4 = 1. / (1.+term3)
qcond = term1*term2*term4
dp = 100. * psfc(i) * del(k)
!c
!c Compute amount of precipitation subject to evaporation
if (rn0 .gt.0. .and. qcond.lt.0.) then
qevap = -qcond * (1.-exp(-.32*sqrt(dt2*rn0)))
if ( qevap .gt. rn0*1000.*g/dp ) then
qevap = rn0*1000.*g/dp
endif
delq2 = delqev + .001 * qevap * dp / g
if (delq2.gt.rntot(i)) then
qevap = 1000.* g * (rntot(i)-delqev) / dp
flg(i) = .false.
endif
endif
!c
!c Adjust t and q profiles plus net rain in response
!c to evaporation of falling precipitation.
if (rn0 .gt.0. .and. qevap.gt.0.) then
qo3(k,i) = qo2(k,i) + qevap
to3(k,i) = to2(k,i) - elocp * qevap
rn0 = rn0 - .001 * qevap * dp / g
delqev = delqev + .001*dp*qevap/g
!c
!c Otherwise, no adjustment to t and q profiles
else
to3(k,i) = to2(k,i)
qo3(k,i) = qo2(k,i)
endif
!c
!c Above cloud top (or not cloud present).
!c Therefore, no changes to t and q.
else
to3(k,i) = to2(k,i)
qo3(k,i) = qo2(k,i)
endif
!c
enddo
!c
!c Save net surface precipitation.
rn1(i) = rn0
enddo
!c
!c
!c Load output arrays.
do i = 1,im
!c
!c If total precipitation is negative, restore initial profiles.
rn11(i) = rn1(i)
if (rn1(i).lt.0.) then
rn1(i) = 0.0
do k = 1,km
t1(k,i) = to(k,i)
q1(k,i) = qo0(k,i)
end do
!c
!c Otherwise, load modified t,q profile to output arrays
else
do k = 1,km
t1(k,i) = to3(k,i)
q1(k,i) = qo3(k,i)
end do
endif
!c
!c Load various diagostic output arrays
ktop(i) = ktcon(i)
kbot(i) = kbcon(i)
kuo(i) = 1
cldwrk(i) = aa1(i)
end do
! Transfer input u,v,cwm to ouptut arrays
do i = 1,im
do k=1,km
u1(k,i) = uo2(k,i)
v1(k,i) = vo2(k,i)
cwm1(k,i)= cwmo2(k,i)
end do
end do
if (.not.adjoint) return
!c
!c
!c
!c=====================================================================
!c
!c
!c Begin adjoint code here
! Adjoint of u,v,cwm transfer
do i = 1,im
do k=1,km
ad_uo2(k,i) = ad_uo2(k,i) + ad_u1(k,i)
ad_vo2(k,i) = ad_vo2(k,i) + ad_v1(k,i)
ad_cwmo2(k,i) = ad_cwmo2(k,i) + ad_cwm1(k,i)
ad_u1(k,i) = 0.0
ad_v1(k,i) = 0.0
ad_cwm1(k,i) = 0.0
end do
end do
!c
!c Adjoint of loading of output arrays
do i = 1,im
!c
!c zero or negative precip case.
if (rn11(i).lt.0.) then
ad_rn1(i) = 0.0
do k = km,1,-1
ad_qo0(k,i)= ad_qo0(k,i)+ ad_q1(k,i)
ad_to(k,i) = ad_to(k,i) + ad_t1(k,i)
ad_q1(k,i) = 0.0
ad_t1(k,i) = 0.0
end do
!c
!c precipitation greater than zero case
else
do k = km,1,-1
ad_qo3(k,i) = ad_qo3(k,i) + ad_q1(k,i)
ad_to3(k,i) = ad_to3(k,i) + ad_t1(k,i)
ad_q1(k,i) = 0.0
ad_t1(k,i) = 0.0
end do
endif
end do
!c
!c Adjoint of final precipitation calculation (including
!c evaporation of falling precipitation)
!c
do i = 1,im
!c
!c Zero local adjoint variables
ad_rn0 = 0.0
ad_delqev = 0.0
!c
!c Adjoint of rn0 to rn1 transfer
ad_rn0 = ad_rn0 + ad_rn1(i)
!c
do k = 1,km
!c
!c Restore saved values from nlm
rn0 = rn0k(k,i)
flg(i) = flgk(k,i)
delqev = delqevk(k,i)
qevap = 0.0
qcond = 0.0
delq2 = 0.0
!c
!c Zero local adjoint variables
ad_qevap = 0.0
ad_qcond = 0.0
!c
!c Adjoint of code for evaporation of falling precipitation
if (flg(i) .and. k.le.ktcon(i)) then
!c
!c Redo nlm
!! if (nint(slimsk(i)) .eq. 1) then
!! evef = 0.07
!! else
!! evef = edt(i)*evfact
!! endif
evef = edt(i)*evfact
term1 = evef
term2 = qo2(k,i) - qeso2(k,i)
term3 = el2orc*qeso2(k,i)/to2(k,i)**2
term4 = 1. / (1.+term3)
qcond = term1*term2*term4
dp = 100. * psfc(i) * del(k)
if (rn0.gt.0 .and. qcond.lt.0.) then
qevap = -qcond * (1.-exp(-.32*sqrt(dt2*rn0)))
qevap0 = qevap
if ( qevap .gt. (rn0*1000.*g/dp) ) then
qevap = rn0*1000.*g/dp
endif
delq2 = delqev + 0.001 * qevap * dp / g
qevap1 = qevap
if (delq2.gt.rntot(i)) then
qevap = 1000.* g * (rntot(i)-delqev) / dp
endif
endif
!c
!c Adjoint of t,q,rn0 adjustments due to evaporation
if ( rn0.gt.0. .and. qevap.gt.0. ) then
ad_qevap = ad_qevap + 0.001*ad_delqev*dp/g
ad_qevap = ad_qevap - 0.001*ad_rn0*dp/g
ad_qevap = ad_qevap - elocp*ad_to3(k,i)
ad_to2(k,i) = ad_to2(k,i) + ad_to3(k,i)
ad_qevap = ad_qevap + ad_qo3(k,i)
ad_qo2(k,i) = ad_qo2(k,i) + ad_qo3(k,i)
!c
!c Otherwise adjoint of no adjustment case
else
ad_qo2(k,i) = ad_qo2(k,i) + ad_qo3(k,i)
ad_to2(k,i) = ad_to2(k,i) + ad_to3(k,i)
endif
!c
!c Adjoint of qevap calculation
if (rn0.gt.0. .and. qcond.lt.0.) then
if (delq2.gt.rntot(i)) then
ad_delqev = ad_delqev - 1000.*g*ad_qevap/dp
ad_rntot(i) = ad_rntot(i) + 1000.*g*ad_qevap/dp
endif
if (qevap0 .gt. rn0*1000.*g/dp) then
ad_rn0 = ad_rn0 + ad_qevap*1000.*g/dp
endif
ad_rn0 = ad_rn0 &
+ qcond * (exp(-0.32*sqrt(dt2*rn0)) * &
(-0.32*0.5/sqrt(dt2*rn0)) * dt2*ad_qevap)
ad_qcond = ad_qcond &
- ad_qevap * (1.-exp(-0.32*sqrt(dt2*rn0)))
endif
!c
!c Adjoint of qcond calculation
ad_term4 = 0.0
ad_term3 = 0.0
ad_term2 = 0.0
ad_term1 = 0.0
ad_evef = 0.0
!c
ad_term4 = ad_term4 + term1*term2*ad_qcond
ad_term2 = ad_term2 + term1*term4*ad_qcond
ad_term1 = ad_term1 + term2*term4*ad_qcond
ad_term3 = ad_term3 - 1./(1.+term3)**2 * ad_term4
ad_to2(k,i) = ad_to2(k,i) &
- 2.*el2orc*qeso2(k,i)/to2(k,i)**3 * ad_term3
ad_qeso2(k,i) = ad_qeso2(k,i) &
+ el2orc*ad_term3/to2(k,i)**2
ad_qeso2(k,i) = ad_qeso2(k,i) - ad_term2
ad_qo2(k,i) = ad_qo2(k,i) + ad_term2
ad_evef = ad_evef + ad_term1
!! if (nint(slimsk(i)) .eq. 1) then
!! ad_evef = 0.0
!! else
!! ad_edt(i) = ad_edt(i) + ad_evef*evfact
!! endif
ad_edt(i) = ad_edt(i) + ad_evef*evfact
!c
!c
!c Adjoint of "2" to "3" transfer. This is the case
!c when above cloud top or no cloud present.
else
ad_to2(k,i) = ad_to2(k,i) + ad_to3(k,i)
ad_qo2(k,i) = ad_qo2(k,i) + ad_qo3(k,i)
endif
!c
!c Adjoint of layer rain calculation
if (k.le.ktcon(i) .and. ktcon(i).ne.1) then
!c
!c Redo nlm
aup = 1.
if (k.le.kb(i)) aup = 0.
adw = 1.
if (k.gt.jmin(i)) adw = 0.
rain = aup*pwo(k,i) + adw*edto1(i)*pwdo(k,i)
!c
!c Adjoint of tlm code.
ad_rain = 0.0
ad_xmb(i) = ad_xmb(i) + rain*ad_rn0*0.001*dt2
ad_rain = ad_rain + ad_rn0*xmb(i)*0.001*dt2
ad_pwdo(k,i) = ad_pwdo(k,i) + adw*edto1(i)*ad_rain
ad_edto1(i) = ad_edto1(i) + adw*ad_rain*pwdo(k,i)
ad_pwo(k,i) = ad_pwo(k,i) + aup*ad_rain
endif
!c
!c Zero local adjoint variables
ad_qevap = 0.0
ad_qcond = 0.0
!c
!c End of loop over vertical coordinate
end do
!c
!c Zero local adjoint variables
ad_rn0 = 0.0
ad_delqev = 0.0
!c
!c End of loop over profiles
end do
!c
!c
!c Adjoint of intial total precipitation calculation
do i = 1,im
do k = 1,kmax
if (k.le.ktcon(i) .and. ktcon(i).ne.1) then
!c
!c Redo nlm
aup = 1.
if (k.le.kb(i)) aup = 0.
adw = 1.
if (k.gt.jmin(i)) adw = 0.
rain = AUP * PWO(k,i) + ADW * EDTO1(I) * PWDO(k,i)
!c
!c Adjoint of tlm code
ad_rain = 0.0
ad_xmb(i) = ad_xmb(i) + rain*0.001*dt2 * ad_rntot(i)
ad_rain = ad_rain + 0.001*xmb(i)*dt2 * ad_rntot(i)
ad_pwdo(k,i) = ad_pwdo(k,i) + adw*edto1(i)*ad_rain
ad_edto1(i) = ad_edto1(i) + adw*pwdo(k,i)*ad_rain
ad_pwo(k,i) = ad_pwo(k,i) + aup*ad_rain
endif
end do
end do
!c
!c Adjoint of to2 and qo2 calculation (prior to evaportation)
do i = 1,im
do k = km,1,-1
ad_dellal1(i) = 0.0
ad_dellal1(i) = ad_dellal1(i) + ad_dellal(i)
if (k.le.ktcon(i) .and. ktcon(i).ne.1) then
!
! Redo nlm
dellat = (dellah(k,i)-hvap*dellaq(k,i))/cp
call adfpvsx(to2(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
!
! Adjoint of tlm code
ad_dellat = 0.0
ad_es = 0.0
ad_es0 = 0.0
if (k.eq.ktcon(i) .and. ncloud.gt.0 .and. cnvflg(i)) then
dp = 100.*psfc(i)*del(k)
ad_dellalk(k,i) = ad_dellalk(k,i) + &
ad_dellal1(i)*xmb(i)*dp/g
ad_xmb(i) = ad_xmb(i) + &
dellalk(k,i)*ad_dellal1(i)*dp/g
ad_dellalk(k,i) = ad_dellalk(k,i) + &
ad_cwmo2(k,i)*xmb(i)*dt2
ad_xmb(i) = ad_xmb(i) + &
dellalk(k,i)*ad_cwmo2(k,i)*dt2
ad_cwmo(k,i) = ad_cwmo(k,i) + ad_cwmo2(k,i)
ad_cwmo2(k,i) = 0.0
else
ad_cwmo(k,i) = ad_cwmo(k,i) + ad_cwmo2(k,i)
ad_dellalk(k,i) = ad_dellalk(k,i) + ad_dellal1(i)
ad_cwmo2(k,i) = 0.0
endif
!c
ad_es = ad_es - eps*es/( epsm1*ad_qeso2(k,i)
ad_es = ad_es + eps*ad_qeso2(k,i)/(
ad_es0 = ad_es0 + 10*ad_es
adt=0.0
call adfpvsx(to2(k,i),es0,adt,ad_es0,adjoint)
ad_to2(k,i) = ad_to2(k,i) + adt
ad_vo(k,i) = ad_vo(k,i) + ad_vo2(k,i)
ad_xmb(i) = ad_xmb(i) + &
ad_vo2(k,i)*dellav(k,i)*dt2
ad_dellav(k,i) = ad_dellav(k,i) + &
ad_vo2(k,i)*xmb(i)*dt2
ad_uo(k,i) = ad_uo(k,i) + ad_uo2(k,i)
ad_xmb(i) = ad_xmb(i) + &
ad_uo2(k,i)*dellau(k,i)*dt2
ad_dellau(k,i) = ad_dellau(k,i) + &
ad_uo2(k,i)*xmb(i)*dt2
ad_qo(k,i) = ad_qo(k,i) + ad_qo2(k,i)
ad_xmb(i) = ad_xmb(i) + &
ad_qo2(k,i)*dellaq(k,i)*dt2
ad_dellaq(k,i) = ad_dellaq(k,i) + &
ad_qo2(k,i)*xmb(i)*dt2
ad_to(k,i) = ad_to(k,i) + ad_to2(k,i)
ad_xmb(i) = ad_xmb(i) + ad_to2(k,i)*dellat*dt2
ad_dellat = ad_dellat + ad_to2(k,i)*xmb(i)*dt2
ad_dellaq(k,i) = ad_dellaq(k,i) - hvap*ad_dellat/cp
ad_dellah(k,i) = ad_dellah(k,i) + ad_dellat/cp
else
!c
!c Redo nlm
call adfpvsx(to2(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
!c
!c Adjoint of tlm code
ad_es = 0.0
ad_es0= 0.0
ad_dellalk(k,i) = ad_dellalk(k,i) + ad_dellal1(i)
ad_es = ad_es - eps*es/( epsm1*ad_qeso2(k,i)
ad_es = ad_es + eps*ad_qeso2(k,i)/(
ad_es0 = ad_es0 + 10.*ad_es
adt = 0.0
call adfpvsx(to2(k,i),es0,adt,ad_es0,adjoint)
ad_to2(k,i) = ad_to2(k,i) + adt
ad_cwmo(k,i) = ad_cwmo(k,i) + ad_cwmo2(k,i)
ad_vo(k,i) = ad_vo(k,i) + ad_vo2(k,i)
ad_uo(k,i) = ad_uo(k,i) + ad_uo2(k,i)
ad_qo(k,i) = ad_qo(k,i) + ad_qo2(k,i)
ad_to(k,i) = ad_to(k,i) + ad_to2(k,i)
endif
ad_dellal(i) = ad_dellal(i) + ad_dellalk(k,i)
ad_dellalk(k,i) = 0.0
ad_dellal1(i) = 0.0
!c
end do
end do
!c
!c Adjoint of kernel (cloud mass flux) calculation
do i = 1,im
!c
!c Redo nlm
if (abs(xk(i)).gt.tiny) then
ratio = -f(i)/xk(i)
else
ratio = 0.0
endif
fixed = tdpscl(i)*tdpcld(i)*tdetrn(i)*tdpdwn(i)
!c
!c Adjoint of tlm
if (xmb0(i).gt.xmbmax(i)) then
!c do nothing since tlm = 0.0
ad_xmb(i) = 0.0
endif
ad_ratio = 0.0
if ( ad_tcwfdn(i) = ad_tcwfdn(i) &
+ fixed*ratio*tcwfup(i)*ad_xmb(i)
ad_tcwfup(i) = ad_tcwfup(i) &
+ fixed*ratio*ad_xmb(i)*tcwfdn(i)
ad_ratio = ad_ratio + fixed*ad_xmb(i)*tcwfup(i)*tcwfdn(i)
else
!c do nothing since tlm = 0.0
ad_xmb(i) = 0.0
endif
if (abs(xk(i)).gt.tiny) then
ad_xk(i) = ad_xk(i) + f(i)/xk(i)**2 * ad_ratio
ad_f(i) = ad_f(i) - ad_ratio/xk(i)
else
!c do nothing since tlm = 0.0
ad_ratio = 0.0
endif
!c
!c Adjoint of xk tlm
ad_aa1(i) = ad_aa1(i) - ad_xk(i)/mbdt
ad_xaa0(i) = ad_xaa0(i) + ad_xk(i)/mbdt
!c
!c Adjoint of f tlm
ad_dtconv(i) = ad_dtconv(i) &
- (aa1(i) - acrt(i)*acrtfct1(i))/dtconv(i)**2 &
* ad_f(i)
ad_acrtfct1(i) = ad_acrtfct1(i) &
- acrt(i)*ad_f(i) / dtconv(i)
ad_aa1(i) = ad_aa1(i) + ad_f(i) / dtconv(i)
!c
end do
!c
!c Adjoint of acrtfct1 and dtconv calculations
!c
do i = 1,im
!c
!c Redo nlm
if (nint(slimsk(i)) .eq. 1) then
w1 = w1l
w2 = w2l
w3 = w3l
w4 = w4l
else
w1 = w1s
w2 = w2s
w3 = w3s
w4 = w4s
endif
!c
!c Adjoint of dtconv calculation
if (dtconv0(i).gt.dtmax) then
!c do nothing since tlm = 0.0
ad_dtconv(i) = 0.0
endif
if (dtconv0(i).lt.dtmin) then
!c do nothing since tlm = 0.0
ad_dtconv(i) = 0.0
endif
ad_pdot(i) = ad_pdot(i) + (1800.-dt2) &
* ad_dtconv(i)/(w1-w2)
!c
!c Adjoint of acrtfct calculation
ad_acrtfct(i) = ad_acrtfct(i) - ad_acrtfct1(i)
if (acrtfct0(i).gt.1.) then
!c do nothing since tlm = 0.0
ad_acrtfct(i) = 0.0
endif
if (acrtfct0(i).lt.-1.) then
!c do nothing since tlm = 0.0
ad_acrtfct(i) = 0.0
endif
if (pdot(i).le.w4) then
ad_pdot(i) = ad_pdot(i) + ad_acrtfct(i)/(w3-w4)
elseif (pdot(i).ge.-w4) then
ad_pdot(i) = ad_pdot(i) - ad_acrtfct(i)/(w4-w3)
else
!c do nothing since tlm = 0.0
ad_acrtfct(i) = 0.0
endif
!c
end do
!c
!c Adjoint of xaa0 downdraft cloudwork function calculation
do i = 1,im
do k = 1,kmax-1
if (dwnflg2(i) .and. k.lt.jmin(i)) then
!c
!c Redo nlm calculations
GAMMA = EL2ORC * XQESOL(k+1,i) / XTOL(k+1,i)**2
DHH = XHCD(I)
DT = XTOL(k+1,i)
DG = GAMMA
DH = XHESO(k+1,i)
dz = xzo(k,i) - xzo(k+1,i)
!c
term1 = edtx1(i)
term2 = dz
term3 = g/(cp*dt)
term4 = (dhh-dh)
term5 = 1./(1.+dg)
term6 = 1. + delta*cp*dg*dt/hvap
dqs = xqesol(k+1,i) - xqol(k+1,i)
!c
!c Adjoint of tlm code
ad_dqs = 0.0
ad_dz = 0.0
ad_dt = 0.0
ad_dg = 0.0
ad_dh = 0.0
ad_dhh = 0.0
if (dqs.gt.0.) then
ad_dqs = ad_dqs + edtx1(i)*dz*g*delta * ad_xaa0(i)
ad_dz = ad_dz + edtx1(i)*g*delta*dqs * ad_xaa0(i)
ad_edtx1(i) = ad_edtx1(i) + dz*g*delta*dqs &
* ad_xaa0(i)
endif
!c
ad_term1 = term2*term3*term4*term5*term6*ad_xaa0(i)
ad_term2 = term1*term3*term4*term5*term6*ad_xaa0(i)
ad_term3 = term1*term2*term4*term5*term6*ad_xaa0(i)
ad_term4 = term1*term2*term3*term5*term6*ad_xaa0(i)
ad_term5 = term1*term2*term3*term4*term6*ad_xaa0(i)
ad_term6 = term1*term2*term3*term4*term5*ad_xaa0(i)
ad_xqol(k+1,i) = ad_xqol(k+1,i) - ad_dqs
ad_xqesol(k+1,i) = ad_xqesol(k+1,i) + ad_dqs
!c
ad_dt = (delta*cp/hvap) * dg*ad_term6
ad_dg = (delta*cp/hvap) * ad_term6*dt
ad_dg = ad_dg - 1./(1.+dg)**2 * ad_term5
ad_dh = ad_dh - ad_term4
ad_dhh = ad_dhh + ad_term4
ad_dt = ad_dt - g/(cp*dt*dt) * ad_term3
ad_dz = ad_dz + ad_term2
ad_edtx1(i) = ad_edtx1(i) + ad_term1
ad_xzo(k,i) = ad_xzo(k,i) + ad_dz
ad_xzo(k+1,i) = ad_xzo(k+1,i) - ad_dz
ad_xheso(k+1,i)= ad_xheso(k+1,i) + ad_dh
ad_gamma = ad_dg
ad_xtol(k+1,i) = ad_xtol(k+1,i) + ad_dt
ad_xhcd(i) = ad_xhcd(i) + ad_dhh
ad_xtol(k+1,i) = ad_xtol(k+1,i) &
- 2.*el2orc*xqesol(k+1,i)/xtol(k+1,i)**3 &
* ad_gamma
ad_xqesol(k+1,i) = ad_xqesol(k+1,i) + &
el2orc*ad_gamma/xtol(k+1,i)**2
endif
end do
end do
!c
!c Adjoint of final edtx1 calculation
do i = 1,im
edtmax = edtmaxl
if (nint(slimsk(i)) .eq. 0) edtmax = edtmaxs
if (dwnflg2(i)) then
if (xpwev(i).ge.0.) then
!c do nothing since tlm is 0.0
ad_edtx1(i) = 0.0
else
if (edtx10(i).gt.edtmax) then
!c do nothing since tlm is 0.0
ad_edtx1(i) = 0.0
else
ad_xpwev(i) = ad_xpwev(i) &
+ edtx(i)*xpwav(i)/xpwev(i)**2 * ad_edtx1(i)
ad_xpwav(i) = ad_xpwav(i) &
- edtx(i)*ad_edtx1(i)/xpwev(i)
ad_edtx(i) = ad_edtx(i) &
- ad_edtx1(i)*xpwav(i)/xpwev(i)
endif
endif
else
!c do nothing since tlm is 0.0
ad_edtx1(i) = 0.0
endif
enddo
!c
!c Adjoint of downdraft xpwev calculation
do i = 1,im
do k = 1,kmax-1
if (dwnflg2(i) .and. k.lt.jmin(i)) then
!c
!c Redo nlm
dq = xqesol(k,i)
dt = xtol(k,i)
gamma = el2orc*dq/dt**2
dh = xhcd(i) - xheso(k,i)
term5 = (1./hvap) * gamma
term6 = 1./(1.+gamma)
detad = etad(k+1,i) - etad(k,i)
term1 = etad(k+1,i)*xqrcd(k+1,i)
term2 = etad(k,i)*xqrcd(k,i)
term3 = detad * 0.5 * (xqrcd(k,i) + xqrcd(k+1,i))
!c
!c Adjoint of tlm code
ad_xpwd = 0.0
ad_detad = 0.0
ad_dh = 0.0
ad_gamma = 0.0
ad_dt = 0.0
ad_dq = 0.0
ad_term1 = 0.0
ad_term2 = 0.0
ad_term3 = 0.0
ad_term5 = 0.0
ad_term6 = 0.0
ad_xpwd = ad_xpwd + ad_xpwev(i)
ad_term1 = ad_term1 + ad_xpwd
ad_term2 = ad_term2 - ad_xpwd
ad_term3 = ad_term3 - ad_xpwd
ad_xqrcd(k+1,i) = ad_xqrcd(k+1,i) + etad(k+1,i)*ad_term1
ad_etad(k+1,i) = ad_etad(k+1,i) + ad_term1*xqrcd(k+1,i)
ad_xqrcd(k,i) = ad_xqrcd(k,i) + etad(k,i)*ad_term2
ad_etad(k,i) = ad_etad(k,i) + ad_term2*xqrcd(k,i)
ad_detad = ad_detad &
+ ad_term3*0.5*(xqrcd(k,i)+xqrcd(k+1,i))
ad_xqrcd(k,i) = ad_xqrcd(k,i) + detad*0.5*ad_term3
ad_xqrcd(k+1,i) = ad_xqrcd(k+1,i) + detad*0.5*ad_term3
ad_etad(k,i) = ad_etad(k,i) - ad_detad
ad_etad(k+1,i) = ad_etad(k+1,i) + ad_detad
ad_dh = ad_dh + term5*term6*ad_xqrcd(k,i)
ad_term6 = ad_term6 + term5*dh*ad_xqrcd(k,i)
ad_term5 = ad_term5 + term6*dh*ad_xqrcd(k,i)
ad_dq = ad_dq + ad_xqrcd(k,i)
ad_gamma = ad_gamma -1./(1.+gamma)**2 * ad_term6
ad_gamma = ad_gamma + (1./hvap)*ad_term5
ad_xheso(k,i) = ad_xheso(k,i) - ad_dh
ad_xhcd(i) = ad_xhcd(i) + ad_dh
ad_dt = ad_dt - 2.*el2orc*dq/dt**3 * ad_gamma
ad_dq = ad_dq + el2orc*ad_gamma/dt**2
ad_xtol(k,i) = ad_xtol(k,i) + ad_dt
ad_xqesol(k,i) = ad_xqesol(k,i) + ad_dq
endif
end do
end do
!c
!c Adjoint of xheo,xqesol --> xhcd,xqrcd transfers
do i = 1,im
if (dwnflg2(i)) then
jmn = jmin(i)
ad_xqesol(jmn,i) = ad_xqesol(jmn,i) + ad_xqrcd(jmn,i)
ad_xheo(jmn,i) = ad_xheo(jmn,i) + ad_xhcd(i)
endif
end do
!c
!c Adjoint of xaa0 calculations
!c
do i = 1,im
do k = kmax-1,2,-1
if (k.gt.kbcon(i).and.k.le.ktcon(i)) then
!c
!c Recalculate nlm
dz1 = xzo(k,i) - xzo(k-1,i)
gamma = el2orc*xqesol(k-1,i)/xtol(k-1,i)**2
rfact = 1. + delta*cp*gamma*xtol(k-1,i)/hvap
xdby = xhcko(k-1,i) - xheso(k-1,i)
term1 = dz1 * (g/cp)
term2 = 1./xtol(k-1,i)
term3 = xdby
term4 = 1./(1.+gamma)
term5 = rfact
dqs = xqesol(k-1,i) - xqol(k-1,i)
!c
!c Zero local adjoint variables
ad_dqs = 0.0
ad_dz1 = 0.0
ad_term1 = 0.0
ad_term2 = 0.0
ad_term3 = 0.0
ad_term4 = 0.0
ad_term5 = 0.0
ad_rfact = 0.0
ad_gamma = 0.0
ad_xdby = 0.0
!c
!c Adjoint of forward model tlm
if (dqs.gt.0.) then
ad_dqs = ad_dqs + g*delta * dz1*ad_xaa0(i)
ad_dz1 = ad_dz1 + g*delta * dqs*ad_xaa0(i)
endif
ad_xqol(k-1,i) = ad_xqol(k-1,i) - ad_dqs
ad_xqesol(k-1,i) = ad_xqesol(k-1,i) + ad_dqs
ad_term5 = ad_term5 + term1*term2*term3*term4*ad_xaa0(i)
ad_term4 = ad_term4 + term1*term2*term3*ad_xaa0(i)*term5
ad_term3 = ad_term3 + term1*term2*ad_xaa0(i)*term4*term5
ad_term2 = ad_term2 + term1*ad_xaa0(i)*term3*term4*term5
ad_term1 = ad_term1 + ad_xaa0(i)*term2*term3*term4*term5
ad_rfact = ad_rfact + ad_term5
ad_gamma = ad_gamma - 1./(1.+gamma)**2 * ad_term4
ad_xdby = ad_xdby + ad_term3
ad_xtol(k-1,i) = ad_xtol(k-1,i) &
- 1./(xtol(k-1,i)**2) * ad_term2
ad_dz1 = ad_dz1 + ad_term1 * (g/cp)
ad_xheso(k-1,i) = ad_xheso(k-1,i) - ad_xdby
ad_xhcko(k-1,i) = ad_xhcko(k-1,i) + ad_xdby
ad_xtol(k-1,i) = ad_xtol(k-1,i) &
+ (delta*cp/hvap)*gamma*ad_rfact
ad_gamma = ad_gamma &
+ (delta*cp/hvap)*ad_rfact*xtol(k-1,i)
ad_xtol(k-1,i) = ad_xtol(k-1,i) &
- 2.*el2orc*xqesol(k-1,i)/(xtol(k-1,i)**3) &
* ad_gamma
ad_xqesol(k-1,i) = ad_xqesol(k-1,i) &
+ el2orc*ad_gamma/xtol(k-1,i)**2
ad_xzo(k-1,i) = ad_xzo(k-1,i) - ad_dz1
ad_xzo(k,i) = ad_xzo(k,i) + ad_dz1
endif
end do
end do
!c
!c Adjoint of xpwav calculation and initial xaa0 calculation
!c
do i = 1,im
do k = kmax-1,2,-1
if (k.gt.kb(i) .and. k.lt.ktcon(i)) then
!c
!c Redo nlm
factor = eta(k-1,i)/eta(k,i)
onemf = 1. - factor
termq = 0.5*(xqol(k,i) + xqol(k+1,i))
gamma = el2orc*xqesol(k,i)/xtol(k,i)**2
xdby = xhcko(k,i) - xheso(k,i)
xtemp = xtemp0(k,i)
xqrch = xqrch0(k,i)
dq = eta(k,i)*(xtemp-xqrch)
!c
!c Adjoint code
ad_xtemp = 0.0
ad_termq = 0.0
ad_onemf = 0.0
ad_factor = 0.0
ad_xqrch = 0.0
ad_gamma = 0.0
ad_xdby = 0.0
ad_xqc = 0.0
ad_qlk = 0.0
ad_xpw = 0.0
ad_dz = 0.0
ad_etah = 0.0
ad_dz1 = 0.0
ad_dq = 0.0
!c
!c Adjoint of actions dependent on sign of dq
if (dq.gt.0.) then
!c
!c Redo nlm calculations
dz = 0.5 * (xzo(k+1,i) - xzo(k-1,i))
dz1 = xzo(k,i) - xzo(k-1,i)
etah = 0.5 * (eta(k,i) + eta(k-1,i))
qlk = dq / (eta(k,i) + etah * c0 * dz)
!c
!c Adjoint code
ad_xpw = ad_xpw + ad_xpwav(i)
ad_qlk = ad_qlk + etah*c0*dz*ad_xpw
ad_dz = ad_dz + etah*c0*ad_xpw*qlk
ad_etah = ad_etah + ad_xpw*c0*dz*qlk
ad_qlk = ad_qlk - dz1*g*ad_xaa0(i)
ad_dz1 = ad_dz1 - ad_xaa0(i)*g*qlk
ad_xqc = ad_xqc + ad_xqcko(k,i)
ad_xqrch = ad_xqrch + ad_xqc
ad_qlk = ad_qlk + ad_xqc
ad_dz = ad_dz &
- (dq/(eta(k,i)+etah*c0*dz)**2) &
* etah*c0*ad_qlk
ad_etah = ad_etah &
- (dq/(eta(k,i)+etah*c0*dz)**2) &
* ad_qlk*c0*dz
ad_eta(k,i) = ad_eta(k,i) &
- (dq/(eta(k,i)+etah*c0*dz)**2) &
* ad_qlk
ad_dq = ad_dq + ad_qlk/(eta(k,i)+etah*c0*dz)
ad_eta(k-1,i) = ad_eta(k-1,i) + 0.5*ad_etah
ad_eta(k,i) = ad_eta(k,i) + 0.5*ad_etah
ad_xzo(k-1,i) = ad_xzo(k-1,i) - ad_dz1
ad_xzo(k,i) = ad_xzo(k,i) + ad_dz1
ad_xzo(k-1,i) = ad_xzo(k-1,i) - 0.5*ad_dz
ad_xzo(k+1,i) = ad_xzo(k+1,i) + 0.5*ad_dz
else
ad_xtemp = ad_xtemp + ad_xqcko(k,i)
endif
!c
!c Adjoint of dq calculation
ad_xqrch = ad_xqrch - eta(k,i)*ad_dq
ad_xtemp = ad_xtemp + eta(k,i)*ad_dq
ad_eta(k,i) = ad_eta(k,i) + ad_dq*(xtemp-xqrch)
!c
!c Adjoint of xqrch calculation
xdby = xhcko(k,i) - xheso(k,i)
if (xdby.gt.0.) then
ad_xqesol(k,i) = ad_xqesol(k,i) + ad_xqrch
ad_gamma = ad_gamma &
- gamma*xdby/((hvap*(1.+gamma))**2) &
* hvap*ad_xqrch
ad_xdby = ad_xdby + gamma*ad_xqrch/(hvap*(1.+gamma))
ad_gamma = ad_gamma + ad_xqrch*xdby/(hvap*(1.+gamma))
ad_xtol(k,i) = ad_xtol(k,i) &
- 2.*el2orc*xqesol(k,i)/xtol(k,i)**3 &
* ad_gamma
ad_xqesol(k,i) = ad_xqesol(k,i) &
+ el2orc*ad_gamma/xtol(k,i)**2
else
ad_xqesol(k,i) = ad_xqesol(k,i) + ad_xqrch
endif
ad_xheso(k,i) = ad_xheso(k,i) - ad_xdby
ad_xhcko(k,i) = ad_xhcko(k,i) + ad_xdby
!c
!c Adjoint of xtemp calculation
ad_termq = ad_termq + onemf*ad_xtemp
ad_onemf = ad_onemf + ad_xtemp*termq
ad_xqcko(k-1,i) = ad_xqcko(k-1,i) + factor*ad_xtemp
ad_factor = ad_factor + ad_xtemp*xqcko(k-1,i)
ad_xqol(k+1,i) = ad_xqol(k+1,i) + 0.5*ad_termq
ad_xqol(k,i) = ad_xqol(k,i) + 0.5*ad_termq
ad_factor = ad_factor - ad_onemf
ad_eta(k,i) = ad_eta(k,i) &
- eta(k-1,i)/eta(k,i)**2 * ad_factor
ad_eta(k-1,i) = ad_eta(k-1,i) + ad_factor/eta(k,i)
!c
endif
end do
end do
!c
!c Adjoint of xhcko calculation
do i = 1,im
do k = kmax-1,2,-1
if (k.gt.kb(i).and.k.le.ktcon(i)) then
!c
!c Recalculate nlm
factor = eta(k-1,i)/eta(k,i)
onemf = 1. - factor
!c
!c Adjoint of tlm code
ad_onemf = 0.0
ad_factor = 0.0
ad_xheo(k+1,i) = ad_xheo(k+1,i) + onemf*0.5*ad_xhcko(k,i)
ad_xheo(k,i) = ad_xheo(k,i) + onemf*0.5*ad_xhcko(k,i)
ad_onemf = ad_onemf + &
ad_xhcko(k,i) * 0.5*(xheo(k,i)+xheo(k+1,i))
ad_xhcko(k-1,i) = ad_xhcko(k-1,i) + factor*ad_xhcko(k,i)
ad_factor = ad_factor + ad_xhcko(k,i)*xhcko(k-1,i)
ad_factor = ad_factor - ad_onemf
ad_eta(k,i) = ad_eta(k,i) - &
(eta(k-1,i)/eta(k,i)**2) * ad_factor
ad_eta(k-1,i) = ad_eta(k-1,i) + ad_factor/eta(k,i)
endif
end do
end do
!c
!c Adjoint of xheo,xqol --> xhcko,xqcko transfer
do i = 1,im
indx = kb(i)
ad_xheo(indx,i) = ad_xheo(indx,i) + ad_xhcko(indx,i)
ad_xqol(indx,i) = ad_xqol(indx,i) + ad_xqcko(indx,i)
end do
!c
!c Adjoint of xheo and xheso calculation at k=kmax
k = kmax
do i = 1,im
ad_xqesol(k,i) = ad_xqesol(k,i) + hvap*ad_xheso(k,i)
ad_xtol(k,i) = ad_xtol(k,i) + cp *ad_xheso(k,i)
ad_xzo(k,i) = ad_xzo(k,i) + g *ad_xheso(k,i)
ad_xqol(k,i) = ad_xqol(k,i) + hvap*ad_xheo(k,i)
ad_xtol(k,i) = ad_xtol(k,i) + cp *ad_xheo(k,i)
ad_xzo(k,i) = ad_xzo(k,i) + g *ad_xheo(k,i)
enddo
!c
!c Adjoint of xheo and xheso calculation below k=kmax
do i = 1,im
do k = kmax-1,1,-1
!c
!c Redo nlm calculation
po = 0.5 * ( call adfpvsx(xtol(k,i),es0,adt,ad_es0,adjoint0)
esl = 10.*es0
!c
!c Adjoint of tlm code
ad_esl = 0.0
ad_es0 = 0.0
ad_xqesol(k,i) = ad_xqesol(k,i) + hvap *ad_xheso(k,i)
ad_xtol(k,i) = ad_xtol(k,i) + cp *ad_xheso(k,i)
ad_xzo(k+1,i) = ad_xzo(k+1,i) + 0.5*g*ad_xheso(k,i)
ad_xzo(k,i) = ad_xzo(k,i) + 0.5*g*ad_xheso(k,i)
ad_xqol(k,i) = ad_xqol(k,i) + hvap *ad_xheo(k,i)
ad_xtol(k,i) = ad_xtol(k,i) + cp *ad_xheo(k,i)
ad_xzo(k+1,i) = ad_xzo(k+1,i) + 0.5*g*ad_xheo(k,i)
ad_xzo(k,i) = ad_xzo(k,i) + 0.5*g*ad_xheo(k,i)
ad_esl = ad_esl &
- eps*esl/((po+epsm1*esl)**2) &
*epsm1*ad_xqesol(k,i)
ad_esl = ad_esl &
+ eps*ad_xqesol(k,i)/(po+epsm1*esl)
ad_es0 = ad_es0 + 10*ad_esl
adt = 0.0
call adfpvsx(xtol(k,i),es0,adt,ad_es0,adjoint)
ad_xtol(k,i) = ad_xtol(k,i) + adt
end do
end do
!c
!c Adjoint of xtol and xqol calculations
do i = 1,im
do k = kmax-1,1,-1
!c
!c Redo nlm calculations
dz = 0.5 * (xzo(k+1,i) - xzo(k,i))
dp = 0.5 * ( call adfpvsx(xto(k+1,i),es0,adt,ad_es0,adjoint0)
es = 10*es0
pprime = p(k+1,i) + epsm1 * es
qs = eps * es / pprime
dqsdp = - qs / pprime
desdt = ES * (FACT1 / xto(k+1,i) + FACT2 / (xto(k+1,i)**2))
dqsdt = qs * p(k+1,i) * desdt / (es * pprime)
gamma = el2orc * xqeso(k+1,i) / (xto(k+1,i)**2)
dt = (g * dz + hvap*dqsdp*dp) / (cp*(1.+gamma))
dq = dqsdt * dt + dqsdp * dp
!c
!c Initialize local ajm variables to 0.0
ad_dq = 0.0
ad_dt = 0.0
ad_dqsdp = 0.0
ad_dqsdt = 0.0
ad_gamma = 0.0
ad_dz = 0.0
ad_pprime = 0.0
ad_es = 0.0
ad_desdt = 0.0
ad_qs = 0.0
ad_es0 = 0.0
!c
!c Adjoint of tlm code
if (xqol0(k,i)>0.) then
ad_xqol0(k,i) = ad_xqol0(k,i) + ad_xqol(k,i)
else
ad_xqol0(k,i) = 0.0
endif
ad_xqo(k+1,i) = ad_xqo(k+1,i) + ad_xqol0(k,i)
ad_dq = ad_dq + ad_xqol0(k,i)
ad_xto(k+1,i) = ad_xto(k+1,i) + ad_xtol(k,i)
ad_dt = ad_dt + ad_xtol(k,i)
ad_dqsdp = ad_dqsdp + ad_dq*dp
ad_dt = ad_dt + dqsdt*ad_dq
ad_dqsdt = ad_dqsdt + ad_dq*dt
ad_gamma = ad_gamma - &
(g*dz+hvap*dqsdp*dp)/((cp*(1.+gamma))**2) * &
cp*ad_dt
ad_dqsdp = ad_dqsdp + hvap*ad_dt*dp/(cp*(1.+gamma))
ad_dz = ad_dz + g*ad_dt/(cp*(1.+gamma))
ad_xto(k+1,i) = ad_xto(k+1,i) - &
2.*el2orc*xqeso(k+1,i)/(xto(k+1,i)**3) * ad_gamma
ad_xqeso(k+1,i) = ad_xqeso(k,i) + &
el2orc*ad_gamma/(xto(k+1,i)**2)
rterm = 1./(es*pprime)
ad_pprime = ad_pprime - (qs*p(k+1,i)*desdt)*rterm**2 * &
es*ad_dqsdt
ad_es = ad_es - (qs*p(k+1,i)*desdt)*rterm**2 * &
ad_dqsdt*pprime
ad_desdt = ad_desdt + p(k+1,i)*qs*ad_dqsdt*rterm
ad_qs = ad_qs + p(k+1,i)*ad_dqsdt*desdt*rterm
ad_es = ad_es + &
ad_desdt*(fact1/xto(k+1,i) + fact2/(xto(k+1,i)**2))
ad_xto(k+1,i) = ad_xto(k+1,i) - ad_desdt * &
(es/(xto(k+1,i)**2)) * (fact1+2.*fact2/xto(k+1,i))
ad_pprime = ad_pprime + (qs/(pprime**2))*ad_dqsdp
ad_qs = ad_qs - ad_dqsdp/pprime
ad_pprime = ad_pprime - (eps*es/(pprime**2))*ad_qs
ad_es = ad_es + eps*ad_qs/pprime
ad_es = ad_es + epsm1*ad_pprime
ad_es0 = ad_es0 + 10*ad_es
adt = 0.0
call adfpvsx(xto(k+1,i),es0,adt,ad_es0,adjoint)
ad_xto(k+1,i) = ad_xto(k+1,i) + adt
ad_xzo(k,i) = ad_xzo(k,i) - 0.5*ad_dz
ad_xzo(k+1,i) = ad_xzo(k+1,i) + 0.5*ad_dz
end do
end do
!c Adjoint of integration of hydrostatic equation to get xzo
do i = 1,im
do k = kmax,2,-1
!c
!c Redo nlm calculations
DLNSIG = LOG(SL(K) / SL(K-1))
term1 = dlnsig * rd / g
!c
!c Adjoint of tlm code
ad_term2 = 0.0
ad_term2 = ad_term2 - term1*ad_xzo(k,i)
ad_xzo(k-1,i) = ad_xzo(k-1,i) + ad_xzo(k,i)
ad_xtvo(k-1,i) = ad_xtvo(k-1,i) + 0.5*ad_term2
ad_xtvo(k,i) = ad_xtvo(k,i) + 0.5*ad_term2
end do
end do
dlnsig = log(sl(1))
do i = 1,im
ad_xtvo(1,i) = ad_xtvo(1,i) - dlnsig*rd/g * ad_xzo(1,i)
end do
!c
!c Adjoint of xqeso and xtvo calculations
do i = 1,im
do k = kmax,1,-1
call adfpvsx(xto(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
ad_es = 0.0
ad_es0 = 0.0
ad_xqo(k,i) = ad_xqo(k,i) + delta*xto(k,i)*ad_xtvo(k,i)
ad_xto(k,i) = ad_xto(k,i) + delta*xqo(k,i)*ad_xtvo(k,i)
ad_xto(k,i) = ad_xto(k,i) + ad_xtvo(k,i)
ad_es = ad_es &
- eps*es/((p(k,i)+epsm1*es)**2) * &
epsm1*ad_xqeso(k,i)
ad_es = ad_es &
+ eps*ad_xqeso(k,i)/(
ad_es0 = ad_es0 + 10*ad_es
adt = 0.0
call adfpvsx(xto(k,i),es0,adt,ad_es0,adjoint)
ad_xto(k,i) = ad_xto(k,i) + adt
end do
end do
!c
!c Adjoint of xto,xqo calculation
do i = 1,im
do k = kmax,1,-1
ad_dellat = 0.0
if ( k.le.ktcon(i) .and. ktcon(i).ne.1 ) then
ad_to(k,i) = ad_to(k,i) + ad_xto(k,i)
ad_dellat = ad_dellat + ad_xto(k,i)*mbdt
ad_dellaq(k,i) = ad_dellaq(k,i) - hvap*ad_dellat/cp
ad_dellah(k,i) = ad_dellah(k,i) + ad_dellat/cp
if (xqo0(k,i)>0.) then
ad_xqo0(k,i) = ad_xqo0(k,i) + ad_xqo(k,i)
else
ad_xqo0(k,i) = 0.0
endif
ad_qo0(k,i) = ad_qo0(k,i) + ad_xqo0(k,i)
ad_dellaq(k,i) = ad_dellaq(k,i) + ad_xqo0(k,i)*mbdt
else
ad_to(k,i) = ad_to(k,i) + ad_xto(k,i)
ad_qo0(k,i)= ad_qo0(k,i)+ ad_xqo(k,i)
endif
end do
end do
!c
!c Adjoint of dellah and dellaq calculations at cloud top
do i = 1,im
if (ktcon(i).ne.1) then
indx = ktcon(i)
!c
!c Redo nlm calculations
dp = 100. * psfc(i)*del(indx)
dv1 = heol(indx-1,i)
dvq1 = qol(indx-1,i)
dvu1 = uol(indx-1,i)
dvv1 = vol(indx-1,i)
!c
!c Adjoint of tlm code
ad_dvv1 = 0.0
ad_dvu1 = 0.0
ad_dvq1 = 0.0
ad_dv1 = 0.0
ad_qlko_ktcon(i) = ad_qlko_ktcon(i) + &
eta(indx-1,i)*ad_dellal(i)*g/dp
ad_eta(indx-1,i) = ad_eta(indx-1,i) + &
ad_dellal(i)*qlko_ktcon(i)*g/dp
ad_dvv1 = ad_dvv1 &
- eta(indx-1,i) * ad_dellav(indx,i) * (g/dp)
ad_vcko(indx-1,i) = ad_vcko(indx-1,i) &
+ eta(indx-1,i) * ad_dellav(indx,i) * (g/dp)
ad_eta(indx-1,i) = ad_eta(indx-1,i) &
+ ad_dellav(indx,i)*(vcko(indx-1,i)-dvv1)*(g/dp)
ad_vol(indx-1,i) = ad_vol(indx-1,i) + ad_dvv1
ad_dvu1 = ad_dvu1 &
- eta(indx-1,i) * ad_dellau(indx,i) * (g/dp)
ad_ucko(indx-1,i) = ad_ucko(indx-1,i) &
+ eta(indx-1,i) * ad_dellau(indx,i) * (g/dp)
ad_eta(indx-1,i) = ad_eta(indx-1,i) &
+ ad_dellau(indx,i)*(ucko(indx-1,i)-dvu1)*(g/dp)
ad_uol(indx-1,i) = ad_uol(indx-1,i) + ad_dvu1
ad_dvq1 = ad_dvq1 &
- eta(indx-1,i) * ad_dellaq(indx,i) * (g/dp)
ad_qcko(indx-1,i) = ad_qcko(indx-1,i) &
+ eta(indx-1,i) * ad_dellaq(indx,i) * (g/dp)
ad_eta(indx-1,i) = ad_eta(indx-1,i) &
+ ad_dellaq(indx,i)*(qcko(indx-1,i)-dvq1)*(g/dp)
ad_qol(indx-1,i) = ad_qol(indx-1,i) + ad_dvq1
ad_dv1 = ad_dv1 &
- eta(indx-1,i) * ad_dellah(indx,i) * (g/dp)
ad_hcko1(indx-1,i) = ad_hcko1(indx-1,i) &
+ eta(indx-1,i) * ad_dellah(indx,i) * (g/dp)
ad_eta(indx-1,i) = ad_eta(indx-1,i) &
+ ad_dellah(indx,i)*(hcko1(indx-1,i)-dv1)*(g/dp)
ad_heol(indx-1,i) = ad_heol(indx-1,i) + ad_dv1
endif
end do
!c
!c Adjoint of dellah and dellaq calculations within cloud
!c
do i = 1,im
do k = kmax-1,2,-1
if ( k.lt.ktcon(i) .and. ktcon(i).ne.1 ) then
!c
!c Redo nlm calculations
aup = 1.
if (k.le.kb(i)) aup = 0.
adw = 1.
if (k.gt.jmin(i)) adw = 0.
!c
dv1 = heol(k,i)
dv2 = 0.5 * (heol(k,i) + heol(k+1,i))
dv3 = heol(k-1,i)
dv1q = qol(k,i)
dv2q = 0.5 * (qol(k,i) + qol(k+1,i))
dv3q = qol(k-1,i)
DV1u = uOL(k,i)
DV2u = 0.5 * (uOL(k,i) + uOL(k+1,i))
DV3u = uOL(k-1,i)
DV1v = vOL(k,i)
DV2v = 0.5 * (vOL(k,i) + vOL(k+1,i))
DV3v = vOL(k-1,i)
dp = 100. * psfc(i) * del(k)
detau = eta(k,i) - eta(k-1,i)
detad = etad(k,i) - etad(k-1,i)
!c
term1 = aup*eta(k,i) - adw*edto1(i)*etad(k,i)
term2 = aup*detau
term3 = aup*eta(k-1,i) - adw*edto1(i)*etad(k-1,i)
term4 = adw*edto1(i)*detad
termq = 0.5*(qrcdo(k,i)+qrcdo(k-1,i))
!c
!c Initialize loop local adjoint variables to 0.0
ad_termq = 0.0
ad_term4 = 0.0
ad_term3 = 0.0
ad_term2 = 0.0
ad_term1 = 0.0
ad_dv3v = 0.0
ad_dv2v = 0.0
ad_dv1v = 0.0
ad_dv3u = 0.0
ad_dv2u = 0.0
ad_dv1u = 0.0
ad_dv3q = 0.0
ad_dv2q = 0.0
ad_dv1q = 0.0
ad_dv3 = 0.0
ad_dv2 = 0.0
ad_dv1 = 0.0
ad_detad = 0.0
ad_detau = 0.0
!c
!c Adjoint of tlm code.
ad_vcdo(i) = ad_vcdo(i) + (g/dp)*term4*ad_dellav(k,i)
ad_term4 = ad_term4 + (g/dp)*ad_dellav(k,i)*vcdo(i)
ad_dv2v = ad_dv2v - (g/dp)*term2*ad_dellav(k,i)
ad_term2 = ad_term2 - (g/dp)*ad_dellav(k,i)*dv2v
ad_dv3v = ad_dv3v - (g/dp)*term3*ad_dellav(k,i)
ad_term3 = ad_term3 - (g/dp)*ad_dellav(k,i)*dv3v
ad_dv1v = ad_dv1v + (g/dp)*term1*ad_dellav(k,i)
ad_term1 = ad_term1 + (g/dp)*ad_dellav(k,i)*dv1v
ad_ucdo(i) = ad_ucdo(i) + (g/dp)*term4*ad_dellau(k,i)
ad_term4 = ad_term4 + (g/dp)*ad_dellau(k,i)*ucdo(i)
ad_dv2u = ad_dv2u - (g/dp)*term2*ad_dellau(k,i)
ad_term2 = ad_term2 - (g/dp)*ad_dellau(k,i)*dv2u
ad_dv3u = ad_dv3u - (g/dp)*term3*ad_dellau(k,i)
ad_term3 = ad_term3 - (g/dp)*ad_dellau(k,i)*dv3u
ad_dv1u = ad_dv1u + (g/dp)*term1*ad_dellau(k,i)
ad_term1 = ad_term1 + (g/dp)*ad_dellau(k,i)*dv1u
ad_termq = ad_termq + (g/dp)*term4*ad_dellaq(k,i)
ad_term4 = ad_term4 + (g/dp)*termq*ad_dellaq(k,i)
ad_dv2q = ad_dv2q - (g/dp)*term2*ad_dellaq(k,i)
ad_term2 = ad_term2 - (g/dp)*dv2q *ad_dellaq(k,i)
ad_dv3q = ad_dv3q - (g/dp)*term3*ad_dellaq(k,i)
ad_term3 = ad_term3 - (g/dp)*dv3q *ad_dellaq(k,i)
ad_dv1q = ad_dv1q + (g/dp)*term1*ad_dellaq(k,i)
ad_term1 = ad_term1 + (g/dp)*dv1q *ad_dellaq(k,i)
!c
ad_hcdo(i) = ad_hcdo(i)+(g/dp)*term4 *ad_dellah(k,i)
ad_term4 = ad_term4 +(g/dp)*hcdo(i)*ad_dellah(k,i)
ad_dv2 = ad_dv2 - (g/dp)*term2*ad_dellah(k,i)
ad_term2 = ad_term2 - (g/dp)*dv2 *ad_dellah(k,i)
ad_dv3 = ad_dv3 - (g/dp)*term3*ad_dellah(k,i)
ad_term3 = ad_term3 - (g/dp)*dv3 *ad_dellah(k,i)
ad_dv1 = ad_dv1 + (g/dp)*term1*ad_dellah(k,i)
ad_term1 = ad_term1 + (g/dp)*dv1 *ad_dellah(k,i)
!c
ad_qrcdo(k-1,i) = ad_qrcdo(k-1,i) + 0.5*ad_termq
ad_qrcdo(k,i) = ad_qrcdo(k,i) + 0.5*ad_termq
ad_detad = ad_detad + adw*edto1(i)*ad_term4
ad_edto1(i) = ad_edto1(i) + adw*ad_term4*detad
ad_etad(k-1,i) = ad_etad(k-1,i) - adw*edto1(i)*ad_term3
ad_edto1(i) = ad_edto1(i) - adw*ad_term3*etad(k-1,i)
ad_eta(k-1,i) = ad_eta(k-1,i) + aup*ad_term3
ad_detau = ad_detau + aup*ad_term2
ad_etad(k,i) = ad_etad(k,i) - adw*edto1(i)*ad_term1
ad_edto1(i) = ad_edto1(i) - adw*ad_term1*etad(k,i)
ad_eta(k,i) = ad_eta(k,i) + aup*ad_term1
ad_etad(k-1,i) = ad_etad(k-1,i) - ad_detad
ad_etad(k,i) = ad_etad(k,i) + ad_detad
ad_eta(k-1,i) = ad_eta(k-1,i) - ad_detau
ad_eta(k,i) = ad_eta(k,i) + ad_detau
ad_uol(k-1,i) = ad_uol(k-1,i) + ad_dv3u
ad_uol(k+1,i) = ad_uol(k+1,i) + 0.5*ad_dv2u
ad_uol(k,i) = ad_uol(k,i) + 0.5*ad_dv2u
ad_uol(k,i) = ad_uol(k,i) + ad_dv1u
ad_vol(k-1,i) = ad_vol(k-1,i) + ad_dv3v
ad_vol(k+1,i) = ad_vol(k+1,i) + 0.5*ad_dv2v
ad_vol(k,i) = ad_vol(k,i) + 0.5*ad_dv2v
ad_vol(k,i) = ad_vol(k,i) + ad_dv1v
ad_qol(k-1,i) = ad_qol(k-1,i) + ad_dv3q
ad_qol(k+1,i) = ad_qol(k+1,i) + 0.5*ad_dv2q
ad_qol(k,i) = ad_qol(k,i) + 0.5*ad_dv2q
ad_qol(k,i) = ad_qol(k,i) + ad_dv1q
ad_heol(k-1,i) = ad_heol(k-1,i) + ad_dv3
ad_heol(k+1,i) = ad_heol(k+1,i) + 0.5*ad_dv2
ad_heol(k,i) = ad_heol(k,i) + 0.5*ad_dv2
ad_heol(k,i) = ad_heol(k,i) + ad_dv1
!c
endif
end do
end do
!c
!c Adjoint of dellah and dellaq calculations at k=1
do i = 1,im
if (ktcon(i).ne.1) then
dp = 100. * psfc(i) * del(1)
ad_vol(1,i) = ad_vol(1,i) &
- edto1(i)*etad(1,i)*ad_dellav(1,i)*g/dp
ad_vcdo(i) = ad_vcdo(i) &
+ edto1(i)*etad(1,i)*ad_dellav(1,i)*g/dp
ad_etad(1,i) = ad_etad(1,i) &
+ edto1(i)*ad_dellav(1,i)*(vcdo(i)-vol(1,i))*g/dp
ad_edto1(i) = ad_edto1(i) &
+ ad_dellav(1,i)*etad(1,i)*(vcdo(i)-vol(1,i))*g/dp
ad_uol(1,i) = ad_uol(1,i) &
- edto1(i)*etad(1,i)*ad_dellau(1,i)*g/dp
ad_ucdo(i) = ad_ucdo(i) &
+ edto1(i)*etad(1,i)*ad_dellau(1,i)*g/dp
ad_etad(1,i) = ad_etad(1,i) &
+ edto1(i)*ad_dellau(1,i)*(ucdo(i)-uol(1,i))*g/dp
ad_edto1(i) = ad_edto1(i) &
+ ad_dellau(1,i)*etad(1,i)*(ucdo(i)-uol(1,i))*g/dp
ad_qol(1,i) = ad_qol(1,i) &
- edto1(i)*etad(1,i)*ad_dellaq(1,i)*g/dp
ad_qcdo(i) = ad_qcdo(i) &
+ edto1(i)*etad(1,i)*ad_dellaq(1,i)*g/dp
ad_etad(1,i) = ad_etad(1,i) &
+ edto1(i)*ad_dellaq(1,i)*(qcdo(i)-qol(1,i))*g/dp
ad_edto1(i) = ad_edto1(i) &
+ ad_dellaq(1,i)*etad(1,i)*(qcdo(i)-qol(1,i))*g/dp
ad_heol(1,i) = ad_heol(1,i) &
- edto1(i)*etad(1,i)*ad_dellah(1,i)*g/dp
ad_hcdo(i) = ad_hcdo(i) &
+ edto1(i)*etad(1,i)*ad_dellah(1,i)*g/dp
ad_etad(1,i) = ad_etad(1,i) &
+ edto1(i)*ad_dellah(1,i)*(hcdo(i)-heol(1,i))*g/dp
ad_edto1(i) = ad_edto1(i) &
+ ad_dellah(1,i)*etad(1,i)*(hcdo(i)-heol(1,i))*g/dp
endif
end do
!c
!c Adjoint of cloud mass flux threshold function based
!c on cloud work function after downdraft
do i = 1,im
ad_term = 0.0
term = 1.e3*aa1(i)
if (abs(term).lt.100.) then
ad_term = dftanh(term)*ad_tcwfdn(i)
else
if (term.lt.-100.) then
ad_term = 0.0
else
ad_term = 0.0
endif
endif
ad_aa1(i) = ad_aa1(i) + 1.e3*ad_term
end do
!c
!c Adjoint of downdraft cloudwork function calculation
do i = 1,im
do k = 1,kmax-1
if (dwnflg2(i) .and. k.lt.jmin(i)) then
!c
!c Redo nlm calculations
GAMMA = EL2ORC * QESOL(k+1,i) / TOL(k+1,i)**2
DHH = HCDO(I)
DT = TOL(k+1,i)
DG = GAMMA
DH = HESOL(k+1,i)
dz = zo(k,i) - zo(k+1,i)
!c
term1 = edto1(i)
term2 = dz
term3 = g/(cp*dt)
term4 = (dhh-dh)
term5 = 1./(1.+dg)
term6 = 1. + delta*cp*dg*dt/hvap
!c
dqs = qesol(k+1,i) - qol(k+1,i)
!c
!c Adjoint of tlm code
ad_dqs = 0.0
ad_dz = 0.0
ad_dt = 0.0
ad_dg = 0.0
ad_dh = 0.0
ad_dhh = 0.0
if (dqs.gt.0.) then
ad_dqs = ad_dqs + edto1(i)*dz*g*delta * ad_aa1(i)
ad_dz = ad_dz + edto1(i)*g*delta*dqs * ad_aa1(i)
ad_edto1(i) = ad_edto1(i) + dz*g*delta*dqs * ad_aa1(i)
endif
ad_qol(k+1,i) = ad_qol(k+1,i) - ad_dqs
ad_qesol(k+1,i) = ad_qesol(k+1,i) + ad_dqs
!c
ad_term1 = term2*term3*term4*term5*term6 * ad_aa1(i)
ad_term2 = term1*term3*term4*term5*term6 * ad_aa1(i)
ad_term3 = term1*term2*term4*term5*term6 * ad_aa1(i)
ad_term4 = term1*term2*term3*term5*term6 * ad_aa1(i)
ad_term5 = term1*term2*term3*term4*term6 * ad_aa1(i)
ad_term6 = term1*term2*term3*term4*term5 * ad_aa1(i)
!c
ad_dt = (delta*cp/hvap) * dg*ad_term6
ad_dg = (delta*cp/hvap) * ad_term6*dt
ad_dg = ad_dg - 1./(1.+dg)**2 * ad_term5
ad_dh = ad_dh - ad_term4
ad_dhh = ad_dhh + ad_term4
ad_dt = ad_dt - 1.*g/(cp*dt*dt) * ad_term3
ad_dz = ad_dz + ad_term2
ad_edto1(i) = ad_edto1(i) + ad_term1
!c
ad_zo(k,i) = ad_zo(k,i) + ad_dz
ad_zo(k+1,i) = ad_zo(k+1,i) - ad_dz
ad_hesol(k+1,i)= ad_hesol(k+1,i) + ad_dh
ad_gamma = ad_dg
ad_tol(k+1,i) = ad_tol(k+1,i) + ad_dt
ad_hcdo(i) = ad_hcdo(i) + ad_dhh
ad_tol(k+1,i) = ad_tol(k+1,i) &
- 2.*el2orc*qesol(k+1,i)/tol(k+1,i)**3 * ad_gamma
ad_qesol(k+1,i)= ad_qesol(k+1,i) + &
el2orc*ad_gamma/tol(k+1,i)**2
endif
end do
end do
!c
!c Adjoint of final edto1 calculation
do i = 1,im
edtmax = edtmaxl
if (nint(slimsk(i)) .eq. 0) edtmax = edtmaxs
if (dwnflg2(i)) then
if (pwevo(i).lt.0.) then
if (edto10(i).gt.edtmax) then
!c do nothing since tlm is 0.0
ad_edto1(i) = 0.0
endif
ad_pwevo(i) = ad_pwevo(i) &
+ edto(i)*pwavo(i)/pwevo(i)**2 * ad_edto1(i)
ad_pwavo(i) = ad_pwavo(i) &
- edto(i)*ad_edto1(i)/pwevo(i)
ad_edto(i) = ad_edto(i) &
- ad_edto1(i)*pwavo(i)/pwevo(i)
else
!c do nothing since tlm is 0.0
ad_edto1(i) = 0.0
endif
else
!c do nothing since tlm is 0.0
ad_edto1(i) = 0.0
endif
enddo
!c
!c Adjoint of qrcdo --> qcdo transfer
do i = 1,im
ad_qrcdo(1,i) = ad_qrcdo(1,i) + ad_qcdo(i)
end do
!c
!c Adjoint of pwevo, pwdo, qrcdo calculation
do i = 1,im
do k = 1,kmax-1
if (k.lt.jmin(i)) then
!c
!c Recalculate part of foward model
dq = qesol(k,i)
dt = tol(k,i)
gamma = el2orc * dq / dt**2
dh = hcdo(i) - hesol(k,i)
term5 = (1./hvap) * gamma
term6 = 1./(1.+gamma)
detad = etad(k+1,i) - etad(k,i)
!c
!c Adjoint of forward tlm
ad_term1 = 0.0
ad_term2 = 0.0
ad_term3 = 0.0
ad_pwdo(k,i) = ad_pwdo(k,i) + ad_pwevo(i)
if (k.lt.jmin(i)-1) then
ad_term1 = ad_term1 + ad_pwdo(k,i)
ad_term2 = ad_term2 - ad_pwdo(k,i)
ad_term3 = ad_term3 - ad_pwdo(k,i)
ad_qrcdo(k+1,i) = ad_qrcdo(k+1,i) + detad*0.5*ad_term3
ad_qrcdo(k,i) = ad_qrcdo(k,i) + detad*0.5*ad_term3
ad_detad = +0.5*(qrcdo(k,i)+qrcdo(k+1,i))*ad_term3
ad_qrcdo(k,i) = ad_qrcdo(k,i) + etad(k,i)*ad_term2
ad_etad(k,i) = ad_etad(k,i) + qrcdo(k,i)*ad_term2
ad_qrcdo(k+1,i) =ad_qrcdo(k+1,i) +ad_term1*etad(k+1,i)
ad_etad(k+1,i) =ad_etad(k+1,i) +ad_term1*qrcdo(k+1,i)
else
ad_term1 = ad_term1 + ad_pwdo(k,i)
ad_term2 = ad_term2 - ad_pwdo(k,i)
ad_term3 = ad_term3 - ad_pwdo(k,i)
ad_qesol(k+1,i) = ad_qesol(k+1,i) + detad*0.5*ad_term3
ad_qrcdo(k,i) = ad_qrcdo(k,i) + detad*0.5*ad_term3
ad_detad = ad_term3*0.5*(qrcdo(k,i)+qesol(k+1,i))
ad_qrcdo(k,i) = ad_qrcdo(k,i) + etad(k,i)*ad_term2
ad_etad(k,i) = ad_etad(k,i) + qrcdo(k,i)*ad_term2
ad_qol(k+1,i) = ad_qol(k+1,i) + etad(k+1,i)*ad_term1
ad_etad(k+1,i) = ad_etad(k+1,i) + qol(k+1,i)*ad_term1
endif
ad_etad(k,i) = ad_etad(k,i) - ad_detad
ad_etad(k+1,i) = ad_etad(k+1,i) + ad_detad
ad_dh = term5 * term6 * ad_qrcdo(k,i)
ad_term6 = term5 * ad_qrcdo(k,i) * dh
ad_term5 = ad_qrcdo(k,i) * term6 * dh
ad_dq = ad_qrcdo(k,i)
ad_gamma = -1./(1.+gamma)**2 * ad_term6
ad_gamma = ad_gamma + (1./hvap) * ad_term5
ad_term5 = 0.0
ad_term6 = 0.0
ad_hesol(k,i) = ad_hesol(k,i) - ad_dh
ad_hcdo(i) = ad_hcdo(i) + ad_dh
ad_dh = 0.0
ad_dt = -2.*el2orc*dq/dt**3 * ad_gamma
ad_dq = ad_dq + el2orc*ad_gamma/dt**2
ad_gamma = 0.0
ad_tol(k,i) = ad_tol(k,i) + ad_dt
ad_qesol(k,i) = ad_qesol(k,i) + ad_dq
ad_dt = 0.0
ad_dq = 0.0
endif
end do
end do
!c
!c Adjoint of hcdo and qrcdo transfer
do i = 1,im
jmn = jmin(i)
ad_heol(jmn,i) = ad_heol(jmn,i) + ad_hcdo(i)
ad_qesol(jmn,i) = ad_qesol(jmn,i) + ad_qrcdo(jmn,i)
ad_uol(jmn,i) = ad_uol(jmn,i) + ad_ucdo(i)
ad_vol(jmn,i) = ad_vol(jmn,i) + ad_vcdo(i)
end do
!c
!c Adjoint of etad (downdraft entrainment) calculation at k=1
k = 1
do i = 1,im
ad_dz = 0.0
if (kbdtr(i).gt.1) then
!c
!c Redo nlm
dz = 0.5 * (zo(2,i)-zo(1,i))
!c
!c Adjoint of tlm code.
ad_dz = ad_dz + &
etad(k+1,i)*exp(xlamd(i)*dz) * xlamd(i)*ad_etad(k,i)
ad_xlamd(i) = ad_xlamd(i) + &
etad(k+1,i)*exp(xlamd(i)*dz) * dz*ad_etad(k,i)
ad_etad(k+1,i) = ad_etad(k+1,i) + &
exp(xlamd(i)*dz) * ad_etad(k,i)
ad_zo(1,i) = ad_zo(1,i) - 0.5*ad_dz
ad_zo(2,i) = ad_zo(2,i) + 0.5*ad_dz
endif
end do
!c
!c Adjoint of etad (downdraft entrainment) calculation above k=1
do i = 1,im
do k = 2,kbmax
ad_dz = 0.0
if (k.lt.kbdtr(i)) then
!c
!c Redo nlm
dz = 0.5 * (zo(k+1,i)-zo(k-1,i))
!c
!c Adjoint of tlm code.
ad_dz = ad_dz + &
etad(k+1,i)*exp(xlamd(i)*dz) * xlamd(i)*ad_etad(k,i)
ad_xlamd(i) = ad_xlamd(i) + &
etad(k+1,i)*exp(xlamd(i)*dz) * dz*ad_etad(k,i)
ad_etad(k+1,i) = ad_etad(k+1,i) + &
exp(xlamd(i)*dz) * ad_etad(k,i)
ad_zo(k-1,i) = ad_zo(k-1,i) - 0.5*ad_dz
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_dz
endif
end do
end do
!c
!c Adjoint of xlamd calculation
do i = 1,im
ad_dz = 0.0
beta = betas
if (nint(slimsk(i)) .eq. 1) beta = betal
if (kbdtr(i).gt.1) then
!c
!c Redo nlm calculations
dz = 0.5 * (zo(kbdtr(i),i) + zo(kbdtr(i)-1,i)) &
- zo(1,i)
!c
!c Adjoint code
ad_dz = ad_dz - log(beta)/(dz*dz)*ad_xlamd(i)
ad_zo(1,i) = ad_zo(1,i) - ad_dz
ad_zo(kbdtr(i)-1,i) = ad_zo(kbdtr(i)-1,i) + 0.5*ad_dz
ad_zo(kbdtr(i),i) = ad_zo(kbdtr(i),i) + 0.5*ad_dz
endif
end do
!c
!c Adjoint of upper and lower bounds on precipitation efficiency.
do i = 1,im
ad_edt(i) = ad_edt(i) + ad_edtx(i)
ad_edt(i) = ad_edt(i) + ad_edto(i)
if (edt0(i).gt.0.9) then
ad_edt(i) = 0.0
elseif (edt0(i).lt.0.0) then
ad_edt(i) = 0.0
endif
end do
!c
!c Adjoint for precipitation efficiency (edt) calculation
do i = 1,im
ad_e1 = 0.0
ad_term = 0.0
ad_dz = 0.0
!c
!c Redo nlm calculations
dz = zo(ktcon(i),i)-zo(kb(i),i)
if (abs(dz).gt.tiny) then
term = 1.e3 * vshear(i)/dz
!c
!c Adjoint code
ad_e1 = ad_e1 - ad_edt(i)
ad_term = ad_term + p3*3.*term**2 * ad_e1
ad_term = ad_term + p2*2.*term * ad_e1
ad_term = ad_term + p1 * ad_e1
ad_dz = ad_dz - 1.e3*vshear(i)/(dz*dz) * ad_term
ad_vshear(i) = ad_vshear(i) + 1.e3*ad_term/dz
else
ad_term = 0.0
ad_dz = 0.0
endif
ad_zo(kb(i),i) = ad_zo(kb(i),i) - ad_dz
ad_zo(ktcon(i),i) = ad_zo(ktcon(i),i) + ad_dz
end do
!c
!c Adjoint of windshear calculation
do i = 1,im
do k = kmax,1,-1
ad_shear = 0.0
ad_termv = 0.0
ad_termu = 0.0
if (k.ge.kb(i) .and. k.le.ktcon(i)) then
!c
!c Adjoint code for vshear
ad_shear = ad_shear + ad_vshear(i)
!c
!c Redo nlm calculations
termu = (uol(k+1,i) - uol(k,i))**2
termv = (vol(k+1,i) - vol(k,i))**2
!c
!c Adjoint code for uol and vol
if (termu+termv.gt.tiny) then
ad_termv = ad_termv + &
rcs(i)*(0.5/sqrt(termu+termv)) * ad_shear
ad_termu = ad_termu + &
rcs(i)*(0.5/sqrt(termu+termv)) * ad_shear
else
ad_shear = 0.0
ad_termu = 0.0
ad_termv = 0.0
endif
ad_vol(k,i) = ad_vol(k,i) &
- 2.*(vol(k+1,i)-vol(k,i))*ad_termv
ad_vol(k+1,i) = ad_vol(k+1,i) &
+ 2.*(vol(k+1,i)-vol(k,i))*ad_termv
ad_uol(k,i) = ad_uol(k,i) &
- 2.*(uol(k+1,i)-uol(k,i))*ad_termu
ad_uol(k+1,i) = ad_uol(k+1,i) &
+ 2.*(uol(k+1,i)-uol(k,i))*ad_termu
endif
end do
end do
!c
!c Adjoint of cloud mass flux threshold function based
!c on updraft cloud work function
do i = 1,im
ad_term = 0.0
term = 1.e3*aa10(i)
if (abs(term).lt.100.) then
ad_term = dftanh(term)*ad_tcwfup(i)
else
if (term.lt.-100.) then
ad_term = 0.0
else
ad_term = 0.0
endif
endif
ad_aa1(i) = ad_aa1(i) + 1.e3*ad_term
end do
!c
!c
!c Adjoint of updraft cwf calculation
do k = kmax,1,-1
do i = 1,im
if (k.gt.kbcon(i).and.k.le.ktcon(i)) then
!c
!c Recalculate part of forward model
dz1 = zo(k,i) - zo(k-1,i)
gamma = el2orc * qesol(k-1,i) / (tol(k-1,i)**2)
rfact = 1. + delta*cp*gamma*tol(k-1,i)/hvap
term1 = dz1*g/cp
term2 = 1./tol(k-1,i)
term3 = dbyo(k-1,i)
term4 = 1./(1.+gamma)
term5 = rfact
dqs = qesol(k-1,i) - qol(k-1,i)
!c
!c Adjoint of forward model tlm
ad_dqs = 0.0
ad_dz1 = 0.0
if (dqs.gt.0.) then
ad_dqs = ad_dqs + dz1*g*delta * ad_aa1(i)
ad_dz1 = ad_dz1 + ad_aa1(i) * g*delta*dqs
endif
ad_qol(k-1,i) = ad_qol(k-1,i) - ad_dqs
ad_qesol(k-1,i) = ad_qesol(k-1,i) + ad_dqs
ad_term5 = term1*term2*term3*term4*ad_aa1(i)
ad_term4 = term1*term2*term3*term5*ad_aa1(i)
ad_term3 = term1*term2*term4*term5*ad_aa1(i)
ad_term2 = term1*term3*term4*term5*ad_aa1(i)
ad_term1 = term2*term3*term4*term5*ad_aa1(i)
ad_rfact = ad_term5
ad_gamma = -1./((1.+gamma)**2) * ad_term4
ad_dbyo(k-1,i) = ad_dbyo(k-1,i) + ad_term3
ad_tol(k-1,i) = ad_tol(k-1,i) - &
(1./(tol(k-1,i)**2)) * ad_term2
ad_dz1 = ad_dz1 + ad_term1*g/cp
ad_tol(k-1,i) = ad_tol(k-1,i) + &
delta*cp*gamma*ad_rfact/hvap
ad_gamma = ad_gamma + &
delta*cp*ad_rfact*tol(k-1,i)/hvap
ad_tol(k-1,i) = ad_tol(k-1,i) - &
2.*el2orc*qesol(k-1,i)/(tol(k-1,i)**3)* &
ad_gamma
ad_qesol(k-1,i) = ad_qesol(k-1,i) + &
el2orc*ad_gamma/(tol(k-1,i)**2)
ad_zo(k,i) = ad_zo(k,i) + ad_dz1
ad_zo(k-1,i) = ad_zo(k-1,i) - ad_dz1
endif
end do
end do
! Adjoint of section is for cloud liquid water
if (ncloud.gt.0) then
do i=1,im
k = ktcon(i)
if (cnvflg(i)) then
! Recompute nlm
gamma = EL2ORC * QESOl(K,i) / (TOl(K,i)**2)
QRCH = QESOl(K,i) &
+ GAMMA * DBYO(K,i) / (HVAP * (1. + GAMMA))
DQ = qcko00(K-1,i) - QRCH
! Adjoint of tangent linear model
ad_dq = 0.0
ad_qrch = 0.0
ad_gamma= 0.0
if (dq.gt.0) then
ad_qrch = ad_qrch + ad_qcko(k-1,i)
ad_dq = ad_dq + ad_qlko_ktcon(i)
endif
ad_qrch = ad_qrch - ad_dq
ad_qcko(k-1,i) = ad_qcko(k-1,i) + ad_dq
ad_gamma = ad_gamma + &
ad_qrch*dbyo(k,i)/(hvap*(1+gamma)*(1+gamma))
ad_dbyo(k,i) = ad_dbyo(k,i) + &
gamma*ad_qrch/(hvap*(1.+gamma))
ad_qesol(k,i) = ad_qesol(k,i) + ad_qrch
ad_tol(k,i) = ad_tol(k,i) - ad_gamma * &
2.*el2orc*qesol(k,i)/(tol(k,i)**3)
ad_qesol(k,i) = ad_qesol(k,i) + &
el2orc*ad_gamma/(tol(k,i)**2)
endif
end do
endif
!c
!c Adjoint of pwavo, pwo, qcko, and aa1 calculation
!c
do i = 1,im
do k = kmax,1,-1
if (k.gt.kb(i).and.k.lt.ktcon(i)) then
!c
!c Recalculate part of forward model
factor = eta(k-1,i)/eta(k,i)
onemf = 1.0 - factor
temp = qcko0(k,i)
gamma = el2orc*qesol(k,i)/(tol(k,i)**2)
qrch = qesol(k,i) + &
gamma*dbyo(k,i)/(hvap*(1.+gamma))
dq = dqk(k,i)
!c
if (dq.gt.0.) then
!c
!c Redo nlm
dz = 0.5*(zo(k+1,i)-zo(k-1,i))
dz1 = zo(k,i) - zo(k-1,i)
etah = 0.5*(eta(k,i)+eta(k-1,i))
qlk = dq/(eta(k,i) + etah*c0*dz)
!c
!c Adjoint of tlm
ad_temp = 0.0
ad_qc = ad_qcko(k,i)
ad_qrch = ad_qc
ad_qlk = ad_qc
ad_pwo(k,i) = ad_pwo(k,i) + ad_pwavo(i)
ad_qlk = ad_qlk + etah*c0*dz*ad_pwo(k,i)
ad_dz = etah*c0*ad_pwo(k,i)*qlk
ad_etah = ad_pwo(k,i)*c0*dz*qlk
ad_qlk = ad_qlk - dz1*g*ad_aa1(i)
ad_dz1 = -1.*ad_aa1(i)*g*qlk
ad_dz = ad_dz - &
(dq/((eta(k,i)+etah*c0*dz)**2))* &
etah*c0*ad_qlk
ad_etah = ad_etah - &
(dq/((eta(k,i)+etah*c0*dz)**2))* &
ad_qlk*c0*dz
ad_eta(k,i) = ad_eta(k,i) - &
(dq/((eta(k,i)+etah*c0*dz)**2))* &
ad_qlk
ad_dq = ad_qlk/(eta(k,i)+etah*c0*dz)
ad_eta(k-1,i) = ad_eta(k-1,i) + 0.5*ad_etah
ad_eta(k,i) = ad_eta(k,i) + 0.5*ad_etah
ad_zo(k-1,i) = ad_zo(k-1,i) - ad_dz1
ad_zo(k,i) = ad_zo(k,i) + ad_dz1
ad_zo(k-1,i) = ad_zo(k-1,i) - 0.5*ad_dz
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_dz
else
ad_qc = 0.0
ad_qrch = 0.0
ad_qlk = 0.0
ad_dz = 0.0
ad_etah = 0.0
ad_dz1 = 0.0
ad_dq = 0.0
ad_temp = ad_qcko(k,i)
endif
ad_qrch = ad_qrch - eta(k,i)*ad_dq
ad_temp = ad_temp + eta(k,i)*ad_dq
ad_eta(k,i) = ad_eta(k,i) + ad_dq*(temp-qrch)
!c
ad_gamma = -gamma*dbyo(k,i)/((hvap*(1.+gamma))**2)* &
hvap*ad_qrch
ad_dbyo(k,i) = ad_dbyo(k,i) + &
gamma*ad_qrch/(hvap*(1+gamma))
ad_gamma = ad_gamma + &
ad_qrch*dbyo(k,i)/(hvap*(1.+gamma))
ad_qesol(k,i) = ad_qesol(k,i) + ad_qrch
ad_tol(k,i) = ad_tol(k,i) - &
el2orc*2*qesol(k,i)/(tol(k,i)**3)*ad_gamma
ad_qesol(k,i) = ad_qesol(k,i) + &
el2orc/(tol(k,i)**2)*ad_gamma
ad_qol(k+1,i) = ad_qol(k+1,i) + onemf*0.5*ad_temp
ad_qol(k,i) = ad_qol(k,i) + onemf*0.5*ad_temp
ad_onemf = ad_temp*0.5*(qol(k,i)+qol(k+1,i))
ad_qcko(k-1,i) = ad_qcko(k-1,i) + factor*ad_temp
ad_factor = ad_temp*qcko(k-1,i)
ad_factor = ad_factor - ad_onemf
ad_eta(k,i) = ad_eta(k,i) - &
(eta(k-1,i)/(eta(k,i)**2))*ad_factor
ad_eta(k-1,i) = ad_eta(k-1,i) + &
(1./eta(k,i))*ad_factor
endif
end do
end do
!c
!c Adjoint of qol to qcko transfer
do i = 1,im
indx = kb(i)
ad_qol(indx,i) = ad_qol(indx,i) + ad_qcko(indx,i)
end do
!c
!c Adjoint of hcko,hckod transfer to hcko1 and dbyo calculation
do i = 1,im
do k = kmax-1,2,-1
if (k.gt.kb(i)) then
ad_hesol(k,i) = ad_hesol(k,i) - ad_dbyo(k,i)
ad_hcko1(k,i) = ad_hcko1(k,i) + ad_dbyo(k,i)
endif
if (k.gt.kb(i).and.k.le.kbcon(i)) then
ad_hcko(k,i) = ad_hcko(k,i) + ad_hcko1(k,i)
elseif (k.gt.kbcon(i).and.k.le.ktcon(i)) then
if (.not.dwnflg(i)) then
ad_hcko(k,i) = ad_hcko(k,i) + ad_hcko1(k,i)
else
ad_hckod(k,i) = ad_hckod(k,i) + ad_hcko1(k,i)
endif
elseif (k.gt.ktcon(i)) then
ad_hcko(k,i) = ad_hcko(k,i) + ad_hcko1(k,i)
endif
end do
end do
!c
!c Adjoint of modification of hcko1 by detrainment process
do i = 1,im
do k = kmax-1,2,-1
if (k.gt.kbcon(i) .and. k.le.ktcon(i)) then
!c
!c Redo nlm calculations
factor = eta(k-1,i)/eta(k,i)
onemf = 1. - factor
fuv = etau(k-1,i)/etau(k,i)
onemfu = 1. - fuv
heol2 = 0.5*(heol(k,i) + heol(k+1,i))
uol2 = 0.5*(uol(k,i)+uol(k+1,i))
vol2 = 0.5*(vol(k,i)+vol(k+1,i))
!c
!c Adjoint code
ad_fuv = 0.0
ad_onemfu = 0.0
ad_onemf = 0.0
ad_factor = 0.0
ad_heol2 = 0.0
ad_vol2 = 0.0
ad_uol2 = 0.0
ad_vol2 = ad_vol2 + onemfu*ad_vckod(k,i)
ad_onemfu = ad_onemfu + ad_vckod(k,i)*vol2
ad_vckod(k-1,i) = ad_vckod(k-1,i) + fuv*ad_vckod(k,i)
ad_fuv = ad_fuv + ad_vckod(k,i)*vckod(k-1,i)
ad_uol2 = ad_uol2 + onemfu*ad_uckod(k,i)
ad_onemfu = ad_onemfu + ad_uckod(k,i)*uol2
ad_uckod(k-1,i) = ad_uckod(k-1,i) + fuv*ad_uckod(k,i)
ad_fuv = ad_fuv + ad_uckod(k,i)*uckod(k-1,i)
ad_heol2 = ad_heol2 + onemf*ad_hckod(k,i)
ad_onemf = ad_onemf + heol2*ad_hckod(k,i)
ad_hckod(k-1,i) = ad_hckod(k-1,i) + factor*ad_hckod(k,i)
ad_factor = ad_factor + hckod(k-1,i)*ad_hckod(k,i)
ad_heol(k,i) = ad_heol(k,i) + 0.5*ad_heol2
ad_heol(k+1,i) = ad_heol(k+1,i) + 0.5*ad_heol2
ad_fuv = ad_fuv - ad_onemfu
ad_etau(k-1,i) = ad_etau(k-1,i) + ad_fuv/etau(k,i)
ad_etau(k,i) = ad_etau(k,i) - ad_fuv * &
etau(k-1,i)/(etau(k,i)**2)
ad_factor = ad_factor - ad_onemf
ad_eta(k,i) = ad_eta(k,i) - &
ad_factor*(eta(k-1,i)/(eta(k,i)**2))
ad_eta(k-1,i) = ad_eta(k-1,i) + ad_factor/eta(k,i)
endif
end do
end do
!c
!c Adjoint of hcko to hckod transfer at kbcon
do i = 1,im
if (dwnflg(i)) then
indx = kbcon(i)
ad_hcko(indx,i) = ad_hcko(indx,i) + ad_hckod(indx,i)
ad_ucko(indx,i) = ad_ucko(indx,i) + ad_uckod(indx,i)
ad_vcko(indx,i) = ad_vcko(indx,i) + ad_vckod(indx,i)
endif
end do
!c
!c Adjoint of eta calculation
do i = 1,im
if (dwnflg(i)) then
do k = kmax-1,2,-1
if (k.gt.jmin(i) .and. k.le.kt2(i)) then
! Redo nlm calculation
dz = 0.5*(zo(k+1,i)-zo(k-1,i))
! Adjoint code
ad_dz = 0.0
ad_dz = ad_dz + etau(k-1,i) * &
(xlamdet(i)+xlambu)*ad_etau(k,i)
ad_xlamdet(i) = ad_xlamdet(i) + &
etau(k-1,i)*ad_etau(k,i)*dz
ad_etau(k-1,i) = ad_etau(k-1,i) + &
ad_etau(k,i)*(1.+(xlamdet(i)+xlambu)*dz)
ad_dz = ad_dz + eta(k-1,i)*xlamdet(i)*ad_eta(k,i)
ad_xlamdet(i) = ad_xlamdet(i) + &
eta(k-1,i)*ad_eta(k,i)*dz
ad_eta(k-1,i) = ad_eta(k-1,i) + &
ad_eta(k,i)*(1.+xlamdet(i)*dz)
ad_zo(k-1,i) = ad_zo(k-1,i) - 0.5*ad_dz
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_dz
endif
end do
else
do k=kmax-1,2,-1
if (k.gt.jmin(i) .and. k.le.ktcon0(i)) then
! Redo nlm calculation
DZ = .5 * (ZO(k+1,i) - ZO(k-1,i))
! Adjoint code
ad_dz = 0.0
ad_dz = ad_dz + etau(k-1,i)*xlambu*ad_etau(k,i)
ad_etau(k-1,i) = ad_etau(k-1,i) + &
ad_etau(k,i)*(1.+xlambu*dz)
endif
end do
endif
end do
!c
!c
!c Adjoint of modification of cloud properties below cloud
!c base via the entrainment process
do i = 1,im
tem1 = hcko(jmin(i),i) - hesol(kt2(i),i)
tem2 = sumz(kt2(i),i) * hesol(kt2(i),i) - sumh(kt2(i),i)
ad_tem1 = 0.0
ad_tem2 = 0.0
if (sumz(kt2(i),i) .gt. 0.000001) then
term = 2.3/sumz(kt2(i),i)
if (xlamdet1(i).lt.term) then
ad_xlamdet1(i) = ad_xlamdet1(i) + ad_xlamdet(i)
else
ad_term = 0.0
ad_term = ad_term + ad_xlamdet(i)
ad_sumz(kt2(i),i) = ad_sumz(kt2(i),i) - ad_term * &
2.3/(sumz(kt2(i),i)**2)
endif
else
ad_xlamdet1(i) = ad_xlamdet1(i) + ad_xlamdet(i)
endif
if (xlamdet0(i).ge.0) then
ad_xlamdet0(i) = ad_xlamdet0(i) + ad_xlamdet1(i)
else
ad_xlamdet0(i) = 0.0
endif
if (abs(tem2) .gt. 0.000001) then
ad_tem1 = 0.0
ad_tem2 = 0.0
ad_tem1 = ad_tem1 + ad_xlamdet0(i)/tem2
ad_tem2 = ad_tem2 - tem1/(tem2**2) * ad_xlamdet0(i)
else
ad_tem1 = 0.0
ad_tem2 = 0.0
endif
ad_sumz(kt2(i),i) =ad_sumz(kt2(i),i) + ad_tem2*hesol(kt2(i),i)
ad_hesol(kt2(i),i)=ad_hesol(kt2(i),i)+ ad_tem2*sumz(kt2(i),i)
ad_sumh(kt2(i),i) =ad_sumh(kt2(i),i) - ad_tem2
ad_hesol(kt2(i),i) = ad_hesol(kt2(i),i) - ad_tem1
ad_hcko(jmin(i),i) = ad_hcko(jmin(i),i) + ad_tem1
end do
do i=1,im
do k=kmax-1,2,-1
if (k.gt.jmin(i) .and. k.le.ktcon0(i)) then
ad_heol(k,i) = ad_heol(k,i) + &
0.5 * (zo(k+1,i) - zo(k-1,i)) * ad_sumh(k,i)
ad_zo(k-1,i) = ad_zo(k-1,i) - 0.5*ad_sumh(k,i)*heol(k,i)
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_sumh(k,i)*heol(k,i)
ad_sumh(k-1,i) = ad_sumh(k-1,i) + ad_sumh(k,i)
ad_zo(k-1,i) = ad_zo(k-1,i) - 0.5*ad_sumz(k,i)
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_sumz(k,i)
ad_sumz(k-1,i) = ad_sumz(k-1,i) + ad_sumz(k,i)
endif
end do
end do
! Adjoint of cloud property below cloud base is modified by the
! entrainment process
do i=1,im
do k=kmax-1,2,-1
!c Adjoint of last if-then block in tlm
if (k.gt.kbcon(i)) then
ad_hcko(kbcon(i),i) = ad_hcko(kbcon(i),i) + ad_hcko(k,i)
ad_ucko(kbcon(i),i) = ad_ucko(kbcon(i),i) + ad_ucko(k,i)
ad_vcko(kbcon(i),i) = ad_vcko(kbcon(i),i) + ad_vcko(k,i)
endif
!c
!c Adjoint of first if-then block in tlm
if (k.gt.kb(i).and.k.le.kbcon(i)) then
!c
!c Redo nlm calculations
factor = eta0(k-1,i)/eta0(k,i)
onemf = 1. - factor
!c
!c Adjoint code
ad_factor = 0.0
ad_onemf = 0.0
ad_vol(k+1,i) = ad_vol(k+1,i) + onemf*0.5*ad_vcko(k,i)
ad_vol(k,i) = ad_vol(k,i) + onemf*0.5*ad_vcko(k,i)
ad_onemf = ad_onemf + ad_vcko(k,i)*0.5* &
(vol(k,i)+vol(k+1,i))
ad_vcko(k-1,i) = ad_vcko(k-1,i) + factor*ad_vcko(k,i)
ad_factor = ad_factor + ad_vcko(k,i)*vcko0(k-1,i)
ad_uol(k+1,i) = ad_uol(k+1,i) + onemf*0.5*ad_ucko(k,i)
ad_uol(k,i) = ad_uol(k,i) + onemf*0.5*ad_ucko(k,i)
ad_onemf = ad_onemf + ad_ucko(k,i)*0.5* &
(uol(k,i)+uol(k+1,i))
ad_ucko(k-1,i) = ad_ucko(k-1,i) + factor*ad_ucko(k,i)
ad_factor = ad_factor + ad_ucko(k,i)*ucko0(k-1,i)
ad_heol(k+1,i) = ad_heol(k+1,i) + onemf*0.5*ad_hcko(k,i)
ad_heol(k,i) = ad_heol(k,i) + onemf*0.5*ad_hcko(k,i)
ad_onemf = ad_onemf + ad_hcko(k,i)*0.5* &
(heol(k,i)+heol(k+1,i))
ad_hcko(k-1,i) = ad_hcko(k-1,i) + factor*ad_hcko(k,i)
ad_factor = ad_factor + ad_hcko(k,i)*hcko0(k-1,i)
ad_factor = ad_factor - ad_onemf
ad_eta(k,i) = ad_eta(k,i) - ad_factor* &
eta0(k-1,i)/(eta0(k,i)*eta0(k,i))
ad_eta(k-1,i) = ad_eta(k-1,i) + ad_factor/eta0(k,i)
endif
end do
end do
!c
!c Adjoint of hkbo --> hcko transfer
do i = 1,im
indx = kb(i)
ad_hkbo(i) = ad_hkbo(i) + ad_hcko(indx,i)
ad_ukbo(i) = ad_ukbo(i) + ad_ucko(indx,i)
ad_vkbo(i) = ad_vkbo(i) + ad_vcko(indx,i)
end do
!c
!c Adjoint of updraft mass flux calculation
do i = 1,im
ad_dz = 0.0
if (kb(i).eq.1 .and. kbcon(i).gt.1) then
dz = 0.5*(zo(2,i)-zo(1,i))
ad_eta(1,i) = ad_eta(1,i) + ad_etau(1,i)
ad_dz = ad_dz - &
eta0(2,i)*exp(-xlamb(i)*dz)*xlamb(i)* &
ad_eta(1,i)
ad_xlamb(i) = ad_xlamb(i) - &
eta0(2,i)*exp(-xlamb(i)*dz)*dz*ad_eta(1,i)
ad_eta(2,i) = ad_eta(2,i) + ad_eta(1,i)* &
exp(-xlamb(i)*dz)
ad_zo(1,i) = ad_zo(1,i) - 0.5*ad_dz
ad_zo(2,i) = ad_zo(2,i) + 0.5*ad_dz
endif
end do
do i = 1,im
do k = 2,kbmax
ad_dz = 0.0
if (k.lt.kbcon(i).and.k.ge.kb(i)) then
dz = 0.5 * (zo(k+1,i) - zo(k-1,i))
ad_eta(k,i) = ad_eta(k,i) + ad_etau(k,i)
ad_dz = ad_dz - eta0(k+1,i)*exp(-xlamb(i)*dz) * &
xlamb(i)*ad_eta(k,i)
ad_xlamb(i) = ad_xlamb(i) - &
eta0(k+1,i)*exp(-xlamb(i)*dz)*dz*ad_eta(k,i)
ad_eta(k+1,i) = ad_eta(k+1,i) + ad_eta(k,i)* &
exp(-xlamb(i)*dz)
ad_zo(k-1,i) = ad_zo(k-1,i) - 0.5*ad_dz
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_dz
endif
end do
end do
!c
!c Adjoint of determination of entrainment rate between kb and kbcon
do i = 1,im
!c
!c Recalculate nlm
alpha = alphas
if (nint(slimsk(i)) .eq. 1) alpha = alphal
if (kb(i).eq.1) then
dz = 0.5 * (zo(kbcon(i),i)+zo(kbcon(i)-1,i)) - zo(1,i)
else
dz = 0.5*(zo(kbcon(i),i)+zo(kbcon(i)-1,i)) - &
0.5*(zo(kb(i),i) +zo(kb(i)-1,i))
endif
!c
!c Adjoint of tlm code.
ad_dz = 0.0
if (kbcon(i).ne.kb(i)) then
ad_dz = ad_dz + 2.*log(alpha)/(dz*dz) * ad_xlamb(i)
else
ad_dz = 0.0
endif
if (kb(i).eq.1) then
ad_zo(1,i) = ad_zo(1,i) - ad_dz
ad_zo(kbcon(i)-1,i) = ad_zo(kbcon(i)-1,i) + 0.5*ad_dz
ad_zo(kbcon(i),i) = ad_zo(kbcon(i),i) + 0.5*ad_dz
else
ad_zo(kb(i)-1,i) = ad_zo(kb(i)-1,i) - 0.5*ad_dz
ad_zo(kb(i),i) = ad_zo(kb(i),i) - 0.5*ad_dz
ad_zo(kbcon(i),i) = ad_zo(kbcon(i),i) + 0.5*ad_dz
ad_zo(kbcon(i)-1,i) = ad_zo(kbcon(i)-1,i) + 0.5*ad_dz
endif
end do
!c
!c Adjoint of dot to pdot transfer at kbcon
do i = 1,im
ad_dot0(kbcon(i),i) = ad_dot0(kbcon(i),i) + 10.*ad_pdot(i)
end do
!c
!c Adjoint of heol --> hkbo transfer
do i = 1,im
indx = kb(i)
ad_heol(indx,i) = ad_heol(indx,i) + ad_hkbo(i)
ad_uol(indx,i) = ad_uol(indx,i) + ad_ukbo(i)
ad_vol(indx,i) = ad_vol(indx,i) + ad_vkbo(i)
end do
!c
!c Adjoint of heol and hesol calculation at kmax
k = kmax
do i = 1,im
ad_heo(k,i) = ad_heo(k,i) + ad_heol(k,i)
ad_heso(k,i) = ad_heso(k,i) + ad_hesol(k,i)
end do
!c
! Adjoint of heol and hesol calculation below kmax
do i = 1,im
do k = kmax-1,1,-1
!c
!c Make necessary nlm calculations
PO = 0.5 * ( call adfpvsx(tol(k,i),es0,adt,ad_es0,adjoint0)
esl = 10.*es0
!c
!c Zero local adjoint variables
ad_es0 = 0.0
ad_esl = 0.0
!c
!c Adjoint of tlm code
ad_vo(k+1,i) = ad_vo(k+1,i) + 0.5*ad_vol(k,i)
ad_vo(k,i) = ad_vo(k,i) + 0.5*ad_vol(k,i)
ad_uo(k+1,i) = ad_uo(k+1,i) + 0.5*ad_uol(k,i)
ad_uo(k,i) = ad_uo(k,i) + 0.5*ad_uol(k,i)
ad_qesol(k,i) = ad_qesol(k,i) + hvap*ad_hesol(k,i)
ad_tol(k,i) = ad_tol(k,i) + cp*ad_hesol(k,i)
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*g*ad_hesol(k,i)
ad_zo(k,i) = ad_zo(k,i) + 0.5*g*ad_hesol(k,i)
ad_qol(k,i) = ad_qol(k,i) + hvap*ad_heol(k,i)
ad_tol(k,i) = ad_tol(k,i) + cp*ad_heol(k,i)
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*g*ad_heol(k,i)
ad_zo(k,i) = ad_zo(k,i) + 0.5*g*ad_heol(k,i)
ad_esl = ad_esl - eps*esl/((po+epsm1*esl)**2) * &
epsm1 * ad_qesol(k,i)
ad_esl = ad_esl + eps*ad_qesol(k,i)/(po+epsm1*esl)
ad_es0 = ad_es0 + 10*ad_esl
adt = 0.0
call adfpvsx(tol(k,i),es0,adt,ad_es0,adjoint)
ad_tol(k,i) = ad_tol(k,i) + adt
end do
end do
!c
!c Adjoint of tol and qol calculation
do i = 1,im
do k = kmax-1,1,-1
!c
!c Make necessary nlm calculations
DZ = 0.5 * (ZO(k+1,i) - ZO(k,i))
DP = 0.5 * (!! ES = 10. * FESB(TO(k+1,i))
call adfpvsx(to(k+1,i),es0,adt,ad_es0,adjoint0)
es = 10*es0
PPRIME= P(k+1,i) + EPSM1 * ES
QS = EPS * ES / PPRIME
DQSDP = - QS / PPRIME
DESDT = ES * (FACT1 / TO(k+1,i) + FACT2 / (TO(k+1,i)**2))
DQSDT = QS * P(k+1,i) * DESDT / (ES * PPRIME)
GAMMA = EL2ORC * QESO(k+1,i) / (TO(k+1,i)**2)
DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
!c
!c Initialize local ajm variables to 0.0
ad_es0 = 0.0
ad_dq = 0.0
ad_dt = 0.0
ad_dqsdp = 0.0
ad_dqsdt = 0.0
ad_gamma = 0.0
ad_dz = 0.0
ad_pprime = 0.0
ad_es = 0.0
ad_desdt = 0.0
ad_qs = 0.0
!c
!c Adjoint of tlm code
if (qol0(k,i)>0.0) then
ad_qol0(k,i) = ad_qol0(k,i) + ad_qol(k,i)
else
ad_qol0(k,i) = 0.0
endif
ad_qo(k+1,i) = ad_qo(k+1,i) + ad_qol0(k,i)
ad_dq = ad_dq + ad_qol0(k,i)
ad_to(k+1,i) = ad_to(k+1,i) + ad_tol(k,i)
ad_dt = ad_dt + ad_tol(k,i)
ad_dqsdp = ad_dqsdp + ad_dq*dp
ad_dt = ad_dt + dqsdt*ad_dq
ad_dqsdt = ad_dqsdt + ad_dq*dt
ad_gamma = ad_gamma - &
(g*dz+hvap*dqsdp*dp)/((cp*(1.+gamma))**2) * &
cp*ad_dt
ad_dqsdp = ad_dqsdp + hvap*ad_dt*dp/(cp*(1.+gamma))
ad_dz = ad_dz + g*ad_dt/(cp*(1.+gamma))
ad_to(k+1,i) = ad_to(k+1,i) - &
2.*el2orc*qeso(k+1,i)/(to(k+1,i)**3) * ad_gamma
ad_qeso(k+1,i) = ad_qeso(k,i) + &
el2orc*ad_gamma/(to(k+1,i)**2)
rterm = 1./(es*pprime)
ad_pprime = ad_pprime - (qs*p(k+1,i)*desdt)*rterm**2 * &
es*ad_dqsdt
ad_es = ad_es - (qs*p(k+1,i)*desdt)*rterm**2 * &
ad_dqsdt*pprime
ad_desdt = ad_desdt + p(k+1,i)*qs*ad_dqsdt*rterm
ad_qs = ad_qs + p(k+1,i)*ad_dqsdt*desdt*rterm
ad_to(k+1,i) = ad_to(k+1,i) + ad_desdt * es * &
(-1.*fact1/(to(k+1,i)**2) - 2.*fact2/(to(k+1,i)**3))
ad_es = ad_es + ad_desdt * &
(fact1/to(k+1,i) + fact2/(to(k+1,i)**2))
ad_pprime = ad_pprime + (qs/(pprime**2))*ad_dqsdp
ad_qs = ad_qs - ad_dqsdp/pprime
ad_pprime = ad_pprime - (eps*es/(pprime**2))*ad_qs
ad_es = ad_es + eps*ad_qs/pprime
ad_es = ad_es + epsm1*ad_pprime
ad_es0 = ad_es0 + 10.*ad_es
adt=0.0
call adfpvsx(to(k+1,i),es0,adt,ad_es0,adjoint0)
ad_to(k+1,i) = ad_to(k+1,i) + adt
ad_zo(k,i) = ad_zo(k,i) - 0.5*ad_dz
ad_zo(k+1,i) = ad_zo(k+1,i) + 0.5*ad_dz
end do
end do
!c
!c
!c Adjoint of moist static energy calculation
do i = 1,im
do k = 1,kmax
ad_qeso(k,i) = ad_qeso(k,i) + hvap*ad_heso(k,i)
ad_to(k,i) = ad_to(k,i) + cp*ad_heso(k,i)
ad_zo(k,i) = ad_zo(k,i) + g*ad_heso(k,i)
ad_qo(k,i) = ad_qo(k,i) + hvap*ad_heo(k,i)
ad_to(k,i) = ad_to(k,i) + cp*ad_heo(k,i)
ad_zo(k,i) = ad_zo(k,i) + g*ad_heo(k,i)
end do
end do
!c
!c
!c Adjoint of zo calculation
do i = 1,im
do k = kmax,2,-1
dlnsig = log(sl(k)/sl(k-1))
term1 = dlnsig * rd / g
ad_zo(k-1,i) = ad_zo(k-1,i) + ad_zo(k,i)
ad_term2 = -1.*term1*ad_zo(k,i)
ad_tvo(k-1,i) = ad_tvo(k-1,i) + 0.5*ad_term2
ad_tvo(k,i) = ad_tvo(k,i) + 0.5*ad_term2
end do
end do
!c
dlnsig = log(sl(1))
do i = 1,im
ad_tvo(1,i) = ad_tvo(1,i) - (dlnsig*rd/g)*ad_zo(1,i)
end do
!c
!c
!c Adjoint of initial qeso and tvo calculation.
do i = 1,im
do k = 1,kmax
call adfpvsx(to(k,i),es0,adt,ad_es0,adjoint0)
es = 10.*es0
ad_qo(k,i) = ad_qo(k,i) + delta*to(k,i)*ad_tvo(k,i)
ad_to(k,i) = ad_to(k,i) + delta*qo(k,i)*ad_tvo(k,i)
ad_to(k,i) = ad_to(k,i) + ad_tvo(k,i)
ad_es = 0.0
ad_es = ad_es + ad_qeso(k,i) * &
eps*p(k,i)/(
ad_es0 = 0.0
ad_es0 = 10.*ad_es
adt=0.0
call adfpvsx(to(k,i),es0,adt,ad_es0,adjoint)
ad_to(k,i) = ad_to(k,i) + adt
end do
end do
!c
!c Adjoint of initial 0 --> o transfer
do i = 1,im
do k = 1,km
ad_q0(k,i) = ad_q0(k,i) + ad_qo0(k,i)
ad_cwm0(k,i) = ad_cwm0(k,i) + ad_cwmo(k,i)
ad_v0(k,i) = ad_v0(k,i) + ad_vo(k,i)
ad_u0(k,i) = ad_u0(k,i) + ad_uo(k,i)
if (q0(k,i)>0.) then
ad_q0(k,i) = ad_q0(k,i) + ad_qo(k,i)
else
ad_q0(k,i) = 0.0
endif
ad_t0(k,i) = ad_t0(k,i) + ad_to(k,i)
end do
end do
!c
!c End of routine
RETURN
END