FV3 Bundle
cloud_ad.F90
Go to the documentation of this file.
1 MODULE cloud_ad
2 
3 USE cloud
4 USE cloud_tl, only: ls_cloud_d
5 
6 IMPLICIT NONE
7 
8 PRIVATE
9 PUBLIC cloud_driver_b
10 
11 CONTAINS
12 ! Generated by TAPENADE (INRIA, Tropics team)
13 ! Tapenade 3.9 (r5096) - 24 Feb 2014 16:53
14 !
15 ! Differentiation of cloud_driver in reverse (adjoint) mode:
16 ! gradient of useful results: th qi_con q cf_ls cf_con ql_ls
17 ! ql_con qi_ls
18 ! with respect to varying inputs: th qi_con cnv_prc3 q cnv_dqldt
19 ! cf_ls cnv_updf cf_con cnv_mfd ql_ls ql_con qi_ls
20 ! RW status of diff variables: th:in-out qi_con:in-out cnv_prc3:out
21 ! q:in-out cnv_dqldt:out cf_ls:in-out cnv_updf:out
22 ! cf_con:in-out cnv_mfd:out ql_ls:in-out ql_con:in-out
23 ! qi_ls:in-out
24 SUBROUTINE cloud_driver_b(dt, im, jm, lm, th, thb, q, qb, ple, cnv_dqldt&
25 & , cnv_dqldtb, cnv_mfd, cnv_mfdb, cnv_prc3, cnv_prc3b, cnv_updf, &
26 & cnv_updfb, qi_ls, qi_lsb, ql_ls, ql_lsb, qi_con, qi_conb, ql_con, &
27 & ql_conb, cf_ls, cf_lsb, cf_con, cf_conb, frland, physparams, estblx, &
28 & khu, khl, cons_runiv, cons_kappa, cons_airmw, cons_h2omw, cons_grav, &
29 & cons_alhl, cons_alhf, cons_pi, cons_rgas, cons_cp, cons_vireps, &
30 & cons_alhs, cons_tice, cons_rvap, cons_p00, do_moist_physics)
31  IMPLICIT NONE
32 !INPUTS
33  INTEGER, INTENT(IN) :: im, jm, lm, do_moist_physics
34  REAL*8, INTENT(IN) :: dt, frland(im, jm), physparams(:)
35  REAL*8, DIMENSION(im, jm, lm), INTENT(IN) :: cnv_dqldt, cnv_mfd, &
36 & cnv_updf, cnv_prc3
37  REAL*8, DIMENSION(im, jm, lm) :: cnv_dqldtb, cnv_mfdb, cnv_updfb, &
38 & cnv_prc3b
39  REAL*8, INTENT(IN) :: estblx(:)
40  REAL*8, DIMENSION(im, jm, 0:lm), INTENT(IN) :: ple
41  INTEGER, DIMENSION(im, jm), INTENT(IN) :: khu, khl
42 !MAPL_CONSTANTS REDEFINED FOR USE IN AUTODIFF TOOL
43  REAL*8, INTENT(IN) :: cons_runiv, cons_kappa, cons_airmw
44  REAL*8, INTENT(IN) :: cons_h2omw, cons_grav, cons_alhl
45  REAL*8, INTENT(IN) :: cons_alhf, cons_pi, cons_rgas
46  REAL*8, INTENT(IN) :: cons_cp, cons_vireps, cons_alhs
47  REAL*8, INTENT(IN) :: cons_tice, cons_rvap, cons_p00
48 !PROGNOSTICS
49  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: th, q
50  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: thb
51  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: qi_ls, ql_ls, qi_con, &
52 & ql_con, cf_con, cf_ls
53  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: qi_lsb
54 !OUTPUTS (DIAGNOSTICS)
55 !LOCALS
56  INTEGER :: i, j, k, l, ktop
57  REAL*8, DIMENSION(im, jm, lm) :: ph, pih, mass, imass, t, dzet, qddf3&
58 & , rh, dp, dm
59  REAL*8, DIMENSION(im, jm, lm) :: tb, dzetb, qddf3b
60  REAL*8, DIMENSION(im, jm, lm + 1) :: zet
61  REAL*8, DIMENSION(im, jm, lm+1) :: zetb
62  REAL*8, DIMENSION(im, jm, 0:lm) :: p, pi
63  REAL*8, DIMENSION(im, jm, lm) :: qs, dqsdt, dqs
64  REAL*8, DIMENSION(im, jm, lm) :: qsb, dqsdtb, dqsb
65  REAL*8, DIMENSION(im, jm) :: vmip
66  REAL*8, DIMENSION(im, jm) :: vmipb
67 !Precip amounts and fall rate
68  REAL*8 :: cf_tot
69  REAL*8 :: cf_totb
70  REAL*8 :: alpha, alhx3, rhcrit
71  REAL*8 :: alhx3b
72 !Microphyiscal constants
73  REAL*8 :: aa, bb
74  REAL*8 :: aab, bbb
75  REAL*8, PARAMETER :: pi_0=4.*atan(1.)
76 !Density of liquid water in kg/m^3
77  REAL*8, PARAMETER :: rho_w=1.0e3
78  REAL*8 :: t_ice_max
79 !PHYSPARAMS constants
80  REAL*8 :: cnv_beta, anv_beta, ls_beta, rh00, c_00, lwcrit, c_acc, &
81 & c_ev_r, c_ev_s, cldvol2frc
82  REAL*8 :: rhsup_ice, shr_evap_fac, min_cld_water, cld_evp_eff, &
83 & ls_sdqv2, ls_sdqv3, ls_sdqvt1
84  REAL*8 :: anv_sdqv2, anv_sdqv3, anv_sdqvt1, anv_to_ls, n_warm, n_ice, &
85 & n_anvil, n_pbl
86  REAL*8 :: anv_icefall_c, ls_icefall_c, revap_off_p, cnvenvfc, wrhodep&
87 & , t_ice_all, cnviceparam
88  REAL*8 :: cnvddrfc, anvddrfc, lsddrfc, minrhcrit, maxrhcrit, &
89 & turnrhcrit, maxrhcritland
90  REAL*8 :: min_rl, min_ri, max_rl, max_ri, ri_anv
91  INTEGER :: nsmax, disable_rad, icefrpwr, tanhrhcrit, fr_ls_wat, &
92 & fr_ls_ice, fr_an_wat, fr_an_ice, pdfflag
93  REAL*8 :: lsenvfc, anvenvfc
94  REAL*8 :: qrn_cu, qsn_cu, qrn_an, qsn_an, qrn_ls, qsn_ls, qrn_cu_1d
95  REAL*8 :: qsn_cub, qrn_anb, qsn_anb, qrn_lsb, qsn_lsb, qrn_cu_1db
96  REAL*8 :: qt_tmpi_1, qt_tmpi_2, qlt_tmp, qit_tmp
97  REAL*8 :: qt_tmpi_1b, qt_tmpi_2b, qlt_tmpb, qit_tmpb
98  REAL*8 :: prn_above_cu_new, prn_above_an_new, prn_above_ls_new
99  REAL*8 :: prn_above_cu_newb, prn_above_an_newb, prn_above_ls_newb
100  REAL*8 :: prn_above_cu_old, prn_above_an_old, prn_above_ls_old
101  REAL*8 :: prn_above_cu_oldb, prn_above_an_oldb, prn_above_ls_oldb
102  REAL*8 :: psn_above_cu_new, psn_above_an_new, psn_above_ls_new
103  REAL*8 :: psn_above_cu_newb, psn_above_an_newb, psn_above_ls_newb
104  REAL*8 :: psn_above_cu_old, psn_above_an_old, psn_above_ls_old
105  REAL*8 :: psn_above_cu_oldb, psn_above_an_oldb, psn_above_ls_oldb
106  REAL*8 :: evap_dd_cu_above_new, evap_dd_an_above_new, &
107 & evap_dd_ls_above_new
108  REAL*8 :: evap_dd_cu_above_newb, evap_dd_an_above_newb, &
109 & evap_dd_ls_above_newb
110  REAL*8 :: evap_dd_cu_above_old, evap_dd_an_above_old, &
111 & evap_dd_ls_above_old
112  REAL*8 :: evap_dd_cu_above_oldb, evap_dd_an_above_oldb, &
113 & evap_dd_ls_above_oldb
114  REAL*8 :: subl_dd_cu_above_new, subl_dd_an_above_new, &
115 & subl_dd_ls_above_new
116  REAL*8 :: subl_dd_cu_above_newb, subl_dd_an_above_newb, &
117 & subl_dd_ls_above_newb
118  REAL*8 :: subl_dd_cu_above_old, subl_dd_an_above_old, &
119 & subl_dd_ls_above_old
120  REAL*8 :: subl_dd_cu_above_oldb, subl_dd_an_above_oldb, &
121 & subl_dd_ls_above_oldb
122  REAL*8 :: area_ls_prc1, area_upd_prc1, area_anv_prc1
123  REAL*8 :: area_ls_prc1b, area_upd_prc1b, area_anv_prc1b
124  REAL*8 :: tot_prec_upd, tot_prec_anv, tot_prec_ls, area_upd_prc, &
125 & area_anv_prc, area_ls_prc
126  REAL*8 :: tot_prec_updb, tot_prec_anvb, tot_prec_lsb, area_upd_prcb, &
127 & area_anv_prcb, area_ls_prcb
128  REAL*8 :: qtmp2
129  REAL*8 :: rhexcess, tpw, negtpw
130  REAL*8 :: tpwb, negtpwb
131  INTEGER :: cloud_pertmod
132  INTRINSIC atan
133  INTRINSIC int
134  INTRINSIC sum
135  INTRINSIC max
136  INTEGER :: branch
137  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: qb
138  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: ql_lsb
139  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: ql_conb
140  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: qi_conb
141  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: cf_conb
142  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: cf_lsb
143  REAL*8 :: temp2
144  REAL*8 :: temp1(im, jm, lm)
145  LOGICAL :: mask0(im, jm, lm)
146  REAL*8 :: temp0
147  REAL*8 :: tempb9
148  REAL*8 :: tempb8
149  REAL*8 :: tempb7
150  REAL*8 :: tempb6
151  REAL*8 :: tempb5
152  REAL*8 :: tempb4
153  REAL*8 :: tempb3
154  REAL*8 :: tempb2
155  REAL*8 :: tempb1
156  REAL*8 :: tempb16
157  REAL*8 :: tempb0
158  REAL*8 :: tempb15
159  REAL*8 :: tempb14(im, jm, lm)
160  REAL*8 :: tempb13
161  REAL*8 :: tempb12
162  REAL*8 :: tempb11
163  REAL*8 :: tempb10
164  REAL*8 :: tempb
165  LOGICAL :: mask(im, jm, lm)
166  REAL*8 :: temp
167 
168  !LS_CLOUD FILTERING
169  integer :: ii
170  real(8) :: xx(8)
171  real(8) :: ttraj, qtraj, qi_lstraj, qi_contraj, ql_lstraj, ql_contraj, cf_lstraj, cf_contraj, phtraj
172  real(8) :: tpert, qpert, qi_lspert, qi_conpert, ql_lspert, ql_conpert, cf_lspert, cf_conpert
173  real(8) :: jacobian(8,8), a(8,8)
174 
175  !Eigenvalue computation
176  integer, parameter :: n1 = 8
177  integer, parameter :: lda = 8
178  REAL(8) :: wr(n1), wi(n1)
179  INTEGER, PARAMETER :: ldvl = n1 !Left vecotrs
180  REAL(8) :: vl(ldvl,n1)
181  INTEGER, PARAMETER :: ldvr = n1 !Right vectors
182  REAL(8) :: vr(ldvr,n1)
183  INTEGER :: info, lwork
184  INTEGER, PARAMETER :: lwmax = 10000
185  REAL(8) :: work(lwmax)
186  EXTERNAL :: dgeev
187  REAL*8 :: maxeval
188 
189  !Total Filtering
190  real(8) :: totfilt_t, totfilt_ql, totfilt_qi
191  real(8) :: t_p_preall, ql_ls_p_preall, ql_con_p_preall, qi_ls_p_preall, qi_con_p_preall
192 
193  !Sink Filtering
194  real(8) :: sinkfilt_ql, sinkfilt_qi, sinkfilt_cf
195  real(8) :: t_p_presink, q_p_presink
196  real(8) :: ql_ls_p_presink, ql_con_p_presink
197  real(8) :: qi_ls_p_presink, qi_con_p_presink
198  real(8) :: cf_con_p_presink
199 
200 
201 !Highest level of calculations
202  ktop = 30
203 !Get Constants from CLOUDPARAMS
204 ! Area factor for convective rain showers (non-dim)
205  cnv_beta = physparams(1)
206 ! Area factor for anvil rain showers (non-dim)
207  anv_beta = physparams(2)
208 ! Area factor for Large Scale rain showers (non-dim)
209  ls_beta = physparams(3)
210 ! Critical relative humidity
211  rh00 = physparams(4)
212  c_00 = physparams(5)
213  lwcrit = physparams(6)
214  c_acc = physparams(7)
215  c_ev_r = physparams(8)
216  c_ev_s = physparams(56)
217  cld_evp_eff = physparams(13)
218  ls_sdqv2 = physparams(15)
219  ls_sdqv3 = physparams(16)
220  ls_sdqvt1 = physparams(17)
221  anv_sdqv2 = physparams(18)
222  anv_sdqv3 = physparams(19)
223  anv_sdqvt1 = physparams(20)
224  anv_icefall_c = physparams(28)
225  ls_icefall_c = physparams(29)
226  revap_off_p = physparams(30)
227  cnvenvfc = physparams(31)
228  wrhodep = physparams(32)
229  t_ice_all = physparams(33) + cons_tice
230  icefrpwr = int(physparams(35) + .001)
231  cnvddrfc = physparams(36)
232  anvddrfc = physparams(37)
233  lsddrfc = physparams(38)
234  minrhcrit = physparams(42)
235  maxrhcrit = physparams(43)
236  turnrhcrit = physparams(45)
237  maxrhcritland = physparams(46)
238  pdfflag = int(physparams(57))
239  t_ice_max = cons_tice
240 !Initialize the saving of downdraft values.
241  prn_above_cu_new = 0.
242  prn_above_an_new = 0.
243  prn_above_ls_new = 0.
244  psn_above_cu_new = 0.
245  psn_above_an_new = 0.
246  psn_above_ls_new = 0.
247  evap_dd_cu_above_new = 0.
248  evap_dd_an_above_new = 0.
249  evap_dd_ls_above_new = 0.
250  subl_dd_cu_above_new = 0.
251  subl_dd_an_above_new = 0.
252  subl_dd_ls_above_new = 0.
253 !Convert to hPa and average pressure to temperature levels
254  p = ple*0.01
255  ph = 0.5*(p(:, :, 0:lm-1)+p(:, :, 1:lm))
256 !Calculate Exner Pressure at temperature levels
257  pi = (p/1000.)**(cons_rgas/cons_cp)
258  pih = (ph/1000.)**(cons_rgas/cons_cp)
259 !Calculate temperature
260  t = th*pih
261 !Compute QS and DQSDT
262  CALL dqsat_bac(dqsdt, qs, t, ph, im, jm, lm, estblx, cons_h2omw, &
263 & cons_airmw)
264 !Relative humidity
265 !Compute layer mass and 1/mass
266  mass = (p(:, :, 1:lm)-p(:, :, 0:lm-1))*100./cons_grav
267  imass = 1/mass
268 !Level thickness
269  dzet(:, :, 1:lm) = th(:, :, 1:lm)*(pi(:, :, 1:lm)-pi(:, :, 0:lm-1))*&
270 & cons_cp/cons_grav
271 !Level heights
272  zet(:, :, lm+1) = 0.0
273  DO k=lm,1,-1
274  zet(:, :, k) = zet(:, :, k+1) + dzet(:, :, k)
275  END DO
276  mask(:, :, 1:lm) = zet(:, :, 1:lm) .LT. 3000.
277  WHERE (mask(:, :, 1:lm))
278  qddf3 = -((zet(:, :, 1:lm)-3000.)*zet(:, :, 1:lm)*mass)
279  ELSEWHERE
280  qddf3 = 0.
281  END WHERE
282  DO i=1,im
283  DO j=1,jm
284  vmip(i, j) = sum(qddf3(i, j, :))
285  END DO
286  END DO
287  DO k=1,lm
288  CALL pushreal8array(qddf3(:, :, k), im*jm)
289  qddf3(:, :, k) = qddf3(:, :, k)/vmip
290  END DO
291 !Pressure and mass thickness for use in cleanup.
292  dp = ple(:, :, 1:lm) - ple(:, :, 0:lm-1)
293  dm = dp*(1./cons_grav)
294 !Begin loop over all grid boxes.
295  DO i=1,im
296  DO j=1,jm
297  DO k=ktop,lm
298  IF (k .EQ. ktop) THEN
299  CALL pushreal8(tot_prec_upd)
300  tot_prec_upd = 0.
301  CALL pushreal8(tot_prec_anv)
302  tot_prec_anv = 0.
303  CALL pushreal8(tot_prec_ls)
304  tot_prec_ls = 0.
305  CALL pushreal8(area_upd_prc)
306  area_upd_prc = 0.
307  CALL pushreal8(area_anv_prc)
308  area_anv_prc = 0.
309  CALL pushreal8(area_ls_prc)
310  area_ls_prc = 0.
311  CALL pushcontrol1b(1)
312  ELSE
313  CALL pushcontrol1b(0)
314  END IF
315 !Initialize precips, except QRN_CU which comes from RAS
316  qrn_ls = 0.
317  qrn_an = 0.
318  qsn_ls = 0.
319  qsn_an = 0.
320  qsn_cu = 0.
321 !Ras Rain
322  qrn_cu_1d = cnv_prc3(i, j, k)
323 !Tidy up where fractions or cloud is too low
324  CALL pushreal8(cf_con(i, j, k))
325  CALL pushreal8(qi_con(i, j, k))
326  CALL pushreal8(ql_con(i, j, k))
327  CALL pushreal8(cf_ls(i, j, k))
328  CALL pushreal8(qi_ls(i, j, k))
329  CALL pushreal8(ql_ls(i, j, k))
330  CALL pushreal8(t(i, j, k))
331  CALL cloud_tidy(q(i, j, k), t(i, j, k), ql_ls(i, j, k), qi_ls(i&
332 & , j, k), cf_ls(i, j, k), ql_con(i, j, k), qi_con(i, j&
333 & , k), cf_con(i, j, k), cons_alhl, cons_alhs, cons_cp)
334 !Phase changes for large scale cloud.
335  CALL pushreal8(qi_ls(i, j, k))
336  CALL pushreal8(ql_ls(i, j, k))
337  CALL pushreal8(t(i, j, k))
338  CALL meltfreeze(dt, t(i, j, k), ql_ls(i, j, k), qi_ls(i, j, k), &
339 & t_ice_all, t_ice_max, icefrpwr, cons_alhl, cons_alhs, &
340 & cons_cp)
341 !Phase changes for convective cloud.
342  CALL pushreal8(qi_con(i, j, k))
343  CALL pushreal8(ql_con(i, j, k))
344  CALL pushreal8(t(i, j, k))
345  CALL meltfreeze(dt, t(i, j, k), ql_con(i, j, k), qi_con(i, j, k)&
346 & , t_ice_all, t_ice_max, icefrpwr, cons_alhl, cons_alhs&
347 & , cons_cp)
348 !STAGE 1 - Compute convective clouds from RAS diagnostics
349  CALL pushreal8(cf_con(i, j, k))
350  CALL pushreal8(qi_con(i, j, k))
351  CALL pushreal8(ql_con(i, j, k))
352  CALL pushreal8(q(i, j, k))
353  CALL pushreal8(t(i, j, k))
354  CALL convec_src(dt, mass(i, j, k), imass(i, j, k), t(i, j, k), q&
355 & (i, j, k), cnv_dqldt(i, j, k), cnv_mfd(i, j, k), &
356 & ql_con(i, j, k), qi_con(i, j, k), cf_con(i, j, k), qs(&
357 & i, j, k), cons_alhs, cons_alhl, cons_cp, t_ice_all, &
358 & t_ice_max, icefrpwr)
359 !STAGE 2a - Get PDF attributes
360  CALL pushreal8(alpha)
361  CALL pdf_width(ph(i, j, k), frland(i, j), maxrhcrit, &
362 & maxrhcritland, turnrhcrit, minrhcrit, pi_0, alpha)
363  IF (alpha .LT. 1.0 - rh00) THEN
364  alpha = 1.0 - rh00
365  ELSE
366  alpha = alpha
367  END IF
368  CALL pushreal8(rhcrit)
369  rhcrit = 1.0 - alpha
370  cloud_pertmod = 1
371 !STAGE 2b - Use PDF to compute large scale cloud effects and diagnostics,
372 ! also update convection clouds
373  CALL pushreal8(cf_con(i, j, k))
374  CALL pushreal8(cf_ls(i, j, k))
375  CALL pushreal8(qi_con(i, j, k))
376  CALL pushreal8(qi_ls(i, j, k))
377  CALL pushreal8(ql_con(i, j, k))
378  CALL pushreal8(ql_ls(i, j, k))
379  CALL pushreal8(q(i, j, k))
380  CALL pushreal8(t(i, j, k))
381  CALL ls_cloud(dt, alpha, pdfflag, ph(i, j, k), t(i, j, k), q(i, &
382 & j, k), ql_ls(i, j, k), ql_con(i, j, k), qi_ls(i, j, k), &
383 & qi_con(i, j, k), cf_ls(i, j, k), cf_con(i, j, k), &
384 & cons_alhl, cons_alhf, cons_alhs, cons_cp, cons_h2omw, &
385 & cons_airmw, t_ice_all, t_ice_max, icefrpwr, estblx, &
386 & cloud_pertmod, do_moist_physics)
387 !Clean up where too much overall cloud.
388  CALL pushreal8(cf_tot)
389  cf_tot = cf_ls(i, j, k) + cf_con(i, j, k)
390  IF (cf_tot .GT. 1.00) THEN
391  CALL pushreal8(cf_ls(i, j, k))
392  cf_ls(i, j, k) = cf_ls(i, j, k)*(1.00/cf_tot)
393  CALL pushreal8(cf_con(i, j, k))
394  cf_con(i, j, k) = cf_con(i, j, k)*(1.00/cf_tot)
395  CALL pushcontrol1b(0)
396  ELSE
397  CALL pushcontrol1b(1)
398  END IF
399 !STAGE 3 - Evap, Sublimation and Autoconversion
400 !Evaporation and sublimation of anvil cloud
401  CALL pushreal8(cf_con(i, j, k))
402  CALL pushreal8(ql_con(i, j, k))
403  CALL pushreal8(q(i, j, k))
404  CALL pushreal8(t(i, j, k))
405  CALL evap_cnv(dt, rhcrit, ph(i, j, k), t(i, j, k), q(i, j, k), &
406 & ql_con(i, j, k), qi_con(i, j, k), cf_con(i, j, k), cf_ls&
407 & (i, j, k), qs(i, j, k), rho_w, cld_evp_eff, cons_h2omw, &
408 & cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, &
409 & cons_cp)
410  CALL pushreal8(cf_con(i, j, k))
411  CALL pushreal8(qi_con(i, j, k))
412  CALL pushreal8(q(i, j, k))
413  CALL pushreal8(t(i, j, k))
414  CALL subl_cnv(dt, rhcrit, ph(i, j, k), t(i, j, k), q(i, j, k), &
415 & ql_con(i, j, k), qi_con(i, j, k), cf_con(i, j, k), cf_ls&
416 & (i, j, k), qs(i, j, k), rho_w, cld_evp_eff, cons_h2omw, &
417 & cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, &
418 & cons_cp, cons_alhs)
419 !Autoconversion
420  CALL pushreal8(cf_ls(i, j, k))
421  CALL pushreal8(ql_ls(i, j, k))
422  CALL autoconversion_ls(dt, ql_ls(i, j, k), qrn_ls, t(i, j, k), &
423 & ph(i, j, k), cf_ls(i, j, k), ls_sdqv2, ls_sdqv3&
424 & , ls_sdqvt1, c_00, lwcrit, dzet(i, j, k))
425  CALL pushreal8(ql_con(i, j, k))
426  CALL autoconversion_cnv(dt, ql_con(i, j, k), qrn_an, t(i, j, k)&
427 & , ph(i, j, k), cf_con(i, j, k), anv_sdqv2, &
428 & anv_sdqv3, anv_sdqvt1, c_00, lwcrit, dzet(i, j&
429 & , k))
430 !STAGE 4 - Fall and Re-evaporation of precip
431  CALL pushreal8(qi_con(i, j, k))
432  CALL ice_settlefall_cnv(wrhodep, qi_con(i, j, k), ph(i, j, k), t&
433 & (i, j, k), cf_con(i, j, k), cons_rgas, khu(i, &
434 & j), khl(i, j), k, dt, dzet(i, j, k), qsn_an, &
435 & anv_icefall_c)
436  CALL pushreal8(cf_ls(i, j, k))
437  CALL pushreal8(qi_ls(i, j, k))
438  CALL ice_settlefall_ls(wrhodep, qi_ls(i, j, k), ph(i, j, k), t(i&
439 & , j, k), cf_ls(i, j, k), cons_rgas, khu(i, j), &
440 & khl(i, j), k, dt, dzet(i, j, k), qsn_ls, &
441 & ls_icefall_c)
442 !"Freeze" out any conv. precip, not done in RAS. This is
443 ! precip w/ large particles, so freezing is strict.
444  IF (t(i, j, k) .LT. cons_tice) THEN
445  qsn_cu = qrn_cu_1d
446  qrn_cu_1d = 0.
447  CALL pushreal8(t(i, j, k))
448  t(i, j, k) = t(i, j, k) + qsn_cu*(cons_alhs-cons_alhl)/cons_cp
449  CALL pushcontrol1b(0)
450  ELSE
451  CALL pushcontrol1b(1)
452  END IF
453 !Area
454  CALL pushreal8(area_ls_prc1)
455  area_ls_prc1 = 0.0
456  CALL pushreal8(area_upd_prc1)
457  area_upd_prc1 = 0.0
458  CALL pushreal8(area_anv_prc1)
459  area_anv_prc1 = 0.0
460  CALL pushreal8(tot_prec_upd)
461  tot_prec_upd = tot_prec_upd + (qrn_cu_1d+qsn_cu)*mass(i, j, k)
462  CALL pushreal8(area_upd_prc)
463  area_upd_prc = area_upd_prc + cnv_updf(i, j, k)*(qrn_cu_1d+&
464 & qsn_cu)*mass(i, j, k)
465  CALL pushreal8(tot_prec_anv)
466  tot_prec_anv = tot_prec_anv + (qrn_an+qsn_an)*mass(i, j, k)
467  CALL pushreal8(area_anv_prc)
468  area_anv_prc = area_anv_prc + cf_con(i, j, k)*(qrn_an+qsn_an)*&
469 & mass(i, j, k)
470  CALL pushreal8(tot_prec_ls)
471  tot_prec_ls = tot_prec_ls + (qrn_ls+qsn_ls)*mass(i, j, k)
472  CALL pushreal8(area_ls_prc)
473  area_ls_prc = area_ls_prc + cf_ls(i, j, k)*(qrn_ls+qsn_ls)*mass(&
474 & i, j, k)
475  IF (tot_prec_anv .GT. 0.0) THEN
476  IF (area_anv_prc/tot_prec_anv .LT. 1.e-6) THEN
477  CALL pushcontrol2b(2)
478  area_anv_prc1 = 1.e-6
479  ELSE
480  area_anv_prc1 = area_anv_prc/tot_prec_anv
481  CALL pushcontrol2b(1)
482  END IF
483  ELSE
484  CALL pushcontrol2b(0)
485  END IF
486  IF (tot_prec_upd .GT. 0.0) THEN
487  IF (area_upd_prc/tot_prec_upd .LT. 1.e-6) THEN
488  CALL pushcontrol2b(2)
489  area_upd_prc1 = 1.e-6
490  ELSE
491  area_upd_prc1 = area_upd_prc/tot_prec_upd
492  CALL pushcontrol2b(1)
493  END IF
494  ELSE
495  CALL pushcontrol2b(0)
496  END IF
497  IF (tot_prec_ls .GT. 0.0) THEN
498  IF (area_ls_prc/tot_prec_ls .LT. 1.e-6) THEN
499  CALL pushcontrol2b(2)
500  area_ls_prc1 = 1.e-6
501  ELSE
502  area_ls_prc1 = area_ls_prc/tot_prec_ls
503  CALL pushcontrol2b(1)
504  END IF
505  ELSE
506  CALL pushcontrol2b(0)
507  END IF
508  area_ls_prc1 = ls_beta*area_ls_prc1
509  area_upd_prc1 = cnv_beta*area_upd_prc1
510  area_anv_prc1 = anv_beta*area_anv_prc1
511  IF (k .EQ. lm) THEN
512 ! We have accumulated over the whole column
513  IF (tot_prec_anv .GT. 0.0) THEN
514  IF (area_anv_prc/tot_prec_anv .LT. 1.e-6) THEN
515  CALL pushreal8(area_anv_prc)
516  area_anv_prc = 1.e-6
517  CALL pushcontrol2b(2)
518  ELSE
519  CALL pushreal8(area_anv_prc)
520  area_anv_prc = area_anv_prc/tot_prec_anv
521  CALL pushcontrol2b(1)
522  END IF
523  ELSE
524  CALL pushcontrol2b(0)
525  END IF
526  IF (tot_prec_upd .GT. 0.0) THEN
527  IF (area_upd_prc/tot_prec_upd .LT. 1.e-6) THEN
528  CALL pushreal8(area_upd_prc)
529  area_upd_prc = 1.e-6
530  CALL pushcontrol2b(2)
531  ELSE
532  CALL pushreal8(area_upd_prc)
533  area_upd_prc = area_upd_prc/tot_prec_upd
534  CALL pushcontrol2b(1)
535  END IF
536  ELSE
537  CALL pushcontrol2b(0)
538  END IF
539  IF (tot_prec_ls .GT. 0.0) THEN
540  IF (area_ls_prc/tot_prec_ls .LT. 1.e-6) THEN
541  CALL pushreal8(area_ls_prc)
542  area_ls_prc = 1.e-6
543  CALL pushcontrol2b(2)
544  ELSE
545  CALL pushreal8(area_ls_prc)
546  area_ls_prc = area_ls_prc/tot_prec_ls
547  CALL pushcontrol2b(1)
548  END IF
549  ELSE
550  CALL pushcontrol2b(0)
551  END IF
552  CALL pushreal8(area_ls_prc)
553  area_ls_prc = ls_beta*area_ls_prc
554  CALL pushreal8(area_upd_prc)
555  area_upd_prc = cnv_beta*area_upd_prc
556  CALL pushreal8(area_anv_prc)
557  area_anv_prc = anv_beta*area_anv_prc
558  CALL pushcontrol1b(0)
559  ELSE
560  CALL pushcontrol1b(1)
561  END IF
562 !Get micro-physical constants
563  CALL pushreal8(alhx3)
564  CALL cons_alhx(t(i, j, k), alhx3, t_ice_max, t_ice_all, &
565 & cons_alhs, cons_alhl)
566  CALL pushreal8(bb)
567  CALL pushreal8(aa)
568  CALL cons_microphys(t(i, j, k), ph(i, j, k), qs(i, j, k), aa, bb&
569 & , cons_h2omw, cons_airmw, cons_rvap, alhx3)
570 !Precip Scheme Expects Total Cloud Liquid
571  CALL pushreal8(qlt_tmp)
572  qlt_tmp = ql_ls(i, j, k) + ql_con(i, j, k)
573  CALL pushreal8(qit_tmp)
574  qit_tmp = qi_ls(i, j, k) + qi_con(i, j, k)
575  CALL pushreal8(prn_above_cu_old)
576  prn_above_cu_old = prn_above_cu_new
577  CALL pushreal8(psn_above_cu_old)
578  psn_above_cu_old = psn_above_cu_new
579  CALL pushreal8(evap_dd_cu_above_old)
580  evap_dd_cu_above_old = evap_dd_cu_above_new
581  CALL pushreal8(subl_dd_cu_above_old)
582  subl_dd_cu_above_old = subl_dd_cu_above_new
583 !Precip and Evap for Convection
584  CALL pushreal8(q(i, j, k))
585  CALL pushreal8(t(i, j, k))
586  CALL pushreal8(qlt_tmp)
587  CALL pushreal8(qsn_cu)
588  CALL pushreal8(qrn_cu_1d)
589  CALL precipandevap(k, ktop, lm, dt, frland(i, j), rhcrit, &
590 & qrn_cu_1d, qsn_cu, qlt_tmp, qit_tmp, t(i, j, k), q(&
591 & i, j, k), mass(i, j, k), imass(i, j, k), ph(i, j, k&
592 & ), dzet(i, j, k), qddf3(i, j, k), aa, bb, &
593 & area_upd_prc1, prn_above_cu_old, prn_above_cu_new, &
594 & psn_above_cu_old, psn_above_cu_new, &
595 & evap_dd_cu_above_old, evap_dd_cu_above_new, &
596 & subl_dd_cu_above_old, subl_dd_cu_above_new, &
597 & cnvenvfc, cnvddrfc, cons_alhf, cons_alhs, cons_alhl&
598 & , cons_cp, cons_tice, cons_h2omw, cons_airmw, &
599 & revap_off_p, c_acc, c_ev_r, c_ev_s, rho_w, estblx)
600  CALL pushreal8(prn_above_an_old)
601  prn_above_an_old = prn_above_an_new
602  CALL pushreal8(psn_above_an_old)
603  psn_above_an_old = psn_above_an_new
604  CALL pushreal8(evap_dd_an_above_old)
605  evap_dd_an_above_old = evap_dd_an_above_new
606  CALL pushreal8(subl_dd_an_above_old)
607  subl_dd_an_above_old = subl_dd_an_above_new
608 !Precip and Evap for Anvil
609  anvenvfc = 1.0
610  CALL pushreal8(q(i, j, k))
611  CALL pushreal8(t(i, j, k))
612  CALL pushreal8(qlt_tmp)
613  CALL pushreal8(qsn_an)
614  CALL pushreal8(qrn_an)
615  CALL precipandevap(k, ktop, lm, dt, frland(i, j), rhcrit, qrn_an&
616 & , qsn_an, qlt_tmp, qit_tmp, t(i, j, k), q(i, j, k)&
617 & , mass(i, j, k), imass(i, j, k), ph(i, j, k), dzet(&
618 & i, j, k), qddf3(i, j, k), aa, bb, area_anv_prc1, &
619 & prn_above_an_old, prn_above_an_new, &
620 & psn_above_an_old, psn_above_an_new, &
621 & evap_dd_an_above_old, evap_dd_an_above_new, &
622 & subl_dd_an_above_old, subl_dd_an_above_new, &
623 & anvenvfc, anvddrfc, cons_alhf, cons_alhs, cons_alhl&
624 & , cons_cp, cons_tice, cons_h2omw, cons_airmw, &
625 & revap_off_p, c_acc, c_ev_r, c_ev_s, rho_w, estblx)
626  CALL pushreal8(prn_above_ls_old)
627  prn_above_ls_old = prn_above_ls_new
628  CALL pushreal8(psn_above_ls_old)
629  psn_above_ls_old = psn_above_ls_new
630  CALL pushreal8(evap_dd_ls_above_old)
631  evap_dd_ls_above_old = evap_dd_ls_above_new
632  CALL pushreal8(subl_dd_ls_above_old)
633  subl_dd_ls_above_old = subl_dd_ls_above_new
634 !Precip and Evap for Large Scale
635  lsenvfc = 1.0
636  CALL pushreal8(q(i, j, k))
637  CALL pushreal8(t(i, j, k))
638  CALL pushreal8(qlt_tmp)
639  CALL pushreal8(qsn_ls)
640  CALL pushreal8(qrn_ls)
641  CALL precipandevap(k, ktop, lm, dt, frland(i, j), rhcrit, qrn_ls&
642 & , qsn_ls, qlt_tmp, qit_tmp, t(i, j, k), q(i, j, k)&
643 & , mass(i, j, k), imass(i, j, k), ph(i, j, k), dzet(&
644 & i, j, k), qddf3(i, j, k), aa, bb, area_ls_prc1, &
645 & prn_above_ls_old, prn_above_ls_new, &
646 & psn_above_ls_old, psn_above_ls_new, &
647 & evap_dd_ls_above_old, evap_dd_ls_above_new, &
648 & subl_dd_ls_above_old, subl_dd_ls_above_new, lsenvfc&
649 & , lsddrfc, cons_alhf, cons_alhs, cons_alhl, cons_cp&
650 & , cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
651 & c_acc, c_ev_r, c_ev_s, rho_w, estblx)
652  IF (ql_ls(i, j, k) + ql_con(i, j, k) .GT. 0.00) THEN
653  CALL pushreal8(qt_tmpi_1)
654  qt_tmpi_1 = 1./(ql_ls(i, j, k)+ql_con(i, j, k))
655  CALL pushcontrol1b(0)
656  ELSE
657  CALL pushreal8(qt_tmpi_1)
658  qt_tmpi_1 = 0.0
659  CALL pushcontrol1b(1)
660  END IF
661  CALL pushreal8(ql_ls(i, j, k))
662  ql_ls(i, j, k) = ql_ls(i, j, k)*qlt_tmp*qt_tmpi_1
663  CALL pushreal8(ql_con(i, j, k))
664  ql_con(i, j, k) = ql_con(i, j, k)*qlt_tmp*qt_tmpi_1
665  IF (qi_ls(i, j, k) + qi_con(i, j, k) .GT. 0.00) THEN
666  CALL pushreal8(qt_tmpi_2)
667  qt_tmpi_2 = 1./(qi_ls(i, j, k)+qi_con(i, j, k))
668  CALL pushcontrol1b(0)
669  ELSE
670  CALL pushreal8(qt_tmpi_2)
671  qt_tmpi_2 = 0.0
672  CALL pushcontrol1b(1)
673  END IF
674  CALL pushreal8(qi_ls(i, j, k))
675  qi_ls(i, j, k) = qi_ls(i, j, k)*qit_tmp*qt_tmpi_2
676  CALL pushreal8(qi_con(i, j, k))
677  qi_con(i, j, k) = qi_con(i, j, k)*qit_tmp*qt_tmpi_2
678  END DO
679  END DO
680  END DO
681 !Clean up of excess relative humidity
682  rhexcess = 1.1
683  CALL pushreal8array(qs, im*jm*lm)
684  CALL dqsat_bac(dqsdt, qs, t, ph, im, jm, lm, estblx, cons_h2omw, &
685 & cons_airmw)
686  mask0 = q .GT. rhexcess*qs
687  WHERE (mask0)
688  dqs = (q-rhexcess*qs)/(1.0+rhexcess*dqsdt*cons_alhl/cons_cp)
689  ELSEWHERE
690  dqs = 0.0
691  END WHERE
692  CALL pushreal8array(q, im*jm*lm)
693  q = q - dqs
694 !Clean up Q<0
695  DO j=1,jm
696  DO i=1,im
697 !Total precipitable water
698  CALL pushreal8(tpw)
699  tpw = sum(q(i, j, :)*dm(i, j, :))
700  CALL pushreal8(negtpw)
701  negtpw = 0.
702  DO l=1,lm
703  IF (q(i, j, l) .LT. 0.0) THEN
704  negtpw = negtpw + q(i, j, l)*dm(i, j, l)
705  q(i, j, l) = 0.0
706  CALL pushcontrol1b(1)
707  ELSE
708  CALL pushcontrol1b(0)
709  END IF
710  END DO
711  DO l=1,lm
712  IF (q(i, j, l) .GE. 0.0) THEN
713  CALL pushcontrol1b(1)
714  ELSE
715  CALL pushcontrol1b(0)
716  END IF
717  END DO
718  END DO
719  END DO
720  tb = 0.0_8
721  tb = thb/pih
722  DO j=jm,1,-1
723  DO i=im,1,-1
724  tpwb = 0.0_8
725  negtpwb = 0.0_8
726  DO l=lm,1,-1
727  CALL popcontrol1b(branch)
728  IF (branch .NE. 0) THEN
729  temp2 = negtpw/(tpw-negtpw)
730  tempb15 = q(i, j, l)*qb(i, j, l)/(tpw-negtpw)
731  tempb16 = -(temp2*tempb15)
732  negtpwb = negtpwb + tempb15 - tempb16
733  tpwb = tpwb + tempb16
734  qb(i, j, l) = (temp2+1.0)*qb(i, j, l)
735  END IF
736  END DO
737  DO l=lm,1,-1
738  CALL popcontrol1b(branch)
739  IF (branch .NE. 0) qb(i, j, l) = dm(i, j, l)*negtpwb
740  END DO
741  CALL popreal8(negtpw)
742  CALL popreal8(tpw)
743  qb(i, j, :) = qb(i, j, :) + dm(i, j, :)*tpwb
744  END DO
745  END DO
746  dqsb = 0.0_8
747  dqsb = cons_alhl*tb/cons_cp - qb
748  CALL popreal8array(q, im*jm*lm)
749  dqsdtb = 0.0_8
750  qsb = 0.0_8
751  temp1 = rhexcess*cons_alhl*dqsdt/cons_cp + 1.0
752  WHERE (.NOT.mask0)
753  dqsb = 0.0_8
754  ELSEWHERE
755  tempb14 = dqsb/temp1
756  qb = qb + tempb14
757  qsb = -(rhexcess*tempb14)
758  dqsdtb = -(rhexcess*cons_alhl*(q-rhexcess*qs)*tempb14/(cons_cp*temp1&
759 & ))
760  END WHERE
761  CALL popreal8array(qs, im*jm*lm)
762  CALL dqsat_bac_b(dqsdt, dqsdtb, qs, qsb, t, tb, ph, im, jm, lm, estblx&
763 & , cons_h2omw, cons_airmw)
764  cnv_prc3b = 0.0_8
765  cnv_dqldtb = 0.0_8
766  cnv_updfb = 0.0_8
767  cnv_mfdb = 0.0_8
768  prn_above_cu_newb = 0.0_8
769  psn_above_cu_newb = 0.0_8
770  evap_dd_an_above_newb = 0.0_8
771  area_upd_prcb = 0.0_8
772  prn_above_ls_newb = 0.0_8
773  evap_dd_ls_above_newb = 0.0_8
774  evap_dd_cu_above_newb = 0.0_8
775  psn_above_ls_newb = 0.0_8
776  tot_prec_lsb = 0.0_8
777  area_ls_prcb = 0.0_8
778  tot_prec_updb = 0.0_8
779  prn_above_an_newb = 0.0_8
780  area_anv_prcb = 0.0_8
781  psn_above_an_newb = 0.0_8
782  alhx3b = 0.0_8
783  subl_dd_an_above_newb = 0.0_8
784  qddf3b = 0.0_8
785  tot_prec_anvb = 0.0_8
786  dzetb = 0.0_8
787  subl_dd_ls_above_newb = 0.0_8
788  subl_dd_cu_above_newb = 0.0_8
789  DO i=im,1,-1
790  DO j=jm,1,-1
791  DO k=lm,ktop,-1
792 
793  !TOTAL FILTERING
794  totfilt_t = 0.25
795 
796  t_p_preall = (1.0-totfilt_t)*tb(i,j,k)
797  tb(i,j,k) = totfilt_t *tb(i,j,k)
798 
799  if (do_moist_physics == 1) then
800  totfilt_ql = 0.75
801  totfilt_qi = 1.0
802  elseif (do_moist_physics == 2) then
803  totfilt_ql = 0.5
804  totfilt_qi = 1.0!0.5
805  endif
806 
807  ql_ls_p_preall = (1.0-totfilt_ql)*ql_lsb(i,j,k)
808  ql_lsb(i,j,k) = totfilt_ql * ql_lsb(i,j,k)
809  ql_con_p_preall = (1.0-totfilt_ql)*ql_conb(i,j,k)
810  ql_conb(i,j,k) = totfilt_ql * ql_conb(i,j,k)
811 
812  qi_ls_p_preall = (1.0-totfilt_qi) * qi_lsb(i,j,k)
813  qi_lsb(i,j,k) = totfilt_qi * qi_lsb(i,j,k)
814  qi_con_p_preall = (1.0-totfilt_qi) * qi_conb(i,j,k)
815  qi_conb(i,j,k) = totfilt_qi * qi_conb(i,j,k)
816 
817  !SINK FILTERING
818  if (do_moist_physics == 1) then
819  sinkfilt_ql = 0.65
820  sinkfilt_qi = 0.65
821  sinkfilt_cf = 1.0
822  elseif (do_moist_physics == 2) then
823  sinkfilt_ql = 0.9
824  sinkfilt_qi = 0.9
825  sinkfilt_cf = 1.0
826  endif
827 
828  qi_ls_p_presink = 0.0
829  qi_con_p_presink = 0.0
830  q_p_presink = 0.0
831  ql_ls_p_presink = 0.0
832  ql_con_p_presink = 0.0
833 
834  if (k < 50) then
835  qi_ls_p_presink = (1.0-sinkfilt_qi) * qi_lsb(i,j,k)
836  qi_lsb(i,j,k) = sinkfilt_qi * qi_lsb(i,j,k)
837  qi_con_p_presink = (1.0-sinkfilt_qi) * qi_conb(i,j,k)
838  qi_conb(i,j,k) = sinkfilt_qi * qi_conb(i,j,k)
839  q_p_presink = (1.0-sinkfilt_qi) * qb(i,j,k)
840  qb(i,j,k) = sinkfilt_qi * qb(i,j,k)
841  endif
842 
843  if ( abs(k-62) .le. 2) then
844  ql_ls_p_presink = (1.0-sinkfilt_ql) * ql_lsb(i,j,k)
845  ql_lsb(i,j,k) = sinkfilt_ql * ql_lsb(i,j,k)
846  ql_con_p_presink = (1.0-sinkfilt_ql) * ql_conb(i,j,k)
847  ql_conb(i,j,k) = sinkfilt_ql * ql_conb(i,j,k)
848  endif
849 
850  !cf_con_p_presink = (1.0-SINKfilt_CF) * cf_conb(i,j,k)
851  !cf_conb(i,j,k) = SINKfilt_CF * cf_conb(i,j,k)
852 
853  cf_con_p_presink = 0.0
854  cf_conb(i,j,k) = 0.0
855 
856  CALL popreal8(qi_con(i, j, k))
857  tempb12 = qi_con(i, j, k)*qi_conb(i, j, k)
858  qi_conb(i, j, k) = qit_tmp*qt_tmpi_2*qi_conb(i, j, k)
859  CALL popreal8(qi_ls(i, j, k))
860  tempb13 = qi_ls(i, j, k)*qi_lsb(i, j, k)
861  qit_tmpb = qt_tmpi_2*tempb13 + qt_tmpi_2*tempb12
862  qt_tmpi_2b = qit_tmp*tempb13 + qit_tmp*tempb12
863  qi_lsb(i, j, k) = qit_tmp*qt_tmpi_2*qi_lsb(i, j, k)
864  CALL popcontrol1b(branch)
865  IF (branch .EQ. 0) THEN
866  CALL popreal8(qt_tmpi_2)
867  temp0 = qi_ls(i, j, k) + qi_con(i, j, k)
868  tempb11 = -(qt_tmpi_2b/temp0**2)
869  qi_lsb(i, j, k) = qi_lsb(i, j, k) + tempb11
870  qi_conb(i, j, k) = qi_conb(i, j, k) + tempb11
871  ELSE
872  CALL popreal8(qt_tmpi_2)
873  END IF
874  CALL popreal8(ql_con(i, j, k))
875  tempb9 = ql_con(i, j, k)*ql_conb(i, j, k)
876  ql_conb(i, j, k) = qlt_tmp*qt_tmpi_1*ql_conb(i, j, k)
877  CALL popreal8(ql_ls(i, j, k))
878  tempb10 = ql_ls(i, j, k)*ql_lsb(i, j, k)
879  qlt_tmpb = qt_tmpi_1*tempb10 + qt_tmpi_1*tempb9
880  qt_tmpi_1b = qlt_tmp*tempb10 + qlt_tmp*tempb9
881  ql_lsb(i, j, k) = qlt_tmp*qt_tmpi_1*ql_lsb(i, j, k)
882  CALL popcontrol1b(branch)
883  IF (branch .EQ. 0) THEN
884  CALL popreal8(qt_tmpi_1)
885  temp = ql_ls(i, j, k) + ql_con(i, j, k)
886  tempb8 = -(qt_tmpi_1b/temp**2)
887  ql_lsb(i, j, k) = ql_lsb(i, j, k) + tempb8
888  ql_conb(i, j, k) = ql_conb(i, j, k) + tempb8
889  ELSE
890  CALL popreal8(qt_tmpi_1)
891  END IF
892  lsenvfc = 1.0
893  CALL popreal8(qrn_ls)
894  CALL popreal8(qsn_ls)
895  CALL popreal8(qlt_tmp)
896  CALL popreal8(t(i, j, k))
897  CALL popreal8(q(i, j, k))
898  aab = 0.0_8
899  bbb = 0.0_8
900  CALL precipandevap_b(k, ktop, lm, dt, frland(i, j), rhcrit, &
901 & qrn_ls, qrn_lsb, qsn_ls, qsn_lsb, qlt_tmp, &
902 & qlt_tmpb, qit_tmp, t(i, j, k), tb(i, j, k), q(i, &
903 & j, k), qb(i, j, k), mass(i, j, k), imass(i, j, k)&
904 & , ph(i, j, k), dzet(i, j, k), dzetb(i, j, k), &
905 & qddf3(i, j, k), qddf3b(i, j, k), aa, aab, bb, bbb&
906 & , area_ls_prc1, area_ls_prc1b, prn_above_ls_old, &
907 & prn_above_ls_oldb, prn_above_ls_new, &
908 & prn_above_ls_newb, psn_above_ls_old, &
909 & psn_above_ls_oldb, psn_above_ls_new, &
910 & psn_above_ls_newb, evap_dd_ls_above_old, &
911 & evap_dd_ls_above_oldb, evap_dd_ls_above_new, &
912 & evap_dd_ls_above_newb, subl_dd_ls_above_old, &
913 & subl_dd_ls_above_oldb, subl_dd_ls_above_new, &
914 & subl_dd_ls_above_newb, lsenvfc, lsddrfc, &
915 & cons_alhf, cons_alhs, cons_alhl, cons_cp, &
916 & cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
917 & c_acc, c_ev_r, c_ev_s, rho_w, estblx)
918  CALL popreal8(subl_dd_ls_above_old)
919  subl_dd_ls_above_newb = subl_dd_ls_above_oldb
920  CALL popreal8(evap_dd_ls_above_old)
921  evap_dd_ls_above_newb = evap_dd_ls_above_oldb
922  CALL popreal8(psn_above_ls_old)
923  psn_above_ls_newb = psn_above_ls_oldb
924  CALL popreal8(prn_above_ls_old)
925  prn_above_ls_newb = prn_above_ls_oldb
926  anvenvfc = 1.0
927  CALL popreal8(qrn_an)
928  CALL popreal8(qsn_an)
929  CALL popreal8(qlt_tmp)
930  CALL popreal8(t(i, j, k))
931  CALL popreal8(q(i, j, k))
932  CALL precipandevap_b(k, ktop, lm, dt, frland(i, j), rhcrit, &
933 & qrn_an, qrn_anb, qsn_an, qsn_anb, qlt_tmp, &
934 & qlt_tmpb, qit_tmp, t(i, j, k), tb(i, j, k), q(i, &
935 & j, k), qb(i, j, k), mass(i, j, k), imass(i, j, k)&
936 & , ph(i, j, k), dzet(i, j, k), dzetb(i, j, k), &
937 & qddf3(i, j, k), qddf3b(i, j, k), aa, aab, bb, bbb&
938 & , area_anv_prc1, area_anv_prc1b, prn_above_an_old&
939 & , prn_above_an_oldb, prn_above_an_new, &
940 & prn_above_an_newb, psn_above_an_old, &
941 & psn_above_an_oldb, psn_above_an_new, &
942 & psn_above_an_newb, evap_dd_an_above_old, &
943 & evap_dd_an_above_oldb, evap_dd_an_above_new, &
944 & evap_dd_an_above_newb, subl_dd_an_above_old, &
945 & subl_dd_an_above_oldb, subl_dd_an_above_new, &
946 & subl_dd_an_above_newb, anvenvfc, anvddrfc, &
947 & cons_alhf, cons_alhs, cons_alhl, cons_cp, &
948 & cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
949 & c_acc, c_ev_r, c_ev_s, rho_w, estblx)
950  CALL popreal8(subl_dd_an_above_old)
951  subl_dd_an_above_newb = subl_dd_an_above_oldb
952  CALL popreal8(evap_dd_an_above_old)
953  evap_dd_an_above_newb = evap_dd_an_above_oldb
954  CALL popreal8(psn_above_an_old)
955  psn_above_an_newb = psn_above_an_oldb
956  CALL popreal8(prn_above_an_old)
957  prn_above_an_newb = prn_above_an_oldb
958  CALL popreal8(qrn_cu_1d)
959  CALL popreal8(qsn_cu)
960  CALL popreal8(qlt_tmp)
961  CALL popreal8(t(i, j, k))
962  CALL popreal8(q(i, j, k))
963  CALL precipandevap_b(k, ktop, lm, dt, frland(i, j), rhcrit, &
964 & qrn_cu_1d, qrn_cu_1db, qsn_cu, qsn_cub, qlt_tmp, &
965 & qlt_tmpb, qit_tmp, t(i, j, k), tb(i, j, k), q(i, &
966 & j, k), qb(i, j, k), mass(i, j, k), imass(i, j, k)&
967 & , ph(i, j, k), dzet(i, j, k), dzetb(i, j, k), &
968 & qddf3(i, j, k), qddf3b(i, j, k), aa, aab, bb, bbb&
969 & , area_upd_prc1, area_upd_prc1b, prn_above_cu_old&
970 & , prn_above_cu_oldb, prn_above_cu_new, &
971 & prn_above_cu_newb, psn_above_cu_old, &
972 & psn_above_cu_oldb, psn_above_cu_new, &
973 & psn_above_cu_newb, evap_dd_cu_above_old, &
974 & evap_dd_cu_above_oldb, evap_dd_cu_above_new, &
975 & evap_dd_cu_above_newb, subl_dd_cu_above_old, &
976 & subl_dd_cu_above_oldb, subl_dd_cu_above_new, &
977 & subl_dd_cu_above_newb, cnvenvfc, cnvddrfc, &
978 & cons_alhf, cons_alhs, cons_alhl, cons_cp, &
979 & cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
980 & c_acc, c_ev_r, c_ev_s, rho_w, estblx)
981  CALL popreal8(subl_dd_cu_above_old)
982  subl_dd_cu_above_newb = subl_dd_cu_above_oldb
983  CALL popreal8(evap_dd_cu_above_old)
984  evap_dd_cu_above_newb = evap_dd_cu_above_oldb
985  CALL popreal8(psn_above_cu_old)
986  psn_above_cu_newb = psn_above_cu_oldb
987  CALL popreal8(prn_above_cu_old)
988  prn_above_cu_newb = prn_above_cu_oldb
989  CALL popreal8(qit_tmp)
990  qi_lsb(i, j, k) = qi_lsb(i, j, k) + qit_tmpb
991  qi_conb(i, j, k) = qi_conb(i, j, k) + qit_tmpb
992  CALL popreal8(qlt_tmp)
993  ql_lsb(i, j, k) = ql_lsb(i, j, k) + qlt_tmpb
994  ql_conb(i, j, k) = ql_conb(i, j, k) + qlt_tmpb
995  CALL popreal8(aa)
996  CALL popreal8(bb)
997  CALL cons_microphys_b(t(i, j, k), tb(i, j, k), ph(i, j, k), qs(i&
998 & , j, k), qsb(i, j, k), aa, aab, bb, bbb, &
999 & cons_h2omw, cons_airmw, cons_rvap, alhx3, alhx3b&
1000 & )
1001  CALL popreal8(alhx3)
1002  CALL cons_alhx_b(t(i, j, k), tb(i, j, k), alhx3, alhx3b, &
1003 & t_ice_max, t_ice_all, cons_alhs, cons_alhl)
1004  CALL popcontrol1b(branch)
1005  IF (branch .EQ. 0) THEN
1006  CALL popreal8(area_anv_prc)
1007  area_anv_prcb = anv_beta*area_anv_prcb
1008  CALL popreal8(area_upd_prc)
1009  area_upd_prcb = cnv_beta*area_upd_prcb
1010  CALL popreal8(area_ls_prc)
1011  area_ls_prcb = ls_beta*area_ls_prcb
1012  CALL popcontrol2b(branch)
1013  IF (branch .NE. 0) THEN
1014  IF (branch .EQ. 1) THEN
1015  CALL popreal8(area_ls_prc)
1016  tot_prec_lsb = tot_prec_lsb - area_ls_prc*area_ls_prcb/&
1017 & tot_prec_ls**2
1018  area_ls_prcb = area_ls_prcb/tot_prec_ls
1019  ELSE
1020  CALL popreal8(area_ls_prc)
1021  area_ls_prcb = 0.0_8
1022  END IF
1023  END IF
1024  CALL popcontrol2b(branch)
1025  IF (branch .NE. 0) THEN
1026  IF (branch .EQ. 1) THEN
1027  CALL popreal8(area_upd_prc)
1028  tot_prec_updb = tot_prec_updb - area_upd_prc*area_upd_prcb&
1029 & /tot_prec_upd**2
1030  area_upd_prcb = area_upd_prcb/tot_prec_upd
1031  ELSE
1032  CALL popreal8(area_upd_prc)
1033  area_upd_prcb = 0.0_8
1034  END IF
1035  END IF
1036  CALL popcontrol2b(branch)
1037  IF (branch .NE. 0) THEN
1038  IF (branch .EQ. 1) THEN
1039  CALL popreal8(area_anv_prc)
1040  tot_prec_anvb = tot_prec_anvb - area_anv_prc*area_anv_prcb&
1041 & /tot_prec_anv**2
1042  area_anv_prcb = area_anv_prcb/tot_prec_anv
1043  ELSE
1044  CALL popreal8(area_anv_prc)
1045  area_anv_prcb = 0.0_8
1046  END IF
1047  END IF
1048  END IF
1049  area_anv_prc1b = anv_beta*area_anv_prc1b
1050  area_upd_prc1b = cnv_beta*area_upd_prc1b
1051  area_ls_prc1b = ls_beta*area_ls_prc1b
1052  CALL popcontrol2b(branch)
1053  IF (branch .NE. 0) THEN
1054  IF (branch .EQ. 1) THEN
1055  area_ls_prcb = area_ls_prcb + area_ls_prc1b/tot_prec_ls
1056  tot_prec_lsb = tot_prec_lsb - area_ls_prc*area_ls_prc1b/&
1057 & tot_prec_ls**2
1058  END IF
1059  END IF
1060  CALL popcontrol2b(branch)
1061  IF (branch .NE. 0) THEN
1062  IF (branch .EQ. 1) THEN
1063  area_upd_prcb = area_upd_prcb + area_upd_prc1b/tot_prec_upd
1064  tot_prec_updb = tot_prec_updb - area_upd_prc*area_upd_prc1b/&
1065 & tot_prec_upd**2
1066  END IF
1067  END IF
1068  CALL popcontrol2b(branch)
1069  IF (branch .NE. 0) THEN
1070  IF (branch .EQ. 1) THEN
1071  area_anv_prcb = area_anv_prcb + area_anv_prc1b/tot_prec_anv
1072  tot_prec_anvb = tot_prec_anvb - area_anv_prc*area_anv_prc1b/&
1073 & tot_prec_anv**2
1074  END IF
1075  END IF
1076  tempb7 = mass(i, j, k)*tot_prec_updb
1077  tempb4 = mass(i, j, k)*tot_prec_anvb
1078  tempb1 = mass(i, j, k)*tot_prec_lsb
1079  CALL popreal8(area_ls_prc)
1080  tempb = mass(i, j, k)*area_ls_prcb
1081  tempb0 = cf_ls(i, j, k)*tempb
1082  cf_lsb(i, j, k) = cf_lsb(i, j, k) + (qrn_ls+qsn_ls)*tempb
1083  qrn_lsb = qrn_lsb + tempb1 + tempb0
1084  qsn_lsb = qsn_lsb + tempb1 + tempb0
1085  CALL popreal8(tot_prec_ls)
1086  CALL popreal8(area_anv_prc)
1087  tempb2 = mass(i, j, k)*area_anv_prcb
1088  tempb3 = cf_con(i, j, k)*tempb2
1089  cf_conb(i, j, k) = cf_conb(i, j, k) + (qrn_an+qsn_an)*tempb2
1090  qrn_anb = qrn_anb + tempb4 + tempb3
1091  qsn_anb = qsn_anb + tempb4 + tempb3
1092  CALL popreal8(tot_prec_anv)
1093  CALL popreal8(area_upd_prc)
1094  tempb5 = mass(i, j, k)*area_upd_prcb
1095  tempb6 = cnv_updf(i, j, k)*tempb5
1096  cnv_updfb(i, j, k) = cnv_updfb(i, j, k) + (qrn_cu_1d+qsn_cu)*&
1097 & tempb5
1098  qrn_cu_1db = qrn_cu_1db + tempb7 + tempb6
1099  qsn_cub = qsn_cub + tempb7 + tempb6
1100  CALL popreal8(tot_prec_upd)
1101  CALL popreal8(area_anv_prc1)
1102  CALL popreal8(area_upd_prc1)
1103  CALL popreal8(area_ls_prc1)
1104  CALL popcontrol1b(branch)
1105  IF (branch .EQ. 0) THEN
1106  CALL popreal8(t(i, j, k))
1107  qsn_cub = qsn_cub + (cons_alhs-cons_alhl)*tb(i, j, k)/cons_cp
1108  qrn_cu_1db = qsn_cub
1109  END IF
1110  CALL popreal8(qi_ls(i, j, k))
1111  CALL popreal8(cf_ls(i, j, k))
1112  CALL ice_settlefall_ls_b(wrhodep, qi_ls(i, j, k), qi_lsb(i, j, k&
1113 & ), ph(i, j, k), t(i, j, k), tb(i, j, k), &
1114 & cf_ls(i, j, k), cf_lsb(i, j, k), cons_rgas, &
1115 & khu(i, j), khl(i, j), k, dt, dzet(i, j, k), &
1116 & dzetb(i, j, k), qsn_ls, qsn_lsb, ls_icefall_c&
1117 & )
1118  CALL popreal8(qi_con(i, j, k))
1119  CALL ice_settlefall_cnv_b(wrhodep, qi_con(i, j, k), qi_conb(i, j&
1120 & , k), ph(i, j, k), t(i, j, k), tb(i, j, k), &
1121 & cf_con(i, j, k), cf_conb(i, j, k), cons_rgas&
1122 & , khu(i, j), khl(i, j), k, dt, dzet(i, j, k)&
1123 & , dzetb(i, j, k), qsn_an, qsn_anb, &
1124 & anv_icefall_c)
1125  CALL popreal8(ql_con(i, j, k))
1126  CALL autoconversion_cnv_b(dt, ql_con(i, j, k), ql_conb(i, j, k)&
1127 & , qrn_an, qrn_anb, t(i, j, k), tb(i, j, k), &
1128 & ph(i, j, k), cf_con(i, j, k), cf_conb(i, j, &
1129 & k), anv_sdqv2, anv_sdqv3, anv_sdqvt1, c_00, &
1130 & lwcrit, dzet(i, j, k))
1131  CALL popreal8(ql_ls(i, j, k))
1132  CALL popreal8(cf_ls(i, j, k))
1133  CALL autoconversion_ls_b(dt, ql_ls(i, j, k), ql_lsb(i, j, k), &
1134 & qrn_ls, qrn_lsb, t(i, j, k), tb(i, j, k), ph(&
1135 & i, j, k), cf_ls(i, j, k), cf_lsb(i, j, k), &
1136 & ls_sdqv2, ls_sdqv3, ls_sdqvt1, c_00, lwcrit, &
1137 & dzet(i, j, k))
1138  CALL popreal8(t(i, j, k))
1139  CALL popreal8(q(i, j, k))
1140  CALL popreal8(qi_con(i, j, k))
1141  CALL popreal8(cf_con(i, j, k))
1142  CALL subl_cnv_b(dt, rhcrit, ph(i, j, k), t(i, j, k), tb(i, j, k)&
1143 & , q(i, j, k), qb(i, j, k), ql_con(i, j, k), ql_conb(i&
1144 & , j, k), qi_con(i, j, k), qi_conb(i, j, k), cf_con(i, &
1145 & j, k), cf_conb(i, j, k), cf_ls(i, j, k), qs(i, j, k), &
1146 & qsb(i, j, k), rho_w, cld_evp_eff, cons_h2omw, &
1147 & cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, &
1148 & cons_cp, cons_alhs)
1149  CALL popreal8(t(i, j, k))
1150  CALL popreal8(q(i, j, k))
1151  CALL popreal8(ql_con(i, j, k))
1152  CALL popreal8(cf_con(i, j, k))
1153  CALL evap_cnv_b(dt, rhcrit, ph(i, j, k), t(i, j, k), tb(i, j, k)&
1154 & , q(i, j, k), qb(i, j, k), ql_con(i, j, k), ql_conb(i&
1155 & , j, k), qi_con(i, j, k), qi_conb(i, j, k), cf_con(i, &
1156 & j, k), cf_conb(i, j, k), cf_ls(i, j, k), qs(i, j, k), &
1157 & qsb(i, j, k), rho_w, cld_evp_eff, cons_h2omw, &
1158 & cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, &
1159 & cons_cp)
1160  CALL popcontrol1b(branch)
1161  IF (branch .EQ. 0) THEN
1162  CALL popreal8(cf_con(i, j, k))
1163  CALL popreal8(cf_ls(i, j, k))
1164  cf_totb = -(cf_ls(i, j, k)*cf_lsb(i, j, k)/cf_tot**2) - cf_con&
1165 & (i, j, k)*cf_conb(i, j, k)/cf_tot**2
1166  cf_conb(i, j, k) = cf_conb(i, j, k)/cf_tot
1167  cf_lsb(i, j, k) = cf_lsb(i, j, k)/cf_tot
1168  ELSE
1169  cf_totb = 0.0_8
1170  END IF
1171  CALL popreal8(cf_tot)
1172  cf_lsb(i, j, k) = cf_lsb(i, j, k) + cf_totb
1173  cf_conb(i, j, k) = cf_conb(i, j, k) + cf_totb
1174 
1175  !ADJOINT OF SAVE PRESINKS INPUTS
1176  qb(i,j,k) = qb(i,j,k) + q_p_presink
1177  qi_lsb(i,j,k) = qi_lsb(i,j,k) + qi_ls_p_presink
1178  qi_conb(i,j,k) = qi_conb(i,j,k) + qi_con_p_presink
1179  ql_lsb(i,j,k) = ql_lsb(i,j,k) + ql_ls_p_presink
1180  ql_conb(i,j,k) = ql_conb(i,j,k) + ql_con_p_presink
1181  cf_conb(i,j,k) = cf_conb(i,j,k) + cf_con_p_presink
1182 
1183  CALL popreal8(t(i, j, k))
1184  CALL popreal8(q(i, j, k))
1185  CALL popreal8(ql_ls(i, j, k))
1186  CALL popreal8(ql_con(i, j, k))
1187  CALL popreal8(qi_ls(i, j, k))
1188  CALL popreal8(qi_con(i, j, k))
1189  CALL popreal8(cf_ls(i, j, k))
1190  CALL popreal8(cf_con(i, j, k))
1191 
1192  if (do_moist_physics == 1) then
1193 
1194  !For data assimilation purposes always employ the perturbation model for cloud fraction
1195  cloud_pertmod = 1
1196 
1197  elseif (do_moist_physics == 2) then
1198 
1199  !For longer forecasts only employ perturabtion model if filtering deems necessary
1200  cloud_pertmod = 0
1201 
1202  !COMPUTE THE JACOBIAN OF LS_CLOUD AND USE RESULTS TO FILTER 'BAD' POINTS
1203  jacobian = 0.0
1204  DO ii = 1,8
1205  xx = 0.0
1206  xx(ii) = 1.0
1207 
1208  !Store seperate trajectory so as not to overwrite
1209  ttraj = t(i,j,k)
1210  qtraj = q(i,j,k)
1211  qi_lstraj = qi_ls(i,j,k)
1212  qi_contraj = qi_con(i,j,k)
1213  ql_lstraj = ql_ls(i,j,k)
1214  ql_contraj = ql_con(i,j,k)
1215  cf_lstraj = cf_ls(i,j,k)
1216  cf_contraj = cf_con(i,j,k)
1217  phtraj = ph(i,j,k)
1218 
1219  tpert = xx(1)
1220  qpert = xx(2)
1221  qi_lspert = xx(3)
1222  qi_conpert = xx(4)
1223  ql_lspert = xx(5)
1224  ql_conpert = xx(6)
1225  cf_lspert = xx(7)
1226  cf_conpert = xx(8)
1227 
1228  CALL ls_cloud_d(dt, alpha, pdfflag, phtraj, ttraj, tpert, qtraj, qpert, &
1229  ql_lstraj, ql_lspert, ql_contraj, ql_conpert, &
1230  qi_lstraj, qi_lspert, qi_contraj, qi_conpert, &
1231  cf_lstraj, cf_lspert, cf_contraj, cf_conpert, &
1232  cons_alhl, cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw, &
1233  t_ice_all, t_ice_max, icefrpwr, estblx, 0, do_moist_physics)
1234 
1235  jacobian(1,ii) = tpert
1236  jacobian(2,ii) = qpert
1237  jacobian(3,ii) = qi_lspert
1238  jacobian(4,ii) = qi_conpert
1239  jacobian(5,ii) = ql_lspert
1240  jacobian(6,ii) = ql_conpert
1241  jacobian(7,ii) = cf_lspert
1242  jacobian(8,ii) = cf_conpert
1243 
1244  endDO
1245 
1246  !Compute eigenvalues of the Jacobian
1247  a = jacobian
1248 
1249  lwork = -1
1250  work = 0.0
1251  info = 0
1252 
1253  CALL dgeev( 'Vectors', 'Vectors', n1, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info )
1254 
1255  lwork = min( lwmax, int( work( 1 ) ) )
1256  CALL dgeev( 'Vectors', 'Vectors', n1, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info )
1257 
1258  maxeval = maxval(abs(wr))
1259 
1260  !Filter based on the eigenvalues
1261  if (maxeval > 1.001) then
1262  cloud_pertmod = 1
1263  endif
1264 
1265  !Filter based on the Jacobian values
1266  if ( ( jacobian(1,1) < 0.6 ) .or. &
1267  ( jacobian(2,1) > 0.75e-4 ) .or. &
1268  ( jacobian(5,1) < -0.75e-4 ) .or. &
1269  ( jacobian(7,1) < -1.10 ) ) then
1270  cloud_pertmod = 1
1271  endif
1272 
1273  endif
1274 
1275  CALL ls_cloud_b(dt, alpha, pdfflag, ph(i, j, k), t(i, j, k), tb(&
1276 & i, j, k), q(i, j, k), qb(i, j, k), ql_ls(i, j, k), &
1277 & ql_lsb(i, j, k), ql_con(i, j, k), ql_conb(i, j, k), &
1278 & qi_ls(i, j, k), qi_lsb(i, j, k), qi_con(i, j, k), &
1279 & qi_conb(i, j, k), cf_ls(i, j, k), cf_lsb(i, j, k), &
1280 & cf_con(i, j, k), cf_conb(i, j, k), cons_alhl, &
1281 & cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw&
1282 & , t_ice_all, t_ice_max, icefrpwr, estblx, &
1283 & cloud_pertmod, do_moist_physics)
1284 
1285  CALL popreal8(rhcrit)
1286  CALL popreal8(alpha)
1287  CALL popreal8(t(i, j, k))
1288  CALL popreal8(q(i, j, k))
1289  CALL popreal8(ql_con(i, j, k))
1290  CALL popreal8(qi_con(i, j, k))
1291  CALL popreal8(cf_con(i, j, k))
1292  CALL convec_src_b(dt, mass(i, j, k), imass(i, j, k), t(i, j, k)&
1293 & , tb(i, j, k), q(i, j, k), qb(i, j, k), cnv_dqldt(i&
1294 & , j, k), cnv_dqldtb(i, j, k), cnv_mfd(i, j, k), &
1295 & cnv_mfdb(i, j, k), ql_con(i, j, k), ql_conb(i, j, k)&
1296 & , qi_con(i, j, k), qi_conb(i, j, k), cf_con(i, j, k)&
1297 & , cf_conb(i, j, k), qs(i, j, k), qsb(i, j, k), &
1298 & cons_alhs, cons_alhl, cons_cp, t_ice_all, t_ice_max&
1299 & , icefrpwr)
1300  CALL popreal8(t(i, j, k))
1301  CALL popreal8(ql_con(i, j, k))
1302  CALL popreal8(qi_con(i, j, k))
1303  CALL meltfreeze_b(dt, t(i, j, k), tb(i, j, k), ql_con(i, j, k), &
1304 & ql_conb(i, j, k), qi_con(i, j, k), qi_conb(i, j, k)&
1305 & , t_ice_all, t_ice_max, icefrpwr, cons_alhl, &
1306 & cons_alhs, cons_cp)
1307  CALL popreal8(t(i, j, k))
1308  CALL popreal8(ql_ls(i, j, k))
1309  CALL popreal8(qi_ls(i, j, k))
1310  CALL meltfreeze_b(dt, t(i, j, k), tb(i, j, k), ql_ls(i, j, k), &
1311 & ql_lsb(i, j, k), qi_ls(i, j, k), qi_lsb(i, j, k), &
1312 & t_ice_all, t_ice_max, icefrpwr, cons_alhl, cons_alhs&
1313 & , cons_cp)
1314  CALL popreal8(t(i, j, k))
1315  CALL popreal8(ql_ls(i, j, k))
1316  CALL popreal8(qi_ls(i, j, k))
1317  CALL popreal8(cf_ls(i, j, k))
1318  CALL popreal8(ql_con(i, j, k))
1319  CALL popreal8(qi_con(i, j, k))
1320  CALL popreal8(cf_con(i, j, k))
1321  CALL cloud_tidy_b(q(i, j, k), qb(i, j, k), t(i, j, k), tb(i, j, &
1322 & k), ql_ls(i, j, k), ql_lsb(i, j, k), qi_ls(i, j, k)&
1323 & , qi_lsb(i, j, k), cf_ls(i, j, k), cf_lsb(i, j, k), &
1324 & ql_con(i, j, k), ql_conb(i, j, k), qi_con(i, j, k), &
1325 & qi_conb(i, j, k), cf_con(i, j, k), cf_conb(i, j, k)&
1326 & , cons_alhl, cons_alhs, cons_cp)
1327  cnv_prc3b(i, j, k) = cnv_prc3b(i, j, k) + qrn_cu_1db
1328  CALL popcontrol1b(branch)
1329  IF (branch .NE. 0) THEN
1330  CALL popreal8(area_ls_prc)
1331  CALL popreal8(area_anv_prc)
1332  CALL popreal8(area_upd_prc)
1333  CALL popreal8(tot_prec_ls)
1334  CALL popreal8(tot_prec_anv)
1335  CALL popreal8(tot_prec_upd)
1336  area_upd_prcb = 0.0_8
1337  tot_prec_lsb = 0.0_8
1338  area_ls_prcb = 0.0_8
1339  tot_prec_updb = 0.0_8
1340  area_anv_prcb = 0.0_8
1341  tot_prec_anvb = 0.0_8
1342  END IF
1343 
1344  !Adjoint save the inputs to the scheme for filtering
1345  tb(i, j, k) = tb(i, j, k) + t_p_preall
1346  ql_lsb(i, j, k) = ql_lsb(i, j, k) + ql_ls_p_preall
1347  ql_conb(i, j, k) = ql_conb(i, j, k) + ql_con_p_preall
1348  qi_lsb(i, j, k) = qi_lsb(i, j, k) + qi_ls_p_preall
1349  qi_conb(i, j, k) = qi_conb(i, j, k) + qi_con_p_preall
1350 
1351  END DO
1352  END DO
1353  END DO
1354  vmipb = 0.0_8
1355  DO k=lm,1,-1
1356  CALL popreal8array(qddf3(:, :, k), im*jm)
1357  vmipb = vmipb - qddf3(:, :, k)*qddf3b(:, :, k)/vmip**2
1358  qddf3b(:, :, k) = qddf3b(:, :, k)/vmip
1359  END DO
1360  DO i=im,1,-1
1361  DO j=jm,1,-1
1362  qddf3b(i, j, :) = qddf3b(i, j, :) + vmipb(i, j)
1363  vmipb(i, j) = 0.0_8
1364  END DO
1365  END DO
1366  WHERE (.NOT.mask(:, :, 1:lm)) qddf3b = 0.0_8
1367  zetb = 0.0_8
1368  WHERE (mask(:, :, 1:lm)) zetb(:, :, 1:lm) = zetb(:, :, 1:lm) - (zet(:&
1369 & , :, 1:lm)-3000.)*mass*qddf3b - zet(:, :, 1:lm)*mass*qddf3b
1370  DO k=1,lm,1
1371  zetb(:, :, k+1) = zetb(:, :, k+1) + zetb(:, :, k)
1372  dzetb(:, :, k) = dzetb(:, :, k) + zetb(:, :, k)
1373  zetb(:, :, k) = 0.0_8
1374  END DO
1375  CALL dqsat_bac_b(dqsdt, dqsdtb, qs, qsb, t, tb, ph, im, jm, lm, estblx&
1376 & , cons_h2omw, cons_airmw)
1377  thb = 0.0_8
1378  thb(:, :, 1:lm) = pih*tb + (pi(:, :, 1:lm)-pi(:, :, 0:lm-1))*cons_cp*&
1379 & dzetb(:, :, 1:lm)/cons_grav
1380 END SUBROUTINE cloud_driver_b
1381 
1382 ! Differentiation of cloud_tidy in reverse (adjoint) mode:
1383 ! gradient of useful results: af qv qla qlc qia qic cf te
1384 ! with respect to varying inputs: af qv qla qlc qia qic cf te
1385 !SUBROUTINES
1386 SUBROUTINE cloud_tidy_b(qv, qvb, te, teb, qlc, qlcb, qic, qicb, cf, cfb&
1387 & , qla, qlab, qia, qiab, af, afb, cons_alhl, cons_alhs, cons_cp)
1388  IMPLICIT NONE
1389  REAL*8, INTENT(INOUT) :: te, qv, qlc, cf, qla, af, qic, qia
1390  REAL*8 :: teb, qvb, qlcb, cfb, qlab, afb, qicb, qiab
1391  REAL*8, INTENT(IN) :: cons_alhl, cons_alhs, cons_cp
1392  INTEGER :: branch
1393 !Fix if Anvil cloud fraction too small
1394  IF (af .LT. 1.e-5) THEN
1395  CALL pushreal8(qla)
1396  qla = 0.
1397  CALL pushreal8(qia)
1398  qia = 0.
1399  CALL pushcontrol1b(0)
1400  ELSE
1401  CALL pushcontrol1b(1)
1402  END IF
1403 !Fix if LS cloud fraction too small
1404 ! if ( CF < 1.E-5 ) then
1405 ! QV = QV + QLC + QIC
1406 ! TE = TE - (CONS_ALHL/CONS_CP)*QLC - (CONS_ALHS/CONS_CP)*QIC
1407 ! CF = 0.
1408 ! QLC = 0.
1409 ! QIC = 0.
1410 ! end if
1411 !LS LIQUID too small
1412  IF (qlc .LT. 1.e-8) THEN
1413  CALL pushreal8(qlc)
1414  qlc = 0.
1415  CALL pushcontrol1b(0)
1416  ELSE
1417  CALL pushcontrol1b(1)
1418  END IF
1419 !LS ICE too small
1420  IF (qic .LT. 1.e-8) THEN
1421  CALL pushreal8(qic)
1422  qic = 0.
1423  CALL pushcontrol1b(0)
1424  ELSE
1425  CALL pushcontrol1b(1)
1426  END IF
1427 !Anvil LIQUID too small
1428  IF (qla .LT. 1.e-8) THEN
1429  CALL pushreal8(qla)
1430  qla = 0.
1431  CALL pushcontrol1b(0)
1432  ELSE
1433  CALL pushcontrol1b(1)
1434  END IF
1435 !Anvil ICE too small
1436  IF (qia .LT. 1.e-8) THEN
1437  CALL pushreal8(qia)
1438  qia = 0.
1439  CALL pushcontrol1b(0)
1440  ELSE
1441  CALL pushcontrol1b(1)
1442  END IF
1443 !Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small
1444  IF (qla + qia .LT. 1.e-8) THEN
1445  CALL pushcontrol1b(0)
1446  ELSE
1447  CALL pushcontrol1b(1)
1448  END IF
1449 !Ditto if LS cloud LIQUID+ICE too small
1450  IF (qlc + qic .LT. 1.e-8) THEN
1451  qlcb = qvb - cons_alhl*teb/cons_cp
1452  qicb = qvb - cons_alhs*teb/cons_cp
1453  cfb = 0.0_8
1454  END IF
1455  CALL popcontrol1b(branch)
1456  IF (branch .EQ. 0) THEN
1457  qlab = qvb - cons_alhl*teb/cons_cp
1458  qiab = qvb - cons_alhs*teb/cons_cp
1459  afb = 0.0_8
1460  END IF
1461  CALL popcontrol1b(branch)
1462  IF (branch .EQ. 0) THEN
1463  CALL popreal8(qia)
1464  qiab = qvb - cons_alhs*teb/cons_cp
1465  END IF
1466  CALL popcontrol1b(branch)
1467  IF (branch .EQ. 0) THEN
1468  CALL popreal8(qla)
1469  qlab = qvb - cons_alhl*teb/cons_cp
1470  END IF
1471  CALL popcontrol1b(branch)
1472  IF (branch .EQ. 0) THEN
1473  CALL popreal8(qic)
1474  qicb = qvb - cons_alhs*teb/cons_cp
1475  END IF
1476  CALL popcontrol1b(branch)
1477  IF (branch .EQ. 0) THEN
1478  CALL popreal8(qlc)
1479  qlcb = qvb - cons_alhl*teb/cons_cp
1480  END IF
1481  CALL popcontrol1b(branch)
1482  IF (branch .EQ. 0) THEN
1483  CALL popreal8(qia)
1484  CALL popreal8(qla)
1485  qlab = qvb - cons_alhl*teb/cons_cp
1486  qiab = qvb - cons_alhs*teb/cons_cp
1487  afb = 0.0_8
1488  END IF
1489 END SUBROUTINE cloud_tidy_b
1490 
1491 ! Differentiation of meltfreeze in reverse (adjoint) mode:
1492 ! gradient of useful results: qi ql te
1493 ! with respect to varying inputs: qi ql te
1494 SUBROUTINE meltfreeze_b(dt, te, teb, ql, qlb, qi, qib, t_ice_all, &
1495 & t_ice_max, icefrpwr, cons_alhl, cons_alhs, cons_cp)
1496  IMPLICIT NONE
1497 !Inputs
1498  REAL*8, INTENT(IN) :: dt, t_ice_all, t_ice_max
1499  INTEGER, INTENT(IN) :: icefrpwr
1500  REAL*8, INTENT(IN) :: cons_alhl, cons_alhs, cons_cp
1501 !Prognostic
1502  REAL*8, INTENT(INOUT) :: te, ql, qi
1503  REAL*8 :: teb, qlb, qib
1504 !Locals
1505  REAL*8 :: fqi, dqil
1506  REAL*8 :: fqib, dqilb
1507  REAL*8, PARAMETER :: taufrz=1000.
1508  INTRINSIC exp
1509  INTRINSIC max
1510  INTRINSIC min
1511  INTEGER :: branch
1512  REAL*8 :: temp
1513  dqil = 0.0
1514  CALL get_ice_fraction(te, t_ice_all, t_ice_max, icefrpwr, fqi)
1515 !Freeze liquid
1516  IF (te .LE. t_ice_max) THEN
1517  dqil = ql*(1.0-exp(-(dt*fqi/taufrz)))
1518  CALL pushcontrol1b(1)
1519  ELSE
1520  CALL pushcontrol1b(0)
1521  END IF
1522  IF (0. .LT. dqil) THEN
1523  CALL pushcontrol1b(0)
1524  dqil = dqil
1525  ELSE
1526  dqil = 0.
1527  CALL pushcontrol1b(1)
1528  END IF
1529  CALL pushreal8(qi)
1530  qi = qi + dqil
1531  CALL pushreal8(te)
1532  te = te + (cons_alhs-cons_alhl)*dqil/cons_cp
1533  dqil = 0.
1534 !Melt ice instantly above 0^C
1535  IF (te .GT. t_ice_max) THEN
1536  dqil = -qi
1537  CALL pushcontrol1b(1)
1538  ELSE
1539  CALL pushcontrol1b(0)
1540  END IF
1541  IF (0. .GT. dqil) THEN
1542  CALL pushcontrol1b(0)
1543  ELSE
1544  CALL pushcontrol1b(1)
1545  END IF
1546  dqilb = qib - qlb + (cons_alhs-cons_alhl)*teb/cons_cp
1547  CALL popcontrol1b(branch)
1548  IF (branch .NE. 0) dqilb = 0.0_8
1549  CALL popcontrol1b(branch)
1550  IF (branch .NE. 0) qib = qib - dqilb
1551  CALL popreal8(te)
1552  dqilb = qib - qlb + (cons_alhs-cons_alhl)*teb/cons_cp
1553  CALL popreal8(qi)
1554  CALL popcontrol1b(branch)
1555  IF (branch .NE. 0) dqilb = 0.0_8
1556  CALL popcontrol1b(branch)
1557  IF (branch .EQ. 0) THEN
1558  fqib = 0.0_8
1559  ELSE
1560  temp = -(dt*fqi/taufrz)
1561  qlb = qlb + (1.0-exp(temp))*dqilb
1562  fqib = dt*exp(temp)*ql*dqilb/taufrz
1563  END IF
1564  CALL get_ice_fraction_b(te, teb, t_ice_all, t_ice_max, icefrpwr, fqi, &
1565 & fqib)
1566 END SUBROUTINE meltfreeze_b
1567 
1568 ! Differentiation of convec_src in reverse (adjoint) mode:
1569 ! gradient of useful results: af qs qv qla qia dcf dmf te
1570 ! with respect to varying inputs: af qs qv qla qia dcf dmf te
1571 SUBROUTINE convec_src_b(dt, mass, imass, te, teb, qv, qvb, dcf, dcfb, &
1572 & dmf, dmfb, qla, qlab, qia, qiab, af, afb, qs, qsb, cons_alhs, &
1573 & cons_alhl, cons_cp, t_ice_all, t_ice_max, icefrpwr)
1574  IMPLICIT NONE
1575 !Inputs
1576  REAL*8, INTENT(IN) :: dt, t_ice_all, t_ice_max
1577  INTEGER, INTENT(IN) :: icefrpwr
1578  REAL*8, INTENT(IN) :: mass, imass, qs
1579  REAL*8 :: qsb
1580  REAL*8, INTENT(IN) :: dmf, dcf
1581  REAL*8 :: dmfb, dcfb
1582  REAL*8, INTENT(IN) :: cons_alhs, cons_alhl, cons_cp
1583 !Prognostic
1584  REAL*8, INTENT(INOUT) :: te, qv
1585  REAL*8 :: teb, qvb
1586  REAL*8, INTENT(INOUT) :: qla, qia, af
1587  REAL*8 :: qlab, qiab, afb
1588 !Locals
1589 !Minimum allowed env RH
1590  REAL*8, PARAMETER :: minrhx=0.001
1591  REAL*8 :: tend, qvx, fqi
1592  REAL*8 :: tendb, fqib
1593  INTRINSIC min
1594  INTEGER :: branch
1595  REAL*8 :: tempb0
1596  REAL*8 :: tempb
1597 !Namelist
1598 !DT - Timestep
1599 !MASS - Level Mass
1600 !iMASS - 1/Level Mass
1601 !TE - Temperature
1602 !QV - Specific Humidity
1603 !DCF - CNV_DQL from RAS
1604 !DMF - CNV_MFD from RAS
1605 !QLA - Convective cloud liquid water
1606 !QIA - Convective cloud liquid ice
1607 !AF - Convective cloud fraction
1608 !QS - Qsat
1609 !Zero out locals
1610 !Addition of condensate from RAS
1611  CALL get_ice_fraction(te, t_ice_all, t_ice_max, icefrpwr, fqi)
1612 !Convective condensation has never frozen so latent heat of fusion
1613 !Compute Tiedtke-style anvil fraction
1614  tend = dmf*imass
1615  CALL pushreal8(af)
1616  af = af + tend*dt
1617  IF (af .GT. 0.99) THEN
1618  af = 0.99
1619  CALL pushcontrol1b(0)
1620  ELSE
1621  CALL pushcontrol1b(1)
1622  af = af
1623  END IF
1624 ! Check for funny (tiny, negative) external QV, resulting from assumed QV=QSAT within anvil.
1625 ! Simply constrain AF assume condensate just gets "packed" in
1626  IF (af .LT. 1.0) THEN
1627  qvx = (qv-qs*af)/(1.-af)
1628  ELSE
1629  qvx = qs
1630  END IF
1631 !If saturated over critial value and there is Anvil
1632  IF (qvx - minrhx*qs .LT. 0.0 .AND. af .GT. 0.) THEN
1633  af = (qv-minrhx*qs)/(qs*(1.0-minrhx))
1634  CALL pushcontrol1b(0)
1635  ELSE
1636  CALL pushcontrol1b(1)
1637  END IF
1638  IF (af .LT. 0.) THEN
1639  qlab = qvb - cons_alhl*teb/cons_cp
1640  qiab = qvb - cons_alhs*teb/cons_cp
1641  afb = 0.0_8
1642  END IF
1643  CALL popcontrol1b(branch)
1644  IF (branch .EQ. 0) THEN
1645  tempb0 = afb/((1.0-minrhx)*qs)
1646  qvb = qvb + tempb0
1647  qsb = qsb + (-((qv-minrhx*qs)/qs)-minrhx)*tempb0
1648  afb = 0.0_8
1649  END IF
1650  CALL popcontrol1b(branch)
1651  IF (branch .EQ. 0) afb = 0.0_8
1652  CALL popreal8(af)
1653  tendb = dt*afb
1654  dmfb = dmfb + imass*tendb
1655  tend = dcf*imass
1656  tempb = (cons_alhs-cons_alhl)*dt*teb/cons_cp
1657  fqib = dt*tend*qiab - tend*dt*qlab + tend*tempb
1658  tendb = dt*fqi*qiab + dt*(1.0-fqi)*qlab + fqi*tempb
1659  CALL get_ice_fraction_b(te, teb, t_ice_all, t_ice_max, icefrpwr, fqi, &
1660 & fqib)
1661  dcfb = dcfb + imass*tendb
1662 END SUBROUTINE convec_src_b
1663 
1664 ! Differentiation of ls_cloud in reverse (adjoint) mode:
1665 ! gradient of useful results: qai qal af qv qci qcl cf te
1666 ! with respect to varying inputs: qai qal af qv qci qcl cf te
1667 SUBROUTINE ls_cloud_b(dt, alpha, pdfshape, pl, te, teb, qv, qvb, qcl, &
1668 & qclb, qal, qalb, qci, qcib, qai, qaib, cf, cfb, af, afb, cons_alhl, &
1669 & cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw, t_ice_all, &
1670 & t_ice_max, icefrpwr, estblx, cloud_pertmod, dmp)
1671  IMPLICIT NONE
1672 !Inputs
1673  REAL*8, INTENT(IN) :: dt, alpha, pl, t_ice_all, t_ice_max
1674  INTEGER, INTENT(IN) :: pdfshape, cloud_pertmod, dmp
1675  INTEGER, INTENT(IN) :: icefrpwr
1676  REAL*8, INTENT(IN) :: cons_alhl, cons_alhf, cons_alhs, cons_cp, &
1677 & cons_h2omw, cons_airmw
1678  REAL*8, INTENT(IN) :: estblx(:)
1679 !Prognostic
1680  REAL*8, INTENT(INOUT) :: te, qv, qcl, qci, qal, qai, cf, af
1681  REAL*8 :: teb, qvb, qclb, qcib, qalb, qaib, cfb, afb
1682 !Locals
1683  INTEGER :: n
1684  REAL*8 :: qco, cfo, qao, qt, qmx, qmn, dq
1685  REAL*8 :: qcob, cfob, qaob, qtb
1686  REAL*8 :: teo, qsx, dqsx, qs, dqs, tmparr
1687  REAL*8 :: teob, qsxb, dqsxb, dqsb, tmparrb
1688  REAL*8 :: qcx, qvx, cfx, qax, qc, qa, fqi, fqi_a, dqai, dqal, dqci, &
1689 & dqcl
1690  REAL*8 :: qcxb, qvxb, cfxb, qaxb, qcb, qab, fqib, dqaib, dqalb, dqcib&
1691 & , dqclb
1692  REAL*8 :: ten, qsp, cfp, qvp, qcp
1693  REAL*8 :: tenb, qcpb
1694  REAL*8 :: tep, qsn, cfn, qvn, qcn
1695  REAL*8 :: tepb, qsnb, cfnb, qcnb
1696  REAL*8 :: alhx, sigmaqt1, sigmaqt2
1697  REAL*8 :: alhxb, sigmaqt1b, sigmaqt2b
1698  REAL*8, DIMENSION(1) :: dqsx1, qsx1, teo1, pl1
1699  INTRINSIC max
1700  INTEGER :: branch
1701  REAL*8 :: temp1
1702  REAL*8 :: temp0
1703  REAL*8 :: tempb5
1704  REAL*8 :: tempb4
1705  REAL*8 :: tempb3
1706  REAL*8 :: tempb2
1707  REAL*8 :: tempb1
1708  REAL*8 :: tempb0
1709  REAL*8 :: tempb
1710  REAL*8 :: temp
1711 !Namelist
1712 !DT - Timestep
1713 !ALPHA - PDF half width
1714 !PL - Pressure (hPa)
1715 !TE - Temperature
1716 !QV - Specific humidity
1717 !QCl - Convective cloud liquid water
1718 !QAl - Large scale cloud liquid water
1719 !QCi - Convective cloud liquid ice
1720 !QAi - Large scale cloud liquid ice
1721 !CF - Large scale cloud fraction
1722 !AF - Convective cloud fraction
1723  qc = qcl + qci
1724  qa = qal + qai
1725  teo = te
1726  CALL dqsats_bac(dqsx, qsx, teo, pl, estblx, cons_h2omw, cons_airmw)
1727  IF (af .LT. 1.0) THEN
1728  IF (dmp .EQ. 1) THEN
1729  IF (1. - af .GT. 0.02) THEN
1730  tmparr = 1./(1.-af)
1731  CALL pushcontrol3b(0)
1732  ELSE
1733  tmparr = 1./(1.-af)
1734  CALL pushcontrol3b(1)
1735  END IF
1736  ELSE IF (dmp .EQ. 2) THEN
1737  tmparr = 1./(1.-af)
1738  CALL pushcontrol3b(2)
1739  ELSE
1740  CALL pushcontrol3b(3)
1741  END IF
1742  ELSE
1743  CALL pushcontrol3b(4)
1744  tmparr = 0.0
1745  END IF
1746  cfx = cf*tmparr
1747  qcx = qc*tmparr
1748  qvx = (qv-qsx*af)*tmparr
1749  IF (af .GE. 1.0) THEN
1750  qvx = qsx*1.e-4
1751  CALL pushcontrol1b(0)
1752  ELSE
1753  CALL pushcontrol1b(1)
1754  END IF
1755  IF (af .GT. 0.) THEN
1756  qax = qa/af
1757  CALL pushcontrol1b(0)
1758  ELSE
1759  CALL pushcontrol1b(1)
1760  qax = 0.
1761  END IF
1762  qt = qcx + qvx
1763  ten = teo
1764  cfn = cfx
1765  qcn = qcx
1766 !Begin iteration
1767 !do n=1,4
1768  n = 1
1769  qcp = qcn
1770 !Dont call again as not looping
1771  dqs = dqsx
1772  qsn = qsx
1773 !call DQSATs_BAC(DQS, QSn, TEn, PL, ESTBLX, CONS_H2OMW, CONS_AIRMW)
1774  tep = ten
1775  CALL get_ice_fraction(tep, t_ice_all, t_ice_max, icefrpwr, fqi)
1776  sigmaqt1 = alpha*qsn
1777  sigmaqt2 = alpha*qsn
1778  IF (pdfshape .EQ. 2) THEN
1779 !For triangular, symmetric: sigmaqt1 = sigmaqt2 = alpha*qsn (alpha is half width)
1780 !For triangular, skewed r : sigmaqt1 < sigmaqt2
1781  sigmaqt1 = alpha*qsn
1782  sigmaqt2 = alpha*qsn
1783  CALL pushcontrol1b(0)
1784  ELSE
1785  CALL pushcontrol1b(1)
1786  END IF
1787 !Compute cloud fraction
1788  IF (cloud_pertmod .EQ. 0) THEN
1789  CALL pdffrac(1, qt, sigmaqt1, sigmaqt2, qsn, cfn)
1790  CALL pushcontrol2b(0)
1791  ELSE IF (cloud_pertmod .EQ. 1) THEN
1792  CALL pdffrac(4, qt, sigmaqt1, sigmaqt2, qsn, cfn)
1793  CALL pushcontrol2b(1)
1794  ELSE
1795  CALL pushcontrol2b(2)
1796  END IF
1797 !Compute cloud condensate
1798  CALL pdfcondensate(pdfshape, qt, sigmaqt1, sigmaqt2, qsn, qcn)
1799 !Adjustments to anvil condensate due to the assumption of a stationary TOTAL
1800 !water PDF subject to a varying QSAT value during the iteration
1801  IF (af .GT. 0.) THEN
1802 ! + QSx - QS
1803  qao = qax
1804  CALL pushcontrol1b(0)
1805  ELSE
1806  qao = 0.
1807  CALL pushcontrol1b(1)
1808  END IF
1809  alhx = (1.0-fqi)*cons_alhl + fqi*cons_alhs
1810  IF (pdfshape .EQ. 1) THEN
1811  CALL pushreal8(qcn)
1812  qcn = qcp + (qcn-qcp)/(1.-(cfn*(alpha-1.)-qcn/qsn)*dqs*alhx/cons_cp)
1813  CALL pushcontrol2b(0)
1814  ELSE IF (pdfshape .EQ. 2) THEN
1815 !This next line needs correcting - need proper d(del qc)/dT derivative for triangular
1816 !for now, just use relaxation of 1/2.
1817  IF (n .NE. 4) THEN
1818  qcn = qcp + (qcn-qcp)*0.5
1819  CALL pushcontrol2b(1)
1820  ELSE
1821  CALL pushcontrol2b(2)
1822  END IF
1823  ELSE
1824  CALL pushcontrol2b(3)
1825  END IF
1826 !enddo ! qsat iteration
1827  cfo = cfn
1828  qco = qcn
1829 ! Update prognostic variables. QCo, QAo become updated grid means.
1830  IF (af .LT. 1.0) THEN
1831  CALL pushreal8(cf)
1832  cf = cfo*(1.-af)
1833  CALL pushreal8(qco)
1834  qco = qco*(1.-af)
1835  CALL pushreal8(qao)
1836  qao = qao*af
1837  CALL pushcontrol2b(0)
1838  ELSE
1839 !Grid box filled with anvil
1840 ! Special case AF=1, i.e., box filled with anvil. Note: no guarantee QV_box > QS_box
1841 ! Remove any other cloud
1842  CALL pushreal8(cf)
1843  cf = 0.
1844 ! Add any LS condensate to anvil type
1845  qao = qa + qc
1846 ! Remove same from LS
1847  qco = 0.
1848 ! Total water
1849  qt = qao + qv
1850  IF (qt - qsx .LT. 0.) THEN
1851  qao = 0.
1852  CALL pushcontrol2b(2)
1853  ELSE
1854  qao = qt - qsx
1855  CALL pushcontrol2b(1)
1856  END IF
1857  END IF
1858 !Partition new condensate into ice and liquid taking care to keep both >=0 separately.
1859 !New condensate can be less than old, so Delta can be < 0.
1860  qcx = qco - qc
1861  dqcl = (1.0-fqi)*qcx
1862  dqci = fqi*qcx
1863 !Large Scale Partition
1864  IF (qcl + dqcl .LT. 0.) THEN
1865  dqci = dqci + (qcl+dqcl)
1866 !== dQCl - (QCl+dQCl)
1867  CALL pushcontrol1b(0)
1868  ELSE
1869  CALL pushcontrol1b(1)
1870  END IF
1871  IF (qci + dqci .LT. 0.) THEN
1872  CALL pushcontrol1b(0)
1873  ELSE
1874  CALL pushcontrol1b(1)
1875  END IF
1876  qax = qao - qa
1877 ! (1.0-fQi)*QAx
1878  dqal = qax
1879 ! fQi * QAx
1880  dqai = 0.
1881 !Convective partition
1882  IF (qal + dqal .LT. 0.) THEN
1883  dqai = dqai + (qal+dqal)
1884  CALL pushcontrol1b(0)
1885  ELSE
1886  CALL pushcontrol1b(1)
1887  END IF
1888  IF (qai + dqai .LT. 0.) THEN
1889  CALL pushcontrol1b(0)
1890  ELSE
1891  CALL pushcontrol1b(1)
1892  END IF
1893 ! Clean-up cloud if fractions are too small
1894  IF (af .LT. 1.e-5) THEN
1895  CALL pushcontrol1b(0)
1896  ELSE
1897  CALL pushcontrol1b(1)
1898  END IF
1899  IF (cf .LT. 1.e-5) THEN
1900  CALL pushcontrol1b(0)
1901  ELSE
1902  CALL pushcontrol1b(1)
1903  END IF
1904 !Update specific humidity
1905 !Update temperature
1906 !Take care of situations where QS moves past QA during QSAT iteration (QAo negative).
1907 !"Evaporate" offending QA
1908  IF (qao .LE. 0.) THEN
1909  qaib = qvb - cons_alhs*teb/cons_cp
1910  qalb = qvb - cons_alhl*teb/cons_cp
1911  afb = 0.0_8
1912  END IF
1913  tempb4 = teb/cons_cp
1914  tempb5 = cons_alhl*tempb4
1915  dqaib = qaib - qvb + cons_alhf*tempb4 + tempb5
1916  dqcib = qcib - qvb + cons_alhf*tempb4 + tempb5
1917  dqalb = qalb - qvb + tempb5
1918  dqclb = qclb - qvb + tempb5
1919  CALL popcontrol1b(branch)
1920  IF (branch .EQ. 0) THEN
1921  qclb = qclb - dqclb
1922  qcib = qcib - dqcib
1923  dqcib = 0.0_8
1924  dqclb = 0.0_8
1925  END IF
1926  CALL popcontrol1b(branch)
1927  IF (branch .EQ. 0) THEN
1928  qalb = qalb - dqalb
1929  qaib = qaib - dqaib
1930  dqaib = 0.0_8
1931  dqalb = 0.0_8
1932  END IF
1933  CALL popcontrol1b(branch)
1934  IF (branch .EQ. 0) THEN
1935  qaib = qaib + dqalb - dqaib
1936  dqaib = dqalb
1937  END IF
1938  CALL popcontrol1b(branch)
1939  IF (branch .EQ. 0) THEN
1940  qalb = qalb + dqaib - dqalb
1941  dqalb = dqaib
1942  END IF
1943  qaxb = dqalb
1944  qaob = qaxb
1945  qab = -qaxb
1946  CALL popcontrol1b(branch)
1947  IF (branch .EQ. 0) THEN
1948  qcib = qcib + dqclb - dqcib
1949  dqcib = dqclb
1950  END IF
1951  CALL popcontrol1b(branch)
1952  IF (branch .EQ. 0) THEN
1953  qclb = qclb + dqcib - dqclb
1954  dqclb = dqcib
1955  END IF
1956  fqib = qcx*dqcib - qcx*dqclb
1957  qcxb = (1.0-fqi)*dqclb + fqi*dqcib
1958  qcob = qcxb
1959  qcb = -qcxb
1960  CALL popcontrol2b(branch)
1961  IF (branch .EQ. 0) THEN
1962  CALL popreal8(qao)
1963  CALL popreal8(qco)
1964  afb = afb + qao*qaob - cfo*cfb - qco*qcob
1965  qaob = af*qaob
1966  qcob = (1.-af)*qcob
1967  CALL popreal8(cf)
1968  cfob = (1.-af)*cfb
1969  qsxb = 0.0_8
1970  ELSE
1971  IF (branch .EQ. 1) THEN
1972  qtb = qaob
1973  qsxb = -qaob
1974  ELSE
1975  qtb = 0.0_8
1976  qsxb = 0.0_8
1977  END IF
1978  qaob = qtb
1979  qvb = qvb + qtb
1980  qab = qab + qaob
1981  qcb = qcb + qaob
1982  CALL popreal8(cf)
1983  qaob = 0.0_8
1984  cfob = 0.0_8
1985  qcob = 0.0_8
1986  END IF
1987  qcnb = qcob
1988  cfnb = cfob
1989  CALL popcontrol2b(branch)
1990  IF (branch .LT. 2) THEN
1991  IF (branch .EQ. 0) THEN
1992  CALL popreal8(qcn)
1993  temp0 = dqs*alhx/cons_cp
1994  temp1 = (alpha-1.)*cfn - qcn/qsn
1995  temp = -(temp1*temp0) + 1.
1996  tempb0 = qcnb/temp
1997  tempb1 = -((qcn-qcp)*tempb0/temp)
1998  tempb2 = temp0*tempb1/qsn
1999  tempb3 = -(temp1*tempb1/cons_cp)
2000  qcpb = qcnb - tempb0
2001  cfnb = cfnb - temp0*(alpha-1.)*tempb1
2002  qsnb = -(qcn*tempb2/qsn)
2003  dqsb = alhx*tempb3
2004  alhxb = dqs*tempb3
2005  qcnb = tempb2 + tempb0
2006  GOTO 100
2007  ELSE
2008  qcpb = 0.5*qcnb
2009  qcnb = 0.5*qcnb
2010  END IF
2011  ELSE IF (branch .EQ. 2) THEN
2012  qcpb = 0.0_8
2013  ELSE
2014  qcpb = 0.0_8
2015  END IF
2016  dqsb = 0.0_8
2017  alhxb = 0.0_8
2018  qsnb = 0.0_8
2019  100 fqib = fqib + (cons_alhs-cons_alhl)*alhxb
2020  CALL popcontrol1b(branch)
2021  IF (branch .EQ. 0) THEN
2022  qaxb = qaob
2023  ELSE
2024  qaxb = 0.0_8
2025  END IF
2026  qcx = qc*tmparr
2027  qt = qcx + qvx
2028  CALL pdfcondensate_b(pdfshape, qt, qtb, sigmaqt1, sigmaqt1b, sigmaqt2&
2029 & , sigmaqt2b, qsn, qsnb, qcn, qcnb)
2030  CALL popcontrol2b(branch)
2031  IF (branch .EQ. 0) THEN
2032  CALL pdffrac_b(1, qt, qtb, sigmaqt1, sigmaqt1b, sigmaqt2, sigmaqt2b&
2033 & , qsn, qsnb, cfn, cfnb)
2034  ELSE IF (branch .EQ. 1) THEN
2035  CALL pdffrac_b(4, qt, qtb, sigmaqt1, sigmaqt1b, sigmaqt2, sigmaqt2b&
2036 & , qsn, qsnb, cfn, cfnb)
2037  END IF
2038  CALL popcontrol1b(branch)
2039  IF (branch .EQ. 0) THEN
2040  qsnb = qsnb + alpha*sigmaqt1b + alpha*sigmaqt2b
2041  sigmaqt1b = 0.0_8
2042  sigmaqt2b = 0.0_8
2043  END IF
2044  qsnb = qsnb + alpha*sigmaqt1b + alpha*sigmaqt2b
2045  tepb = 0.0_8
2046  CALL get_ice_fraction_b(tep, tepb, t_ice_all, t_ice_max, icefrpwr, fqi&
2047 & , fqib)
2048  tenb = tepb
2049  qsxb = qsxb + qsnb
2050  dqsxb = dqsb
2051  qcnb = qcpb
2052  qcxb = qtb + qcnb
2053  cfxb = cfnb
2054  teob = tenb
2055  qvxb = qtb
2056  CALL popcontrol1b(branch)
2057  IF (branch .EQ. 0) THEN
2058  qab = qab + qaxb/af
2059  afb = afb - qa*qaxb/af**2
2060  END IF
2061  CALL popcontrol1b(branch)
2062  IF (branch .EQ. 0) THEN
2063  qsxb = qsxb + 1.e-4*qvxb
2064  qvxb = 0.0_8
2065  END IF
2066  tempb = tmparr*qvxb
2067  qvb = qvb + tempb
2068  qsxb = qsxb - af*tempb
2069  afb = afb - qsx*tempb
2070  tmparrb = qc*qcxb + cf*cfxb + (qv-qsx*af)*qvxb
2071  qcb = qcb + tmparr*qcxb
2072  cfb = tmparr*cfxb
2073  CALL popcontrol3b(branch)
2074  IF (branch .LT. 2) THEN
2075  IF (branch .EQ. 0) THEN
2076  afb = afb + tmparrb/(1.-af)**2
2077  ELSE
2078  afb = afb + tmparrb/(0.02)**2
2079  END IF
2080  ELSE IF (branch .EQ. 2) THEN
2081  afb = afb + tmparrb/(1.-af)**2
2082  END IF
2083  CALL dqsats_bac_b(dqsx, dqsxb, qsx, qsxb, teo, teob, pl, estblx, &
2084 & cons_h2omw, cons_airmw)
2085  teb = teb + teob
2086  qalb = qalb + qab
2087  qaib = qaib + qab
2088  qclb = qclb + qcb
2089  qcib = qcib + qcb
2090 END SUBROUTINE ls_cloud_b
2091 
2092 ! Differentiation of pdffrac in reverse (adjoint) mode:
2093 ! gradient of useful results: qtmean sigmaqt1 sigmaqt2 qstar
2094 ! clfrac
2095 ! with respect to varying inputs: qtmean sigmaqt1 sigmaqt2 qstar
2096 ! clfrac
2097 SUBROUTINE pdffrac_b(flag, qtmean, qtmeanb, sigmaqt1, sigmaqt1b, &
2098 & sigmaqt2, sigmaqt2b, qstar, qstarb, clfrac, clfracb)
2099  IMPLICIT NONE
2100 !Regularization
2101 !clfracd = 0.2*clfracd
2102 !Inputs
2103  INTEGER, INTENT(IN) :: flag
2104  REAL*8, INTENT(IN) :: qtmean, sigmaqt1, sigmaqt2, qstar
2105  REAL*8 :: qtmeanb, sigmaqt1b, sigmaqt2b, qstarb
2106 !Prognostic
2107  REAL*8, INTENT(INOUT) :: clfrac
2108  REAL*8 :: clfracb
2109 !LOCALS
2110  REAL*8 :: qtmode, qtmin, qtmax
2111  REAL*8 :: qtmodeb, qtminb, qtmaxb
2112  REAL*8 :: rh, rhd, q1, q2
2113  REAL*8 :: rhb
2114  INTRINSIC min
2115  INTRINSIC tanh
2116  INTEGER :: branch
2117  REAL*8 :: temp0
2118  REAL*8 :: min1
2119  REAL*8 :: tempb5
2120  REAL*8 :: tempb4
2121  REAL*8 :: tempb3
2122  REAL*8 :: tempb2
2123  REAL*8 :: tempb1
2124  REAL*8 :: tempb0
2125  REAL*8 :: min1b
2126  REAL*8 :: tempb
2127  REAL*8 :: temp
2128  IF (flag .EQ. 1) THEN
2129 !Tophat PDF
2130  IF (qtmean + sigmaqt1 .GE. qstar) THEN
2131  IF (sigmaqt1 .GT. 0.) THEN
2132  IF (qtmean + sigmaqt1 - qstar .GT. 2.*sigmaqt1) THEN
2133  min1 = 2.*sigmaqt1
2134  CALL pushcontrol1b(0)
2135  ELSE
2136  min1 = qtmean + sigmaqt1 - qstar
2137  CALL pushcontrol1b(1)
2138  END IF
2139  tempb = clfracb/(2.*sigmaqt1)
2140  min1b = tempb
2141  sigmaqt1b = sigmaqt1b - min1*tempb/sigmaqt1
2142  CALL popcontrol1b(branch)
2143  IF (branch .EQ. 0) THEN
2144  sigmaqt1b = sigmaqt1b + 2.*min1b
2145  ELSE
2146  qtmeanb = qtmeanb + min1b
2147  sigmaqt1b = sigmaqt1b + min1b
2148  qstarb = qstarb - min1b
2149  END IF
2150  END IF
2151  END IF
2152  clfracb = 0.0_8
2153  ELSE IF (flag .EQ. 2) THEN
2154 !Triangular PDF
2155  qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.
2156  IF (qtmode - sigmaqt1 .GT. 0.) THEN
2157  CALL pushcontrol1b(0)
2158  qtmin = 0.
2159  ELSE
2160  qtmin = qtmode - sigmaqt1
2161  CALL pushcontrol1b(1)
2162  END IF
2163  qtmax = qtmode + sigmaqt2
2164  IF (qtmax .LT. qstar) THEN
2165  clfracb = 0.0_8
2166  qtmaxb = 0.0_8
2167  qtminb = 0.0_8
2168  qtmodeb = 0.0_8
2169  ELSE IF (qtmode .LE. qstar .AND. qstar .LT. qtmax) THEN
2170  temp = (qtmax-qtmin)*(qtmax-qtmode)
2171  tempb0 = clfracb/temp
2172  tempb1 = 2*(qtmax-qstar)*tempb0
2173  tempb2 = -((qtmax-qstar)**2*tempb0/temp)
2174  qtmaxb = (2*qtmax-qtmin-qtmode)*tempb2 + tempb1
2175  qstarb = qstarb - tempb1
2176  qtminb = -((qtmax-qtmode)*tempb2)
2177  qtmodeb = -((qtmax-qtmin)*tempb2)
2178  clfracb = 0.0_8
2179  ELSE IF (qtmin .LE. qstar .AND. qstar .LT. qtmode) THEN
2180  temp0 = (qtmax-qtmin)*(qtmode-qtmin)
2181  tempb3 = -(clfracb/temp0)
2182  tempb4 = 2*(qstar-qtmin)*tempb3
2183  tempb5 = -((qstar-qtmin)**2*tempb3/temp0)
2184  qstarb = qstarb + tempb4
2185  qtminb = (2*qtmin-qtmax-qtmode)*tempb5 - tempb4
2186  qtmaxb = (qtmode-qtmin)*tempb5
2187  qtmodeb = (qtmax-qtmin)*tempb5
2188  clfracb = 0.0_8
2189  ELSE
2190  IF (qstar .LE. qtmin) clfracb = 0.0_8
2191  qtmaxb = 0.0_8
2192  qtminb = 0.0_8
2193  qtmodeb = 0.0_8
2194  END IF
2195  qtmodeb = qtmodeb + qtmaxb
2196  sigmaqt2b = sigmaqt2b + qtmaxb
2197  CALL popcontrol1b(branch)
2198  IF (branch .NE. 0) THEN
2199  qtmodeb = qtmodeb + qtminb
2200  sigmaqt1b = sigmaqt1b - qtminb
2201  END IF
2202  qtmeanb = qtmeanb + qtmodeb
2203  sigmaqt1b = sigmaqt1b + qtmodeb/3.
2204  sigmaqt2b = sigmaqt2b - qtmodeb/3.
2205  ELSE IF (flag .EQ. 3) THEN
2206 
2207  ! (REGULARIZATION) * (GRADIENT) * (PERTURBATION)
2208  clfracb = (0.66*( cosh(10*(rh-1.0))**-2)) * clfracb
2209 
2210  rh = qtmean/qstar
2211  q1 = 22.556
2212  rhb = (1.0-tanh(q1*(rh-1.0))**2)*0.5*q1*clfracb
2213  qtmeanb = qtmeanb + rhb/qstar
2214  qstarb = qstarb - qtmean*rhb/qstar**2
2215  clfracb = 0.0_8
2216 
2217  ELSE IF (flag .EQ. 4) THEN
2218 
2219  !Linear for the perturbation part
2220 
2221  !Regularization
2222  clfracb = clfracb * 0.2
2223 
2224  rh = qtmean/qstar
2225  q1 = 0.9335
2226  q2 = 1.0665
2227  IF (rh .LT. q1) THEN
2228  rhb = 0.0_8
2229  ELSE IF (rh .GE. q1 .AND. rh .LT. q2) THEN
2230  rhb = clfracb/((q2/q1-1)*q1)
2231  ELSE
2232  rhb = 0.0_8
2233  END IF
2234  qtmeanb = qtmeanb + rhb/qstar
2235  qstarb = qstarb - qtmean*rhb/qstar**2
2236  clfracb = 0.0_8
2237  END IF
2238 
2239 END SUBROUTINE pdffrac_b
2240 
2241 ! Differentiation of pdfcondensate in reverse (adjoint) mode:
2242 ! gradient of useful results: qstar4 condensate4
2243 ! with respect to varying inputs: qtmean4 qstar4 sigmaqt14 sigmaqt24
2244 SUBROUTINE pdfcondensate_b(flag, qtmean4, qtmean4b, sigmaqt14, &
2245 & sigmaqt14b, sigmaqt24, sigmaqt24b, qstar4, qstar4b, condensate4, &
2246 & condensate4b)
2247  IMPLICIT NONE
2248 !Inputs
2249  INTEGER, INTENT(IN) :: flag
2250  REAL*8, INTENT(IN) :: qtmean4, sigmaqt14, sigmaqt24, qstar4
2251  REAL*8 :: qtmean4b, sigmaqt14b, sigmaqt24b, qstar4b
2252 !Prognostic
2253  REAL*8, INTENT(INOUT) :: condensate4
2254  REAL*8 :: condensate4b
2255 !Locals
2256  REAL*8 :: qtmode, qtmin, qtmax, consta, constb, cloudf
2257  REAL*8 :: qtmodeb, qtminb, qtmaxb, constab, constbb, cloudfb
2258  REAL*8 :: term1, term2, term3
2259  REAL*8 :: term1b, term2b, term3b
2260  REAL*8 :: qtmean, sigmaqt1, sigmaqt2, qstar, condensate
2261  REAL*8 :: qtmeanb, sigmaqt1b, sigmaqt2b, qstarb, condensateb
2262  INTRINSIC min
2263  INTEGER :: branch
2264  REAL*8 :: temp2
2265  REAL*8 :: temp1
2266  REAL*8 :: temp0
2267  REAL*8 :: tempb9
2268  REAL*8 :: tempb8
2269  REAL*8 :: tempb7
2270  REAL*8 :: tempb6
2271  REAL*8 :: min1
2272  REAL*8 :: tempb5
2273  REAL*8 :: tempb4
2274  REAL*8 :: tempb3
2275  REAL*8 :: tempb2
2276  REAL*8 :: tempb1
2277  REAL*8 :: tempb0
2278  REAL*8 :: min1b
2279  REAL*8 :: tempb
2280  REAL*8 :: temp
2281  qtmean = qtmean4
2282  sigmaqt1 = sigmaqt14
2283  sigmaqt2 = sigmaqt24
2284  qstar = qstar4
2285  IF (flag .EQ. 1) THEN
2286  IF (qtmean + sigmaqt1 .LT. qstar) THEN
2287  CALL pushcontrol4b(0)
2288  ELSE IF (qstar .GT. qtmean - sigmaqt1) THEN
2289  IF (sigmaqt1 .GT. 0.d0) THEN
2290  IF (qtmean + sigmaqt1 - qstar .GT. 2.d0*sigmaqt1) THEN
2291  min1 = 2.d0*sigmaqt1
2292  CALL pushcontrol1b(0)
2293  ELSE
2294  min1 = qtmean + sigmaqt1 - qstar
2295  CALL pushcontrol1b(1)
2296  END IF
2297  CALL pushcontrol4b(1)
2298  ELSE
2299  CALL pushcontrol4b(2)
2300  END IF
2301  ELSE
2302  CALL pushcontrol4b(3)
2303  END IF
2304  ELSE IF (flag .EQ. 2) THEN
2305  qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.d0
2306  IF (qtmode - sigmaqt1 .GT. 0.d0) THEN
2307  qtmin = 0.d0
2308  CALL pushcontrol1b(0)
2309  ELSE
2310  qtmin = qtmode - sigmaqt1
2311  CALL pushcontrol1b(1)
2312  END IF
2313  qtmax = qtmode + sigmaqt2
2314  IF (qtmax .LT. qstar) THEN
2315  CALL pushcontrol4b(4)
2316  ELSE IF (qtmode .LE. qstar .AND. qstar .LT. qtmax) THEN
2317  constb = 2.d0/((qtmax-qtmin)*(qtmax-qtmode))
2318  cloudf = (qtmax-qstar)*(qtmax-qstar)/((qtmax-qtmin)*(qtmax-qtmode)&
2319 & )
2320  term1 = qstar*qstar*qstar/3.d0
2321  term2 = qtmax*qstar*qstar/2.d0
2322  term3 = qtmax*qtmax*qtmax/6.d0
2323  CALL pushcontrol4b(5)
2324  ELSE IF (qtmin .LE. qstar .AND. qstar .LT. qtmode) THEN
2325  consta = 2.d0/((qtmax-qtmin)*(qtmode-qtmin))
2326  cloudf = 1.d0 - (qstar-qtmin)*(qstar-qtmin)/((qtmax-qtmin)*(qtmode&
2327 & -qtmin))
2328  term1 = qstar*qstar*qstar/3.d0
2329  term2 = qtmin*qstar*qstar/2.d0
2330  term3 = qtmin*qtmin*qtmin/6.d0
2331  CALL pushcontrol4b(6)
2332  ELSE IF (qstar .LE. qtmin) THEN
2333  CALL pushcontrol4b(7)
2334  ELSE
2335  CALL pushcontrol4b(8)
2336  END IF
2337  ELSE IF (flag .EQ. 3) THEN
2338 !Reference part done normally
2339 !IF (qtmean + sigmaqt1 .LT. qstar) THEN
2340 ! condensate = 0.d0
2341 !ELSE IF (qstar .GT. qtmean - sigmaqt1) THEN
2342 ! IF (sigmaqt1 .GT. 0.d0) THEN
2343 ! IF (qtmean + sigmaqt1 - qstar .GT. 2.d0*sigmaqt1) THEN
2344 ! min1 = 2.d0*sigmaqt1
2345 ! ELSE
2346 ! min1 = qtmean + sigmaqt1 - qstar
2347 ! END IF
2348 ! condensate = min1**2/(4.d0*sigmaqt1)
2349 ! ELSE
2350 ! condensate = qtmean - qstar
2351 ! END IF
2352 !ELSE
2353 ! condensate = qtmean - qstar
2354 !END IF
2355 !Perturbation part from linear
2356  IF (qtmean - qstar .GT. -0.5e-3) THEN
2357  CALL pushcontrol4b(9)
2358  ELSE
2359  CALL pushcontrol4b(10)
2360  END IF
2361  ELSE
2362  CALL pushcontrol4b(11)
2363  END IF
2364  condensateb = condensate4b
2365  CALL popcontrol4b(branch)
2366  IF (branch .LT. 6) THEN
2367  IF (branch .LT. 3) THEN
2368  IF (branch .EQ. 0) THEN
2369  qtmeanb = 0.0_8
2370  sigmaqt1b = 0.0_8
2371  qstarb = 0.0_8
2372  ELSE IF (branch .EQ. 1) THEN
2373  tempb = condensateb/(4.d0*sigmaqt1)
2374  min1b = 2*min1*tempb
2375  sigmaqt1b = -(min1**2*tempb/sigmaqt1)
2376  CALL popcontrol1b(branch)
2377  IF (branch .EQ. 0) THEN
2378  sigmaqt1b = sigmaqt1b + 2.d0*min1b
2379  qtmeanb = 0.0_8
2380  qstarb = 0.0_8
2381  ELSE
2382  qtmeanb = min1b
2383  sigmaqt1b = sigmaqt1b + min1b
2384  qstarb = -min1b
2385  END IF
2386  ELSE
2387  qtmeanb = condensateb
2388  qstarb = -condensateb
2389  sigmaqt1b = 0.0_8
2390  END IF
2391  ELSE IF (branch .EQ. 3) THEN
2392  qtmeanb = condensateb
2393  qstarb = -condensateb
2394  sigmaqt1b = 0.0_8
2395  ELSE
2396  IF (branch .EQ. 4) THEN
2397  qtmeanb = 0.0_8
2398  qtmaxb = 0.0_8
2399  qstarb = 0.0_8
2400  qtminb = 0.0_8
2401  qtmodeb = 0.0_8
2402  ELSE
2403  cloudfb = -(qstar*condensateb)
2404  temp0 = (qtmax-qtmin)*(qtmax-qtmode)
2405  tempb4 = cloudfb/temp0
2406  tempb1 = 2*(qtmax-qstar)*tempb4
2407  tempb0 = constb*condensateb
2408  constbb = (term1-term2+term3)*condensateb
2409  term1b = tempb0
2410  term2b = -tempb0
2411  term3b = tempb0
2412  qstarb = qtmax*2*qstar*term2b/2.d0 - tempb1 + 3*qstar**2*term1b/&
2413 & 3.d0 - cloudf*condensateb
2414  tempb3 = -((qtmax-qstar)**2*tempb4/temp0)
2415  temp = (qtmax-qtmin)*(qtmax-qtmode)
2416  tempb2 = -(2.d0*constbb/temp**2)
2417  qtmaxb = qstar**2*term2b/2.d0 + (2*qtmax-qtmin-qtmode)*tempb2 + &
2418 & (2*qtmax-qtmin-qtmode)*tempb3 + tempb1 + 3*qtmax**2*term3b/&
2419 & 6.d0
2420  qtminb = -((qtmax-qtmode)*tempb2) - (qtmax-qtmode)*tempb3
2421  qtmodeb = -((qtmax-qtmin)*tempb2) - (qtmax-qtmin)*tempb3
2422  qtmeanb = 0.0_8
2423  END IF
2424  GOTO 100
2425  END IF
2426  sigmaqt2b = 0.0_8
2427  GOTO 110
2428  ELSE IF (branch .LT. 9) THEN
2429  IF (branch .EQ. 6) THEN
2430  cloudfb = -(qstar*condensateb)
2431  temp2 = (qtmax-qtmin)*(qtmode-qtmin)
2432  tempb9 = -(cloudfb/temp2)
2433  tempb6 = 2*(qstar-qtmin)*tempb9
2434  tempb5 = -(consta*condensateb)
2435  qtmeanb = condensateb
2436  constab = -((term1-term2+term3)*condensateb)
2437  term1b = tempb5
2438  term2b = -tempb5
2439  term3b = tempb5
2440  qstarb = qtmin*2*qstar*term2b/2.d0 + tempb6 + 3*qstar**2*term1b/&
2441 & 3.d0 - cloudf*condensateb
2442  tempb8 = -((qstar-qtmin)**2*tempb9/temp2)
2443  temp1 = (qtmax-qtmin)*(qtmode-qtmin)
2444  tempb7 = -(2.d0*constab/temp1**2)
2445  qtminb = qstar**2*term2b/2.d0 + (2*qtmin-qtmax-qtmode)*tempb7 + (2&
2446 & *qtmin-qtmax-qtmode)*tempb8 - tempb6 + 3*qtmin**2*term3b/6.d0
2447  qtmaxb = (qtmode-qtmin)*tempb7 + (qtmode-qtmin)*tempb8
2448  qtmodeb = (qtmax-qtmin)*tempb7 + (qtmax-qtmin)*tempb8
2449  ELSE
2450  IF (branch .EQ. 7) THEN
2451  qtmeanb = condensateb
2452  qstarb = -condensateb
2453  ELSE
2454  qtmeanb = 0.0_8
2455  qstarb = 0.0_8
2456  END IF
2457  qtmaxb = 0.0_8
2458  qtminb = 0.0_8
2459  qtmodeb = 0.0_8
2460  END IF
2461  ELSE
2462  IF (branch .EQ. 9) THEN
2463  qtmeanb = condensateb
2464  qstarb = -condensateb
2465  ELSE IF (branch .EQ. 10) THEN
2466  qtmeanb = 0.0_8
2467  qstarb = 0.0_8
2468  ELSE
2469  qtmeanb = 0.0_8
2470  qstarb = 0.0_8
2471  END IF
2472  sigmaqt1b = 0.0_8
2473  sigmaqt2b = 0.0_8
2474  GOTO 110
2475  END IF
2476  100 qtmodeb = qtmodeb + qtmaxb
2477  sigmaqt2b = qtmaxb
2478  CALL popcontrol1b(branch)
2479  IF (branch .EQ. 0) THEN
2480  sigmaqt1b = 0.0_8
2481  ELSE
2482  qtmodeb = qtmodeb + qtminb
2483  sigmaqt1b = -qtminb
2484  END IF
2485  qtmeanb = qtmeanb + qtmodeb
2486  sigmaqt1b = sigmaqt1b + qtmodeb/3.d0
2487  sigmaqt2b = sigmaqt2b - qtmodeb/3.d0
2488  110 qstar4b = qstar4b + qstarb
2489  sigmaqt24b = sigmaqt2b
2490  sigmaqt14b = sigmaqt1b
2491  qtmean4b = qtmeanb
2492 END SUBROUTINE pdfcondensate_b
2493 
2494 ! Differentiation of evap_cnv in reverse (adjoint) mode:
2495 ! gradient of useful results: f qi ql qs qv te
2496 ! with respect to varying inputs: f qi ql qs qv te
2497 SUBROUTINE evap_cnv_b(dt, rhcr, pl, te, teb, qv, qvb, ql, qlb, qi, qib, &
2498 & f, fb, xf, qs, qsb, rho_w, cld_evp_eff, cons_h2omw, cons_airmw, &
2499 & cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp)
2500  IMPLICIT NONE
2501 !Inputs
2502  REAL*8, INTENT(IN) :: dt, rhcr, pl, xf, qs, rho_w, cld_evp_eff
2503  REAL*8 :: qsb
2504  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_alhl, cons_rvap, &
2505 & cons_rgas, cons_pi, cons_cp
2506 !Prognostics
2507  REAL*8, INTENT(INOUT) :: te, qv, ql, qi, f
2508  REAL*8 :: teb, qvb, qlb, qib, fb
2509 !Locals
2510  REAL*8 :: es, radius, k1, k2, teff, qcm, evap, rhx, qc, a_eff, epsilon
2511  REAL*8 :: esb, radiusb, k1b, k2b, teffb, qcmb, evapb, rhxb, qcb
2512  REAL*8, PARAMETER :: k_cond=2.4e-2
2513  REAL*8, PARAMETER :: diffu=2.2e-5
2514  REAL*8, PARAMETER :: nn=50.*1.0e6
2515  INTRINSIC min
2516  INTEGER :: branch
2517  REAL*8 :: temp3
2518  REAL*8 :: temp2
2519  REAL*8 :: temp1
2520  REAL*8 :: temp0
2521  REAL*8 :: tempb4
2522  REAL*8 :: tempb3
2523  REAL*8 :: tempb2
2524  REAL*8 :: tempb1
2525  REAL*8 :: tempb0
2526  REAL*8 :: tempb
2527  REAL*8 :: temp
2528  epsilon = cons_h2omw/cons_airmw
2529  a_eff = cld_evp_eff
2530 !EVAPORATION OF CLOUD WATER.
2531 ! (100 <-^ convert from mbar to Pa)
2532  es = 100.*pl*qs/(epsilon+(1.0-epsilon)*qs)
2533  IF (qv/qs .GT. 1.00) THEN
2534  CALL pushcontrol1b(0)
2535  rhx = 1.00
2536  ELSE
2537  rhx = qv/qs
2538  CALL pushcontrol1b(1)
2539  END IF
2540  k1 = cons_alhl**2*rho_w/(k_cond*cons_rvap*te**2)
2541  k2 = cons_rvap*te*rho_w/(diffu*(1000./pl)*es)
2542 !Here DIFFU is given for 1000 mb so 1000./PR accounts for increased diffusivity at lower pressure.
2543  IF (f .GT. 0. .AND. ql .GT. 0.) THEN
2544  qcm = ql/f
2545  CALL pushcontrol1b(0)
2546  ELSE
2547  CALL pushcontrol1b(1)
2548  qcm = 0.
2549  END IF
2550  CALL ldradius(pl, te, qcm, nn, rho_w, radius, cons_rgas, cons_pi)
2551  IF (rhx .LT. rhcr .AND. radius .GT. 0.0) THEN
2552 ! / (1.00 - RHx)
2553  teff = (rhcr-rhx)/((k1+k2)*radius**2)
2554  CALL pushcontrol1b(1)
2555  ELSE
2556 ! -999.
2557  teff = 0.0
2558  CALL pushcontrol1b(0)
2559  END IF
2560  evap = a_eff*ql*dt*teff
2561  IF (evap .GT. ql) THEN
2562  evap = ql
2563  CALL pushcontrol1b(0)
2564  ELSE
2565  CALL pushcontrol1b(1)
2566  evap = evap
2567  END IF
2568  qc = ql + qi
2569  IF (qc .GT. 0.) THEN
2570  CALL pushcontrol1b(0)
2571  ELSE
2572  CALL pushcontrol1b(1)
2573  END IF
2574  evapb = qvb - qlb - cons_alhl*teb/cons_cp
2575  CALL popcontrol1b(branch)
2576  IF (branch .EQ. 0) THEN
2577  temp3 = f/qc
2578  tempb4 = (qc-evap)*fb/qc
2579  qcb = temp3*fb - temp3*tempb4
2580  evapb = evapb - temp3*fb
2581  fb = tempb4
2582  ELSE
2583  qcb = 0.0_8
2584  END IF
2585  qlb = qlb + qcb
2586  qib = qib + qcb
2587  CALL popcontrol1b(branch)
2588  IF (branch .EQ. 0) THEN
2589  qlb = qlb + evapb
2590  evapb = 0.0_8
2591  END IF
2592  tempb3 = a_eff*dt*evapb
2593  qlb = qlb + teff*tempb3
2594  teffb = ql*tempb3
2595  CALL popcontrol1b(branch)
2596  IF (branch .EQ. 0) THEN
2597  radiusb = 0.0_8
2598  k1b = 0.0_8
2599  k2b = 0.0_8
2600  rhxb = 0.0_8
2601  ELSE
2602  temp2 = (k1+k2)*radius**2
2603  tempb1 = -((rhcr-rhx)*teffb/temp2**2)
2604  tempb2 = radius**2*tempb1
2605  rhxb = -(teffb/temp2)
2606  k1b = tempb2
2607  k2b = tempb2
2608  radiusb = (k1+k2)*2*radius*tempb1
2609  END IF
2610  CALL ldradius_b(pl, te, teb, qcm, qcmb, nn, rho_w, radius, radiusb, &
2611 & cons_rgas, cons_pi)
2612 qcmb = 0.0_8
2613  CALL popcontrol1b(branch)
2614  IF (branch .EQ. 0) THEN
2615  qlb = qlb + qcmb/f
2616  fb = fb - ql*qcmb/f**2
2617  END IF
2618  temp0 = k_cond*cons_rvap*te**2
2619  temp1 = 1000.*diffu*es
2620  tempb0 = cons_rvap*rho_w*pl*k2b/temp1
2621  teb = teb + tempb0 - k_cond*cons_rvap*cons_alhl**2*rho_w*2*te*k1b/&
2622 & temp0**2
2623  esb = -(te*diffu*1000.*tempb0/temp1)
2624  CALL popcontrol1b(branch)
2625  IF (branch .NE. 0) THEN
2626  qvb = qvb + rhxb/qs
2627  qsb = qsb - qv*rhxb/qs**2
2628  END IF
2629  temp = epsilon + (-epsilon+1.0)*qs
2630  tempb = pl*100.*esb/temp
2631  qsb = qsb + (1.0_8-qs*(1.0-epsilon)/temp)*tempb
2632 END SUBROUTINE evap_cnv_b
2633 
2634 ! Differentiation of subl_cnv in reverse (adjoint) mode:
2635 ! gradient of useful results: f qi ql qs qv te
2636 ! with respect to varying inputs: f qi ql qs qv te
2637 SUBROUTINE subl_cnv_b(dt, rhcr, pl, te, teb, qv, qvb, ql, qlb, qi, qib, &
2638 & f, fb, xf, qs, qsb, rho_w, cld_evp_eff, cons_h2omw, cons_airmw, &
2639 & cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp, cons_alhs)
2640  IMPLICIT NONE
2641 !INPUTS
2642  REAL*8, INTENT(IN) :: dt, rhcr, pl, xf, qs, rho_w, cld_evp_eff
2643  REAL*8 :: qsb
2644  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_alhl, cons_rvap, &
2645 & cons_rgas, cons_pi, cons_cp, cons_alhs
2646 !PROGNOSTIC
2647  REAL*8, INTENT(INOUT) :: te, qv, ql, qi, f
2648  REAL*8 :: teb, qvb, qlb, qib, fb
2649 !LOCALS
2650  REAL*8 :: es, radius, k1, k2, teff, qcm, subl, rhx, qc, a_eff, nn, &
2651 & epsilon
2652  REAL*8 :: esb, radiusb, k1b, k2b, teffb, qcmb, sublb, rhxb, qcb
2653  REAL*8, PARAMETER :: k_cond=2.4e-2
2654  REAL*8, PARAMETER :: diffu=2.2e-5
2655  INTRINSIC min
2656  INTEGER :: branch
2657  REAL*8 :: temp3
2658  REAL*8 :: temp2
2659  REAL*8 :: temp1
2660  REAL*8 :: temp0
2661  REAL*8 :: tempb4
2662  REAL*8 :: tempb3
2663  REAL*8 :: tempb2
2664  REAL*8 :: tempb1
2665  REAL*8 :: tempb0
2666  REAL*8 :: tempb
2667  REAL*8 :: temp
2668  epsilon = cons_h2omw/cons_airmw
2669  a_eff = cld_evp_eff
2670  nn = 5.*1.0e6
2671 ! (100 s <-^ convert from mbar to Pa)
2672  es = 100.*pl*qs/(epsilon+(1.0-epsilon)*qs)
2673  IF (qv/qs .GT. 1.00) THEN
2674  CALL pushcontrol1b(0)
2675  rhx = 1.00
2676  ELSE
2677  rhx = qv/qs
2678  CALL pushcontrol1b(1)
2679  END IF
2680  k1 = cons_alhl**2*rho_w/(k_cond*cons_rvap*te**2)
2681  k2 = cons_rvap*te*rho_w/(diffu*(1000./pl)*es)
2682 !Here DIFFU is given for 1000 mb so 1000./PR accounts for increased diffusivity at lower pressure.
2683  IF (f .GT. 0. .AND. qi .GT. 0.) THEN
2684  qcm = qi/f
2685  CALL pushcontrol1b(0)
2686  ELSE
2687  CALL pushcontrol1b(1)
2688  qcm = 0.
2689  END IF
2690  CALL ldradius(pl, te, qcm, nn, rho_w, radius, cons_rgas, cons_pi)
2691  IF (rhx .LT. rhcr .AND. radius .GT. 0.0) THEN
2692 ! / (1.00 - RHx)
2693  teff = (rhcr-rhx)/((k1+k2)*radius**2)
2694  CALL pushcontrol1b(1)
2695  ELSE
2696 ! -999.
2697  teff = 0.0
2698  CALL pushcontrol1b(0)
2699  END IF
2700  subl = a_eff*qi*dt*teff
2701  IF (subl .GT. qi) THEN
2702  subl = qi
2703  CALL pushcontrol1b(0)
2704  ELSE
2705  CALL pushcontrol1b(1)
2706  subl = subl
2707  END IF
2708  qc = ql + qi
2709  IF (qc .GT. 0.) THEN
2710  CALL pushcontrol1b(0)
2711  ELSE
2712  CALL pushcontrol1b(1)
2713  END IF
2714  sublb = qvb - qib - cons_alhs*teb/cons_cp
2715  CALL popcontrol1b(branch)
2716  IF (branch .EQ. 0) THEN
2717  temp3 = f/qc
2718  tempb4 = (qc-subl)*fb/qc
2719  qcb = temp3*fb - temp3*tempb4
2720  sublb = sublb - temp3*fb
2721  fb = tempb4
2722  ELSE
2723  qcb = 0.0_8
2724  END IF
2725  qlb = qlb + qcb
2726  qib = qib + qcb
2727  CALL popcontrol1b(branch)
2728  IF (branch .EQ. 0) THEN
2729  qib = qib + sublb
2730  sublb = 0.0_8
2731  END IF
2732  tempb3 = a_eff*dt*sublb
2733  qib = qib + teff*tempb3
2734  teffb = qi*tempb3
2735  CALL popcontrol1b(branch)
2736  IF (branch .EQ. 0) THEN
2737  radiusb = 0.0_8
2738  k1b = 0.0_8
2739  k2b = 0.0_8
2740  rhxb = 0.0_8
2741  ELSE
2742  temp2 = (k1+k2)*radius**2
2743  tempb1 = -((rhcr-rhx)*teffb/temp2**2)
2744  tempb2 = radius**2*tempb1
2745  rhxb = -(teffb/temp2)
2746  k1b = tempb2
2747  k2b = tempb2
2748  radiusb = (k1+k2)*2*radius*tempb1
2749  END IF
2750  CALL ldradius_b(pl, te, teb, qcm, qcmb, nn, rho_w, radius, radiusb, &
2751 & cons_rgas, cons_pi)
2752 qcmb = 0.0
2753  CALL popcontrol1b(branch)
2754  IF (branch .EQ. 0) THEN
2755  qib = qib + qcmb/f
2756  fb = fb - qi*qcmb/f**2
2757  END IF
2758  temp0 = k_cond*cons_rvap*te**2
2759  temp1 = 1000.*diffu*es
2760  tempb0 = cons_rvap*rho_w*pl*k2b/temp1
2761  teb = teb + tempb0 - k_cond*cons_rvap*cons_alhl**2*rho_w*2*te*k1b/&
2762 & temp0**2
2763  esb = -(te*diffu*1000.*tempb0/temp1)
2764  CALL popcontrol1b(branch)
2765  IF (branch .NE. 0) THEN
2766  qvb = qvb + rhxb/qs
2767  qsb = qsb - qv*rhxb/qs**2
2768  END IF
2769  temp = epsilon + (-epsilon+1.0)*qs
2770  tempb = pl*100.*esb/temp
2771  qsb = qsb + (1.0_8-qs*(1.0-epsilon)/temp)*tempb
2772 END SUBROUTINE subl_cnv_b
2773 
2774 ! Differentiation of ldradius in reverse (adjoint) mode:
2775 ! gradient of useful results: radius te
2776 ! with respect to varying inputs: qcl te
2777 SUBROUTINE ldradius_b(pl, te, teb, qcl, qclb, nn, rho_w, radius, radiusb&
2778 & , cons_rgas, cons_pi)
2779  IMPLICIT NONE
2780 !Inputs
2781  REAL*8, INTENT(IN) :: te, pl, nn, qcl, rho_w
2782  REAL*8 :: teb, qclb
2783  REAL*8, INTENT(IN) :: cons_rgas, cons_pi
2784 !Outputs
2785  REAL*8 :: radius
2786  REAL*8 :: radiusb
2787 !Equiv. Spherical Cloud Particle Radius in m
2788  REAL*8 :: temp2
2789  REAL*8 :: temp1
2790  REAL*8 :: temp0
2791  REAL*8 :: tempb
2792  REAL*8 :: temp
2793  temp2 = 4.*cons_rgas*cons_pi*nn*rho_w
2794  temp1 = temp2*te
2795  temp = qcl/temp1
2796  temp0 = 3.*pl*100.
2797  IF (temp0*temp .LE. 0.0 .AND. (1.0/3. .EQ. 0.0 .OR. 1.0/3. .NE. int(&
2798 & 1.0/3.))) THEN
2799  tempb = 0.0
2800  ELSE
2801  tempb = temp0*(temp0*temp)**(1.0/3.-1)*radiusb/(3.*temp1)
2802  END IF
2803  qclb = tempb
2804  teb = teb - temp*temp2*tempb
2805 END SUBROUTINE ldradius_b
2806 
2807 ! Differentiation of autoconversion_ls in reverse (adjoint) mode:
2808 ! gradient of useful results: f qc qp te
2809 ! with respect to varying inputs: f qc te
2810 SUBROUTINE autoconversion_ls_b(dt, qc, qcb, qp, qpb, te, teb, pl, f, fb&
2811 & , sundqv2, sundqv3, sundqt1, c_00, lwcrit, dzet)
2812  IMPLICIT NONE
2813 !Inputs
2814  REAL*8, INTENT(IN) :: dt, te, pl, dzet, sundqv2, sundqv3, sundqt1, &
2815 & c_00, lwcrit
2816  REAL*8 :: teb
2817 !Prognostic
2818  REAL*8, INTENT(INOUT) :: qc, qp, f
2819  REAL*8 :: qcb, qpb, fb
2820 !Locals
2821  REAL*8 :: acf0, acf, c00x, iqccrx, f2, f3, rate, dqp, qcm, dqfac
2822  REAL*8 :: c00xb, iqccrxb, f2b, f3b, rateb, dqpb, qcmb, dqfacb
2823  INTRINSIC exp
2824  INTRINSIC min
2825  INTRINSIC max
2826  INTEGER :: branch
2827  REAL*8 :: tempb1
2828  REAL*8 :: tempb0
2829  REAL*8 :: x4
2830  REAL*8 :: x3
2831  REAL*8 :: x2
2832  REAL*8 :: x1
2833  REAL*8 :: x2b
2834  REAL*8 :: tempb
2835  REAL*8 :: x1b
2836  REAL*8 :: x4b
2837  REAL*8 :: temp
2838 !Zero Locals
2839  f2 = 0.0
2840  f3 = 0.0
2841  CALL pushreal8(f2)
2842  CALL cons_sundq3(te, sundqv2, sundqv3, sundqt1, f2, f3)
2843  c00x = c_00*f2*f3
2844  iqccrx = f2*f3/lwcrit
2845  IF (f .GT. 0. .AND. qc .GT. 0.) THEN
2846  qcm = qc/f
2847  CALL pushcontrol1b(0)
2848  ELSE
2849  CALL pushcontrol1b(1)
2850  qcm = 0.
2851  END IF
2852  rate = c00x*(1.0-exp(-((qcm*iqccrx)**2)))
2853 !Temporary kluge until we can figure a better to make thicker low clouds.
2854  CALL pushreal8(f3)
2855  f3 = 1.0
2856 !Implement ramps for gradual change in autoconv
2857 !Thicken low high lat clouds
2858  IF (pl .GE. 775. .AND. te .LE. 275.) f3 = 0.2
2859 !F3 = max(-0.016 * PL + 13.4, 0.2)
2860  IF (pl .GE. 825. .AND. te .LE. 282.) f3 = 0.2
2861 !F3 = max(0.11 * TE - 30.02, 0.2)
2862  IF (pl .GE. 775. .AND. pl .LT. 825. .AND. te .LE. 282. .AND. te .GT. &
2863 & 275.) f3 = 0.2
2864 !F3 = min(max(-0.016*PL + 0.11 * TE - 16.85, 0.2),1.)
2865  IF (pl .GE. 825. .AND. te .LE. 275.) f3 = 0.2
2866  IF (pl .LE. 775. .OR. te .GT. 282.) f3 = 1.
2867 !Thin-out low tropical clouds
2868  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
2869  IF (0.2*te - 56 .GT. 2.) THEN
2870  CALL pushcontrol2b(2)
2871  f3 = 2.
2872  ELSE
2873  f3 = 0.2*te - 56
2874  CALL pushcontrol2b(1)
2875  END IF
2876  ELSE
2877  CALL pushcontrol2b(0)
2878  END IF
2879  IF (pl .GE. 925. .AND. te .GE. 290.) THEN
2880  IF (0.04*pl - 36. .GT. 2.) THEN
2881  CALL pushcontrol1b(1)
2882  f3 = 2.
2883  ELSE
2884  CALL pushcontrol1b(1)
2885  f3 = 0.04*pl - 36.
2886  END IF
2887  ELSE
2888  CALL pushcontrol1b(0)
2889  END IF
2890  IF (pl .GE. 925. .AND. pl .LT. 950. .AND. te .GT. 285. .AND. te .LT. &
2891 & 290.) THEN
2892  IF (0.04*pl + 0.2*te - 94. .GT. 2.) THEN
2893  CALL pushcontrol1b(0)
2894  x1 = 2.
2895  ELSE
2896  x1 = 0.04*pl + 0.2*te - 94.
2897  CALL pushcontrol1b(1)
2898  END IF
2899  IF (x1 .LT. 1.) THEN
2900  f3 = 1.
2901  CALL pushcontrol2b(2)
2902  ELSE
2903  f3 = x1
2904  CALL pushcontrol2b(1)
2905  END IF
2906  ELSE
2907  CALL pushcontrol2b(0)
2908  END IF
2909  IF (pl .GE. 950. .AND. te .GE. 290.) THEN
2910  f3 = 2.
2911  CALL pushcontrol1b(1)
2912  ELSE
2913  CALL pushcontrol1b(0)
2914  END IF
2915  IF (f3 .LT. 0.1) THEN
2916  f3 = 0.1
2917  CALL pushcontrol1b(0)
2918  ELSE
2919  CALL pushcontrol1b(1)
2920  f3 = f3
2921  END IF
2922  CALL pushreal8(rate)
2923  rate = f3*rate
2924  dqp = qc*(1.0-exp(-(rate*dt)))
2925  IF (dqp .LT. 0.0) THEN
2926  dqp = 0.0
2927  CALL pushcontrol1b(0)
2928  ELSE
2929  CALL pushcontrol1b(1)
2930  dqp = dqp
2931  END IF
2932 !Wipe-out warm fogs
2933  dqfac = 0.
2934  IF (pl .GE. 975. .AND. te .GE. 280.) THEN
2935  IF (0.2*te - 56. .GT. 1.) THEN
2936  CALL pushcontrol1b(0)
2937  x2 = 1.
2938  ELSE
2939  x2 = 0.2*te - 56.
2940  CALL pushcontrol1b(1)
2941  END IF
2942  IF (x2 .LT. 0.) THEN
2943  dqfac = 0.
2944  CALL pushcontrol2b(2)
2945  ELSE
2946  dqfac = x2
2947  CALL pushcontrol2b(1)
2948  END IF
2949  ELSE
2950  CALL pushcontrol2b(0)
2951  END IF
2952  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
2953  IF (0.04*pl - 38. .GT. 1.) THEN
2954  CALL pushcontrol1b(1)
2955  x3 = 1.
2956  ELSE
2957  CALL pushcontrol1b(1)
2958  x3 = 0.04*pl - 38.
2959  END IF
2960  IF (x3 .LT. 0.) THEN
2961  dqfac = 0.
2962  ELSE
2963  dqfac = x3
2964  END IF
2965  ELSE
2966  CALL pushcontrol1b(0)
2967  END IF
2968  IF (pl .GE. 950. .AND. pl .LT. 975. .AND. te .GT. 280. .AND. te .LT. &
2969 & 285.) THEN
2970  IF (0.04*pl + 0.2*te - 95. .GT. 1.) THEN
2971  CALL pushcontrol1b(0)
2972  x4 = 1.
2973  ELSE
2974  x4 = 0.04*pl + 0.2*te - 95.
2975  CALL pushcontrol1b(1)
2976  END IF
2977  IF (x4 .LT. 0.) THEN
2978  dqfac = 0.
2979  CALL pushcontrol2b(2)
2980  ELSE
2981  dqfac = x4
2982  CALL pushcontrol2b(1)
2983  END IF
2984  ELSE
2985  CALL pushcontrol2b(0)
2986  END IF
2987  IF (pl .GE. 975. .AND. te .GE. 285.) THEN
2988  dqfac = 1.
2989  CALL pushcontrol1b(1)
2990  ELSE
2991  CALL pushcontrol1b(0)
2992  END IF
2993  IF (dqp .LT. dqfac*qc) THEN
2994  dqp = dqfac*qc
2995  CALL pushcontrol1b(0)
2996  ELSE
2997  dqp = dqp
2998  CALL pushcontrol1b(1)
2999  END IF
3000  CALL pushreal8(qc)
3001  qc = qc - dqp
3002 !IF LARGE SCALE THEN
3003  IF (qc + dqp .GT. 0.) THEN
3004  tempb0 = fb/(qc+dqp)
3005  tempb1 = -(qc*f*tempb0/(qc+dqp))
3006  qcb = qcb + tempb1 + f*tempb0
3007  dqpb = tempb1
3008  fb = qc*tempb0
3009  ELSE
3010  dqpb = 0.0_8
3011  END IF
3012  dqpb = dqpb + qpb - qcb
3013  CALL popreal8(qc)
3014  CALL popcontrol1b(branch)
3015  IF (branch .EQ. 0) THEN
3016  dqfacb = qc*dqpb
3017  qcb = qcb + dqfac*dqpb
3018  dqpb = 0.0_8
3019  ELSE
3020  dqfacb = 0.0_8
3021  END IF
3022  CALL popcontrol1b(branch)
3023  IF (branch .NE. 0) dqfacb = 0.0_8
3024  CALL popcontrol2b(branch)
3025  IF (branch .NE. 0) THEN
3026  IF (branch .EQ. 1) THEN
3027  x4b = dqfacb
3028  ELSE
3029  x4b = 0.0_8
3030  END IF
3031  CALL popcontrol1b(branch)
3032  IF (branch .NE. 0) teb = teb + 0.2*x4b
3033  dqfacb = 0.0_8
3034  END IF
3035  CALL popcontrol1b(branch)
3036  IF (branch .NE. 0) dqfacb = 0.0_8
3037  CALL popcontrol2b(branch)
3038  IF (branch .NE. 0) THEN
3039  IF (branch .EQ. 1) THEN
3040  x2b = dqfacb
3041  ELSE
3042  x2b = 0.0_8
3043  END IF
3044  CALL popcontrol1b(branch)
3045  IF (branch .NE. 0) teb = teb + 0.2*x2b
3046  END IF
3047  CALL popcontrol1b(branch)
3048  IF (branch .EQ. 0) dqpb = 0.0_8
3049  qcb = qcb + (1.0-exp(-(dt*rate)))*dqpb
3050  rateb = exp(-(dt*rate))*qc*dt*dqpb
3051  CALL popreal8(rate)
3052  f3b = rate*rateb
3053  rateb = f3*rateb
3054  CALL popcontrol1b(branch)
3055  IF (branch .EQ. 0) f3b = 0.0_8
3056  CALL popcontrol1b(branch)
3057  IF (branch .NE. 0) f3b = 0.0_8
3058  CALL popcontrol2b(branch)
3059  IF (branch .NE. 0) THEN
3060  IF (branch .EQ. 1) THEN
3061  x1b = f3b
3062  ELSE
3063  x1b = 0.0_8
3064  END IF
3065  CALL popcontrol1b(branch)
3066  IF (branch .NE. 0) teb = teb + 0.2*x1b
3067  f3b = 0.0_8
3068  END IF
3069  CALL popcontrol1b(branch)
3070  IF (branch .NE. 0) f3b = 0.0_8
3071  CALL popcontrol2b(branch)
3072  IF (branch .NE. 0) THEN
3073  IF (branch .EQ. 1) teb = teb + 0.2*f3b
3074  END IF
3075  CALL popreal8(f3)
3076  temp = qcm*iqccrx
3077  tempb = exp(-(temp**2))*c00x*2*temp*rateb
3078  c00xb = (1.0-exp(-(temp**2)))*rateb
3079  qcmb = iqccrx*tempb
3080  iqccrxb = qcm*tempb
3081  CALL popcontrol1b(branch)
3082  IF (branch .EQ. 0) THEN
3083  qcb = qcb + qcmb/f
3084  fb = fb - qc*qcmb/f**2
3085  END IF
3086  f2b = c_00*f3*c00xb + f3*iqccrxb/lwcrit
3087  f2b = 0.5*f2b
3088  CALL popreal8(f2)
3089  CALL cons_sundq3_b(te, teb, sundqv2, sundqv3, sundqt1, f2, f2b, f3)
3090 END SUBROUTINE autoconversion_ls_b
3091 
3092 ! Differentiation of autoconversion_cnv in reverse (adjoint) mode:
3093 ! gradient of useful results: f qc qp te
3094 ! with respect to varying inputs: f qc te
3095 SUBROUTINE autoconversion_cnv_b(dt, qc, qcb, qp, qpb, te, teb, pl, f, fb&
3096 & , sundqv2, sundqv3, sundqt1, c_00, lwcrit, dzet)
3097  IMPLICIT NONE
3098 !Inputs
3099  REAL*8, INTENT(IN) :: dt, te, pl, dzet, sundqv2, sundqv3, sundqt1, &
3100 & c_00, lwcrit
3101  REAL*8 :: teb
3102 !Prognostic
3103  REAL*8, INTENT(INOUT) :: qc, qp, f
3104  REAL*8 :: qcb, qpb, fb
3105 !Locals
3106  REAL*8 :: acf0, acf, c00x, iqccrx, f2, f3, rate, dqp, qcm, dqfac
3107  REAL*8 :: c00xb, iqccrxb, f2b, f3b, rateb, dqpb, qcmb, dqfacb
3108  INTRINSIC exp
3109  INTRINSIC min
3110  INTRINSIC max
3111  INTEGER :: branch
3112  REAL*8 :: x4
3113  REAL*8 :: x3
3114  REAL*8 :: x2
3115  REAL*8 :: x1
3116  REAL*8 :: x2b
3117  REAL*8 :: tempb
3118  REAL*8 :: x1b
3119  REAL*8 :: x4b
3120  REAL*8 :: temp
3121 !Zero Locals
3122  f2 = 0.0
3123  f3 = 0.0
3124  CALL pushreal8(f2)
3125  CALL cons_sundq3(te, sundqv2, sundqv3, sundqt1, f2, f3)
3126  c00x = c_00*f2*f3
3127  iqccrx = f2*f3/lwcrit
3128  IF (f .GT. 0. .AND. qc .GT. 0.) THEN
3129  qcm = qc/f
3130  CALL pushcontrol1b(0)
3131  ELSE
3132  CALL pushcontrol1b(1)
3133  qcm = 0.
3134  END IF
3135  rate = c00x*(1.0-exp(-((qcm*iqccrx)**2)))
3136 !Temporary kluge until we can figure a better to make thicker low clouds.
3137  CALL pushreal8(f3)
3138  f3 = 1.0
3139 !Implement ramps for gradual change in autoconv
3140 !Thicken low high lat clouds
3141  IF (pl .GE. 775. .AND. te .LE. 275.) f3 = 0.2
3142 !F3 = max(-0.016 * PL + 13.4, 0.2)
3143  IF (pl .GE. 825. .AND. te .LE. 282.) f3 = 0.2
3144 !F3 = max(0.11 * TE - 30.02, 0.2)
3145  IF (pl .GE. 775. .AND. pl .LT. 825. .AND. te .LE. 282. .AND. te .GT. &
3146 & 275.) f3 = 0.2
3147 !F3 = min(max(-0.016*PL + 0.11 * TE - 16.85, 0.2),1.)
3148  IF (pl .GE. 825. .AND. te .LE. 275.) f3 = 0.2
3149  IF (pl .LE. 775. .OR. te .GT. 282.) f3 = 1.
3150 !Thin-out low tropical clouds
3151  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
3152  IF (0.2*te - 56 .GT. 2.) THEN
3153  CALL pushcontrol2b(2)
3154  f3 = 2.
3155  ELSE
3156  f3 = 0.2*te - 56
3157  CALL pushcontrol2b(1)
3158  END IF
3159  ELSE
3160  CALL pushcontrol2b(0)
3161  END IF
3162  IF (pl .GE. 925. .AND. te .GE. 290.) THEN
3163  IF (0.04*pl - 36. .GT. 2.) THEN
3164  CALL pushcontrol1b(1)
3165  f3 = 2.
3166  ELSE
3167  CALL pushcontrol1b(1)
3168  f3 = 0.04*pl - 36.
3169  END IF
3170  ELSE
3171  CALL pushcontrol1b(0)
3172  END IF
3173  IF (pl .GE. 925. .AND. pl .LT. 950. .AND. te .GT. 285. .AND. te .LT. &
3174 & 290.) THEN
3175  IF (0.04*pl + 0.2*te - 94. .GT. 2.) THEN
3176  CALL pushcontrol1b(0)
3177  x1 = 2.
3178  ELSE
3179  x1 = 0.04*pl + 0.2*te - 94.
3180  CALL pushcontrol1b(1)
3181  END IF
3182  IF (x1 .LT. 1.) THEN
3183  f3 = 1.
3184  CALL pushcontrol2b(2)
3185  ELSE
3186  f3 = x1
3187  CALL pushcontrol2b(1)
3188  END IF
3189  ELSE
3190  CALL pushcontrol2b(0)
3191  END IF
3192  IF (pl .GE. 950. .AND. te .GE. 290.) THEN
3193  f3 = 2.
3194  CALL pushcontrol1b(1)
3195  ELSE
3196  CALL pushcontrol1b(0)
3197  END IF
3198  IF (f3 .LT. 0.1) THEN
3199  f3 = 0.1
3200  CALL pushcontrol1b(0)
3201  ELSE
3202  CALL pushcontrol1b(1)
3203  f3 = f3
3204  END IF
3205  CALL pushreal8(rate)
3206  rate = f3*rate
3207  dqp = qc*(1.0-exp(-(rate*dt)))
3208  IF (dqp .LT. 0.0) THEN
3209  dqp = 0.0
3210  CALL pushcontrol1b(0)
3211  ELSE
3212  CALL pushcontrol1b(1)
3213  dqp = dqp
3214  END IF
3215 !Wipe-out warm fogs
3216  dqfac = 0.
3217  IF (pl .GE. 975. .AND. te .GE. 280.) THEN
3218  IF (0.2*te - 56. .GT. 1.) THEN
3219  CALL pushcontrol1b(0)
3220  x2 = 1.
3221  ELSE
3222  x2 = 0.2*te - 56.
3223  CALL pushcontrol1b(1)
3224  END IF
3225  IF (x2 .LT. 0.) THEN
3226  dqfac = 0.
3227  CALL pushcontrol2b(2)
3228  ELSE
3229  dqfac = x2
3230  CALL pushcontrol2b(1)
3231  END IF
3232  ELSE
3233  CALL pushcontrol2b(0)
3234  END IF
3235  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
3236  IF (0.04*pl - 38. .GT. 1.) THEN
3237  CALL pushcontrol1b(1)
3238  x3 = 1.
3239  ELSE
3240  CALL pushcontrol1b(1)
3241  x3 = 0.04*pl - 38.
3242  END IF
3243  IF (x3 .LT. 0.) THEN
3244  dqfac = 0.
3245  ELSE
3246  dqfac = x3
3247  END IF
3248  ELSE
3249  CALL pushcontrol1b(0)
3250  END IF
3251  IF (pl .GE. 950. .AND. pl .LT. 975. .AND. te .GT. 280. .AND. te .LT. &
3252 & 285.) THEN
3253  IF (0.04*pl + 0.2*te - 95. .GT. 1.) THEN
3254  CALL pushcontrol1b(0)
3255  x4 = 1.
3256  ELSE
3257  x4 = 0.04*pl + 0.2*te - 95.
3258  CALL pushcontrol1b(1)
3259  END IF
3260  IF (x4 .LT. 0.) THEN
3261  dqfac = 0.
3262  CALL pushcontrol2b(2)
3263  ELSE
3264  dqfac = x4
3265  CALL pushcontrol2b(1)
3266  END IF
3267  ELSE
3268  CALL pushcontrol2b(0)
3269  END IF
3270  IF (pl .GE. 975. .AND. te .GE. 285.) THEN
3271  dqfac = 1.
3272  CALL pushcontrol1b(1)
3273  ELSE
3274  CALL pushcontrol1b(0)
3275  END IF
3276  IF (dqp .LT. dqfac*qc) THEN
3277  CALL pushcontrol1b(0)
3278  ELSE
3279  CALL pushcontrol1b(1)
3280  END IF
3281  dqpb = qpb - qcb
3282  CALL popcontrol1b(branch)
3283  IF (branch .EQ. 0) THEN
3284  dqfacb = qc*dqpb
3285  qcb = qcb + dqfac*dqpb
3286  dqpb = 0.0_8
3287  ELSE
3288  dqfacb = 0.0_8
3289  END IF
3290  CALL popcontrol1b(branch)
3291  IF (branch .NE. 0) dqfacb = 0.0_8
3292  CALL popcontrol2b(branch)
3293  IF (branch .NE. 0) THEN
3294  IF (branch .EQ. 1) THEN
3295  x4b = dqfacb
3296  ELSE
3297  x4b = 0.0_8
3298  END IF
3299  CALL popcontrol1b(branch)
3300  IF (branch .NE. 0) teb = teb + 0.2*x4b
3301  dqfacb = 0.0_8
3302  END IF
3303  CALL popcontrol1b(branch)
3304  IF (branch .NE. 0) dqfacb = 0.0_8
3305  CALL popcontrol2b(branch)
3306  IF (branch .NE. 0) THEN
3307  IF (branch .EQ. 1) THEN
3308  x2b = dqfacb
3309  ELSE
3310  x2b = 0.0_8
3311  END IF
3312  CALL popcontrol1b(branch)
3313  IF (branch .NE. 0) teb = teb + 0.2*x2b
3314  END IF
3315  CALL popcontrol1b(branch)
3316  IF (branch .EQ. 0) dqpb = 0.0_8
3317  qcb = qcb + (1.0-exp(-(dt*rate)))*dqpb
3318  rateb = exp(-(dt*rate))*qc*dt*dqpb
3319  CALL popreal8(rate)
3320  f3b = rate*rateb
3321  rateb = f3*rateb
3322  CALL popcontrol1b(branch)
3323  IF (branch .EQ. 0) f3b = 0.0_8
3324  CALL popcontrol1b(branch)
3325  IF (branch .NE. 0) f3b = 0.0_8
3326  CALL popcontrol2b(branch)
3327  IF (branch .NE. 0) THEN
3328  IF (branch .EQ. 1) THEN
3329  x1b = f3b
3330  ELSE
3331  x1b = 0.0_8
3332  END IF
3333  CALL popcontrol1b(branch)
3334  IF (branch .NE. 0) teb = teb + 0.2*x1b
3335  f3b = 0.0_8
3336  END IF
3337  CALL popcontrol1b(branch)
3338  IF (branch .NE. 0) f3b = 0.0_8
3339  CALL popcontrol2b(branch)
3340  IF (branch .NE. 0) THEN
3341  IF (branch .EQ. 1) teb = teb + 0.2*f3b
3342  END IF
3343  CALL popreal8(f3)
3344  temp = qcm*iqccrx
3345  tempb = exp(-(temp**2))*c00x*2*temp*rateb
3346  c00xb = (1.0-exp(-(temp**2)))*rateb
3347  qcmb = iqccrx*tempb
3348  iqccrxb = qcm*tempb
3349  CALL popcontrol1b(branch)
3350  IF (branch .EQ. 0) THEN
3351  qcb = qcb + qcmb/f
3352  fb = fb - qc*qcmb/f**2
3353  END IF
3354  f2b = c_00*f3*c00xb + f3*iqccrxb/lwcrit
3355  f2b = 0.5*f2b
3356  CALL popreal8(f2)
3357  CALL cons_sundq3_b(te, teb, sundqv2, sundqv3, sundqt1, f2, f2b, f3)
3358 END SUBROUTINE autoconversion_cnv_b
3359 
3360 ! Differentiation of get_ice_fraction in reverse (adjoint) mode:
3361 ! gradient of useful results: temp icefrct
3362 ! with respect to varying inputs: temp
3363 SUBROUTINE get_ice_fraction_b(temp, tempb, t_ice_all, t_ice_max, &
3364 & icefrpwr, icefrct, icefrctb)
3365  IMPLICIT NONE
3366 !Inputs
3367  REAL*8, INTENT(IN) :: temp, t_ice_all, t_ice_max
3368  REAL*8 :: tempb
3369  INTEGER, INTENT(IN) :: icefrpwr
3370 !Outputs
3371  REAL*8 :: icefrct
3372  REAL*8 :: icefrctb
3373  INTRINSIC min
3374  INTRINSIC max
3375  INTEGER :: branch
3376  icefrct = 0.00
3377  IF (temp .LE. t_ice_all) THEN
3378  CALL pushcontrol2b(2)
3379  icefrct = 1.000
3380  ELSE IF (temp .GT. t_ice_all .AND. temp .LE. t_ice_max) THEN
3381  icefrct = 1.00 - (temp-t_ice_all)/(t_ice_max-t_ice_all)
3382  CALL pushcontrol2b(1)
3383  ELSE
3384  CALL pushcontrol2b(0)
3385  END IF
3386  IF (icefrct .GT. 1.00) THEN
3387  icefrct = 1.00
3388  CALL pushcontrol1b(0)
3389  ELSE
3390  CALL pushcontrol1b(1)
3391  icefrct = icefrct
3392  END IF
3393  IF (icefrct .LT. 0.00) THEN
3394  icefrct = 0.00
3395  CALL pushcontrol1b(0)
3396  ELSE
3397  CALL pushcontrol1b(1)
3398  icefrct = icefrct
3399  END IF
3400  IF (icefrct .LE. 0.0 .AND. (icefrpwr .EQ. 0.0 .OR. icefrpwr .NE. int(&
3401 & icefrpwr))) THEN
3402  icefrctb = 0.0
3403  ELSE
3404  icefrctb = icefrpwr*icefrct**(icefrpwr-1)*icefrctb
3405  END IF
3406  CALL popcontrol1b(branch)
3407  IF (branch .EQ. 0) icefrctb = 0.0_8
3408  CALL popcontrol1b(branch)
3409  IF (branch .EQ. 0) icefrctb = 0.0_8
3410  CALL popcontrol2b(branch)
3411  IF (branch .NE. 0) THEN
3412  IF (branch .EQ. 1) tempb = tempb - icefrctb/(t_ice_max-t_ice_all)
3413  END IF
3414 END SUBROUTINE get_ice_fraction_b
3415 
3416 ! Differentiation of cons_sundq3 in reverse (adjoint) mode:
3417 ! gradient of useful results: temp f2
3418 ! with respect to varying inputs: temp
3419 SUBROUTINE cons_sundq3_b(temp, tempb, rate2, rate3, te1, f2, f2b, f3)
3420  IMPLICIT NONE
3421 !Inputs
3422  REAL*8, INTENT(IN) :: rate2, rate3, te1, temp
3423  REAL*8 :: tempb
3424 !Outputs
3425  REAL*8 :: f2, f3
3426  REAL*8 :: f2b
3427 !Locals
3428  REAL*8, PARAMETER :: te0=273.
3429  REAL*8, PARAMETER :: te2=200.
3430  REAL*8 :: jump1
3431  INTRINSIC abs
3432  INTRINSIC min
3433  INTEGER :: branch
3434  REAL*8 :: abs0
3435  jump1 = (rate2-1.0)/(te0-te1)**0.333
3436 !Ice - phase treatment
3437  IF (temp .GE. te0) f2 = 1.0
3438  IF (temp .GE. te1 .AND. temp .LT. te0) THEN
3439  IF (te0 - temp .GE. 0.) THEN
3440  abs0 = te0 - temp
3441  ELSE
3442  abs0 = -(te0-temp)
3443  END IF
3444  IF (abs0 .GT. 0.0) THEN
3445 !Linearisation security
3446  f2 = 1.0 + jump1*(te0-temp)**0.3333
3447  CALL pushcontrol2b(0)
3448  ELSE
3449  CALL pushcontrol2b(1)
3450  f2 = 1.0
3451  END IF
3452  ELSE
3453  CALL pushcontrol2b(2)
3454  END IF
3455  IF (temp .LT. te1) THEN
3456  f2 = rate2 + (rate3-rate2)*(te1-temp)/(te1-te2)
3457  CALL pushcontrol1b(1)
3458  ELSE
3459  CALL pushcontrol1b(0)
3460  END IF
3461  IF (f2 .GT. 27.0) f2b = 0.0_8
3462  CALL popcontrol1b(branch)
3463  IF (branch .NE. 0) THEN
3464  tempb = tempb - (rate3-rate2)*f2b/(te1-te2)
3465  f2b = 0.0_8
3466  END IF
3467  CALL popcontrol2b(branch)
3468  IF (branch .EQ. 0) tempb = tempb - 0.3333*(te0-temp)**(-0.6667)*jump1*&
3469 & f2b
3470 END SUBROUTINE cons_sundq3_b
3471 
3472 ! Differentiation of cons_microphys in reverse (adjoint) mode:
3473 ! gradient of useful results: aa temp bb q_sat alhx3
3474 ! with respect to varying inputs: temp q_sat alhx3
3475 SUBROUTINE cons_microphys_b(temp, tempb1, pr, q_sat, q_satb, aa, aab, bb&
3476 & , bbb, cons_h2omw, cons_airmw, cons_rvap, alhx3, alhx3b)
3477  IMPLICIT NONE
3478 !Inputs
3479  REAL*8, INTENT(IN) :: temp, q_sat, pr, alhx3
3480  REAL*8 :: tempb1, q_satb, alhx3b
3481  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_rvap
3482 !Outputs
3483  REAL*8 :: aa, bb
3484  REAL*8 :: aab, bbb
3485 !Locals
3486  REAL*8, PARAMETER :: k_cond=2.4e-2
3487  REAL*8, PARAMETER :: diffu=2.2e-5
3488  REAL*8 :: e_sat, epsi
3489  REAL*8 :: e_satb
3490  REAL*8 :: temp2
3491  REAL*8 :: temp1
3492  REAL*8 :: temp0
3493  REAL*8 :: tempb0
3494  REAL*8 :: tempb
3495  epsi = cons_h2omw/cons_airmw
3496 ! (100 converts from mbar to Pa)
3497  e_sat = 100.*pr*q_sat/(epsi+(1.0-epsi)*q_sat)
3498  temp1 = k_cond*cons_rvap*temp**2
3499  temp2 = 1000.*diffu*e_sat
3500  tempb = cons_rvap*pr*bbb/temp2
3501  tempb1 = tempb1 + tempb - k_cond*cons_rvap*alhx3**2*2*temp*aab/temp1**&
3502 & 2
3503  e_satb = -(temp*diffu*1000.*tempb/temp2)
3504  alhx3b = alhx3b + 2*alhx3*aab/temp1
3505  temp0 = epsi + (-epsi+1.0)*q_sat
3506  tempb0 = pr*100.*e_satb/temp0
3507  q_satb = q_satb + (1.0_8-q_sat*(1.0-epsi)/temp0)*tempb0
3508 END SUBROUTINE cons_microphys_b
3509 
3510 ! Differentiation of cons_alhx in reverse (adjoint) mode:
3511 ! gradient of useful results: t alhx3
3512 ! with respect to varying inputs: t alhx3
3513 SUBROUTINE cons_alhx_b(t, tb, alhx3, alhx3b, t_ice_max, t_ice_all, &
3514 & cons_alhs, cons_alhl)
3515  IMPLICIT NONE
3516 !Inputs
3517  REAL*8, INTENT(IN) :: t, t_ice_max, t_ice_all
3518  REAL*8 :: tb
3519  REAL*8, INTENT(IN) :: cons_alhs, cons_alhl
3520 !Outputs
3521  REAL*8 :: alhx3
3522  REAL*8 :: alhx3b
3523  INTEGER :: branch
3524  IF (t .LT. t_ice_all) THEN
3525  CALL pushcontrol1b(0)
3526  ELSE
3527  CALL pushcontrol1b(1)
3528  END IF
3529  IF (t .GT. t_ice_max) THEN
3530  CALL pushcontrol1b(0)
3531  ELSE
3532  CALL pushcontrol1b(1)
3533  END IF
3534  IF (t .LE. t_ice_max .AND. t .GE. t_ice_all) THEN
3535  tb = tb + (cons_alhl-cons_alhs)*alhx3b/(t_ice_max-t_ice_all)
3536  alhx3b = 0.0_8
3537  END IF
3538  CALL popcontrol1b(branch)
3539  IF (branch .EQ. 0) alhx3b = 0.0_8
3540  CALL popcontrol1b(branch)
3541  IF (branch .EQ. 0) alhx3b = 0.0_8
3542 END SUBROUTINE cons_alhx_b
3543 
3544 ! Differentiation of ice_settlefall_cnv in reverse (adjoint) mode:
3545 ! gradient of useful results: f qi qp dz te
3546 ! with respect to varying inputs: f qi dz te
3547 SUBROUTINE ice_settlefall_cnv_b(wxr, qi, qib, pl, te, teb, f, fb, &
3548 & cons_rgas, khu, khl, k, dt, dz, dzb, qp, qpb, anv_icefall_c)
3549  IMPLICIT NONE
3550 !Inputs
3551  REAL*8, INTENT(IN) :: wxr, pl, te, dz, dt, anv_icefall_c
3552  REAL*8 :: teb, dzb
3553  REAL*8, INTENT(IN) :: cons_rgas
3554  INTEGER, INTENT(IN) :: khu, khl, k
3555  REAL*8, INTENT(INOUT) :: qi, f, qp
3556  REAL*8 :: qib, fb, qpb
3557 !Locals
3558  REAL*8 :: rho, xim, lxim, qixp, vf
3559  REAL*8 :: rhob, ximb, lximb, qixpb, vfb
3560  INTRINSIC log10
3561  INTRINSIC max
3562  INTRINSIC min
3563  INTEGER :: branch
3564  REAL*8 :: tempb0
3565  REAL*8 :: tempb
3566  REAL*8 :: max1
3567 ! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
3568  rho = 1000.*100.*pl/(cons_rgas*te)
3569  IF (f .GT. 0. .AND. qi .GT. 0.) THEN
3570  xim = qi/f*rho
3571  CALL pushcontrol1b(0)
3572  ELSE
3573  xim = 0.
3574  CALL pushcontrol1b(1)
3575  END IF
3576  IF (xim .GT. 0.) THEN
3577  lxim = log10(xim)
3578  CALL pushcontrol1b(0)
3579  ELSE
3580  lxim = 0.0
3581  CALL pushcontrol1b(1)
3582  END IF
3583  vf = 128.6 + 53.2*lxim + 5.5*lxim**2
3584 !VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! )
3585 ! Assume unmodified they represent situation at 100 mb
3586  IF (wxr .GT. 0.) THEN
3587  IF (pl .LT. 10.) THEN
3588  max1 = 10.
3589  ELSE
3590  max1 = pl
3591  END IF
3592  vf = vf*(100./max1)**wxr
3593  CALL pushcontrol1b(0)
3594  ELSE
3595  CALL pushcontrol1b(1)
3596  END IF
3597  vf = vf/100.
3598  IF (khu .GT. 0 .AND. khl .GT. 0) THEN
3599  IF (k - 1 .GE. khu .AND. k - 1 .LE. khl) THEN
3600  vf = 0.01*vf
3601  CALL pushcontrol2b(2)
3602  ELSE
3603  CALL pushcontrol2b(1)
3604  END IF
3605  ELSE
3606  CALL pushcontrol2b(0)
3607  END IF
3608  vf = anv_icefall_c*vf
3609  qixp = qi*(vf*dt/dz)
3610  IF (qixp .GT. qi) THEN
3611  qixp = qi
3612  CALL pushcontrol1b(0)
3613  ELSE
3614  CALL pushcontrol1b(1)
3615  qixp = qixp
3616  END IF
3617  IF (qixp .LT. 0.0) THEN
3618  CALL pushcontrol1b(0)
3619  ELSE
3620  CALL pushcontrol1b(1)
3621  END IF
3622  qixpb = qpb - qib
3623  CALL popcontrol1b(branch)
3624  IF (branch .EQ. 0) qixpb = 0.0_8
3625  CALL popcontrol1b(branch)
3626  IF (branch .EQ. 0) THEN
3627  qib = qib + qixpb
3628  qixpb = 0.0_8
3629  END IF
3630  tempb0 = dt*qixpb/dz
3631  qib = qib + vf*tempb0
3632  vfb = qi*tempb0
3633  dzb = dzb - qi*vf*tempb0/dz
3634  vfb = anv_icefall_c*vfb
3635  CALL popcontrol2b(branch)
3636  IF (branch .NE. 0) THEN
3637  IF (branch .NE. 1) vfb = 0.01*vfb
3638  END IF
3639  vfb = vfb/100.
3640  CALL popcontrol1b(branch)
3641  IF (branch .EQ. 0) vfb = (100./max1)**wxr*vfb
3642  lximb = (5.5*2*lxim+53.2)*vfb
3643  CALL popcontrol1b(branch)
3644  IF (branch .EQ. 0) THEN
3645  ximb = lximb/(xim*log(10.0))
3646  ELSE
3647  ximb = 0.0_8
3648  END IF
3649  CALL popcontrol1b(branch)
3650  IF (branch .EQ. 0) THEN
3651  tempb = ximb/f
3652  qib = qib + rho*tempb
3653  rhob = qi*tempb
3654  fb = fb - qi*rho*tempb/f
3655  ELSE
3656  rhob = 0.0_8
3657  END IF
3658  teb = teb - pl*100.*1000.*rhob/(cons_rgas*te**2)
3659 END SUBROUTINE ice_settlefall_cnv_b
3660 
3661 ! Differentiation of ice_settlefall_ls in reverse (adjoint) mode:
3662 ! gradient of useful results: f qi qp dz te
3663 ! with respect to varying inputs: f qi dz te
3664 SUBROUTINE ice_settlefall_ls_b(wxr, qi, qib, pl, te, teb, f, fb, &
3665 & cons_rgas, khu, khl, k, dt, dz, dzb, qp, qpb, ls_icefall_c)
3666  IMPLICIT NONE
3667 !Inputs
3668  REAL*8, INTENT(IN) :: wxr, pl, te, dz, dt, ls_icefall_c
3669  REAL*8 :: teb, dzb
3670  REAL*8, INTENT(IN) :: cons_rgas
3671  INTEGER, INTENT(IN) :: khu, khl, k
3672  REAL*8, INTENT(INOUT) :: qi, f, qp
3673  REAL*8 :: qib, fb, qpb
3674 !Locals
3675  REAL*8 :: rho, xim, lxim, qixp, vf
3676  REAL*8 :: rhob, ximb, qixpb, vfb
3677  INTRINSIC log10
3678  INTRINSIC abs
3679  INTRINSIC max
3680  INTRINSIC min
3681  INTEGER :: branch
3682  REAL*8 :: tempb2
3683  REAL*8 :: tempb1
3684  REAL*8 :: tempb0
3685  REAL*8 :: tempb
3686  REAL*8 :: abs0
3687  REAL*8 :: max1
3688 ! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
3689  rho = 1000.*100.*pl/(cons_rgas*te)
3690  IF (f .GT. 0. .AND. qi .GT. 0.) THEN
3691  xim = qi/f*rho
3692  CALL pushcontrol1b(0)
3693  ELSE
3694  xim = 0.
3695  CALL pushcontrol1b(1)
3696  END IF
3697  IF (xim .GE. 0.) THEN
3698  abs0 = xim
3699  ELSE
3700  abs0 = -xim
3701  END IF
3702  IF (abs0 .GT. 0.0) THEN
3703 !Linearisation security
3704  vf = 109.0*xim**0.16
3705  CALL pushcontrol1b(0)
3706  ELSE
3707  vf = 0.0
3708  CALL pushcontrol1b(1)
3709  END IF
3710 !VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! )
3711 ! Assume unmodified they represent situation at 100 mb
3712  IF (wxr .GT. 0.) THEN
3713  IF (pl .LT. 10.) THEN
3714  max1 = 10.
3715  ELSE
3716  max1 = pl
3717  END IF
3718  vf = vf*(100./max1)**wxr
3719  CALL pushcontrol1b(0)
3720  ELSE
3721  CALL pushcontrol1b(1)
3722  END IF
3723  vf = vf/100.
3724  IF (khu .GT. 0 .AND. khl .GT. 0) THEN
3725  IF (k - 1 .GE. khu .AND. k - 1 .LE. khl) THEN
3726  vf = 0.01*vf
3727  CALL pushcontrol2b(2)
3728  ELSE
3729  CALL pushcontrol2b(1)
3730  END IF
3731  ELSE
3732  CALL pushcontrol2b(0)
3733  END IF
3734  vf = ls_icefall_c*vf
3735  qixp = qi*(vf*dt/dz)
3736  IF (qixp .GT. qi) THEN
3737  qixp = qi
3738  CALL pushcontrol1b(0)
3739  ELSE
3740  CALL pushcontrol1b(1)
3741  qixp = qixp
3742  END IF
3743  IF (qixp .LT. 0.0) THEN
3744  qixp = 0.0
3745  CALL pushcontrol1b(0)
3746  ELSE
3747  CALL pushcontrol1b(1)
3748  qixp = qixp
3749  END IF
3750  CALL pushreal8(qi)
3751  qi = qi - qixp
3752  IF (qi + qixp .GT. 0.) THEN
3753  tempb1 = fb/(qi+qixp)
3754  tempb2 = -(qi*f*tempb1/(qi+qixp))
3755  qib = qib + tempb2 + f*tempb1
3756  qixpb = tempb2
3757  fb = qi*tempb1
3758  ELSE
3759  qixpb = 0.0_8
3760  END IF
3761  CALL popreal8(qi)
3762  qixpb = qixpb + qpb - qib
3763  CALL popcontrol1b(branch)
3764  IF (branch .EQ. 0) qixpb = 0.0_8
3765  CALL popcontrol1b(branch)
3766  IF (branch .EQ. 0) THEN
3767  qib = qib + qixpb
3768  qixpb = 0.0_8
3769  END IF
3770  tempb0 = dt*qixpb/dz
3771  qib = qib + vf*tempb0
3772  vfb = qi*tempb0
3773  dzb = dzb - qi*vf*tempb0/dz
3774  vfb = ls_icefall_c*vfb
3775  CALL popcontrol2b(branch)
3776  IF (branch .NE. 0) THEN
3777  IF (branch .NE. 1) vfb = 0.01*vfb
3778  END IF
3779  vfb = vfb/100.
3780  CALL popcontrol1b(branch)
3781  IF (branch .EQ. 0) vfb = (100./max1)**wxr*vfb
3782  CALL popcontrol1b(branch)
3783  IF (branch .EQ. 0) THEN
3784  ximb = 109.0*0.16*xim**(-0.84)*vfb
3785  ELSE
3786  ximb = 0.0_8
3787  END IF
3788  CALL popcontrol1b(branch)
3789  IF (branch .EQ. 0) THEN
3790  tempb = ximb/f
3791  qib = qib + rho*tempb
3792  rhob = qi*tempb
3793  fb = fb - qi*rho*tempb/f
3794  ELSE
3795  rhob = 0.0_8
3796  END IF
3797  teb = teb - pl*100.*1000.*rhob/(cons_rgas*te**2)
3798 END SUBROUTINE ice_settlefall_ls_b
3799 
3800 ! Differentiation of precipandevap in reverse (adjoint) mode:
3801 ! gradient of useful results: evap_dd_above_out aa qv bb
3802 ! subl_dd_above_out qcl pfi_above_out dze qddf3
3803 ! pfl_above_out te
3804 ! with respect to varying inputs: aa area pfl_above_in qv pfi_above_in
3805 ! bb evap_dd_above_in qcl qpi qpl subl_dd_above_in
3806 ! dze qddf3 te
3807 SUBROUTINE precipandevap_b(k, ktop, lm, dt, frland, rhcr3, qpl, qplb, &
3808 & qpi, qpib, qcl, qclb, qci, te, teb, qv, qvb, mass, imass, pl, dze, &
3809 & dzeb, qddf3, qddf3b, aa, aab, bb, bbb, area, areab, pfl_above_in, &
3810 & pfl_above_inb, pfl_above_out, pfl_above_outb, pfi_above_in, &
3811 & pfi_above_inb, pfi_above_out, pfi_above_outb, evap_dd_above_in, &
3812 & evap_dd_above_inb, evap_dd_above_out, evap_dd_above_outb, &
3813 & subl_dd_above_in, subl_dd_above_inb, subl_dd_above_out, &
3814 & subl_dd_above_outb, envfc, ddrfc, cons_alhf, cons_alhs, cons_alhl, &
3815 & cons_cp, cons_tice, cons_h2omw, cons_airmw, revap_off_p, c_acc, c_ev_r&
3816 & , c_ev_s, rho_w, estblx)
3817  IMPLICIT NONE
3818 !Inputs
3819  INTEGER, INTENT(IN) :: k, lm, ktop
3820  REAL*8, INTENT(IN) :: dt, mass, imass, pl, aa, bb, rhcr3, dze, qddf3, &
3821 & area, frland, envfc, ddrfc
3822  REAL*8 :: aab, bbb, dzeb, qddf3b, areab
3823  REAL*8, INTENT(IN) :: cons_alhf, cons_alhs, cons_alhl, cons_cp, &
3824 & cons_tice, cons_h2omw, cons_airmw
3825  REAL*8, INTENT(IN) :: revap_off_p
3826  REAL*8, INTENT(IN) :: c_acc, c_ev_r, c_ev_s, rho_w
3827  REAL*8, INTENT(IN) :: estblx(:)
3828 !Prognostics
3829  REAL*8, INTENT(INOUT) :: qv, qpl, qpi, qcl, qci, te
3830  REAL*8 :: qvb, qplb, qpib, qclb, teb
3831  REAL*8, INTENT(INOUT) :: pfl_above_in, pfl_above_out, pfi_above_in, &
3832 & pfi_above_out
3833  REAL*8 :: pfl_above_inb, pfl_above_outb, pfi_above_inb, pfi_above_outb
3834  REAL*8, INTENT(INOUT) :: evap_dd_above_in, evap_dd_above_out, &
3835 & subl_dd_above_in, subl_dd_above_out
3836  REAL*8 :: evap_dd_above_inb, evap_dd_above_outb, subl_dd_above_inb, &
3837 & subl_dd_above_outb
3838 !Locals
3839  INTEGER :: ns, nsmx, itr, l
3840  REAL*8 :: pfi, pfl, qs, dqs, envfrac, tko, qko, qstko, dqstko, rh_box&
3841 & , t_ed, qplko, qpiko
3842  REAL*8 :: pfib, pflb, qsb, dqsb, tkob, qkob, qstkob, dqstkob, rh_boxb&
3843 & , t_edb
3844  REAL*8 :: ifactor, rainrat0, snowrat0, fallrn, fallsn, vesn, vern, &
3845 & nrain, nsnow, efactor
3846  REAL*8 :: ifactorb, rainrat0b, snowrat0b, fallrnb, fallsnb, vesnb, &
3847 & vernb, efactorb
3848  REAL*8 :: tinlayerrn, diamrn, droprad, tinlayersn, diamsn, flakrad
3849  REAL*8 :: tinlayerrnb, diamrnb, dropradb, tinlayersnb, diamsnb, &
3850 & flakradb
3851  REAL*8 :: evap, subl, accr, mltfrz, evapx, sublx, evap_dd, subl_dd, &
3852 & ddfract, landseaf
3853  REAL*8 :: evapb, sublb, accrb, mltfrzb, evap_ddb, subl_ddb
3854  REAL*8 :: tau_frz, tau_mlt
3855 !m/s
3856  REAL*8, PARAMETER :: trmv_l=1.0
3857  LOGICAL, PARAMETER :: taneff=.false.
3858 !Fraction of precip falling through "environment" vs through cloud
3859  REAL*8, PARAMETER :: b_sub=1.00
3860  INTRINSIC max
3861  INTRINSIC min
3862  INTRINSIC exp
3863  INTEGER :: branch
3864  REAL*8 :: temp3
3865  REAL*8 :: temp2
3866  REAL*8 :: temp1
3867  REAL*8 :: temp0
3868  REAL*8 :: tempb8
3869  REAL*8 :: tempb7
3870  REAL*8 :: tempb6
3871  REAL*8 :: tempb5
3872  REAL*8 :: tempb4
3873  REAL*8 :: tempb3
3874  REAL*8 :: tempb2
3875  REAL*8 :: tempb1
3876  REAL*8 :: tempb0
3877  REAL*8 :: tempb
3878  REAL*8 :: temp
3879  REAL*8 :: temp4
3880  envfrac = envfc
3881  IF (area .GT. 0.) THEN
3882  ifactor = 1./area
3883  CALL pushcontrol1b(1)
3884  ELSE
3885  ifactor = 1.00
3886  CALL pushcontrol1b(0)
3887  END IF
3888  IF (ifactor .LT. 1.) THEN
3889  ifactor = 1.
3890  CALL pushcontrol1b(0)
3891  ELSE
3892  CALL pushcontrol1b(1)
3893  ifactor = ifactor
3894  END IF
3895 !Start at top of precip column:
3896 !
3897 ! a) Accrete
3898 ! b) Evaporate/Sublimate
3899 ! c) Rain/Snow-out to next level down
3900 ! d) return to (a)
3901 !Update saturated humidity
3902  CALL dqsats_bac(dqs, qs, te, pl, estblx, cons_h2omw, cons_airmw)
3903  ddfract = ddrfc
3904  IF (k .EQ. ktop) THEN
3905  evap_dd = 0.
3906  subl_dd = 0.
3907  CALL pushcontrol1b(0)
3908  ELSE
3909  CALL pushreal8(qpl)
3910  qpl = qpl + pfl_above_in*imass
3911  CALL pushreal8(qpi)
3912  qpi = qpi + pfi_above_in*imass
3913  accr = b_sub*c_acc*(qpl*mass)*qcl
3914  IF (accr .GT. qcl) THEN
3915  accr = qcl
3916  CALL pushcontrol1b(0)
3917  ELSE
3918  CALL pushcontrol1b(1)
3919  accr = accr
3920  END IF
3921  CALL pushreal8(qpl)
3922  qpl = qpl + accr
3923  CALL pushreal8(qcl)
3924  qcl = qcl - accr
3925 !Accretion of liquid condensate by falling ice/snow
3926  accr = b_sub*c_acc*(qpi*mass)*qcl
3927  IF (accr .GT. qcl) THEN
3928  accr = qcl
3929  CALL pushcontrol1b(0)
3930  ELSE
3931  CALL pushcontrol1b(1)
3932  accr = accr
3933  END IF
3934  CALL pushreal8(qpi)
3935  qpi = qpi + accr
3936 !! Liquid freezes when accreted by snow
3937  CALL pushreal8(te)
3938  te = te + cons_alhf*accr/cons_cp
3939  rainrat0 = ifactor*qpl*mass/dt
3940  snowrat0 = ifactor*qpi*mass/dt
3941  CALL pushreal8(diamrn)
3942  CALL marshpalm(rainrat0, pl, diamrn, nrain, fallrn, vern)
3943  CALL pushreal8(diamsn)
3944  CALL marshpalm(snowrat0, pl, diamsn, nsnow, fallsn, vesn)
3945  tinlayerrn = dze/(fallrn+0.01)
3946  tinlayersn = dze/(fallsn+0.01)
3947 !Melting of Frozen precipitation
3948 ! time scale for freezing (s).
3949  tau_frz = 5000.
3950  IF (te .GT. cons_tice .AND. te .LE. cons_tice + 5.) THEN
3951  mltfrz = tinlayersn*qpi*(te-cons_tice)/tau_frz
3952  IF (qpi .GT. mltfrz) THEN
3953  CALL pushcontrol1b(0)
3954  mltfrz = mltfrz
3955  ELSE
3956  mltfrz = qpi
3957  CALL pushcontrol1b(1)
3958  END IF
3959  CALL pushreal8(te)
3960  te = te - cons_alhf*mltfrz/cons_cp
3961  CALL pushreal8(qpl)
3962  qpl = qpl + mltfrz
3963  CALL pushreal8(qpi)
3964  qpi = qpi - mltfrz
3965  CALL pushcontrol1b(0)
3966  ELSE
3967  CALL pushcontrol1b(1)
3968  END IF
3969  IF (te .GT. cons_tice + 5.) THEN
3970 ! Go Ahead and melt any snow/hail left above 5 C
3971  mltfrz = qpi
3972  te = te - cons_alhf*mltfrz/cons_cp
3973  CALL pushreal8(qpl)
3974  qpl = qpl + mltfrz
3975  CALL pushreal8(qpi)
3976  qpi = qpi - mltfrz
3977  CALL pushcontrol1b(0)
3978  ELSE
3979  CALL pushcontrol1b(1)
3980  END IF
3981  IF (k .GE. lm - 1) THEN
3982  IF (te .GT. cons_tice + 0.) THEN
3983 ! Go Ahead and melt any snow/hail left above 0 C in lowest layers
3984  mltfrz = qpi
3985  te = te - cons_alhf*mltfrz/cons_cp
3986  CALL pushreal8(qpl)
3987  qpl = qpl + mltfrz
3988  CALL pushreal8(qpi)
3989  qpi = qpi - mltfrz
3990  CALL pushcontrol2b(0)
3991  ELSE
3992  CALL pushcontrol2b(1)
3993  END IF
3994  ELSE
3995  CALL pushcontrol2b(2)
3996  END IF
3997 !Freezing of liquid precipitation
3998  IF (te .LE. cons_tice) THEN
3999  te = te + cons_alhf*qpl/cons_cp
4000  CALL pushreal8(qpi)
4001  qpi = qpl + qpi
4002  CALL pushreal8(qpl)
4003  qpl = 0.
4004  CALL pushcontrol1b(1)
4005  ELSE
4006  CALL pushcontrol1b(0)
4007  END IF
4008 !In the exp below, evaporation time scale is determined "microphysically" from temp,
4009 !press, and drop size. In this context C_EV becomes a dimensionless fudge-fraction.
4010 !Also remember that these microphysics are still only for liquid.
4011  qko = qv
4012  tko = te
4013 !do itr = 1,1
4014  dqstko = dqs
4015  qstko = qs + dqstko*(tko-te)
4016  IF (qstko .LT. 1.0e-7) THEN
4017  qstko = 1.0e-7
4018  CALL pushcontrol1b(0)
4019  ELSE
4020  CALL pushcontrol1b(1)
4021  qstko = qstko
4022  END IF
4023  rh_box = qko/qstko
4024  IF (rh_box .LT. rhcr3) THEN
4025  efactor = rho_w*(aa+bb)/(rhcr3-rh_box)
4026  CALL pushcontrol1b(0)
4027  ELSE
4028  efactor = 9.99e9
4029  CALL pushcontrol1b(1)
4030  END IF
4031  landseaf = 1.00
4032 !Rain falling
4033  IF (rh_box .LT. rhcr3 .AND. diamrn .GT. 0.00 .AND. pl .GT. 100. &
4034 & .AND. pl .LT. revap_off_p) THEN
4035  droprad = 0.5*diamrn
4036  t_ed = efactor*droprad**2
4037  CALL pushreal8(t_ed)
4038  t_ed = t_ed*(1.0+dqstko*cons_alhl/cons_cp)
4039  evap = qpl*(1.0-exp(-(c_ev_r*vern*landseaf*envfrac*tinlayerrn/t_ed&
4040 & )))
4041  CALL pushcontrol1b(0)
4042  ELSE
4043  evap = 0.0
4044  CALL pushcontrol1b(1)
4045  END IF
4046 !Snow falling
4047  IF (rh_box .LT. rhcr3 .AND. diamsn .GT. 0.00 .AND. pl .GT. 100. &
4048 & .AND. pl .LT. revap_off_p) THEN
4049  flakrad = 0.5*diamsn
4050  CALL pushreal8(t_ed)
4051  t_ed = efactor*flakrad**2
4052  CALL pushreal8(t_ed)
4053  t_ed = t_ed*(1.0+dqstko*cons_alhs/cons_cp)
4054  subl = qpi*(1.0-exp(-(c_ev_s*vesn*landseaf*envfrac*tinlayersn/t_ed&
4055 & )))
4056  CALL pushcontrol1b(0)
4057  ELSE
4058  subl = 0.0
4059  CALL pushcontrol1b(1)
4060  END IF
4061 !if (itr == 1) then
4062 ! EVAPx = EVAP
4063 ! SUBLx = SUBL
4064 !else
4065 ! EVAP = (EVAP+EVAPx) /2.0
4066 ! SUBL = (SUBL+SUBLx) /2.0
4067 !endif
4068 !enddo
4069 !Put some re-evap/re-subl precip in to a \quote{downdraft} to be applied later
4070  evap_dd = evap_dd_above_in + ddfract*evap*mass
4071  subl_dd = subl_dd_above_in + ddfract*subl*mass
4072  CALL pushcontrol1b(1)
4073  END IF
4074  sublb = qvb - cons_alhs*teb/cons_cp
4075  evapb = qvb - cons_alhl*teb/cons_cp
4076  subl_ddb = qddf3*sublb/mass + subl_dd_above_outb
4077  evap_ddb = qddf3*evapb/mass + evap_dd_above_outb
4078  pfib = pfi_above_outb
4079  pflb = pfl_above_outb
4080  qddf3b = qddf3b + evap_dd*evapb/mass + subl_dd*sublb/mass
4081  CALL popcontrol1b(branch)
4082  IF (branch .EQ. 0) THEN
4083  qpib = mass*pfib
4084  qplb = mass*pflb
4085  pfl_above_inb = 0.0_8
4086  pfi_above_inb = 0.0_8
4087  evap_dd_above_inb = 0.0_8
4088  subl_dd_above_inb = 0.0_8
4089  dqsb = 0.0_8
4090  qsb = 0.0_8
4091  ifactorb = 0.0_8
4092  ELSE
4093  qpib = mass*pfib
4094  qplb = mass*pflb
4095  evapb = qvb - cons_alhl*teb/cons_cp
4096  sublb = qvb - cons_alhs*teb/cons_cp
4097  sublb = ddfract*mass*subl_ddb - qpib + (1.0_8-ddfract)*sublb
4098  subl_dd_above_inb = subl_ddb
4099  evapb = ddfract*mass*evap_ddb - qplb + (1.0_8-ddfract)*evapb
4100  evap_dd_above_inb = evap_ddb
4101  CALL popcontrol1b(branch)
4102  IF (branch .EQ. 0) THEN
4103  temp2 = vesn*tinlayersn/t_ed
4104  temp4 = c_ev_s*landseaf*envfrac
4105  temp3 = -(temp4*temp2)
4106  tempb8 = temp4*exp(temp3)*qpi*sublb/t_ed
4107  qpib = qpib + (1.0-exp(temp3))*sublb
4108  vesnb = tinlayersn*tempb8
4109  tinlayersnb = vesn*tempb8
4110  t_edb = -(temp2*tempb8)
4111  CALL popreal8(t_ed)
4112  dqstkob = cons_alhs*t_ed*t_edb/cons_cp
4113  t_edb = (cons_alhs*(dqstko/cons_cp)+1.0)*t_edb
4114  CALL popreal8(t_ed)
4115  efactorb = flakrad**2*t_edb
4116  flakradb = efactor*2*flakrad*t_edb
4117  diamsnb = 0.5*flakradb
4118  ELSE
4119  dqstkob = 0.0_8
4120  efactorb = 0.0_8
4121  diamsnb = 0.0_8
4122  vesnb = 0.0_8
4123  tinlayersnb = 0.0_8
4124  END IF
4125  CALL popcontrol1b(branch)
4126  IF (branch .EQ. 0) THEN
4127  temp = vern*tinlayerrn/t_ed
4128  temp1 = c_ev_r*landseaf*envfrac
4129  temp0 = -(temp1*temp)
4130  tempb7 = temp1*exp(temp0)*qpl*evapb/t_ed
4131  qplb = qplb + (1.0-exp(temp0))*evapb
4132  vernb = tinlayerrn*tempb7
4133  tinlayerrnb = vern*tempb7
4134  t_edb = -(temp*tempb7)
4135  CALL popreal8(t_ed)
4136  dqstkob = dqstkob + cons_alhl*t_ed*t_edb/cons_cp
4137  t_edb = (cons_alhl*(dqstko/cons_cp)+1.0)*t_edb
4138  efactorb = efactorb + droprad**2*t_edb
4139  dropradb = efactor*2*droprad*t_edb
4140  diamrnb = 0.5*dropradb
4141  ELSE
4142  diamrnb = 0.0_8
4143  vernb = 0.0_8
4144  tinlayerrnb = 0.0_8
4145  END IF
4146  CALL popcontrol1b(branch)
4147  IF (branch .EQ. 0) THEN
4148  tempb6 = rho_w*efactorb/(rhcr3-rh_box)
4149  aab = aab + tempb6
4150  bbb = bbb + tempb6
4151  rh_boxb = (aa+bb)*tempb6/(rhcr3-rh_box)
4152  ELSE
4153  rh_boxb = 0.0_8
4154  END IF
4155  qkob = rh_boxb/qstko
4156  qstkob = -(qko*rh_boxb/qstko**2)
4157  CALL popcontrol1b(branch)
4158  IF (branch .EQ. 0) qstkob = 0.0_8
4159  qsb = qstkob
4160  dqstkob = dqstkob + (tko-te)*qstkob
4161  tkob = dqstko*qstkob
4162  teb = teb + tkob - dqstko*qstkob
4163  dqsb = dqstkob
4164  qvb = qvb + qkob
4165  CALL popcontrol1b(branch)
4166  IF (branch .NE. 0) THEN
4167  CALL popreal8(qpl)
4168  CALL popreal8(qpi)
4169  qplb = cons_alhf*teb/cons_cp + qpib
4170  END IF
4171  CALL popcontrol2b(branch)
4172  IF (branch .EQ. 0) THEN
4173  CALL popreal8(qpi)
4174  mltfrzb = qplb - cons_alhf*teb/cons_cp - qpib
4175  CALL popreal8(qpl)
4176  qpib = qpib + mltfrzb
4177  END IF
4178  CALL popcontrol1b(branch)
4179  IF (branch .EQ. 0) THEN
4180  CALL popreal8(qpi)
4181  mltfrzb = qplb - cons_alhf*teb/cons_cp - qpib
4182  CALL popreal8(qpl)
4183  qpib = qpib + mltfrzb
4184  END IF
4185  CALL popcontrol1b(branch)
4186  IF (branch .EQ. 0) THEN
4187  CALL popreal8(qpi)
4188  mltfrzb = qplb - cons_alhf*teb/cons_cp - qpib
4189  CALL popreal8(qpl)
4190  CALL popreal8(te)
4191  CALL popcontrol1b(branch)
4192  IF (branch .NE. 0) THEN
4193  qpib = qpib + mltfrzb
4194  mltfrzb = 0.0_8
4195  END IF
4196  tempb5 = (te-cons_tice)*mltfrzb/tau_frz
4197  tinlayersnb = tinlayersnb + qpi*tempb5
4198  qpib = qpib + tinlayersn*tempb5
4199  teb = teb + tinlayersn*qpi*mltfrzb/tau_frz
4200  END IF
4201  tempb2 = tinlayerrnb/(fallrn+0.01)
4202  tempb1 = tinlayersnb/(fallsn+0.01)
4203  dzeb = dzeb + tempb2 + tempb1
4204  fallsnb = -(dze*tempb1/(fallsn+0.01))
4205  fallrnb = -(dze*tempb2/(fallrn+0.01))
4206  CALL popreal8(diamsn)
4207  CALL marshpalm_b(snowrat0, snowrat0b, pl, diamsn, diamsnb, nsnow, &
4208 & fallsn, fallsnb, vesn, vesnb)
4209  CALL popreal8(diamrn)
4210  CALL marshpalm_b(rainrat0, rainrat0b, pl, diamrn, diamrnb, nrain, &
4211 & fallrn, fallrnb, vern, vernb)
4212  tempb3 = mass*snowrat0b/dt
4213  qpib = qpib + ifactor*tempb3
4214  tempb4 = mass*rainrat0b/dt
4215  ifactorb = qpl*tempb4 + qpi*tempb3
4216  qplb = qplb + ifactor*tempb4
4217  CALL popreal8(te)
4218  accrb = qpib - qclb + cons_alhf*teb/cons_cp
4219  CALL popreal8(qpi)
4220  CALL popcontrol1b(branch)
4221  IF (branch .EQ. 0) THEN
4222  qclb = qclb + accrb
4223  accrb = 0.0_8
4224  END IF
4225  tempb0 = b_sub*c_acc*mass*accrb
4226  qpib = qpib + qcl*tempb0
4227  qclb = qclb + qpi*tempb0
4228  CALL popreal8(qcl)
4229  accrb = qplb - qclb
4230  CALL popreal8(qpl)
4231  CALL popcontrol1b(branch)
4232  IF (branch .EQ. 0) THEN
4233  qclb = qclb + accrb
4234  accrb = 0.0_8
4235  END IF
4236  tempb = b_sub*c_acc*mass*accrb
4237  qplb = qplb + qcl*tempb
4238  qclb = qclb + qpl*tempb
4239  CALL popreal8(qpi)
4240  pfi_above_inb = imass*qpib
4241  CALL popreal8(qpl)
4242  pfl_above_inb = imass*qplb
4243  END IF
4244  CALL dqsats_bac_b(dqs, dqsb, qs, qsb, te, teb, pl, estblx, cons_h2omw&
4245 & , cons_airmw)
4246  CALL popcontrol1b(branch)
4247  IF (branch .EQ. 0) ifactorb = 0.0_8
4248  CALL popcontrol1b(branch)
4249  IF (branch .EQ. 0) THEN
4250  areab = 0.0_8
4251  ELSE
4252  areab = -(ifactorb/area**2)
4253  END IF
4254 END SUBROUTINE precipandevap_b
4255 
4256 ! Differentiation of marshpalm in reverse (adjoint) mode:
4257 ! gradient of useful results: diam3 w ve
4258 ! with respect to varying inputs: rain
4259 SUBROUTINE marshpalm_b(rain, rainb, pr, diam3, diam3b, ntotal, w, wb, ve&
4260 & , veb)
4261  IMPLICIT NONE
4262 !Inputs
4263 ! in kg m^-2 s^-1, mbar
4264  REAL*8, INTENT(IN) :: rain, pr
4265  REAL*8 :: rainb
4266 !Outputs
4267  REAL*8 :: diam3, ntotal, w, ve
4268  REAL*8 :: diam3b, wb, veb
4269 !Locals
4270  INTEGER :: iqd
4271 !cm^-3
4272  REAL*8, PARAMETER :: n0=0.08
4273  REAL*8 :: rain_day, slopr, diam1
4274  REAL*8 :: rain_dayb
4275  REAL*8 :: rx(8), d3x(8)
4276  INTRINSIC sqrt
4277  INTRINSIC max
4278  INTEGER :: branch
4279 !Marshall-Palmer sizes at different rain-rates: avg(D^3)
4280 !RX = (/ 0. , 5. , 20. , 80. , 320. , 1280., 5120., 20480. /) ! rain per in mm/day
4281  rx(1) = 0.
4282  rx(2) = 5.
4283  rx(3) = 20.
4284  rx(4) = 80.
4285  rx(5) = 320.
4286  rx(6) = 1280.
4287  rx(7) = 5120.
4288  rx(8) = 20480.
4289 !D3X= (/ 0.019, 0.032, 0.043, 0.057, 0.076, 0.102, 0.137, 0.183 /)
4290  d3x(1) = 0.019
4291  d3x(2) = 0.032
4292  d3x(3) = 0.043
4293  d3x(4) = 0.057
4294  d3x(5) = 0.076
4295  d3x(6) = 0.102
4296  d3x(7) = 0.137
4297  d3x(8) = 0.183
4298  rain_day = rain*3600.*24.
4299  IF (rain_day .LE. 0.00) diam3 = 0.00
4300  DO iqd=1,7
4301  IF (rain_day .LE. rx(iqd+1) .AND. rain_day .GT. rx(iqd)) THEN
4302  CALL pushreal8(slopr)
4303  slopr = (d3x(iqd+1)-d3x(iqd))/(rx(iqd+1)-rx(iqd))
4304  diam3 = d3x(iqd) + (rain_day-rx(iqd))*slopr
4305  CALL pushcontrol1b(1)
4306  ELSE
4307  CALL pushcontrol1b(0)
4308  END IF
4309  END DO
4310  IF (rain_day .GE. rx(8)) THEN
4311  diam3 = d3x(8)
4312  CALL pushcontrol1b(1)
4313  ELSE
4314  CALL pushcontrol1b(0)
4315  END IF
4316  diam3 = 0.664*diam3
4317  w = (2483.8*diam3+80.)*sqrt(1000./pr)
4318  IF (0.99*w/100. .LT. 1.000) THEN
4319  CALL pushcontrol1b(0)
4320  ELSE
4321  CALL pushcontrol1b(1)
4322  END IF
4323  wb = wb/100.
4324  diam3b = diam3b/100.
4325  CALL popcontrol1b(branch)
4326  IF (branch .NE. 0) wb = wb + 0.99*veb/100.
4327  diam3b = diam3b + sqrt(1000./pr)*2483.8*wb
4328  diam3b = 0.664*diam3b
4329  CALL popcontrol1b(branch)
4330  IF (branch .NE. 0) diam3b = 0.0_8
4331  rain_dayb = 0.0_8
4332  DO iqd=7,1,-1
4333  CALL popcontrol1b(branch)
4334  IF (branch .NE. 0) THEN
4335  rain_dayb = rain_dayb + slopr*diam3b
4336  CALL popreal8(slopr)
4337  diam3b = 0.0_8
4338  END IF
4339  END DO
4340  rainb = 24.*3600.*rain_dayb
4341 END SUBROUTINE marshpalm_b
4342 
4343 ! Differentiation of dqsat_bac in reverse (adjoint) mode:
4344 ! gradient of useful results: temp dqsi qssi
4345 ! with respect to varying inputs: temp dqsi qssi
4346 SUBROUTINE dqsat_bac_b(dqsi, dqsib, qssi, qssib, temp, tempb, plo, im, &
4347 & jm, lm, estblx, cons_h2omw, cons_airmw)
4348  IMPLICIT NONE
4349 !Inputs
4350  INTEGER :: im, jm, lm
4351  REAL*8, DIMENSION(im, jm, lm) :: temp, plo
4352  REAL*8, DIMENSION(im, jm, lm) :: tempb
4353  REAL*8 :: estblx(:)
4354  REAL*8 :: cons_h2omw, cons_airmw
4355 !Outputs
4356  REAL*8, DIMENSION(im, jm, lm) :: dqsi, qssi
4357  REAL*8, DIMENSION(im, jm, lm) :: dqsib, qssib
4358 !Locals
4359  REAL*8, PARAMETER :: max_mixing_ratio=1.0
4360  REAL*8 :: esfac
4361  INTEGER :: i, j, k
4362  REAL*8 :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
4363  REAL*8 :: tlb, ttb, tib, dqsatb, qsatb, qqb, ddb
4364  INTEGER :: it
4365  INTEGER, PARAMETER :: degsubs=100
4366  REAL*8, PARAMETER :: tmintbl=150.0, tmaxtbl=333.0
4367  INTEGER, PARAMETER :: tablesize=nint(tmaxtbl-tmintbl)*degsubs+1
4368  INTRINSIC nint
4369  INTRINSIC int
4370  INTEGER :: branch
4371  REAL*8 :: temp0
4372  esfac = cons_h2omw/cons_airmw
4373  DO k=1,lm
4374  DO j=1,jm
4375  DO i=1,im
4376  tl = temp(i, j, k)
4377  pl = plo(i, j, k)
4378  pp = pl*100.0
4379  IF (tl .LE. tmintbl) THEN
4380  ti = tmintbl
4381  CALL pushcontrol2b(0)
4382  ELSE IF (tl .GE. tmaxtbl - .001) THEN
4383  ti = tmaxtbl - .001
4384  CALL pushcontrol2b(1)
4385  ELSE
4386  ti = tl
4387  CALL pushcontrol2b(2)
4388  END IF
4389  tt = (ti-tmintbl)*degsubs + 1
4390  it = int(tt)
4391  CALL pushreal8(dqq)
4392  dqq = estblx(it+1) - estblx(it)
4393  CALL pushreal8(qq)
4394  qq = (tt-it)*dqq + estblx(it)
4395  IF (pp .LE. qq) THEN
4396  CALL pushcontrol1b(0)
4397  ELSE
4398  CALL pushcontrol1b(1)
4399  END IF
4400  END DO
4401  END DO
4402  END DO
4403  DO k=lm,1,-1
4404  DO j=jm,1,-1
4405  DO i=im,1,-1
4406  qsatb = qssib(i, j, k)
4407  qssib(i, j, k) = 0.0_8
4408  dqsatb = dqsib(i, j, k)
4409  dqsib(i, j, k) = 0.0_8
4410  CALL popcontrol1b(branch)
4411  IF (branch .EQ. 0) THEN
4412  qqb = 0.0_8
4413  ELSE
4414  pl = plo(i, j, k)
4415  pp = pl*100.0
4416  dd = 1.0/(pp-(1.0-esfac)*qq)
4417  ddb = esfac*qq*qsatb + esfac*degsubs*dqq*pp*2*dd*dqsatb
4418  temp0 = pp - (-esfac+1.0)*qq
4419  qqb = (1.0-esfac)*ddb/temp0**2 + esfac*dd*qsatb
4420  END IF
4421  CALL popreal8(qq)
4422  ttb = dqq*qqb
4423  CALL popreal8(dqq)
4424  tib = degsubs*ttb
4425  CALL popcontrol2b(branch)
4426  IF (branch .EQ. 0) THEN
4427  tlb = 0.0_8
4428  ELSE IF (branch .EQ. 1) THEN
4429  tlb = 0.0_8
4430  ELSE
4431  tlb = tib
4432  END IF
4433  tempb(i, j, k) = tempb(i, j, k) + tlb
4434  END DO
4435  END DO
4436  END DO
4437 END SUBROUTINE dqsat_bac_b
4438 
4439 ! Differentiation of dqsats_bac in reverse (adjoint) mode:
4440 ! gradient of useful results: temp dqsi qssi
4441 ! with respect to varying inputs: temp
4442 SUBROUTINE dqsats_bac_b(dqsi, dqsib, qssi, qssib, temp, tempb, plo, &
4443 & estblx, cons_h2omw, cons_airmw)
4444  IMPLICIT NONE
4445 !Inputs
4446  REAL*8 :: temp, plo
4447  REAL*8 :: tempb
4448  REAL*8 :: estblx(:)
4449  REAL*8 :: cons_h2omw, cons_airmw
4450 !Outputs
4451  REAL*8 :: dqsi, qssi
4452  REAL*8 :: dqsib, qssib
4453 !Locals
4454  REAL*8, PARAMETER :: max_mixing_ratio=1.0
4455  REAL*8 :: esfac
4456  REAL*8 :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
4457  REAL*8 :: tlb, ttb, tib, dqsatb, qsatb, qqb, ddb
4458  INTEGER :: it
4459  INTEGER, PARAMETER :: degsubs=100
4460  REAL*8, PARAMETER :: tmintbl=150.0, tmaxtbl=333.0
4461  INTEGER, PARAMETER :: tablesize=nint(tmaxtbl-tmintbl)*degsubs+1
4462  INTRINSIC nint
4463  INTRINSIC int
4464  INTEGER :: branch
4465  REAL*8 :: temp0
4466  esfac = cons_h2omw/cons_airmw
4467  tl = temp
4468  pl = plo
4469  pp = pl*100.0
4470  IF (tl .LE. tmintbl) THEN
4471  ti = tmintbl
4472  CALL pushcontrol2b(0)
4473  ELSE IF (tl .GE. tmaxtbl - .001) THEN
4474  ti = tmaxtbl - .001
4475  CALL pushcontrol2b(1)
4476  ELSE
4477  ti = tl
4478  CALL pushcontrol2b(2)
4479  END IF
4480  tt = (ti-tmintbl)*degsubs + 1
4481  it = int(tt)
4482  dqq = estblx(it+1) - estblx(it)
4483  qq = (tt-it)*dqq + estblx(it)
4484  IF (pp .LE. qq) THEN
4485  CALL pushcontrol1b(0)
4486  ELSE
4487  dd = 1.0/(pp-(1.0-esfac)*qq)
4488  CALL pushcontrol1b(1)
4489  END IF
4490  qsatb = qssib
4491  dqsatb = dqsib
4492  CALL popcontrol1b(branch)
4493  IF (branch .EQ. 0) THEN
4494  qqb = 0.0_8
4495  ELSE
4496  temp0 = pp - (-esfac+1.0)*qq
4497  ddb = esfac*qq*qsatb + esfac*degsubs*dqq*pp*2*dd*dqsatb
4498  qqb = (1.0-esfac)*ddb/temp0**2 + esfac*dd*qsatb
4499  END IF
4500  ttb = dqq*qqb
4501  tib = degsubs*ttb
4502  CALL popcontrol2b(branch)
4503  IF (branch .EQ. 0) THEN
4504  tlb = 0.0_8
4505  ELSE IF (branch .EQ. 1) THEN
4506  tlb = 0.0_8
4507  ELSE
4508  tlb = tib
4509  END IF
4510  tempb = tempb + tlb
4511 END SUBROUTINE dqsats_bac_b
4512 
4513 END MODULE cloud_ad
cloud_ad::convec_src_b
subroutine convec_src_b(dt, mass, imass, te, teb, qv, qvb, dcf, dcfb, dmf, dmfb, qla, qlab, qia, qiab, af, afb, qs, qsb, cons_alhs, cons_alhl, cons_cp, t_ice_all, t_ice_max, icefrpwr)
Definition: cloud_ad.F90:1574
cloud_ad::autoconversion_ls_b
subroutine autoconversion_ls_b(dt, qc, qcb, qp, qpb, te, teb, pl, f, fb, sundqv2, sundqv3, sundqt1, c_00, lwcrit, dzet)
Definition: cloud_ad.F90:2812
popcontrol2b
subroutine popcontrol2b(cc)
Definition: adBuffer.f:146
cloud_ad::cloud_tidy_b
subroutine cloud_tidy_b(qv, qvb, te, teb, qlc, qlcb, qic, qicb, cf, cfb, qla, qlab, qia, qiab, af, afb, cons_alhl, cons_alhs, cons_cp)
Definition: cloud_ad.F90:1388
cloud_ad::ldradius_b
subroutine ldradius_b(pl, te, teb, qcl, qclb, nn, rho_w, radius, radiusb, cons_rgas, cons_pi)
Definition: cloud_ad.F90:2779
cloud_ad::ice_settlefall_cnv_b
subroutine ice_settlefall_cnv_b(wxr, qi, qib, pl, te, teb, f, fb, cons_rgas, khu, khl, k, dt, dz, dzb, qp, qpb, anv_icefall_c)
Definition: cloud_ad.F90:3549
cloud::autoconversion_cnv
subroutine, public autoconversion_cnv(DT, QC, QP, TE, PL, F, SUNDQV2, SUNDQV3, SUNDQT1, C_00, LWCRIT, DZET)
Definition: cloud.F90:1774
popreal8array
void popreal8array(double *x, int n)
Definition: adStack.c:375
cloud_ad::cloud_driver_b
subroutine, public cloud_driver_b(dt, im, jm, lm, th, thb, q, qb, ple, cnv_dqldt, cnv_dqldtb, cnv_mfd, cnv_mfdb, cnv_prc3, cnv_prc3b, cnv_updf, cnv_updfb, qi_ls, qi_lsb, ql_ls, ql_lsb, qi_con, qi_conb, ql_con, ql_conb, cf_ls, cf_lsb, cf_con, cf_conb, frland, physparams, estblx, khu, khl, cons_runiv, cons_kappa, cons_airmw, cons_h2omw, cons_grav, cons_alhl, cons_alhf, cons_pi, cons_rgas, cons_cp, cons_vireps, cons_alhs, cons_tice, cons_rvap, cons_p00, do_moist_physics)
Definition: cloud_ad.F90:31
cloud::ldradius
subroutine, public ldradius(PL, TE, QCL, NN, RHO_W, RADIUS, CONS_RGAS, CONS_PI)
Definition: cloud.F90:1643
cloud::pdfcondensate
subroutine, public pdfcondensate(flag, qtmean4, sigmaqt14, sigmaqt24, qstar4, condensate4)
Definition: cloud.F90:1428
cloud::dqsat_bac
subroutine, public dqsat_bac(DQSi, QSSi, TEMP, PLO, im, jm, lm, ESTBLX, CONS_H2OMW, CONS_AIRMW)
Definition: cloud.F90:2457
cloud::precipandevap
subroutine, public precipandevap(K, KTOP, LM, DT, FRLAND, RHCR3, QPl, QPi, QCl, QCi, TE, QV, mass, imass, PL, dZE, QDDF3, AA, BB, AREA, PFl_above_in, PFl_above_out, PFi_above_in, PFi_above_out, EVAP_DD_above_in, EVAP_DD_above_out, SUBL_DD_above_in, SUBL_DD_above_out, ENVFC, DDRFC, CONS_ALHF, CONS_ALHS, CONS_ALHL, CONS_CP, CONS_TICE, CONS_H2OMW, CONS_AIRMW, REVAP_OFF_P, C_ACC, C_EV_R, C_EV_S, RHO_W, ESTBLX)
Definition: cloud.F90:2219
cloud::evap_cnv
subroutine, public evap_cnv(DT, RHCR, PL, TE, QV, QL, QI, F, XF, QS, RHO_W, CLD_EVP_EFF, CONS_H2OMW, CONS_AIRMW, CONS_ALHL, CONS_RVAP, CONS_RGAS, CONS_PI, CONS_CP)
Definition: cloud.F90:1525
cloud_ad
Definition: cloud_ad.F90:1
cloud_ad::pdfcondensate_b
subroutine pdfcondensate_b(flag, qtmean4, qtmean4b, sigmaqt14, sigmaqt14b, sigmaqt24, sigmaqt24b, qstar4, qstar4b, condensate4, condensate4b)
Definition: cloud_ad.F90:2247
cloud::ice_settlefall_ls
subroutine, public ice_settlefall_ls(WXR, QI, PL, TE, F, CONS_RGAS, KHu, KHl, k, DT, DZ, QP, LS_ICEFALL_C)
Definition: cloud.F90:2145
cloud::ice_settlefall_cnv
subroutine, public ice_settlefall_cnv(WXR, QI, PL, TE, F, CONS_RGAS, KHu, KHl, k, DT, DZ, QP, ANV_ICEFALL_C)
Definition: cloud.F90:2082
cloud_ad::cons_sundq3_b
subroutine cons_sundq3_b(temp, tempb, rate2, rate3, te1, f2, f2b, f3)
Definition: cloud_ad.F90:3420
pushcontrol1b
subroutine pushcontrol1b(cc)
Definition: adBuffer.f:115
cloud_tl::ls_cloud_d
subroutine, public ls_cloud_d(dt, alpha, pdfshape, pl, te, ted, qv, qvd, qcl, qcld, qal, qald, qci, qcid, qai, qaid, cf, cfd, af, afd, cons_alhl, cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw, t_ice_all, t_ice_max, icefrpwr, estblx, cloud_pertmod, dmp)
Definition: cloud_tl.F90:1157
cloud_ad::ice_settlefall_ls_b
subroutine ice_settlefall_ls_b(wxr, qi, qib, pl, te, teb, f, fb, cons_rgas, khu, khl, k, dt, dz, dzb, qp, qpb, ls_icefall_c)
Definition: cloud_ad.F90:3666
cloud::pdffrac
subroutine, public pdffrac(flag, qtmean, sigmaqt1, sigmaqt2, qstar, clfrac)
Definition: cloud.F90:1353
cloud
Definition: cloud.F90:1
cloud::cons_alhx
subroutine, public cons_alhx(T, ALHX3, T_ICE_MAX, T_ICE_ALL, CONS_ALHS, CONS_ALHL)
Definition: cloud.F90:1978
pushcontrol2b
subroutine pushcontrol2b(cc)
Definition: adBuffer.f:140
cloud_ad::dqsat_bac_b
subroutine dqsat_bac_b(dqsi, dqsib, qssi, qssib, temp, tempb, plo, im, jm, lm, estblx, cons_h2omw, cons_airmw)
Definition: cloud_ad.F90:4348
pushreal8array
void pushreal8array(double *x, int n)
Definition: adStack.c:372
cloud_ad::evap_cnv_b
subroutine evap_cnv_b(dt, rhcr, pl, te, teb, qv, qvb, ql, qlb, qi, qib, f, fb, xf, qs, qsb, rho_w, cld_evp_eff, cons_h2omw, cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp)
Definition: cloud_ad.F90:2500
popreal8
subroutine popreal8(x)
Definition: adBuffer.f:820
cloud::pdf_width
subroutine, public pdf_width(PP, FRLAND, maxrhcrit, maxrhcritland, turnrhcrit, minrhcrit, pi_0, ALPHA)
Definition: cloud.F90:1046
cloud_tl
Definition: cloud_tl.F90:1
cloud_ad::meltfreeze_b
subroutine meltfreeze_b(dt, te, teb, ql, qlb, qi, qib, t_ice_all, t_ice_max, icefrpwr, cons_alhl, cons_alhs, cons_cp)
Definition: cloud_ad.F90:1496
cloud::cons_microphys
subroutine, public cons_microphys(TEMP, PR, Q_SAT, AA, BB, CONS_H2OMW, CONS_AIRMW, CONS_RVAP, ALHX3)
Definition: cloud.F90:1951
cloud_ad::get_ice_fraction_b
subroutine get_ice_fraction_b(temp, tempb, t_ice_all, t_ice_max, icefrpwr, icefrct, icefrctb)
Definition: cloud_ad.F90:3365
cloud::meltfreeze
subroutine, public meltfreeze(DT, TE, QL, QI, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, CONS_ALHL, CONS_ALHS, CONS_CP)
Definition: cloud.F90:923
cloud_ad::dqsats_bac_b
subroutine dqsats_bac_b(dqsi, dqsib, qssi, qssib, temp, tempb, plo, estblx, cons_h2omw, cons_airmw)
Definition: cloud_ad.F90:4444
cloud::convec_src
subroutine, public convec_src(DT, MASS, iMASS, TE, QV, DCF, DMF, QLA, QIA, AF, QS, CONS_ALHS, CONS_ALHL, CONS_CP, T_ICE_ALL, T_ICE_MAX, ICEFRPWR)
Definition: cloud.F90:969
cloud::cloud_tidy
subroutine, public cloud_tidy(QV, TE, QLC, QIC, CF, QLA, QIA, AF, CONS_ALHL, CONS_ALHS, CONS_CP)
Definition: cloud.F90:852
cloud::cons_sundq3
subroutine, public cons_sundq3(TEMP, RATE2, RATE3, TE1, F2, F3)
Definition: cloud.F90:1909
cloud::subl_cnv
subroutine, public subl_cnv(DT, RHCR, PL, TE, QV, QL, QI, F, XF, QS, RHO_W, CLD_EVP_EFF, CONS_H2OMW, CONS_AIRMW, CONS_ALHL, CONS_RVAP, CONS_RGAS, CONS_PI, CONS_CP, CONS_ALHS)
Definition: cloud.F90:1584
popcontrol3b
subroutine popcontrol3b(cc)
Definition: adBuffer.f:175
cloud::dqsats_bac
subroutine, public dqsats_bac(DQSi, QSSi, TEMP, PLO, ESTBLX, CONS_H2OMW, CONS_AIRMW)
Definition: cloud.F90:2529
pushreal8
subroutine pushreal8(x)
Definition: adBuffer.f:763
popcontrol1b
subroutine popcontrol1b(cc)
Definition: adBuffer.f:120
cloud_ad::subl_cnv_b
subroutine subl_cnv_b(dt, rhcr, pl, te, teb, qv, qvb, ql, qlb, qi, qib, f, fb, xf, qs, qsb, rho_w, cld_evp_eff, cons_h2omw, cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp, cons_alhs)
Definition: cloud_ad.F90:2640
cloud::autoconversion_ls
subroutine, public autoconversion_ls(DT, QC, QP, TE, PL, F, SUNDQV2, SUNDQV3, SUNDQT1, C_00, LWCRIT, DZET)
Definition: cloud.F90:1660
cloud_ad::pdffrac_b
subroutine pdffrac_b(flag, qtmean, qtmeanb, sigmaqt1, sigmaqt1b, sigmaqt2, sigmaqt2b, qstar, qstarb, clfrac, clfracb)
Definition: cloud_ad.F90:2099
max
#define max(a, b)
Definition: mosaic_util.h:33
popcontrol4b
subroutine popcontrol4b(cc)
Definition: adBuffer.f:207
cloud_ad::marshpalm_b
subroutine marshpalm_b(rain, rainb, pr, diam3, diam3b, ntotal, w, wb, ve, veb)
Definition: cloud_ad.F90:4261
cloud::marshpalm
subroutine, public marshpalm(RAIN, PR, DIAM3, NTOTAL, W, VE)
Definition: cloud.F90:2003
pushcontrol4b
subroutine pushcontrol4b(cc)
Definition: adBuffer.f:199
cloud::get_ice_fraction
subroutine, public get_ice_fraction(TEMP, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, ICEFRCT)
Definition: cloud.F90:1881
pushcontrol3b
subroutine pushcontrol3b(cc)
Definition: adBuffer.f:168
cloud_ad::autoconversion_cnv_b
subroutine autoconversion_cnv_b(dt, qc, qcb, qp, qpb, te, teb, pl, f, fb, sundqv2, sundqv3, sundqt1, c_00, lwcrit, dzet)
Definition: cloud_ad.F90:3097
min
#define min(a, b)
Definition: mosaic_util.h:32
cloud_ad::ls_cloud_b
subroutine ls_cloud_b(dt, alpha, pdfshape, pl, te, teb, qv, qvb, qcl, qclb, qal, qalb, qci, qcib, qai, qaib, cf, cfb, af, afb, cons_alhl, cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw, t_ice_all, t_ice_max, icefrpwr, estblx, cloud_pertmod, dmp)
Definition: cloud_ad.F90:1671
cloud::ls_cloud
subroutine, public ls_cloud(DT, ALPHA, PDFSHAPE, PL, TE, QV, QCl, QAl, QCi, QAi, CF, AF, CONS_ALHL, CONS_ALHF, CONS_ALHS, CONS_CP, CONS_H2OMW, CONS_AIRMW, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, ESTBLX, cloud_pertmod, dmp)
Definition: cloud.F90:1107
cloud_ad::precipandevap_b
subroutine precipandevap_b(k, ktop, lm, dt, frland, rhcr3, qpl, qplb, qpi, qpib, qcl, qclb, qci, te, teb, qv, qvb, mass, imass, pl, dze, dzeb, qddf3, qddf3b, aa, aab, bb, bbb, area, areab, pfl_above_in, pfl_above_inb, pfl_above_out, pfl_above_outb, pfi_above_in, pfi_above_inb, pfi_above_out, pfi_above_outb, evap_dd_above_in, evap_dd_above_inb, evap_dd_above_out, evap_dd_above_outb, subl_dd_above_in, subl_dd_above_inb, subl_dd_above_out, subl_dd_above_outb, envfc, ddrfc, cons_alhf, cons_alhs, cons_alhl, cons_cp, cons_tice, cons_h2omw, cons_airmw, revap_off_p, c_acc, c_ev_r, c_ev_s, rho_w, estblx)
Definition: cloud_ad.F90:3817
cloud_ad::cons_microphys_b
subroutine cons_microphys_b(temp, tempb1, pr, q_sat, q_satb, aa, aab, bb, bbb, cons_h2omw, cons_airmw, cons_rvap, alhx3, alhx3b)
Definition: cloud_ad.F90:3477
cloud_ad::cons_alhx_b
subroutine cons_alhx_b(t, tb, alhx3, alhx3b, t_ice_max, t_ice_all, cons_alhs, cons_alhl)
Definition: cloud_ad.F90:3515