FV3 Bundle
irrad_ad.F90
Go to the documentation of this file.
1 module irrad_ad
2 
3 USE irradmod
4 
5 IMPLICIT NONE
6 
7 PRIVATE
8 PUBLIC :: irrad_b
9 
10 contains
11 
12 SUBROUTINE irrad_b(np, ple_dev, ta_dev, ta_devb, wa_dev, wa_devb, oa_dev&
13 & , oa_devb, tb_dev, tb_devb, co2, trace, n2o_dev, ch4_dev, cfc11_dev, &
14 & cfc12_dev, cfc22_dev, cwc_dev, cwc_devb, fcld_dev, fcld_devb, ict, icb&
15 & , reff_dev, reff_devb, ns, fs_dev, tg_dev, eg_dev, tv_dev, ev_dev, &
16 & rv_dev, na, nb, taua_dev, taua_devb, ssaa_dev, ssaa_devb, asya_dev, &
17 & asya_devb, flxu_dev, flxu_devb, flxd_dev, flxd_devb, dfdts_dev, &
18 & dfdts_devb, aib_ir, awb_ir, aiw_ir, aww_ir, aig_ir, awg_ir, xkw, xke, &
19 & mw, aw, bw, pm, fkw, gkw, cb, dcb, w11, w12, w13, p11, p12, p13, dwe, &
20 & dpe, c1, c2, c3, oo1, oo2, oo3, h11, h12, h13, h21, h22, h23, h81, h82&
21 & , h83)
22  IMPLICIT NONE
23 !Radiation constants, these need to be inputs for the autodiff tool
24  REAL*8, INTENT(IN) :: aib_ir(3, 10), awb_ir(4, 10), aiw_ir(4, 10)
25  REAL*8, INTENT(IN) :: aww_ir(4, 10), aig_ir(4, 10), awg_ir(4, 10)
26  INTEGER, INTENT(IN) :: mw(9)
27  REAL*8, INTENT(IN) :: xkw(9), xke(9), aw(9), bw(9), pm(9)
28  REAL*8, INTENT(IN) :: fkw(6, 9), gkw(6, 3), cb(6, 10), dcb(5, 10)
29  REAL*8, INTENT(IN) :: w11, w12, w13, p11, p12
30  REAL*8, INTENT(IN) :: p13, dwe, dpe
31  REAL*8, INTENT(IN) :: c1(26, 30), c2(26, 30), c3(26, 30)
32  REAL*8, INTENT(IN) :: oo1(26, 21), oo2(26, 21), oo3(26, 21)
33  REAL*8, INTENT(IN) :: h11(26, 31), h12(26, 31), h13(26, 31)
34  REAL*8, INTENT(IN) :: h21(26, 31), h22(26, 31), h23(26, 31)
35  REAL*8, INTENT(IN) :: h81(26, 31), h82(26, 31), h83(26, 31)
36 !----- INPUTS -----
37  INTEGER, INTENT(IN) :: np, ict, icb, ns, na, nb
38  LOGICAL, INTENT(IN) :: trace
39  REAL*8, INTENT(IN) :: co2
40 !Rank 2 inputs
41  REAL*8, INTENT(IN) :: tb_dev
42  REAL*8 :: tb_devb
43 !Rank 3 (Prognostic variables and tracers)
44  REAL*8, DIMENSION(np), INTENT(IN) :: ta_dev, wa_dev, oa_dev, fcld_dev
45  REAL*8, DIMENSION(np) :: ta_devb, wa_devb, oa_devb, fcld_devb
46  REAL*8, DIMENSION(np), INTENT(IN) :: n2o_dev, ch4_dev, cfc11_dev, &
47 & cfc12_dev, cfc22_dev
48  REAL*8, DIMENSION(np + 1), INTENT(IN) :: ple_dev
49 !Rank 3 (surface types)
50  REAL*8, DIMENSION(ns), INTENT(IN) :: fs_dev, tg_dev, tv_dev
51  REAL*8, DIMENSION(ns, 10), INTENT(IN) :: eg_dev, ev_dev, rv_dev
52 !Rank 3 (diagnostic cloud parts)
53  REAL*8, DIMENSION(np, 4), INTENT(IN) :: cwc_dev, reff_dev
54  REAL*8, DIMENSION(np, 4) :: cwc_devb, reff_devb
55 !Rank 3 (aerosols)
56  REAL*8, DIMENSION(np, nb), INTENT(INOUT) :: taua_dev, ssaa_dev, &
57 & asya_dev
58  REAL*8, DIMENSION(np, nb), INTENT(INOUT) :: taua_devb
59 !----- OUPUTS -----
60  REAL*8, DIMENSION(np + 1) :: flxu_dev
61  REAL*8, DIMENSION(np+1) :: flxu_devb
62  REAL*8, DIMENSION(np + 1) :: flxd_dev
63  REAL*8, DIMENSION(np+1) :: flxd_devb
64  REAL*8, DIMENSION(np + 1) :: dfdts_dev
65  REAL*8, DIMENSION(np+1) :: dfdts_devb
66 !----- LOCALS -----
67  REAL*8, PARAMETER :: cons_grav=9.80665
68  INTEGER, PARAMETER :: nx1=26
69  INTEGER, PARAMETER :: no1=21
70  INTEGER, PARAMETER :: nc1=30
71  INTEGER, PARAMETER :: nh1=31
72 !Temporary arrays
73  REAL*8 :: pa(0:np), dt(0:np)
74  REAL*8 :: dtb(0:np)
75  REAL*8 :: x1, x2, x3
76  REAL*8 :: x1b, x2b, x3b
77  REAL*8 :: dh2o(0:np), dcont(0:np), dco2(0:np), do3(0:np)
78  REAL*8 :: dh2ob(0:np), dcontb(0:np), dco2b(0:np), do3b(0:np)
79  REAL*8 :: dn2o(0:np), dch4(0:np)
80  REAL*8 :: df11(0:np), df12(0:np), df22(0:np)
81  REAL*8 :: th2o(6), tcon(3), tco2(6)
82  REAL*8 :: th2ob(6), tconb(3), tco2b(6)
83  REAL*8 :: tn2o(4), tch4(4), tcom(6)
84  REAL*8 :: tn2ob(4), tch4b(4), tcomb(6)
85  REAL*8 :: tf11, tf12, tf22
86  REAL*8 :: tf11b, tf12b, tf22b
87  REAL*8 :: blayer(0:np+1), blevel(0:np+1)
88  REAL*8 :: blayerb(0:np+1), blevelb(0:np+1)
89  REAL*8 :: bd(0:np+1), bu(0:np+1)
90  REAL*8 :: bdb(0:np+1), bub(0:np+1)
91  REAL*8 :: bs, dbs, rflxs
92  REAL*8 :: dp(0:np)
93  REAL*8 :: trant, tranal
94  REAL*8 :: trantb, tranalb
95  REAL*8 :: transfc(0:np+1)
96  REAL*8 :: transfcb(0:np+1)
97  REAL*8 :: flxu(0:np+1), flxd(0:np+1)
98  REAL*8 :: flxub(0:np+1), flxdb(0:np+1)
99  REAL*8 :: taerlyr(0:np)
100  REAL*8 :: taerlyrb(0:np)
101 !OVERCAST
102  INTEGER :: ncld(3)
103  INTEGER :: icx(0:np)
104 !OVERCAST
105  INTEGER :: idx, rc
106  INTEGER :: k, l, ip, iw, ibn, ik, iq, isb, k1, k2, ne
107  REAL*8 :: enn(0:np)
108  REAL*8 :: ennb(0:np)
109  REAL*8 :: cldhi, cldmd, cldlw, tcldlyr(0:np), fclr
110  REAL*8 :: cldhib, cldmdb, cldlwb, tcldlyrb(0:np), fclrb
111  REAL*8 :: x, xx, yy, p1, a1, b1, fk1, a2, b2, fk2
112  REAL*8 :: xxb, yyb
113  REAL*8 :: w1, ff
114  REAL*8 :: ffb
115  LOGICAL :: oznbnd, co2bnd, h2otable, conbnd, n2obnd
116  LOGICAL :: ch4bnd, combnd, f11bnd, f12bnd, f22bnd, b10bnd
117  LOGICAL :: do_aerosol
118 !Temp arrays and variables for consolidation of tables
119  INTEGER, PARAMETER :: max_num_tables=17
120  REAL*8 :: exptbl(0:np, max_num_tables)
121  REAL*8 :: exptblb(0:np, max_num_tables)
122  TYPE band_table
123  INTEGER :: start
124  INTEGER :: end
125  END TYPE band_table
126  TYPE(band_table) :: h2oexp
127  TYPE(band_table) :: conexp
128  TYPE(band_table) :: co2exp
129  TYPE(band_table) :: n2oexp
130  TYPE(band_table) :: ch4exp
131  TYPE(band_table) :: comexp
132  TYPE(band_table) :: f11exp
133  TYPE(band_table) :: f12exp
134  TYPE(band_table) :: f22exp
135 !Variables for new getirtau routine
136  REAL*8 :: dp_pa(np)
137  REAL*8 :: fcld_col(np)
138  REAL*8 :: fcld_colb(np)
139  REAL*8 :: reff_col(np, 4)
140  REAL*8 :: reff_colb(np, 4)
141  REAL*8 :: cwc_col(np, 4)
142  REAL*8 :: cwc_colb(np, 4)
143  REAL*8 :: h2oexp_tmp(0:np, 5), conexp_tmp(0:np), co2exp_tmp(0:np, 6), &
144 & n2oexp_tmp(0:np, 2)
145  REAL*8 :: h2oexp_tmpb(0:np, 5), co2exp_tmpb(0:np, 6), n2oexp_tmpb(0:np&
146 & , 2)
147  INTRINSIC max
148  INTRINSIC exp
149  INTRINSIC min
150  INTRINSIC log
151  INTEGER :: arg1
152  INTEGER :: branch
153  INTEGER :: ad_from
154  INTEGER :: ad_count
155  INTEGER :: i
156  REAL*8, DIMENSION(np, nb), INTENT(INOUT) :: asya_devb
157  REAL*8, DIMENSION(np, nb), INTENT(INOUT) :: ssaa_devb
158  REAL*8 :: temp2
159  REAL*8 :: temp1
160  REAL*8 :: temp0
161  REAL*8 :: tempb6
162  REAL*8 :: tempb5
163  REAL*8 :: tempb4
164  REAL*8 :: tempb3
165  REAL*8 :: tempb2
166  REAL*8 :: tempb1
167  REAL*8 :: tempb0
168  REAL*8 :: tempb
169  REAL*8 :: temp
170 !BEGIN CALCULATIONS ...
171 !compute layer pressure (pa) and layer temperature minus 250K (dt)
172  DO k=1,np
173  pa(k) = 0.5*(ple_dev(k+1)+ple_dev(k))*0.01
174  dp(k) = (ple_dev(k+1)-ple_dev(k))*0.01
175 ! dp in Pascals for getirtau
176  dp_pa(k) = ple_dev(k+1) - ple_dev(k)
177  dt(k) = ta_dev(k) - 250.0
178 !compute layer absorber amount
179 !dh2o : water vapor amount (g/cm**2)
180 !dcont: scaled water vapor amount for continuum absorption
181 ! (g/cm**2)
182 !dco2 : co2 amount (cm-atm)stp
183 !do3 : o3 amount (cm-atm)stp
184 !dn2o : n2o amount (cm-atm)stp
185 !dch4 : ch4 amount (cm-atm)stp
186 !df11 : cfc11 amount (cm-atm)stp
187 !df12 : cfc12 amount (cm-atm)stp
188 !df22 : cfc22 amount (cm-atm)stp
189 !the factor 1.02 is equal to 1000/980
190 !factors 789 and 476 are for unit conversion
191 !the factor 0.001618 is equal to 1.02/(.622*1013.25)
192 !the factor 6.081 is equal to 1800/296
193  dh2o(k) = 1.02*wa_dev(k)*dp(k)
194  do3(k) = 476.*oa_dev(k)*dp(k)
195  dco2(k) = 789.*co2*dp(k)
196  dch4(k) = 789.*ch4_dev(k)*dp(k)
197  dn2o(k) = 789.*n2o_dev(k)*dp(k)
198  df11(k) = 789.*cfc11_dev(k)*dp(k)
199  df12(k) = 789.*cfc12_dev(k)*dp(k)
200  df22(k) = 789.*cfc22_dev(k)*dp(k)
201  IF (dh2o(k) .LT. 1.e-10) THEN
202  dh2o(k) = 1.e-10
203  CALL pushcontrol1b(0)
204  ELSE
205  CALL pushcontrol1b(1)
206  dh2o(k) = dh2o(k)
207  END IF
208  IF (do3(k) .LT. 1.e-6) THEN
209  do3(k) = 1.e-6
210  CALL pushcontrol1b(0)
211  ELSE
212  CALL pushcontrol1b(1)
213  do3(k) = do3(k)
214  END IF
215  IF (dco2(k) .LT. 1.e-4) THEN
216  dco2(k) = 1.e-4
217  ELSE
218  dco2(k) = dco2(k)
219  END IF
220 !Compute scaled water vapor amount for h2o continuum absorption
221 !following eq. (4.21).
222  xx = pa(k)*0.001618*wa_dev(k)*wa_dev(k)*dp(k)
223  dcont(k) = xx*exp(1800./ta_dev(k)-6.081)
224 !Fill the reff, cwc, and fcld for the column
225  fcld_col(k) = fcld_dev(k)
226  DO l=1,4
227  reff_col(k, l) = reff_dev(k, l)
228  cwc_col(k, l) = cwc_dev(k, l)
229  END DO
230  END DO
231  IF (ple_dev(1)*0.01 .LT. 0.005) THEN
232  CALL pushreal8(dp(0))
233  dp(0) = 0.005
234  CALL pushcontrol1b(0)
235  ELSE
236  CALL pushreal8(dp(0))
237  dp(0) = ple_dev(1)*0.01
238  CALL pushcontrol1b(1)
239  END IF
240  CALL pushreal8(pa(0))
241  pa(0) = 0.5*dp(0)
242  dt(0) = ta_dev(1) - 250.0
243  dh2o(0) = 1.02*wa_dev(1)*dp(0)
244  do3(0) = 476.*oa_dev(1)*dp(0)
245  dco2(0) = 789.*co2*dp(0)
246  dch4(0) = 789.*ch4_dev(1)*dp(0)
247  dn2o(0) = 789.*n2o_dev(1)*dp(0)
248  df11(0) = 789.*cfc11_dev(1)*dp(0)
249  df12(0) = 789.*cfc12_dev(1)*dp(0)
250  df22(0) = 789.*cfc22_dev(1)*dp(0)
251  IF (dh2o(0) .LT. 1.e-10) THEN
252  dh2o(0) = 1.e-10
253  CALL pushcontrol1b(0)
254  ELSE
255  CALL pushcontrol1b(1)
256  dh2o(0) = dh2o(0)
257  END IF
258  IF (do3(0) .LT. 1.e-6) THEN
259  do3(0) = 1.e-6
260  CALL pushcontrol1b(0)
261  ELSE
262  CALL pushcontrol1b(1)
263  do3(0) = do3(0)
264  END IF
265  IF (dco2(0) .LT. 1.e-4) THEN
266  dco2(0) = 1.e-4
267  ELSE
268  dco2(0) = dco2(0)
269  END IF
270  xx = pa(0)*0.001618*wa_dev(1)*wa_dev(1)*dp(0)
271  dcont(0) = xx*exp(1800./ta_dev(1)-6.081)
272 !The surface (np+1) is treated as a layer filled with black clouds.
273 !transfc is the transmittance between the surface and a pressure
274 !level.
275  transfc(np+1) = 1.0
276  CALL pushinteger4(ibn)
277  ad_count = 1
278 !Integration over spectral bands
279  DO ibn=1,10
280  IF (ibn .EQ. 10 .AND. (.NOT.trace)) THEN
281  GOTO 100
282  ELSE
283 !if h2otable, compute h2o (line) transmittance using table look-up.
284 !if conbnd, compute h2o (continuum) transmittance in bands 2-7.
285 !if co2bnd, compute co2 transmittance in band 3.
286 !if oznbnd, compute o3 transmittance in band 5.
287 !if n2obnd, compute n2o transmittance in bands 6 and 7.
288 !if ch4bnd, compute ch4 transmittance in bands 6 and 7.
289 !if combnd, compute co2-minor transmittance in bands 4 and 5.
290 !if f11bnd, compute cfc11 transmittance in bands 4 and 5.
291 !if f12bnd, compute cfc12 transmittance in bands 4 and 6.
292 !if f22bnd, compute cfc22 transmittance in bands 4 and 6.
293 !if b10bnd, compute flux reduction due to n2o in band 10.
294  h2otable = (ibn .EQ. 1 .OR. ibn .EQ. 2) .OR. ibn .EQ. 8
295  conbnd = ibn .GE. 2 .AND. ibn .LE. 7
296  co2bnd = ibn .EQ. 3
297  oznbnd = ibn .EQ. 5
298  n2obnd = ibn .EQ. 6 .OR. ibn .EQ. 7
299  ch4bnd = ibn .EQ. 6 .OR. ibn .EQ. 7
300  combnd = ibn .EQ. 4 .OR. ibn .EQ. 5
301  f11bnd = ibn .EQ. 4 .OR. ibn .EQ. 5
302  f12bnd = ibn .EQ. 4 .OR. ibn .EQ. 6
303  f22bnd = ibn .EQ. 4 .OR. ibn .EQ. 6
304  b10bnd = ibn .EQ. 10
305  do_aerosol = na .GT. 0
306  CALL pushreal8array(exptbl, (np+1)*17)
307  exptbl = 0.0
308 !Control packing of the new exponential tables by band
309  SELECT CASE (ibn)
310  CASE (2)
311  CALL pushinteger4(conexp%start)
312  conexp%start = 1
313  CALL pushinteger4(conexp%end)
314  conexp%end = 1
315  CALL pushcontrol4b(7)
316  CASE (3)
317  CALL pushinteger4(h2oexp%start)
318  h2oexp%start = 1
319  CALL pushinteger4(h2oexp%end)
320  h2oexp%end = 6
321  CALL pushinteger4(conexp%start)
322  conexp%start = 7
323  CALL pushinteger4(conexp%end)
324  conexp%end = 9
325  CALL pushcontrol4b(6)
326  CASE (4)
327  CALL pushinteger4(h2oexp%start)
328  h2oexp%start = 1
329  CALL pushinteger4(h2oexp%end)
330  h2oexp%end = 6
331  CALL pushinteger4(conexp%start)
332  conexp%start = 7
333  CALL pushinteger4(conexp%end)
334  conexp%end = 7
335  CALL pushinteger4(comexp%start)
336  comexp%start = 8
337  CALL pushinteger4(comexp%end)
338  comexp%end = 13
339  CALL pushinteger4(f11exp%start)
340  f11exp%start = 14
341  CALL pushinteger4(f11exp%end)
342  f11exp%end = 14
343  CALL pushinteger4(f12exp%start)
344  f12exp%start = 15
345  CALL pushinteger4(f12exp%end)
346  f12exp%end = 15
347  CALL pushinteger4(f22exp%start)
348  f22exp%start = 16
349  CALL pushinteger4(f22exp%end)
350  f22exp%end = 16
351  CALL pushcontrol4b(5)
352  CASE (5)
353  CALL pushinteger4(h2oexp%start)
354  h2oexp%start = 1
355  CALL pushinteger4(h2oexp%end)
356  h2oexp%end = 6
357  CALL pushinteger4(conexp%start)
358  conexp%start = 7
359  CALL pushinteger4(conexp%end)
360  conexp%end = 7
361  CALL pushinteger4(comexp%start)
362  comexp%start = 8
363  CALL pushinteger4(comexp%end)
364  comexp%end = 13
365  CALL pushinteger4(f11exp%start)
366  f11exp%start = 14
367  CALL pushinteger4(f11exp%end)
368  f11exp%end = 14
369  CALL pushcontrol4b(4)
370  CASE (6)
371  CALL pushinteger4(h2oexp%start)
372  h2oexp%start = 1
373  CALL pushinteger4(h2oexp%end)
374  h2oexp%end = 6
375  CALL pushinteger4(conexp%start)
376  conexp%start = 7
377  CALL pushinteger4(conexp%end)
378  conexp%end = 7
379  CALL pushinteger4(n2oexp%start)
380  n2oexp%start = 8
381  CALL pushinteger4(n2oexp%end)
382  n2oexp%end = 11
383  CALL pushinteger4(ch4exp%start)
384  ch4exp%start = 12
385  CALL pushinteger4(ch4exp%end)
386  ch4exp%end = 15
387  CALL pushinteger4(f12exp%start)
388  f12exp%start = 16
389  CALL pushinteger4(f12exp%end)
390  f12exp%end = 16
391  CALL pushinteger4(f22exp%start)
392  f22exp%start = 17
393  CALL pushinteger4(f22exp%end)
394  f22exp%end = 17
395  CALL pushcontrol4b(3)
396  CASE (7)
397  CALL pushinteger4(h2oexp%start)
398  h2oexp%start = 1
399  CALL pushinteger4(h2oexp%end)
400  h2oexp%end = 6
401  CALL pushinteger4(conexp%start)
402  conexp%start = 7
403  CALL pushinteger4(conexp%end)
404  conexp%end = 7
405  CALL pushinteger4(n2oexp%start)
406  n2oexp%start = 8
407  CALL pushinteger4(n2oexp%end)
408  n2oexp%end = 11
409  CALL pushinteger4(ch4exp%start)
410  ch4exp%start = 12
411  CALL pushinteger4(ch4exp%end)
412  ch4exp%end = 15
413  CALL pushcontrol4b(2)
414  CASE (9)
415  CALL pushinteger4(h2oexp%start)
416  h2oexp%start = 1
417  CALL pushinteger4(h2oexp%end)
418  h2oexp%end = 6
419  CALL pushcontrol4b(1)
420  CASE (10)
421  CALL pushinteger4(h2oexp%start)
422  h2oexp%start = 1
423  CALL pushinteger4(h2oexp%end)
424  h2oexp%end = 5
425  CALL pushinteger4(conexp%start)
426  conexp%start = 6
427  CALL pushinteger4(conexp%end)
428  conexp%end = 6
429  CALL pushinteger4(co2exp%start)
430  co2exp%start = 7
431  CALL pushinteger4(co2exp%end)
432  co2exp%end = 12
433  CALL pushinteger4(n2oexp%start)
434  n2oexp%start = 13
435  CALL pushinteger4(n2oexp%end)
436  n2oexp%end = 14
437  CALL pushcontrol4b(0)
438  CASE DEFAULT
439  CALL pushcontrol4b(8)
440  END SELECT
441 !blayer is the spectrally integrated planck flux of the mean layer
442 !temperature derived from eq. (3.11)
443 !The fitting for the planck flux is valid for the range 160-345 K.
444  DO k=1,np
445  CALL planck(ibn, cb, ta_dev(k), blayer(k))
446  END DO
447 !Index "0" is the layer above the top of the atmosphere.
448  blayer(0) = blayer(1)
449  CALL pushreal8(blevel(0))
450  blevel(0) = blayer(1)
451 !Surface emission and reflectivity. See Section 9.
452 !bs and dbs include the effect of surface emissivity.
453  CALL pushreal8(rflxs)
454  CALL pushreal8(dbs)
455  CALL sfcflux(ibn, cb, dcb, ns, fs_dev, tg_dev, eg_dev, tv_dev, &
456 & ev_dev, rv_dev, bs, dbs, rflxs)
457  blayer(np+1) = bs
458 !interpolate Planck function at model levels (linear in p)
459  DO k=2,np
460  CALL pushreal8(blevel(k))
461  blevel(k) = (blayer(k-1)*dp(k)+blayer(k)*dp(k-1))/(dp(k-1)+dp(k)&
462 & )
463  END DO
464 !Extrapolate blevel(1) from blayer(2) and blayer(1)
465  CALL pushreal8(blevel(1))
466  blevel(1) = blayer(1) + (blayer(1)-blayer(2))*dp(1)/(dp(1)+dp(2))
467  CALL pushreal8(blevel(0))
468  blevel(0) = blevel(1)
469 !If the surface air temperature tb is known, compute blevel(np+1)
470  CALL pushreal8(blevel(np+1))
471  CALL planck(ibn, cb, tb_dev, blevel(np+1))
472 !Compute cloud optical thickness following Eqs. (6.4a,b) and (6.7)
473 ! NOTE: dp_pa is only dims(1:np) as the 0'th level isn't needed in getirtau.
474 ! Plus, the pressures in getirtau *MUST* be in Pascals.
475 ! Slots for reff, hydrometeors and tauall are as follows:
476 ! 1 Cloud Ice
477 ! 2 Cloud Liquid
478 ! 3 Falling Liquid (Rain)
479 ! 4 Falling Ice (Snow)
480  CALL pushreal8array(enn, np + 1)
481  CALL pushreal8array(tcldlyr, np + 1)
482  CALL getirtau1(ibn, np, dp_pa, fcld_col, reff_col, cwc_col, &
483 & tcldlyr, enn, aib_ir, awb_ir, aiw_ir, aww_ir, aig_ir, &
484 & awg_ir, cons_grav)
485  DO k=0,np
486  CALL pushinteger4(icx(k))
487  icx(k) = k
488  END DO
489  CALL pushinteger4array(ncld, 3)
490  CALL pushinteger4array(icx, np + 1)
491  CALL mkicx(np, ict, icb, enn, icx, ncld)
492 !Compute optical thickness, single-scattering albedo and asymmetry
493 !factor for a mixture of "na" aerosol types. Eqs. (7.1)-(7.3)
494  IF (do_aerosol) THEN
495  CALL pushreal8(taerlyr(0))
496  taerlyr(0) = 1.0
497  DO k=1,np
498 !-----taerlyr is the aerosol diffuse transmittance
499  CALL pushreal8(taerlyr(k))
500  taerlyr(k) = 1.0
501  IF (taua_dev(k, ibn) .GT. 0.001) THEN
502  IF (ssaa_dev(k, ibn) .GT. 0.001) THEN
503  CALL pushreal8(asya_dev(k, ibn))
504  asya_dev(k, ibn) = asya_dev(k, ibn)/ssaa_dev(k, ibn)
505  CALL pushreal8(ssaa_dev(k, ibn))
506  ssaa_dev(k, ibn) = ssaa_dev(k, ibn)/taua_dev(k, ibn)
507 !Parameterization of aerosol scattering following
508  ff = .5 + (.3739+(0.0076+0.1185*asya_dev(k, ibn))*asya_dev&
509 & (k, ibn))*asya_dev(k, ibn)
510  CALL pushreal8(taua_dev(k, ibn))
511  taua_dev(k, ibn) = taua_dev(k, ibn)*(1.-ssaa_dev(k, ibn)*&
512 & ff)
513  CALL pushcontrol1b(0)
514  ELSE
515  CALL pushcontrol1b(1)
516  END IF
517  taerlyr(k) = exp(-(1.66*taua_dev(k, ibn)))
518  CALL pushcontrol1b(1)
519  ELSE
520  CALL pushcontrol1b(0)
521  END IF
522  END DO
523  CALL pushcontrol1b(0)
524  ELSE
525  CALL pushcontrol1b(1)
526  END IF
527 !Compute the exponential terms (Eq. 8.21) at each layer due to
528 !water vapor line absorption when k-distribution is used
529  IF (.NOT.h2otable .AND. (.NOT.b10bnd)) THEN
530  CALL pushreal8array(exptbl(:, h2oexp%start:h2oexp%end), (np+1)*(&
531 & h2oexp%end-h2oexp%start+1))
532  CALL h2oexps(ibn, np, dh2o, pa, dt, xkw, aw, bw, pm, mw, exptbl(&
533 & :, h2oexp%start:h2oexp%end))
534  CALL pushcontrol1b(0)
535  ELSE
536  CALL pushcontrol1b(1)
537  END IF
538 !compute the exponential terms (Eq. 4.24) at each layer due to
539 !water vapor continuum absorption.
540  CALL pushinteger4(ne)
541  ne = 0
542  IF (conbnd) THEN
543  ne = 1
544  IF (ibn .EQ. 3) ne = 3
545  CALL pushreal8array(exptbl(:, conexp%start:conexp%end), (np+1)*(&
546 & conexp%end-conexp%start+1))
547  CALL conexps(ibn, np, dcont, xke, exptbl(:, conexp%start:conexp%&
548 & end))
549  CALL pushcontrol1b(0)
550  ELSE
551  CALL pushcontrol1b(1)
552  END IF
553  IF (trace) THEN
554 !compute the exponential terms at each layer due to n2o absorption
555  IF (n2obnd) THEN
556  CALL pushreal8array(exptbl(:, n2oexp%start:n2oexp%end), (np+1)&
557 & *(n2oexp%end-n2oexp%start+1))
558  CALL n2oexps(ibn, np, dn2o, pa, dt, exptbl(:, n2oexp%start:&
559 & n2oexp%end))
560  CALL pushcontrol1b(0)
561  ELSE
562  CALL pushcontrol1b(1)
563  END IF
564 !compute the exponential terms at each layer due to ch4 absorption
565  IF (ch4bnd) THEN
566  CALL pushreal8array(exptbl(:, ch4exp%start:ch4exp%end), (np+1)&
567 & *(ch4exp%end-ch4exp%start+1))
568  CALL ch4exps(ibn, np, dch4, pa, dt, exptbl(:, ch4exp%start:&
569 & ch4exp%end))
570  CALL pushcontrol1b(0)
571  ELSE
572  CALL pushcontrol1b(1)
573  END IF
574 !Compute the exponential terms due to co2 minor absorption
575  IF (combnd) THEN
576  CALL pushreal8array(exptbl(:, comexp%start:comexp%end), (np+1)&
577 & *(comexp%end-comexp%start+1))
578  CALL comexps(ibn, np, dco2, dt, exptbl(:, comexp%start:comexp%&
579 & end))
580  CALL pushcontrol1b(0)
581  ELSE
582  CALL pushcontrol1b(1)
583  END IF
584 !Compute the exponential terms due to cfc11 absorption.
585 !The values of the parameters are given in Table 7.
586  IF (f11bnd) THEN
587  a1 = 1.26610e-3
588  b1 = 3.55940e-6
589  fk1 = 1.89736e+1
590  a2 = 8.19370e-4
591  b2 = 4.67810e-6
592  fk2 = 1.01487e+1
593  CALL cfcexps(ibn, np, a1, b1, fk1, a2, b2, fk2, df11, dt, &
594 & exptbl(:, f11exp%start:f11exp%end))
595  CALL pushcontrol1b(0)
596  ELSE
597  CALL pushcontrol1b(1)
598  END IF
599 !Compute the exponential terms due to cfc12 absorption.
600  IF (f12bnd) THEN
601  a1 = 8.77370e-4
602  b1 = -5.88440e-6
603  fk1 = 1.58104e+1
604  a2 = 8.62000e-4
605  b2 = -4.22500e-6
606  fk2 = 3.70107e+1
607  CALL cfcexps(ibn, np, a1, b1, fk1, a2, b2, fk2, df12, dt, &
608 & exptbl(:, f12exp%start:f12exp%end))
609  CALL pushcontrol1b(0)
610  ELSE
611  CALL pushcontrol1b(1)
612  END IF
613 !Compute the exponential terms due to cfc22 absorption.
614  IF (f22bnd) THEN
615  a1 = 9.65130e-4
616  b1 = 1.31280e-5
617  fk1 = 6.18536e+0
618  a2 = -3.00010e-5
619  b2 = 5.25010e-7
620  fk2 = 3.27912e+1
621  CALL cfcexps(ibn, np, a1, b1, fk1, a2, b2, fk2, df22, dt, &
622 & exptbl(:, f22exp%start:f22exp%end))
623  CALL pushcontrol1b(0)
624  ELSE
625  CALL pushcontrol1b(1)
626  END IF
627 !Compute the exponential terms at each layer in band 10 due to
628 !h2o line and continuum, co2, and n2o absorption
629  IF (b10bnd) THEN
630  CALL pushreal8array(n2oexp_tmp, (np+1)*2)
631  CALL pushreal8array(co2exp_tmp, (np+1)*6)
632  CALL pushreal8array(h2oexp_tmp, (np+1)*5)
633  CALL b10exps(np, dh2o, dcont, dco2, dn2o, pa, dt, h2oexp_tmp, &
634 & exptbl(:, conexp%start:conexp%end), co2exp_tmp, &
635 & n2oexp_tmp)
636  exptbl(:, h2oexp%start:h2oexp%end) = h2oexp_tmp
637  exptbl(:, co2exp%start:co2exp%end) = co2exp_tmp
638  exptbl(:, n2oexp%start:n2oexp%end) = n2oexp_tmp
639  CALL pushcontrol2b(0)
640  ELSE
641  CALL pushcontrol2b(1)
642  END IF
643  ELSE
644  CALL pushcontrol2b(2)
645  END IF
646 !blayer(np+1) includes the effect of surface emissivity.
647 ! ALT: this was undefined, check with Max if 0.0 is good value
648  CALL pushreal8(bu(0))
649  bu(0) = 0.0
650  CALL pushreal8(bd(0))
651  bd(0) = blayer(1)
652  CALL pushreal8(bu(np+1))
653  bu(np+1) = blayer(np+1)
654 !do-loop 1500 is for computing upward (bu) and downward (bd)
655 !Here, trant is the transmittance of the layer k2-1.
656  DO k2=1,np+1
657 !for h2o line transmission
658  IF (.NOT.h2otable) THEN
659  th2o = 1.0
660  CALL pushcontrol1b(0)
661  ELSE
662  CALL pushcontrol1b(1)
663  END IF
664 !for h2o continuum transmission
665  CALL pushreal8array(tcon, 3)
666  tcon = 1.0
667  x1 = 0.0
668  x2 = 0.0
669  x3 = 0.0
670  CALL pushreal8(trant)
671  trant = 1.0
672  IF (h2otable) THEN
673 !Compute water vapor transmittance using table look-up.
674  IF (ibn .EQ. 1) THEN
675  CALL pushreal8(trant)
676  CALL pushreal8(x3)
677  CALL pushreal8(x2)
678  CALL pushreal8(x1)
679  CALL tablup(nx1, nh1, dh2o(k2-1), pa(k2-1), dt(k2-1), x1, x2&
680 & , x3, w11, p11, dwe, dpe, h11, h12, h13, trant)
681  CALL pushcontrol1b(0)
682  ELSE
683  CALL pushcontrol1b(1)
684  END IF
685  IF (ibn .EQ. 2) THEN
686  CALL pushreal8(trant)
687  CALL pushreal8(x3)
688  CALL pushreal8(x2)
689  CALL pushreal8(x1)
690  CALL tablup(nx1, nh1, dh2o(k2-1), pa(k2-1), dt(k2-1), x1, x2&
691 & , x3, w11, p11, dwe, dpe, h21, h22, h23, trant)
692  CALL pushcontrol1b(0)
693  ELSE
694  CALL pushcontrol1b(1)
695  END IF
696  IF (ibn .EQ. 8) THEN
697  CALL pushreal8(trant)
698  CALL pushreal8(x3)
699  CALL pushreal8(x2)
700  CALL pushreal8(x1)
701  CALL tablup(nx1, nh1, dh2o(k2-1), pa(k2-1), dt(k2-1), x1, x2&
702 & , x3, w11, p11, dwe, dpe, h81, h82, h83, trant)
703  CALL pushcontrol1b(0)
704  ELSE
705  CALL pushcontrol1b(1)
706  END IF
707 !for water vapor continuum absorption
708  IF (conbnd) THEN
709 ! Only the first exp
710  CALL pushreal8(tcon(1))
711  tcon(1) = tcon(1)*exptbl(k2-1, conexp%start)
712  CALL pushreal8(trant)
713  trant = trant*tcon(1)
714  CALL pushcontrol2b(0)
715  ELSE
716  CALL pushcontrol2b(1)
717  END IF
718  ELSE IF (.NOT.b10bnd) THEN
719 !compute water vapor transmittance using k-distribution
720  arg1 = k2 - 1
721  CALL pushreal8(trant)
722  CALL pushreal8array(tcon, 3)
723  CALL pushreal8array(th2o, 6)
724  CALL h2okdis(ibn, np, arg1, fkw, gkw, ne, exptbl(:, h2oexp%&
725 & start:h2oexp%end), exptbl(:, conexp%start:conexp%end), &
726 & th2o, tcon, trant)
727  CALL pushcontrol2b(2)
728  ELSE
729  CALL pushcontrol2b(3)
730  END IF
731  IF (co2bnd) THEN
732 !Compute co2 transmittance using table look-up method
733  CALL pushreal8(trant)
734  CALL pushreal8(x3)
735  CALL pushreal8(x2)
736  CALL pushreal8(x1)
737  CALL tablup(nx1, nc1, dco2(k2-1), pa(k2-1), dt(k2-1), x1, x2, &
738 & x3, w12, p12, dwe, dpe, c1, c2, c3, trant)
739  CALL pushcontrol1b(0)
740  ELSE
741  CALL pushcontrol1b(1)
742  END IF
743 !Always use table look-up to compute o3 transmittance.
744  IF (oznbnd) THEN
745  CALL pushreal8(trant)
746  CALL pushreal8(x3)
747  CALL pushreal8(x2)
748  CALL pushreal8(x1)
749  CALL tablup(nx1, no1, do3(k2-1), pa(k2-1), dt(k2-1), x1, x2, &
750 & x3, w13, p13, dwe, dpe, oo1, oo2, oo3, trant)
751  CALL pushcontrol1b(0)
752  ELSE
753  CALL pushcontrol1b(1)
754  END IF
755 !include aerosol effect
756  IF (do_aerosol) THEN
757  CALL pushreal8(trant)
758  trant = trant*taerlyr(k2-1)
759  CALL pushcontrol1b(1)
760  ELSE
761  CALL pushcontrol1b(0)
762  END IF
763 !Compute upward (bu) and downward (bd) emission of the layer k2-1
764  CALL pushreal8(xx)
765  xx = (1.-enn(k2-1))*trant
766  IF (0.9999 .GT. xx) THEN
767  CALL pushreal8(yy)
768  yy = xx
769  CALL pushcontrol1b(0)
770  ELSE
771  CALL pushreal8(yy)
772  yy = 0.9999
773  CALL pushcontrol1b(1)
774  END IF
775  IF (0.00001 .LT. yy) THEN
776  CALL pushcontrol1b(0)
777  yy = yy
778  ELSE
779  yy = 0.00001
780  CALL pushcontrol1b(1)
781  END IF
782  xx = (blevel(k2-1)-blevel(k2))/log(yy)
783  CALL pushreal8(bd(k2-1))
784  bd(k2-1) = (blevel(k2)-blevel(k2-1)*yy)/(1.0-yy) - xx
785  CALL pushreal8(bu(k2-1))
786  bu(k2-1) = blevel(k2-1) + blevel(k2) - bd(k2-1)
787  END DO
788 !Initialize fluxes
789  CALL pushreal8array(flxd, np + 2)
790  flxd = 0.0
791 !Compute upward and downward fluxes for each spectral band, ibn.
792  DO k1=0,np
793 !initialization
794  CALL pushreal8(cldlw)
795  cldlw = 0.0
796  CALL pushreal8(cldmd)
797  cldmd = 0.0
798  CALL pushreal8(cldhi)
799  cldhi = 0.0
800  CALL pushreal8(tranal)
801  tranal = 1.0
802 !for h2o line transmission
803  IF (.NOT.h2otable) THEN
804  th2o = 1.0
805  CALL pushcontrol1b(0)
806  ELSE
807  CALL pushcontrol1b(1)
808  END IF
809 !for h2o continuum transmission
810  CALL pushreal8array(tcon, 3)
811  tcon = 1.0
812  IF (trace) THEN
813 !Add trace gases contribution
814  IF (n2obnd) THEN
815 !n2o
816  tn2o = 1.0
817  CALL pushcontrol1b(0)
818  ELSE
819  CALL pushcontrol1b(1)
820  END IF
821  IF (ch4bnd) THEN
822 !ch4
823  tch4 = 1.0
824  CALL pushcontrol1b(0)
825  ELSE
826  CALL pushcontrol1b(1)
827  END IF
828  IF (combnd) THEN
829 !co2-minor
830  tcom = 1.0
831  CALL pushcontrol1b(0)
832  ELSE
833  CALL pushcontrol1b(1)
834  END IF
835  IF (f11bnd) THEN
836 !cfc-11
837  tf11 = 1.0
838  CALL pushcontrol1b(0)
839  ELSE
840  CALL pushcontrol1b(1)
841  END IF
842  IF (f12bnd) THEN
843 !cfc-12
844  tf12 = 1.0
845  CALL pushcontrol1b(0)
846  ELSE
847  CALL pushcontrol1b(1)
848  END IF
849  IF (f22bnd) THEN
850 !cfc-22
851  tf22 = 1.0
852  CALL pushcontrol1b(0)
853  ELSE
854  CALL pushcontrol1b(1)
855  END IF
856  IF (b10bnd) THEN
857 !
858  th2o = 1.0
859  tco2 = 1.0
860  tcon(1) = 1.0
861  tn2o = 1.0
862  CALL pushcontrol2b(0)
863  ELSE
864  CALL pushcontrol2b(1)
865  END IF
866  ELSE
867  CALL pushcontrol2b(2)
868  END IF
869  x1 = 0.0
870  x2 = 0.0
871  x3 = 0.0
872  ad_from = k1 + 1
873  DO k2=ad_from,np+1
874  CALL pushreal8(trant)
875  trant = 1.0
876  IF (h2otable) THEN
877 !Compute water vapor transmittance using table look-up.
878  IF (ibn .EQ. 1) THEN
879  CALL pushreal8(trant)
880  CALL pushreal8(x3)
881  CALL pushreal8(x2)
882  CALL pushreal8(x1)
883  CALL tablup(nx1, nh1, dh2o(k2-1), pa(k2-1), dt(k2-1), x1, &
884 & x2, x3, w11, p11, dwe, dpe, h11, h12, h13, trant)
885  CALL pushcontrol1b(0)
886  ELSE
887  CALL pushcontrol1b(1)
888  END IF
889  IF (ibn .EQ. 2) THEN
890  CALL pushreal8(trant)
891  CALL pushreal8(x3)
892  CALL pushreal8(x2)
893  CALL pushreal8(x1)
894  CALL tablup(nx1, nh1, dh2o(k2-1), pa(k2-1), dt(k2-1), x1, &
895 & x2, x3, w11, p11, dwe, dpe, h21, h22, h23, trant)
896  CALL pushcontrol1b(0)
897  ELSE
898  CALL pushcontrol1b(1)
899  END IF
900  IF (ibn .EQ. 8) THEN
901  CALL pushreal8(trant)
902  CALL pushreal8(x3)
903  CALL pushreal8(x2)
904  CALL pushreal8(x1)
905  CALL tablup(nx1, nh1, dh2o(k2-1), pa(k2-1), dt(k2-1), x1, &
906 & x2, x3, w11, p11, dwe, dpe, h81, h82, h83, trant)
907  CALL pushcontrol1b(0)
908  ELSE
909  CALL pushcontrol1b(1)
910  END IF
911  IF (conbnd) THEN
912 ! Only the first exp
913  CALL pushreal8(tcon(1))
914  tcon(1) = tcon(1)*exptbl(k2-1, conexp%start)
915  CALL pushreal8(trant)
916  trant = trant*tcon(1)
917  CALL pushcontrol2b(0)
918  ELSE
919  CALL pushcontrol2b(1)
920  END IF
921  ELSE IF (.NOT.b10bnd) THEN
922 !Compute water vapor transmittance using k-distribution
923  arg1 = k2 - 1
924  CALL pushreal8(trant)
925  CALL pushreal8array(tcon, 3)
926  CALL pushreal8array(th2o, 6)
927  CALL h2okdis(ibn, np, arg1, fkw, gkw, ne, exptbl(:, h2oexp%&
928 & start:h2oexp%end), exptbl(:, conexp%start:conexp%end)&
929 & , th2o, tcon, trant)
930  CALL pushcontrol2b(2)
931  ELSE
932  CALL pushcontrol2b(3)
933  END IF
934  IF (co2bnd) THEN
935 !Compute co2 transmittance using table look-up method.
936  CALL pushreal8(trant)
937  CALL pushreal8(x3)
938  CALL pushreal8(x2)
939  CALL pushreal8(x1)
940  CALL tablup(nx1, nc1, dco2(k2-1), pa(k2-1), dt(k2-1), x1, x2&
941 & , x3, w12, p12, dwe, dpe, c1, c2, c3, trant)
942  CALL pushcontrol1b(0)
943  ELSE
944  CALL pushcontrol1b(1)
945  END IF
946  IF (oznbnd) THEN
947 !Always use table look-up to compute o3 transmittanc
948  CALL pushreal8(trant)
949  CALL pushreal8(x3)
950  CALL pushreal8(x2)
951  CALL pushreal8(x1)
952  CALL tablup(nx1, no1, do3(k2-1), pa(k2-1), dt(k2-1), x1, x2&
953 & , x3, w13, p13, dwe, dpe, oo1, oo2, oo3, trant)
954  CALL pushcontrol1b(0)
955  ELSE
956  CALL pushcontrol1b(1)
957  END IF
958  IF (trace) THEN
959 !Add trace gas effects
960  IF (n2obnd) THEN
961 !n2o
962  arg1 = k2 - 1
963  CALL pushreal8(trant)
964  CALL pushreal8array(tn2o, 4)
965  CALL n2okdis(ibn, np, arg1, exptbl(:, n2oexp%start:n2oexp%&
966 & end), tn2o, trant)
967  CALL pushcontrol1b(0)
968  ELSE
969  CALL pushcontrol1b(1)
970  END IF
971  IF (ch4bnd) THEN
972 !ch4
973  arg1 = k2 - 1
974  CALL pushreal8(trant)
975  CALL pushreal8array(tch4, 4)
976  CALL ch4kdis(ibn, np, arg1, exptbl(:, ch4exp%start:ch4exp%&
977 & end), tch4, trant)
978  CALL pushcontrol1b(0)
979  ELSE
980  CALL pushcontrol1b(1)
981  END IF
982  IF (combnd) THEN
983 !co2-minor
984  arg1 = k2 - 1
985  CALL pushreal8(trant)
986  CALL pushreal8array(tcom, 6)
987  CALL comkdis(ibn, np, arg1, exptbl(:, comexp%start:comexp%&
988 & end), tcom, trant)
989  CALL pushcontrol1b(0)
990  ELSE
991  CALL pushcontrol1b(1)
992  END IF
993  IF (f11bnd) THEN
994 !cfc11
995  arg1 = k2 - 1
996  CALL pushreal8(trant)
997  CALL pushreal8(tf11)
998  CALL cfckdis(np, arg1, exptbl(:, f11exp%start:f11exp%end)&
999 & , tf11, trant)
1000  CALL pushcontrol1b(0)
1001  ELSE
1002  CALL pushcontrol1b(1)
1003  END IF
1004  IF (f12bnd) THEN
1005 !cfc12
1006  arg1 = k2 - 1
1007  CALL pushreal8(trant)
1008  CALL pushreal8(tf12)
1009  CALL cfckdis(np, arg1, exptbl(:, f12exp%start:f12exp%end)&
1010 & , tf12, trant)
1011  CALL pushcontrol1b(0)
1012  ELSE
1013  CALL pushcontrol1b(1)
1014  END IF
1015  IF (f22bnd) THEN
1016 !CFC22
1017  arg1 = k2 - 1
1018  CALL pushreal8(trant)
1019  CALL pushreal8(tf22)
1020  CALL cfckdis(np, arg1, exptbl(:, f22exp%start:f22exp%end)&
1021 & , tf22, trant)
1022  CALL pushcontrol1b(0)
1023  ELSE
1024  CALL pushcontrol1b(1)
1025  END IF
1026  IF (b10bnd) THEN
1027  arg1 = k2 - 1
1028  CALL pushreal8array(tn2o, 4)
1029  CALL pushreal8array(tco2, 6)
1030  CALL pushreal8array(tcon, 3)
1031  CALL pushreal8array(th2o, 6)
1032  CALL b10kdis(np, arg1, exptbl(:, h2oexp%start:h2oexp%end)&
1033 & , exptbl(:, conexp%start:conexp%end), exptbl(:, &
1034 & co2exp%start:co2exp%end), exptbl(:, n2oexp%start:&
1035 & n2oexp%end), th2o, tcon, tco2, tn2o, trant)
1036  CALL pushcontrol2b(0)
1037  ELSE
1038  CALL pushcontrol2b(1)
1039  END IF
1040  ELSE
1041  CALL pushcontrol2b(2)
1042  END IF
1043  IF (do_aerosol) THEN
1044  CALL pushreal8(tranal)
1045  tranal = tranal*taerlyr(k2-1)
1046  CALL pushreal8(trant)
1047  trant = trant*tranal
1048  CALL pushcontrol1b(0)
1049  ELSE
1050  CALL pushcontrol1b(1)
1051  END IF
1052  IF (enn(k2-1) .GE. 0.001) THEN
1053  CALL pushreal8(cldlw)
1054  CALL pushreal8(cldmd)
1055  CALL pushreal8(cldhi)
1056  CALL cldovlp(np, k1, k2, ict, icb, icx, ncld, enn, tcldlyr, &
1057 & cldhi, cldmd, cldlw)
1058  CALL pushcontrol1b(0)
1059  ELSE
1060  CALL pushcontrol1b(1)
1061  END IF
1062  CALL pushreal8(fclr)
1063  fclr = (1.0-cldhi)*(1.0-cldmd)*(1.0-cldlw)
1064  IF (k2 .EQ. k1 + 1 .AND. ibn .NE. 10) THEN
1065  flxd(k2) = flxd(k2) + bd(k1)
1066  CALL pushcontrol1b(0)
1067  ELSE
1068  CALL pushcontrol1b(1)
1069  END IF
1070  CALL pushreal8(xx)
1071  IF (k1 .EQ. 0) THEN
1072 !mjs bd(-1) is not defined
1073  xx = -(trant*bd(k1))
1074  CALL pushcontrol1b(0)
1075  ELSE
1076  xx = trant*(bd(k1-1)-bd(k1))
1077  CALL pushcontrol1b(1)
1078  END IF
1079  flxd(k2) = flxd(k2) + xx*fclr
1080  END DO
1081  CALL pushinteger4(ad_from)
1082  CALL pushreal8(transfc(k1))
1083  transfc(k1) = trant*fclr
1084  IF (k1 .GT. 0) THEN
1085  CALL pushcontrol1b(1)
1086  ELSE
1087  CALL pushcontrol1b(0)
1088  END IF
1089  END DO
1090  IF (.NOT.b10bnd) THEN
1091  CALL pushcontrol1b(1)
1092  ELSE
1093  CALL pushcontrol1b(0)
1094  END IF
1095  CALL pushinteger4(ibn)
1096  ad_count = ad_count + 1
1097  END IF
1098  END DO
1099  CALL pushcontrol1b(0)
1100  CALL pushinteger4(ad_count)
1101  GOTO 110
1102  100 CALL pushcontrol1b(1)
1103  CALL pushinteger4(ad_count)
1104  110 CALL popinteger4(ad_count)
1105  DO i=1,ad_count
1106  IF (i .EQ. 1) THEN
1107  CALL popcontrol1b(branch)
1108  IF (branch .EQ. 0) THEN
1109  dh2ob = 0.0_8
1110  dcontb = 0.0_8
1111  n2oexp_tmpb = 0.0_8
1112  dtb = 0.0_8
1113  reff_colb = 0.0_8
1114  blayerb = 0.0_8
1115  tn2ob = 0.0_8
1116  transfcb = 0.0_8
1117  bdb = 0.0_8
1118  h2oexp_tmpb = 0.0_8
1119  tcomb = 0.0_8
1120  tf11b = 0.0_8
1121  tf12b = 0.0_8
1122  trantb = 0.0_8
1123  bub = 0.0_8
1124  th2ob = 0.0_8
1125  tcldlyrb = 0.0_8
1126  tch4b = 0.0_8
1127  tf22b = 0.0_8
1128  ennb = 0.0_8
1129  fcld_colb = 0.0_8
1130  fclrb = 0.0_8
1131  cwc_colb = 0.0_8
1132  tco2b = 0.0_8
1133  co2exp_tmpb = 0.0_8
1134  taerlyrb = 0.0_8
1135  do3b = 0.0_8
1136  blevelb = 0.0_8
1137  ELSE
1138  dh2ob = 0.0_8
1139  dcontb = 0.0_8
1140  n2oexp_tmpb = 0.0_8
1141  dtb = 0.0_8
1142  reff_colb = 0.0_8
1143  blayerb = 0.0_8
1144  tn2ob = 0.0_8
1145  transfcb = 0.0_8
1146  bdb = 0.0_8
1147  h2oexp_tmpb = 0.0_8
1148  tcomb = 0.0_8
1149  tf11b = 0.0_8
1150  tf12b = 0.0_8
1151  trantb = 0.0_8
1152  bub = 0.0_8
1153  th2ob = 0.0_8
1154  tcldlyrb = 0.0_8
1155  tch4b = 0.0_8
1156  tf22b = 0.0_8
1157  ennb = 0.0_8
1158  fcld_colb = 0.0_8
1159  fclrb = 0.0_8
1160  cwc_colb = 0.0_8
1161  tco2b = 0.0_8
1162  co2exp_tmpb = 0.0_8
1163  taerlyrb = 0.0_8
1164  do3b = 0.0_8
1165  blevelb = 0.0_8
1166  END IF
1167  ELSE
1168  flxub = 0.0_8
1169  flxdb = 0.0_8
1170  DO k=np+1,1,-1
1171  flxdb(k) = flxdb(k) + flxd_devb(k)
1172  flxub(k) = flxub(k) + flxu_devb(k)
1173  END DO
1174  CALL popcontrol1b(branch)
1175  IF (branch .NE. 0) THEN
1176  DO k=np+1,1,-1
1177  flxdb(np+1) = flxdb(np+1) - rflxs*transfc(k)*flxub(k)
1178  transfcb(k) = transfcb(k) - rflxs*flxd(np+1)*flxub(k)
1179  END DO
1180  blayerb(np+1) = blayerb(np+1) - flxub(np+1)
1181  flxub(np+1) = 0.0_8
1182  END IF
1183  exptblb = 0.0_8
1184  DO k1=np,0,-1
1185  CALL popcontrol1b(branch)
1186  IF (branch .NE. 0) transfcb(k1) = transfcb(k1) - dbs*dfdts_devb(&
1187 & k1)
1188  CALL popreal8(transfc(k1))
1189  trantb = trantb + fclr*transfcb(k1)
1190  fclrb = fclrb + trant*transfcb(k1)
1191  transfcb(k1) = 0.0_8
1192  cldhib = 0.0_8
1193  tconb = 0.0_8
1194  tranalb = 0.0_8
1195  x1b = 0.0_8
1196  x2b = 0.0_8
1197  x3b = 0.0_8
1198  cldlwb = 0.0_8
1199  cldmdb = 0.0_8
1200  CALL popinteger4(ad_from)
1201  DO k2=np+1,ad_from,-1
1202  xxb = fclr*flxdb(k2)
1203  fclrb = fclrb + xx*flxdb(k2)
1204  CALL popcontrol1b(branch)
1205  IF (branch .EQ. 0) THEN
1206  trantb = trantb - bd(k1)*xxb
1207  bdb(k1) = bdb(k1) - trant*xxb
1208  ELSE
1209  trantb = trantb + (bd(k1-1)-bd(k1))*xxb
1210  bdb(k1-1) = bdb(k1-1) + trant*xxb
1211  bdb(k1) = bdb(k1) - trant*xxb
1212  END IF
1213  xx = trant*(bu(k2-1)-bu(k2))
1214  xxb = fclr*flxub(k1)
1215  fclrb = fclrb + xx*flxub(k1)
1216  CALL popreal8(xx)
1217  trantb = trantb + (bu(k2-1)-bu(k2))*xxb
1218  bub(k2-1) = bub(k2-1) + trant*xxb
1219  bub(k2) = bub(k2) - trant*xxb
1220  CALL popcontrol1b(branch)
1221  IF (branch .EQ. 0) THEN
1222  bdb(k1) = bdb(k1) + flxdb(k2)
1223  bub(k1) = bub(k1) - flxub(k1)
1224  END IF
1225  CALL popreal8(fclr)
1226  cldhib = cldhib - (1.0-cldmd)*(1.0-cldlw)*fclrb
1227  cldmdb = cldmdb - (1.0-cldhi)*(1.0-cldlw)*fclrb
1228  cldlwb = cldlwb - (1.0-cldhi)*(1.0-cldmd)*fclrb
1229  CALL popcontrol1b(branch)
1230  IF (branch .EQ. 0) THEN
1231  CALL popreal8(cldhi)
1232  CALL popreal8(cldmd)
1233  CALL popreal8(cldlw)
1234  CALL cldovlp_b(np, k1, k2, ict, icb, icx, ncld, enn, ennb, &
1235 & tcldlyr, tcldlyrb, cldhi, cldhib, cldmd, cldmdb, &
1236 & cldlw, cldlwb)
1237  END IF
1238  CALL popcontrol1b(branch)
1239  IF (branch .EQ. 0) THEN
1240  CALL popreal8(trant)
1241  tranalb = tranalb + trant*trantb
1242  trantb = tranal*trantb
1243  CALL popreal8(tranal)
1244  taerlyrb(k2-1) = taerlyrb(k2-1) + tranal*tranalb
1245  tranalb = taerlyr(k2-1)*tranalb
1246  END IF
1247  CALL popcontrol2b(branch)
1248  IF (branch .EQ. 0) THEN
1249  arg1 = k2 - 1
1250  CALL popreal8array(th2o, 6)
1251  CALL popreal8array(tcon, 3)
1252  CALL popreal8array(tco2, 6)
1253  CALL popreal8array(tn2o, 4)
1254  CALL b10kdis_b(np, arg1, exptbl(:, h2oexp%start:h2oexp%end)&
1255 & , exptblb(:, h2oexp%start:h2oexp%end), exptbl(:, &
1256 & conexp%start:conexp%end), exptblb(:, conexp%start:&
1257 & conexp%end), exptbl(:, co2exp%start:co2exp%end), &
1258 & exptblb(:, co2exp%start:co2exp%end), exptbl(:, &
1259 & n2oexp%start:n2oexp%end), exptblb(:, n2oexp%start:&
1260 & n2oexp%end), th2o, th2ob, tcon, tconb, tco2, tco2b&
1261 & , tn2o, tn2ob, trant, trantb)
1262  trantb = 0.0_8
1263  ELSE IF (branch .NE. 1) THEN
1264  GOTO 120
1265  END IF
1266  CALL popcontrol1b(branch)
1267  IF (branch .EQ. 0) THEN
1268  arg1 = k2 - 1
1269  CALL popreal8(tf22)
1270  CALL popreal8(trant)
1271  CALL cfckdis_b(np, arg1, exptbl(:, f22exp%start:f22exp%end)&
1272 & , exptblb(:, f22exp%start:f22exp%end), tf22, tf22b&
1273 & , trant, trantb)
1274  END IF
1275  CALL popcontrol1b(branch)
1276  IF (branch .EQ. 0) THEN
1277  arg1 = k2 - 1
1278  CALL popreal8(tf12)
1279  CALL popreal8(trant)
1280  CALL cfckdis_b(np, arg1, exptbl(:, f12exp%start:f12exp%end)&
1281 & , exptblb(:, f12exp%start:f12exp%end), tf12, tf12b&
1282 & , trant, trantb)
1283  END IF
1284  CALL popcontrol1b(branch)
1285  IF (branch .EQ. 0) THEN
1286  arg1 = k2 - 1
1287  CALL popreal8(tf11)
1288  CALL popreal8(trant)
1289  CALL cfckdis_b(np, arg1, exptbl(:, f11exp%start:f11exp%end)&
1290 & , exptblb(:, f11exp%start:f11exp%end), tf11, tf11b&
1291 & , trant, trantb)
1292  END IF
1293  CALL popcontrol1b(branch)
1294  IF (branch .EQ. 0) THEN
1295  arg1 = k2 - 1
1296  CALL popreal8array(tcom, 6)
1297  CALL popreal8(trant)
1298  CALL comkdis_b(ibn, np, arg1, exptbl(:, comexp%start:comexp%&
1299 & end), exptblb(:, comexp%start:comexp%end), tcom, &
1300 & tcomb, trant, trantb)
1301  END IF
1302  CALL popcontrol1b(branch)
1303  IF (branch .EQ. 0) THEN
1304  arg1 = k2 - 1
1305  CALL popreal8array(tch4, 4)
1306  CALL popreal8(trant)
1307  CALL ch4kdis_b(ibn, np, arg1, exptbl(:, ch4exp%start:ch4exp%&
1308 & end), exptblb(:, ch4exp%start:ch4exp%end), tch4, &
1309 & tch4b, trant, trantb)
1310  END IF
1311  CALL popcontrol1b(branch)
1312  IF (branch .EQ. 0) THEN
1313  arg1 = k2 - 1
1314  CALL popreal8array(tn2o, 4)
1315  CALL popreal8(trant)
1316  CALL n2okdis_b(ibn, np, arg1, exptbl(:, n2oexp%start:n2oexp%&
1317 & end), exptblb(:, n2oexp%start:n2oexp%end), tn2o, &
1318 & tn2ob, trant, trantb)
1319  END IF
1320  120 CALL popcontrol1b(branch)
1321  IF (branch .EQ. 0) THEN
1322  CALL popreal8(x1)
1323  CALL popreal8(x2)
1324  CALL popreal8(x3)
1325  CALL popreal8(trant)
1326  CALL tablup_b(nx1, no1, do3(k2-1), do3b(k2-1), pa(k2-1), dt(&
1327 & k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w13, &
1328 & p13, dwe, dpe, oo1, oo2, oo3, trant, trantb)
1329  END IF
1330  CALL popcontrol1b(branch)
1331  IF (branch .EQ. 0) THEN
1332  CALL popreal8(x1)
1333  CALL popreal8(x2)
1334  CALL popreal8(x3)
1335  CALL popreal8(trant)
1336  dco2b = 0.0_8
1337  CALL tablup_b(nx1, nc1, dco2(k2-1), dco2b(k2-1), pa(k2-1), &
1338 & dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w12&
1339 & , p12, dwe, dpe, c1, c2, c3, trant, trantb)
1340  END IF
1341  CALL popcontrol2b(branch)
1342  IF (branch .LT. 2) THEN
1343  IF (branch .EQ. 0) THEN
1344  CALL popreal8(trant)
1345  tconb(1) = tconb(1) + trant*trantb
1346  trantb = tcon(1)*trantb
1347  CALL popreal8(tcon(1))
1348  exptblb(k2-1, conexp%start) = exptblb(k2-1, conexp%start) &
1349 & + tcon(1)*tconb(1)
1350  tconb(1) = exptbl(k2-1, conexp%start)*tconb(1)
1351  END IF
1352  CALL popcontrol1b(branch)
1353  IF (branch .EQ. 0) THEN
1354  CALL popreal8(x1)
1355  CALL popreal8(x2)
1356  CALL popreal8(x3)
1357  CALL popreal8(trant)
1358  CALL tablup_b(nx1, nh1, dh2o(k2-1), dh2ob(k2-1), pa(k2-1)&
1359 & , dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, &
1360 & w11, p11, dwe, dpe, h81, h82, h83, trant, trantb)
1361  END IF
1362  CALL popcontrol1b(branch)
1363  IF (branch .EQ. 0) THEN
1364  CALL popreal8(x1)
1365  CALL popreal8(x2)
1366  CALL popreal8(x3)
1367  CALL popreal8(trant)
1368  CALL tablup_b(nx1, nh1, dh2o(k2-1), dh2ob(k2-1), pa(k2-1)&
1369 & , dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, &
1370 & w11, p11, dwe, dpe, h21, h22, h23, trant, trantb)
1371  END IF
1372  CALL popcontrol1b(branch)
1373  IF (branch .EQ. 0) THEN
1374  CALL popreal8(x1)
1375  CALL popreal8(x2)
1376  CALL popreal8(x3)
1377  CALL popreal8(trant)
1378  CALL tablup_b(nx1, nh1, dh2o(k2-1), dh2ob(k2-1), pa(k2-1)&
1379 & , dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, &
1380 & w11, p11, dwe, dpe, h11, h12, h13, trant, trantb)
1381  END IF
1382  ELSE IF (branch .EQ. 2) THEN
1383  arg1 = k2 - 1
1384  CALL popreal8array(th2o, 6)
1385  CALL popreal8array(tcon, 3)
1386  CALL popreal8(trant)
1387  CALL h2okdis_b(ibn, np, arg1, fkw, gkw, ne, exptbl(:, h2oexp&
1388 & %start:h2oexp%end), exptblb(:, h2oexp%start:h2oexp%&
1389 & end), exptbl(:, conexp%start:conexp%end), exptblb(:&
1390 & , conexp%start:conexp%end), th2o, th2ob, tcon, &
1391 & tconb, trant, trantb)
1392  END IF
1393  CALL popreal8(trant)
1394  trantb = 0.0_8
1395  fclrb = 0.0_8
1396  END DO
1397  CALL popcontrol2b(branch)
1398  IF (branch .EQ. 0) THEN
1399  tn2ob = 0.0_8
1400  th2ob = 0.0_8
1401  tco2b = 0.0_8
1402  ELSE IF (branch .NE. 1) THEN
1403  GOTO 130
1404  END IF
1405  CALL popcontrol1b(branch)
1406  IF (branch .EQ. 0) tf22b = 0.0_8
1407  CALL popcontrol1b(branch)
1408  IF (branch .EQ. 0) tf12b = 0.0_8
1409  CALL popcontrol1b(branch)
1410  IF (branch .EQ. 0) tf11b = 0.0_8
1411  CALL popcontrol1b(branch)
1412  IF (branch .EQ. 0) tcomb = 0.0_8
1413  CALL popcontrol1b(branch)
1414  IF (branch .EQ. 0) tch4b = 0.0_8
1415  CALL popcontrol1b(branch)
1416  IF (branch .EQ. 0) tn2ob = 0.0_8
1417  130 CALL popreal8array(tcon, 3)
1418  CALL popcontrol1b(branch)
1419  IF (branch .EQ. 0) th2ob = 0.0_8
1420  CALL popreal8(tranal)
1421  CALL popreal8(cldhi)
1422  CALL popreal8(cldmd)
1423  CALL popreal8(cldlw)
1424  END DO
1425  CALL popreal8array(flxd, np + 2)
1426  DO k2=np+1,1,-1
1427  bdb(k2-1) = bdb(k2-1) - bub(k2-1)
1428  tempb5 = bdb(k2-1)/(1.0-yy)
1429  xxb = -bdb(k2-1)
1430  temp2 = log(yy)
1431  tempb6 = xxb/temp2
1432  CALL popreal8(bu(k2-1))
1433  blevelb(k2-1) = blevelb(k2-1) + bub(k2-1)
1434  blevelb(k2) = blevelb(k2) + tempb5 + bub(k2-1)
1435  bub(k2-1) = 0.0_8
1436  CALL popreal8(bd(k2-1))
1437  blevelb(k2-1) = blevelb(k2-1) + tempb6 - yy*tempb5
1438  yyb = ((blevel(k2)-blevel(k2-1)*yy)/(1.0-yy)-blevel(k2-1))*&
1439 & tempb5 - (blevel(k2-1)-blevel(k2))*tempb6/(temp2*yy)
1440  bdb(k2-1) = 0.0_8
1441  blevelb(k2) = blevelb(k2) - tempb6
1442  CALL popcontrol1b(branch)
1443  IF (branch .NE. 0) yyb = 0.0_8
1444  CALL popcontrol1b(branch)
1445  IF (branch .EQ. 0) THEN
1446  CALL popreal8(yy)
1447  xxb = yyb
1448  ELSE
1449  CALL popreal8(yy)
1450  xxb = 0.0_8
1451  END IF
1452  CALL popreal8(xx)
1453  ennb(k2-1) = ennb(k2-1) - trant*xxb
1454  trantb = trantb + (1.-enn(k2-1))*xxb
1455  CALL popcontrol1b(branch)
1456  IF (branch .NE. 0) THEN
1457  CALL popreal8(trant)
1458  taerlyrb(k2-1) = taerlyrb(k2-1) + trant*trantb
1459  trantb = taerlyr(k2-1)*trantb
1460  END IF
1461  CALL popcontrol1b(branch)
1462  IF (branch .EQ. 0) THEN
1463  CALL popreal8(x1)
1464  CALL popreal8(x2)
1465  CALL popreal8(x3)
1466  CALL popreal8(trant)
1467  x1b = 0.0_8
1468  x2b = 0.0_8
1469  x3b = 0.0_8
1470  CALL tablup_b(nx1, no1, do3(k2-1), do3b(k2-1), pa(k2-1), dt(k2&
1471 & -1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w13, p13, &
1472 & dwe, dpe, oo1, oo2, oo3, trant, trantb)
1473  ELSE
1474  x1b = 0.0_8
1475  x2b = 0.0_8
1476  x3b = 0.0_8
1477  END IF
1478  CALL popcontrol1b(branch)
1479  IF (branch .EQ. 0) THEN
1480  CALL popreal8(x1)
1481  CALL popreal8(x2)
1482  CALL popreal8(x3)
1483  CALL popreal8(trant)
1484  dco2b = 0.0_8
1485  CALL tablup_b(nx1, nc1, dco2(k2-1), dco2b(k2-1), pa(k2-1), dt(&
1486 & k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w12, p12&
1487 & , dwe, dpe, c1, c2, c3, trant, trantb)
1488  END IF
1489  CALL popcontrol2b(branch)
1490  IF (branch .LT. 2) THEN
1491  IF (branch .EQ. 0) THEN
1492  tconb = 0.0_8
1493  CALL popreal8(trant)
1494  tconb(1) = tconb(1) + trant*trantb
1495  trantb = tcon(1)*trantb
1496  CALL popreal8(tcon(1))
1497  exptblb(k2-1, conexp%start) = exptblb(k2-1, conexp%start) + &
1498 & tcon(1)*tconb(1)
1499  END IF
1500  CALL popcontrol1b(branch)
1501  IF (branch .EQ. 0) THEN
1502  CALL popreal8(x1)
1503  CALL popreal8(x2)
1504  CALL popreal8(x3)
1505  CALL popreal8(trant)
1506  CALL tablup_b(nx1, nh1, dh2o(k2-1), dh2ob(k2-1), pa(k2-1), &
1507 & dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w11&
1508 & , p11, dwe, dpe, h81, h82, h83, trant, trantb)
1509  END IF
1510  CALL popcontrol1b(branch)
1511  IF (branch .EQ. 0) THEN
1512  CALL popreal8(x1)
1513  CALL popreal8(x2)
1514  CALL popreal8(x3)
1515  CALL popreal8(trant)
1516  CALL tablup_b(nx1, nh1, dh2o(k2-1), dh2ob(k2-1), pa(k2-1), &
1517 & dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w11&
1518 & , p11, dwe, dpe, h21, h22, h23, trant, trantb)
1519  END IF
1520  CALL popcontrol1b(branch)
1521  IF (branch .EQ. 0) THEN
1522  CALL popreal8(x1)
1523  CALL popreal8(x2)
1524  CALL popreal8(x3)
1525  CALL popreal8(trant)
1526  CALL tablup_b(nx1, nh1, dh2o(k2-1), dh2ob(k2-1), pa(k2-1), &
1527 & dt(k2-1), dtb(k2-1), x1, x1b, x2, x2b, x3, x3b, w11&
1528 & , p11, dwe, dpe, h11, h12, h13, trant, trantb)
1529  END IF
1530  ELSE IF (branch .EQ. 2) THEN
1531  arg1 = k2 - 1
1532  CALL popreal8array(th2o, 6)
1533  CALL popreal8array(tcon, 3)
1534  CALL popreal8(trant)
1535  tconb = 0.0_8
1536  CALL h2okdis_b(ibn, np, arg1, fkw, gkw, ne, exptbl(:, h2oexp%&
1537 & start:h2oexp%end), exptblb(:, h2oexp%start:h2oexp%end&
1538 & ), exptbl(:, conexp%start:conexp%end), exptblb(:, &
1539 & conexp%start:conexp%end), th2o, th2ob, tcon, tconb, &
1540 & trant, trantb)
1541  END IF
1542  CALL popreal8(trant)
1543  CALL popreal8array(tcon, 3)
1544  CALL popcontrol1b(branch)
1545  IF (branch .EQ. 0) th2ob = 0.0_8
1546  trantb = 0.0_8
1547  END DO
1548  CALL popreal8(bu(np+1))
1549  blayerb(np+1) = blayerb(np+1) + bub(np+1)
1550  bub(np+1) = 0.0_8
1551  CALL popreal8(bd(0))
1552  blayerb(1) = blayerb(1) + bdb(0)
1553  bdb(0) = 0.0_8
1554  CALL popreal8(bu(0))
1555  bub(0) = 0.0_8
1556  CALL popcontrol2b(branch)
1557  IF (branch .EQ. 0) THEN
1558  n2oexp_tmpb = n2oexp_tmpb + exptblb(:, n2oexp%start:n2oexp%end)
1559  exptblb(:, n2oexp%start:n2oexp%end) = 0.0_8
1560  co2exp_tmpb = co2exp_tmpb + exptblb(:, co2exp%start:co2exp%end)
1561  exptblb(:, co2exp%start:co2exp%end) = 0.0_8
1562  h2oexp_tmpb = h2oexp_tmpb + exptblb(:, h2oexp%start:h2oexp%end)
1563  exptblb(:, h2oexp%start:h2oexp%end) = 0.0_8
1564  CALL popreal8array(h2oexp_tmp, (np+1)*5)
1565  CALL popreal8array(co2exp_tmp, (np+1)*6)
1566  CALL popreal8array(n2oexp_tmp, (np+1)*2)
1567  CALL b10exps_b(np, dh2o, dh2ob, dcont, dcontb, dco2, dn2o, pa, &
1568 & dt, dtb, h2oexp_tmp, h2oexp_tmpb, exptbl(:, conexp%&
1569 & start:conexp%end), exptblb(:, conexp%start:conexp%end)&
1570 & , co2exp_tmp, co2exp_tmpb, n2oexp_tmp, n2oexp_tmpb)
1571  ELSE IF (branch .NE. 1) THEN
1572  GOTO 140
1573  END IF
1574  CALL popcontrol1b(branch)
1575  IF (branch .EQ. 0) THEN
1576  a1 = 9.65130e-4
1577  a2 = -3.00010e-5
1578  b1 = 1.31280e-5
1579  b2 = 5.25010e-7
1580  fk1 = 6.18536e+0
1581  fk2 = 3.27912e+1
1582  CALL cfcexps_b(ibn, np, a1, b1, fk1, a2, b2, fk2, df22, dt, dtb&
1583 & , exptbl(:, f22exp%start:f22exp%end), exptblb(:, f22exp&
1584 & %start:f22exp%end))
1585  END IF
1586  CALL popcontrol1b(branch)
1587  IF (branch .EQ. 0) THEN
1588  a1 = 8.77370e-4
1589  a2 = 8.62000e-4
1590  b1 = -5.88440e-6
1591  b2 = -4.22500e-6
1592  fk1 = 1.58104e+1
1593  fk2 = 3.70107e+1
1594  CALL cfcexps_b(ibn, np, a1, b1, fk1, a2, b2, fk2, df12, dt, dtb&
1595 & , exptbl(:, f12exp%start:f12exp%end), exptblb(:, f12exp&
1596 & %start:f12exp%end))
1597  END IF
1598  CALL popcontrol1b(branch)
1599  IF (branch .EQ. 0) THEN
1600  a1 = 1.26610e-3
1601  a2 = 8.19370e-4
1602  b1 = 3.55940e-6
1603  b2 = 4.67810e-6
1604  fk1 = 1.89736e+1
1605  fk2 = 1.01487e+1
1606  CALL cfcexps_b(ibn, np, a1, b1, fk1, a2, b2, fk2, df11, dt, dtb&
1607 & , exptbl(:, f11exp%start:f11exp%end), exptblb(:, f11exp&
1608 & %start:f11exp%end))
1609  END IF
1610  CALL popcontrol1b(branch)
1611  IF (branch .EQ. 0) THEN
1612  CALL popreal8array(exptbl(:, comexp%start:comexp%end), (np+1)*(&
1613 & comexp%end-comexp%start+1))
1614  CALL comexps_b(ibn, np, dco2, dt, dtb, exptbl(:, comexp%start:&
1615 & comexp%end), exptblb(:, comexp%start:comexp%end))
1616  END IF
1617  CALL popcontrol1b(branch)
1618  IF (branch .EQ. 0) THEN
1619  CALL popreal8array(exptbl(:, ch4exp%start:ch4exp%end), (np+1)*(&
1620 & ch4exp%end-ch4exp%start+1))
1621  CALL ch4exps_b(ibn, np, dch4, pa, dt, dtb, exptbl(:, ch4exp%&
1622 & start:ch4exp%end), exptblb(:, ch4exp%start:ch4exp%end))
1623  END IF
1624  CALL popcontrol1b(branch)
1625  IF (branch .EQ. 0) THEN
1626  CALL popreal8array(exptbl(:, n2oexp%start:n2oexp%end), (np+1)*(&
1627 & n2oexp%end-n2oexp%start+1))
1628  CALL n2oexps_b(ibn, np, dn2o, pa, dt, dtb, exptbl(:, n2oexp%&
1629 & start:n2oexp%end), exptblb(:, n2oexp%start:n2oexp%end))
1630  END IF
1631  140 CALL popcontrol1b(branch)
1632  IF (branch .EQ. 0) THEN
1633  CALL popreal8array(exptbl(:, conexp%start:conexp%end), (np+1)*(&
1634 & conexp%end-conexp%start+1))
1635  CALL conexps_b(ibn, np, dcont, dcontb, xke, exptbl(:, conexp%&
1636 & start:conexp%end), exptblb(:, conexp%start:conexp%end))
1637  END IF
1638  CALL popinteger4(ne)
1639  CALL popcontrol1b(branch)
1640  IF (branch .EQ. 0) THEN
1641  CALL popreal8array(exptbl(:, h2oexp%start:h2oexp%end), (np+1)*(&
1642 & h2oexp%end-h2oexp%start+1))
1643  CALL h2oexps_b(ibn, np, dh2o, dh2ob, pa, dt, dtb, xkw, aw, bw, &
1644 & pm, mw, exptbl(:, h2oexp%start:h2oexp%end), exptblb(:, &
1645 & h2oexp%start:h2oexp%end))
1646  exptblb(:, h2oexp%start:h2oexp%end) = 0.0_8
1647  END IF
1648  CALL popcontrol1b(branch)
1649  IF (branch .EQ. 0) THEN
1650  DO k=np,1,-1
1651  CALL popcontrol1b(branch)
1652  IF (branch .NE. 0) THEN
1653  taua_devb(k, ibn) = taua_devb(k, ibn) - exp(-(1.66*taua_dev(&
1654 & k, ibn)))*1.66*taerlyrb(k)
1655  taerlyrb(k) = 0.0_8
1656  CALL popcontrol1b(branch)
1657  IF (branch .EQ. 0) THEN
1658  ff = .5 + (.3739+(0.0076+0.1185*asya_dev(k, ibn))*asya_dev&
1659 & (k, ibn))*asya_dev(k, ibn)
1660  CALL popreal8(taua_dev(k, ibn))
1661  tempb1 = taua_dev(k, ibn)*taua_devb(k, ibn)
1662  ssaa_devb(k, ibn) = ssaa_devb(k, ibn) - ff*tempb1
1663  ffb = -(ssaa_dev(k, ibn)*tempb1)
1664  taua_devb(k, ibn) = (1.-ssaa_dev(k, ibn)*ff)*taua_devb(k, &
1665 & ibn)
1666  tempb2 = asya_dev(k, ibn)*ffb
1667  temp1 = 0.1185*asya_dev(k, ibn) + 0.0076
1668  asya_devb(k, ibn) = asya_devb(k, ibn) + (temp1*asya_dev(k&
1669 & , ibn)+.3739)*ffb + (temp1+asya_dev(k, ibn)*0.1185)*&
1670 & tempb2
1671  CALL popreal8(ssaa_dev(k, ibn))
1672  tempb3 = ssaa_devb(k, ibn)/taua_dev(k, ibn)
1673  taua_devb(k, ibn) = taua_devb(k, ibn) - ssaa_dev(k, ibn)*&
1674 & tempb3/taua_dev(k, ibn)
1675  CALL popreal8(asya_dev(k, ibn))
1676  tempb4 = asya_devb(k, ibn)/ssaa_dev(k, ibn)
1677  ssaa_devb(k, ibn) = tempb3 - asya_dev(k, ibn)*tempb4/&
1678 & ssaa_dev(k, ibn)
1679  asya_devb(k, ibn) = tempb4
1680  END IF
1681  END IF
1682  CALL popreal8(taerlyr(k))
1683  taerlyrb(k) = 0.0_8
1684  END DO
1685  CALL popreal8(taerlyr(0))
1686  taerlyrb(0) = 0.0_8
1687  END IF
1688  CALL popinteger4array(icx, np + 1)
1689  CALL popinteger4array(ncld, 3)
1690  DO k=np,0,-1
1691  CALL popinteger4(icx(k))
1692  END DO
1693  CALL popreal8array(tcldlyr, np + 1)
1694  CALL popreal8array(enn, np + 1)
1695  CALL getirtau1_b(ibn, np, dp_pa, fcld_col, fcld_colb, reff_col, &
1696 & reff_colb, cwc_col, cwc_colb, tcldlyr, tcldlyrb, enn, &
1697 & ennb, aib_ir, awb_ir, aiw_ir, aww_ir, aig_ir, awg_ir, &
1698 & cons_grav)
1699  CALL popreal8(blevel(np+1))
1700  CALL planck_b(ibn, cb, tb_dev, tb_devb, blevel(np+1), blevelb(np+1&
1701 & ))
1702  blevelb(np+1) = 0.0_8
1703  CALL popreal8(blevel(0))
1704  blevelb(1) = blevelb(1) + blevelb(0)
1705  blevelb(0) = 0.0_8
1706  CALL popreal8(blevel(1))
1707  tempb0 = dp(1)*blevelb(1)/(dp(1)+dp(2))
1708  blayerb(1) = blayerb(1) + tempb0 + blevelb(1)
1709  blayerb(2) = blayerb(2) - tempb0
1710  blevelb(1) = 0.0_8
1711  DO k=np,2,-1
1712  CALL popreal8(blevel(k))
1713  tempb = blevelb(k)/(dp(k-1)+dp(k))
1714  blayerb(k-1) = blayerb(k-1) + dp(k)*tempb
1715  blayerb(k) = blayerb(k) + dp(k-1)*tempb
1716  blevelb(k) = 0.0_8
1717  END DO
1718  blayerb(np+1) = 0.0_8
1719  CALL popreal8(dbs)
1720  CALL popreal8(rflxs)
1721  CALL popreal8(blevel(0))
1722  blayerb(1) = blayerb(1) + blayerb(0) + blevelb(0)
1723  blevelb(0) = 0.0_8
1724  blayerb(0) = 0.0_8
1725  DO k=np,1,-1
1726  CALL planck_b(ibn, cb, ta_dev(k), ta_devb(k), blayer(k), blayerb&
1727 & (k))
1728  blayerb(k) = 0.0_8
1729  END DO
1730  CALL popcontrol4b(branch)
1731  IF (branch .LT. 4) THEN
1732  IF (branch .LT. 2) THEN
1733  IF (branch .EQ. 0) THEN
1734  CALL popinteger4(n2oexp%end)
1735  CALL popinteger4(n2oexp%start)
1736  CALL popinteger4(co2exp%end)
1737  CALL popinteger4(co2exp%start)
1738  CALL popinteger4(conexp%end)
1739  CALL popinteger4(conexp%start)
1740  CALL popinteger4(h2oexp%end)
1741  CALL popinteger4(h2oexp%start)
1742  ELSE
1743  CALL popinteger4(h2oexp%end)
1744  CALL popinteger4(h2oexp%start)
1745  END IF
1746  ELSE IF (branch .EQ. 2) THEN
1747  CALL popinteger4(ch4exp%end)
1748  CALL popinteger4(ch4exp%start)
1749  CALL popinteger4(n2oexp%end)
1750  CALL popinteger4(n2oexp%start)
1751  CALL popinteger4(conexp%end)
1752  CALL popinteger4(conexp%start)
1753  CALL popinteger4(h2oexp%end)
1754  CALL popinteger4(h2oexp%start)
1755  ELSE
1756  CALL popinteger4(f22exp%end)
1757  CALL popinteger4(f22exp%start)
1758  CALL popinteger4(f12exp%end)
1759  CALL popinteger4(f12exp%start)
1760  CALL popinteger4(ch4exp%end)
1761  CALL popinteger4(ch4exp%start)
1762  CALL popinteger4(n2oexp%end)
1763  CALL popinteger4(n2oexp%start)
1764  CALL popinteger4(conexp%end)
1765  CALL popinteger4(conexp%start)
1766  CALL popinteger4(h2oexp%end)
1767  CALL popinteger4(h2oexp%start)
1768  END IF
1769  ELSE IF (branch .LT. 6) THEN
1770  IF (branch .EQ. 4) THEN
1771  CALL popinteger4(f11exp%end)
1772  CALL popinteger4(f11exp%start)
1773  CALL popinteger4(comexp%end)
1774  CALL popinteger4(comexp%start)
1775  CALL popinteger4(conexp%end)
1776  CALL popinteger4(conexp%start)
1777  CALL popinteger4(h2oexp%end)
1778  CALL popinteger4(h2oexp%start)
1779  ELSE
1780  CALL popinteger4(f22exp%end)
1781  CALL popinteger4(f22exp%start)
1782  CALL popinteger4(f12exp%end)
1783  CALL popinteger4(f12exp%start)
1784  CALL popinteger4(f11exp%end)
1785  CALL popinteger4(f11exp%start)
1786  CALL popinteger4(comexp%end)
1787  CALL popinteger4(comexp%start)
1788  CALL popinteger4(conexp%end)
1789  CALL popinteger4(conexp%start)
1790  CALL popinteger4(h2oexp%end)
1791  CALL popinteger4(h2oexp%start)
1792  END IF
1793  ELSE IF (branch .EQ. 6) THEN
1794  CALL popinteger4(conexp%end)
1795  CALL popinteger4(conexp%start)
1796  CALL popinteger4(h2oexp%end)
1797  CALL popinteger4(h2oexp%start)
1798  ELSE IF (branch .EQ. 7) THEN
1799  CALL popinteger4(conexp%end)
1800  CALL popinteger4(conexp%start)
1801  END IF
1802  CALL popreal8array(exptbl, (np+1)*17)
1803  END IF
1804  CALL popinteger4(ibn)
1805  END DO
1806  DO k=np+1,1,-1
1807  dfdts_devb(k) = 0.0_8
1808  flxd_devb(k) = 0.0_8
1809  flxu_devb(k) = 0.0_8
1810  END DO
1811  xx = pa(0)*0.001618*wa_dev(1)*wa_dev(1)*dp(0)
1812  temp0 = 1800./ta_dev(1)
1813  xxb = exp(temp0-6.081)*dcontb(0)
1814  ta_devb(1) = ta_devb(1) - exp(temp0-6.081)*xx*temp0*dcontb(0)/ta_dev(1&
1815 & )
1816  dcontb(0) = 0.0_8
1817  wa_devb(1) = wa_devb(1) + pa(0)*0.001618*dp(0)*2*wa_dev(1)*xxb
1818  CALL popcontrol1b(branch)
1819  IF (branch .EQ. 0) do3b(0) = 0.0_8
1820  CALL popcontrol1b(branch)
1821  IF (branch .EQ. 0) dh2ob(0) = 0.0_8
1822  oa_devb(1) = oa_devb(1) + dp(0)*476.*do3b(0)
1823  do3b(0) = 0.0_8
1824  wa_devb(1) = wa_devb(1) + dp(0)*1.02*dh2ob(0)
1825  dh2ob(0) = 0.0_8
1826  ta_devb(1) = ta_devb(1) + dtb(0)
1827  dtb(0) = 0.0_8
1828  CALL popreal8(pa(0))
1829  CALL popcontrol1b(branch)
1830  IF (branch .EQ. 0) THEN
1831  CALL popreal8(dp(0))
1832  ELSE
1833  CALL popreal8(dp(0))
1834  END IF
1835  DO k=np,1,-1
1836  DO l=4,1,-1
1837  cwc_devb(k, l) = cwc_devb(k, l) + cwc_colb(k, l)
1838  cwc_colb(k, l) = 0.0_8
1839  reff_devb(k, l) = reff_devb(k, l) + reff_colb(k, l)
1840  reff_colb(k, l) = 0.0_8
1841  END DO
1842  fcld_devb(k) = fcld_devb(k) + fcld_colb(k)
1843  fcld_colb(k) = 0.0_8
1844  xx = pa(k)*0.001618*wa_dev(k)*wa_dev(k)*dp(k)
1845  temp = 1800./ta_dev(k)
1846  xxb = exp(temp-6.081)*dcontb(k)
1847  ta_devb(k) = ta_devb(k) - exp(temp-6.081)*xx*temp*dcontb(k)/ta_dev(k&
1848 & )
1849  dcontb(k) = 0.0_8
1850  wa_devb(k) = wa_devb(k) + pa(k)*0.001618*dp(k)*2*wa_dev(k)*xxb
1851  CALL popcontrol1b(branch)
1852  IF (branch .EQ. 0) do3b(k) = 0.0_8
1853  CALL popcontrol1b(branch)
1854  IF (branch .EQ. 0) dh2ob(k) = 0.0_8
1855  oa_devb(k) = oa_devb(k) + dp(k)*476.*do3b(k)
1856  do3b(k) = 0.0_8
1857  wa_devb(k) = wa_devb(k) + dp(k)*1.02*dh2ob(k)
1858  dh2ob(k) = 0.0_8
1859  ta_devb(k) = ta_devb(k) + dtb(k)
1860  dtb(k) = 0.0_8
1861  END DO
1862 END SUBROUTINE irrad_b
1863 
1864 ! Differentiation of planck in reverse (adjoint) mode:
1865 ! gradient of useful results: t xlayer
1866 ! with respect to varying inputs: t
1867 !***********************************************************************
1868 SUBROUTINE planck_b(ibn, cb, t, tb, xlayer, xlayerb)
1869  IMPLICIT NONE
1870 ! spectral band index
1871  INTEGER, INTENT(IN) :: ibn
1872 ! Planck table coefficients
1873  REAL*8, INTENT(IN) :: cb(6, 10)
1874 ! temperature (K)
1875  REAL*8, INTENT(IN) :: t
1876  REAL*8 :: tb
1877 ! planck flux (w/m2)
1878  REAL*8 :: xlayer
1879  REAL*8 :: xlayerb
1880  REAL*8 :: temp1
1881  REAL*8 :: temp0
1882  REAL*8 :: tempb
1883  REAL*8 :: temp
1884  temp1 = cb(5, ibn) + cb(6, ibn)*t
1885  temp0 = cb(4, ibn) + t*temp1
1886  temp = cb(3, ibn) + t*temp0
1887  tempb = t**2*xlayerb
1888  tb = tb + (t**2*cb(6, ibn)+t*temp1+temp0)*tempb + (2*(t*temp)+cb(2, &
1889 & ibn))*xlayerb
1890 END SUBROUTINE planck_b
1891 
1892 ! Differentiation of h2oexps in reverse (adjoint) mode:
1893 ! gradient of useful results: dh2o dt h2oexp
1894 ! with respect to varying inputs: dh2o dt
1895 !**********************************************************************
1896 SUBROUTINE h2oexps_b(ib, np, dh2o, dh2ob, pa, dt, dtb, xkw, aw, bw, pm, &
1897 & mw, h2oexp, h2oexpb)
1898  IMPLICIT NONE
1899  INTEGER :: ib, np, ik, k
1900 !---- input parameters ------
1901  REAL*8 :: dh2o(0:np), pa(0:np), dt(0:np)
1902  REAL*8 :: dh2ob(0:np), dtb(0:np)
1903 !---- output parameters -----
1904  REAL*8 :: h2oexp(0:np, 6)
1905  REAL*8 :: h2oexpb(0:np, 6)
1906 !---- static data -----
1907  INTEGER :: mw(9)
1908  REAL*8 :: xkw(9), aw(9), bw(9), pm(9)
1909 !---- temporary arrays -----
1910  REAL*8 :: xh
1911  REAL*8 :: xhb
1912  INTRINSIC exp
1913  INTEGER :: branch
1914  REAL*8 :: tempb
1915  REAL*8 :: temp
1916 !**********************************************************************
1917 ! note that the 3 sub-bands in band 3 use the same set of xkw, aw,
1918 ! and bw, therefore, h2oexp for these sub-bands are identical.
1919 !**********************************************************************
1920  DO k=0,np
1921 !-----xh is the scaled water vapor amount for line absorption
1922 ! computed from Eq. (4.4).
1923  CALL pushreal8(xh)
1924  xh = dh2o(k)*(pa(k)/500.)**pm(ib)*(1.+(aw(ib)+bw(ib)*dt(k))*dt(k))
1925 !-----h2oexp is the water vapor transmittance of the layer k
1926 ! due to line absorption
1927  h2oexp(k, 1) = exp(-(xh*xkw(ib)))
1928 !-----compute transmittances from Eq. (8.22)
1929  DO ik=2,6
1930  IF (mw(ib) .EQ. 6) THEN
1931  CALL pushreal8(xh)
1932  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)
1933  CALL pushreal8(h2oexp(k, ik))
1934  h2oexp(k, ik) = xh*xh*xh
1935  CALL pushcontrol2b(3)
1936  ELSE IF (mw(ib) .EQ. 8) THEN
1937  CALL pushreal8(xh)
1938  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)
1939  CALL pushreal8(xh)
1940  xh = xh*xh
1941  CALL pushreal8(h2oexp(k, ik))
1942  h2oexp(k, ik) = xh*xh
1943  CALL pushcontrol2b(2)
1944  ELSE IF (mw(ib) .EQ. 9) THEN
1945  CALL pushreal8(xh)
1946  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)*h2oexp(k, ik-1)
1947  CALL pushreal8(h2oexp(k, ik))
1948  h2oexp(k, ik) = xh*xh*xh
1949  CALL pushcontrol2b(1)
1950  ELSE
1951  CALL pushreal8(xh)
1952  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)
1953  CALL pushreal8(xh)
1954  xh = xh*xh
1955  CALL pushreal8(xh)
1956  xh = xh*xh
1957  CALL pushreal8(h2oexp(k, ik))
1958  h2oexp(k, ik) = xh*xh
1959  CALL pushcontrol2b(0)
1960  END IF
1961  END DO
1962  END DO
1963  DO k=np,0,-1
1964  DO ik=6,2,-1
1965  CALL popcontrol2b(branch)
1966  IF (branch .LT. 2) THEN
1967  IF (branch .EQ. 0) THEN
1968  CALL popreal8(h2oexp(k, ik))
1969  xhb = 2*xh*h2oexpb(k, ik)
1970  h2oexpb(k, ik) = 0.0_8
1971  CALL popreal8(xh)
1972  xhb = 2*xh*xhb
1973  CALL popreal8(xh)
1974  xhb = 2*xh*xhb
1975  CALL popreal8(xh)
1976  h2oexpb(k, ik-1) = h2oexpb(k, ik-1) + 2*h2oexp(k, ik-1)*xhb
1977  ELSE
1978  CALL popreal8(h2oexp(k, ik))
1979  xhb = 3*xh**2*h2oexpb(k, ik)
1980  h2oexpb(k, ik) = 0.0_8
1981  CALL popreal8(xh)
1982  h2oexpb(k, ik-1) = h2oexpb(k, ik-1) + 3*h2oexp(k, ik-1)**2*xhb
1983  END IF
1984  ELSE IF (branch .EQ. 2) THEN
1985  CALL popreal8(h2oexp(k, ik))
1986  xhb = 2*xh*h2oexpb(k, ik)
1987  h2oexpb(k, ik) = 0.0_8
1988  CALL popreal8(xh)
1989  xhb = 2*xh*xhb
1990  CALL popreal8(xh)
1991  h2oexpb(k, ik-1) = h2oexpb(k, ik-1) + 2*h2oexp(k, ik-1)*xhb
1992  ELSE
1993  CALL popreal8(h2oexp(k, ik))
1994  xhb = 3*xh**2*h2oexpb(k, ik)
1995  h2oexpb(k, ik) = 0.0_8
1996  CALL popreal8(xh)
1997  h2oexpb(k, ik-1) = h2oexpb(k, ik-1) + 2*h2oexp(k, ik-1)*xhb
1998  END IF
1999  END DO
2000  xhb = -(exp(-(xkw(ib)*xh))*xkw(ib)*h2oexpb(k, 1))
2001  h2oexpb(k, 1) = 0.0_8
2002  CALL popreal8(xh)
2003  temp = aw(ib) + bw(ib)*dt(k)
2004  tempb = (pa(k)/500.)**pm(ib)*xhb
2005  dh2ob(k) = dh2ob(k) + (temp*dt(k)+1.)*tempb
2006  dtb(k) = dtb(k) + (dh2o(k)*temp+dt(k)*dh2o(k)*bw(ib))*tempb
2007  END DO
2008 END SUBROUTINE h2oexps_b
2009 
2010 ! Differentiation of conexps in reverse (adjoint) mode:
2011 ! gradient of useful results: dcont conexp
2012 ! with respect to varying inputs: dcont conexp
2013 !**********************************************************************
2014 SUBROUTINE conexps_b(ib, np, dcont, dcontb, xke, conexp, conexpb)
2015  IMPLICIT NONE
2016  INTEGER :: ib, np, k
2017 !---- input parameters ------
2018  REAL*8 :: dcont(0:np)
2019  REAL*8 :: dcontb(0:np)
2020 !---- updated parameters -----
2021  REAL*8 :: conexp(0:np, 3)
2022  REAL*8 :: conexpb(0:np, 3)
2023 !---- static data -----
2024  REAL*8 :: xke(9)
2025  INTRINSIC exp
2026  INTEGER :: branch
2027 !****************************************************************
2028  DO k=0,np
2029  conexp(k, 1) = exp(-(dcont(k)*xke(ib)))
2030 !-----The absorption coefficients for sub-bands 3b and 3a are, respectively,
2031 ! two and four times the absorption coefficient for sub-band 3c (Table 9).
2032 ! Note that conexp(3) is for sub-band 3a.
2033  IF (ib .EQ. 3) THEN
2034  CALL pushreal8(conexp(k, 2))
2035  conexp(k, 2) = conexp(k, 1)*conexp(k, 1)
2036  CALL pushcontrol1b(1)
2037  ELSE
2038  CALL pushcontrol1b(0)
2039  END IF
2040  END DO
2041  DO k=np,0,-1
2042  CALL popcontrol1b(branch)
2043  IF (branch .NE. 0) THEN
2044  conexpb(k, 2) = conexpb(k, 2) + 2*conexp(k, 2)*conexpb(k, 3)
2045  conexpb(k, 3) = 0.0_8
2046  CALL popreal8(conexp(k, 2))
2047  conexpb(k, 1) = conexpb(k, 1) + 2*conexp(k, 1)*conexpb(k, 2)
2048  conexpb(k, 2) = 0.0_8
2049  END IF
2050  dcontb(k) = dcontb(k) - exp(-(xke(ib)*dcont(k)))*xke(ib)*conexpb(k, &
2051 & 1)
2052  conexpb(k, 1) = 0.0_8
2053  END DO
2054 END SUBROUTINE conexps_b
2055 
2056 ! Differentiation of n2oexps in reverse (adjoint) mode:
2057 ! gradient of useful results: dt n2oexp
2058 ! with respect to varying inputs: dt n2oexp
2059 !**********************************************************************
2060 SUBROUTINE n2oexps_b(ib, np, dn2o, pa, dt, dtb, n2oexp, n2oexpb)
2061  IMPLICIT NONE
2062  INTEGER :: ib, k, np
2063 !---- input parameters -----
2064  REAL*8 :: dn2o(0:np), pa(0:np), dt(0:np)
2065  REAL*8 :: dtb(0:np)
2066 !---- output parameters -----
2067  REAL*8 :: n2oexp(0:np, 4)
2068  REAL*8 :: n2oexpb(0:np, 4)
2069 !---- temporary arrays -----
2070  REAL*8 :: xc, xc1, xc2
2071  REAL*8 :: xcb, xc1b, xc2b
2072  INTRINSIC exp
2073  INTEGER :: branch
2074  REAL*8 :: tempb
2075 !-----Scaling and absorption data are given in Table 5.
2076 ! Transmittances are computed using Eqs. (8.21) and (8.22).
2077  DO k=0,np
2078 !-----four exponential by powers of 21 for band 6.
2079  IF (ib .EQ. 6) THEN
2080  CALL pushreal8(xc)
2081  xc = dn2o(k)*(1.+(1.9297e-3+4.3750e-6*dt(k))*dt(k))
2082  n2oexp(k, 1) = exp(-(xc*6.31582e-2))
2083  xc = n2oexp(k, 1)*n2oexp(k, 1)*n2oexp(k, 1)
2084 !-----four exponential by powers of 8 for band 7
2085  CALL pushcontrol1b(1)
2086  ELSE
2087  CALL pushreal8(xc)
2088  xc = dn2o(k)*(pa(k)/500.0)**0.48*(1.+(1.3804e-3+7.4838e-6*dt(k))*&
2089 & dt(k))
2090  n2oexp(k, 1) = exp(-(xc*5.35779e-2))
2091  xc = n2oexp(k, 1)*n2oexp(k, 1)
2092  CALL pushreal8(xc)
2093  xc = xc*xc
2094  CALL pushreal8(n2oexp(k, 2))
2095  n2oexp(k, 2) = xc*xc
2096  CALL pushreal8(xc)
2097  xc = n2oexp(k, 2)*n2oexp(k, 2)
2098  CALL pushreal8(xc)
2099  xc = xc*xc
2100  CALL pushreal8(n2oexp(k, 3))
2101  n2oexp(k, 3) = xc*xc
2102  CALL pushreal8(xc)
2103  xc = n2oexp(k, 3)*n2oexp(k, 3)
2104  CALL pushreal8(xc)
2105  xc = xc*xc
2106  CALL pushcontrol1b(0)
2107  END IF
2108  END DO
2109  DO k=np,0,-1
2110  CALL popcontrol1b(branch)
2111  IF (branch .EQ. 0) THEN
2112  xcb = 2*xc*n2oexpb(k, 4)
2113  n2oexpb(k, 4) = 0.0_8
2114  CALL popreal8(xc)
2115  xcb = 2*xc*xcb
2116  CALL popreal8(xc)
2117  n2oexpb(k, 3) = n2oexpb(k, 3) + 2*n2oexp(k, 3)*xcb
2118  CALL popreal8(n2oexp(k, 3))
2119  xcb = 2*xc*n2oexpb(k, 3)
2120  n2oexpb(k, 3) = 0.0_8
2121  CALL popreal8(xc)
2122  xcb = 2*xc*xcb
2123  CALL popreal8(xc)
2124  n2oexpb(k, 2) = n2oexpb(k, 2) + 2*n2oexp(k, 2)*xcb
2125  CALL popreal8(n2oexp(k, 2))
2126  xcb = 2*xc*n2oexpb(k, 2)
2127  n2oexpb(k, 2) = 0.0_8
2128  CALL popreal8(xc)
2129  xcb = 2*xc*xcb
2130  n2oexpb(k, 1) = n2oexpb(k, 1) + 2*n2oexp(k, 1)*xcb
2131  xc = dn2o(k)*(pa(k)/500.0)**0.48*(1.+(1.3804e-3+7.4838e-6*dt(k))*&
2132 & dt(k))
2133  xcb = -(exp(-(5.35779e-2*xc))*5.35779e-2*n2oexpb(k, 1))
2134  n2oexpb(k, 1) = 0.0_8
2135  CALL popreal8(xc)
2136  tempb = (pa(k)/500.0)**0.48*dn2o(k)*xcb
2137  dtb(k) = dtb(k) + (7.4838e-6*dt(k)+dt(k)*7.4838e-6+1.3804e-3)*&
2138 & tempb
2139  ELSE
2140  xc1 = xc*xc
2141  xc2 = xc1*xc1
2142  xc2b = xc*xc1*n2oexpb(k, 2)
2143  xc1b = 2*xc1*xc2b + xc2*xc*n2oexpb(k, 2)
2144  xcb = 2*xc*xc1b + xc2*xc1*n2oexpb(k, 2)
2145  n2oexpb(k, 2) = 0.0_8
2146  n2oexpb(k, 1) = n2oexpb(k, 1) + 3*n2oexp(k, 1)**2*xcb
2147  xc = dn2o(k)*(1.+(1.9297e-3+4.3750e-6*dt(k))*dt(k))
2148  xcb = -(exp(-(6.31582e-2*xc))*6.31582e-2*n2oexpb(k, 1))
2149  n2oexpb(k, 1) = 0.0_8
2150  CALL popreal8(xc)
2151  dtb(k) = dtb(k) + (dn2o(k)*(4.3750e-6*dt(k)+1.9297e-3)+dt(k)*dn2o(&
2152 & k)*4.3750e-6)*xcb
2153  END IF
2154  END DO
2155 END SUBROUTINE n2oexps_b
2156 
2157 ! Differentiation of ch4exps in reverse (adjoint) mode:
2158 ! gradient of useful results: dt ch4exp
2159 ! with respect to varying inputs: dt ch4exp
2160 !**********************************************************************
2161 SUBROUTINE ch4exps_b(ib, np, dch4, pa, dt, dtb, ch4exp, ch4expb)
2162  IMPLICIT NONE
2163  INTEGER :: ib, np, k
2164 !---- input parameters -----
2165  REAL*8 :: dch4(0:np), pa(0:np), dt(0:np)
2166  REAL*8 :: dtb(0:np)
2167 !---- output parameters -----
2168  REAL*8 :: ch4exp(0:np, 4)
2169  REAL*8 :: ch4expb(0:np, 4)
2170 !---- temporary arrays -----
2171  REAL*8 :: xc
2172  REAL*8 :: xcb
2173  INTRINSIC exp
2174  INTEGER :: branch
2175  REAL*8 :: tempb
2176 !***** Scaling and absorption data are given in Table 5 *****
2177  DO k=0,np
2178 !-----four exponentials for band 6
2179  IF (ib .EQ. 6) THEN
2180  CALL pushreal8(xc)
2181 !-----four exponentials by powers of 12 for band 7
2182  CALL pushcontrol1b(1)
2183  ELSE
2184  CALL pushreal8(xc)
2185  xc = dch4(k)*(pa(k)/500.0)**0.65*(1.+(5.9590e-4-2.2931e-6*dt(k))*&
2186 & dt(k))
2187  ch4exp(k, 1) = exp(-(xc*6.29247e-2))
2188  xc = ch4exp(k, 1)*ch4exp(k, 1)*ch4exp(k, 1)
2189  CALL pushreal8(xc)
2190  xc = xc*xc
2191  CALL pushreal8(ch4exp(k, 2))
2192  ch4exp(k, 2) = xc*xc
2193  CALL pushreal8(xc)
2194  xc = ch4exp(k, 2)*ch4exp(k, 2)*ch4exp(k, 2)
2195  CALL pushreal8(xc)
2196  xc = xc*xc
2197  CALL pushreal8(ch4exp(k, 3))
2198  ch4exp(k, 3) = xc*xc
2199  CALL pushreal8(xc)
2200  xc = ch4exp(k, 3)*ch4exp(k, 3)*ch4exp(k, 3)
2201  CALL pushreal8(xc)
2202  xc = xc*xc
2203  CALL pushcontrol1b(0)
2204  END IF
2205  END DO
2206  DO k=np,0,-1
2207  CALL popcontrol1b(branch)
2208  IF (branch .EQ. 0) THEN
2209  xcb = 2*xc*ch4expb(k, 4)
2210  ch4expb(k, 4) = 0.0_8
2211  CALL popreal8(xc)
2212  xcb = 2*xc*xcb
2213  CALL popreal8(xc)
2214  ch4expb(k, 3) = ch4expb(k, 3) + 3*ch4exp(k, 3)**2*xcb
2215  CALL popreal8(ch4exp(k, 3))
2216  xcb = 2*xc*ch4expb(k, 3)
2217  ch4expb(k, 3) = 0.0_8
2218  CALL popreal8(xc)
2219  xcb = 2*xc*xcb
2220  CALL popreal8(xc)
2221  ch4expb(k, 2) = ch4expb(k, 2) + 3*ch4exp(k, 2)**2*xcb
2222  CALL popreal8(ch4exp(k, 2))
2223  xcb = 2*xc*ch4expb(k, 2)
2224  ch4expb(k, 2) = 0.0_8
2225  CALL popreal8(xc)
2226  xcb = 2*xc*xcb
2227  ch4expb(k, 1) = ch4expb(k, 1) + 3*ch4exp(k, 1)**2*xcb
2228  xc = dch4(k)*(pa(k)/500.0)**0.65*(1.+(5.9590e-4-2.2931e-6*dt(k))*&
2229 & dt(k))
2230  xcb = -(exp(-(6.29247e-2*xc))*6.29247e-2*ch4expb(k, 1))
2231  ch4expb(k, 1) = 0.0_8
2232  CALL popreal8(xc)
2233  tempb = (pa(k)/500.0)**0.65*dch4(k)*xcb
2234  dtb(k) = dtb(k) + (5.9590e-4-dt(k)*2.2931e-6-2.2931e-6*dt(k))*&
2235 & tempb
2236  ELSE
2237  xc = dch4(k)*(1.+(1.7007e-2+1.5826e-4*dt(k))*dt(k))
2238  xcb = -(exp(-(5.80708e-3*xc))*5.80708e-3*ch4expb(k, 1))
2239  ch4expb(k, 1) = 0.0_8
2240  CALL popreal8(xc)
2241  dtb(k) = dtb(k) + (dch4(k)*(1.5826e-4*dt(k)+1.7007e-2)+dt(k)*dch4(&
2242 & k)*1.5826e-4)*xcb
2243  END IF
2244  END DO
2245 END SUBROUTINE ch4exps_b
2246 
2247 ! Differentiation of comexps in reverse (adjoint) mode:
2248 ! gradient of useful results: dt comexp
2249 ! with respect to varying inputs: dt comexp
2250 !**********************************************************************
2251 SUBROUTINE comexps_b(ib, np, dcom, dt, dtb, comexp, comexpb)
2252  IMPLICIT NONE
2253  INTEGER :: ib, ik, np, k
2254 !---- input parameters -----
2255  REAL*8 :: dcom(0:np), dt(0:np)
2256  REAL*8 :: dtb(0:np)
2257 !---- output parameters -----
2258  REAL*8 :: comexp(0:np, 6)
2259  REAL*8 :: comexpb(0:np, 6)
2260 !---- temporary arrays -----
2261  REAL*8 :: xc
2262  REAL*8 :: xcb
2263  INTRINSIC exp
2264  INTEGER :: branch
2265 !***** Scaling and absorpton data are given in Table 6 *****
2266  DO k=0,np
2267  IF (ib .EQ. 4) THEN
2268  CALL pushreal8(xc)
2269  xc = dcom(k)*(1.+(3.5775e-2+4.0447e-4*dt(k))*dt(k))
2270  CALL pushcontrol1b(0)
2271  ELSE
2272  CALL pushcontrol1b(1)
2273  END IF
2274  IF (ib .EQ. 5) THEN
2275  CALL pushreal8(xc)
2276  xc = dcom(k)*(1.+(3.4268e-2+3.7401e-4*dt(k))*dt(k))
2277  CALL pushcontrol1b(0)
2278  ELSE
2279  CALL pushcontrol1b(1)
2280  END IF
2281  comexp(k, 1) = exp(-(xc*1.922e-7))
2282  DO ik=2,6
2283  CALL pushreal8(xc)
2284  xc = comexp(k, ik-1)*comexp(k, ik-1)
2285  CALL pushreal8(xc)
2286  xc = xc*xc
2287  CALL pushreal8(comexp(k, ik))
2288  comexp(k, ik) = xc*comexp(k, ik-1)
2289  END DO
2290  END DO
2291  xcb = 0.0_8
2292  DO k=np,0,-1
2293  DO ik=6,2,-1
2294  CALL popreal8(comexp(k, ik))
2295  xcb = xcb + comexp(k, ik-1)*comexpb(k, ik)
2296  comexpb(k, ik-1) = comexpb(k, ik-1) + xc*comexpb(k, ik)
2297  comexpb(k, ik) = 0.0_8
2298  CALL popreal8(xc)
2299  xcb = 2*xc*xcb
2300  CALL popreal8(xc)
2301  comexpb(k, ik-1) = comexpb(k, ik-1) + 2*comexp(k, ik-1)*xcb
2302  xcb = 0.0_8
2303  END DO
2304  xcb = xcb - exp(-(1.922e-7*xc))*1.922e-7*comexpb(k, 1)
2305  comexpb(k, 1) = 0.0_8
2306  CALL popcontrol1b(branch)
2307  IF (branch .EQ. 0) THEN
2308  CALL popreal8(xc)
2309  dtb(k) = dtb(k) + (dcom(k)*(3.7401e-4*dt(k)+3.4268e-2)+dt(k)*dcom(&
2310 & k)*3.7401e-4)*xcb
2311  xcb = 0.0_8
2312  END IF
2313  CALL popcontrol1b(branch)
2314  IF (branch .EQ. 0) THEN
2315  CALL popreal8(xc)
2316  dtb(k) = dtb(k) + (dcom(k)*(4.0447e-4*dt(k)+3.5775e-2)+dt(k)*dcom(&
2317 & k)*4.0447e-4)*xcb
2318  xcb = 0.0_8
2319  END IF
2320  END DO
2321 END SUBROUTINE comexps_b
2322 
2323 ! Differentiation of cfcexps in reverse (adjoint) mode:
2324 ! gradient of useful results: dt cfcexp
2325 ! with respect to varying inputs: dt cfcexp
2326 !**********************************************************************
2327 SUBROUTINE cfcexps_b(ib, np, a1, b1, fk1, a2, b2, fk2, dcfc, dt, dtb, &
2328 & cfcexp, cfcexpb)
2329  IMPLICIT NONE
2330  INTEGER :: ib, np, k
2331 !---- input parameters -----
2332  REAL*8 :: dcfc(0:np), dt(0:np)
2333  REAL*8 :: dtb(0:np)
2334 !---- output parameters -----
2335  REAL*8 :: cfcexp(0:np)
2336  REAL*8 :: cfcexpb(0:np)
2337 !---- static data -----
2338  REAL*8 :: a1, b1, fk1, a2, b2, fk2
2339 !---- temporary arrays -----
2340  REAL*8 :: xf
2341  REAL*8 :: xfb
2342  INTRINSIC exp
2343  INTEGER :: branch
2344 !**********************************************************************
2345  DO k=0,np
2346 !-----compute the scaled cfc amount (xf) and exponential (cfcexp)
2347  IF (ib .EQ. 4) THEN
2348  CALL pushcontrol1b(1)
2349  ELSE
2350  CALL pushcontrol1b(0)
2351  END IF
2352  END DO
2353  DO k=np,0,-1
2354  CALL popcontrol1b(branch)
2355  IF (branch .EQ. 0) THEN
2356  xf = dcfc(k)*(1.+(a2+b2*dt(k))*dt(k))
2357  xfb = -(exp(-(fk2*xf))*fk2*cfcexpb(k))
2358  cfcexpb(k) = 0.0_8
2359  dtb(k) = dtb(k) + (dcfc(k)*(a2+b2*dt(k))+dt(k)*dcfc(k)*b2)*xfb
2360  ELSE
2361  xf = dcfc(k)*(1.+(a1+b1*dt(k))*dt(k))
2362  xfb = -(exp(-(fk1*xf))*fk1*cfcexpb(k))
2363  cfcexpb(k) = 0.0_8
2364  dtb(k) = dtb(k) + (dcfc(k)*(a1+b1*dt(k))+dt(k)*dcfc(k)*b1)*xfb
2365  END IF
2366  END DO
2367 END SUBROUTINE cfcexps_b
2368 
2369 ! Differentiation of b10exps in reverse (adjoint) mode:
2370 ! gradient of useful results: dh2o dcont dt co2exp h2oexp
2371 ! n2oexp conexp
2372 ! with respect to varying inputs: dh2o dcont dt co2exp h2oexp
2373 ! n2oexp conexp
2374 !**********************************************************************
2375 SUBROUTINE b10exps_b(np, dh2o, dh2ob, dcont, dcontb, dco2, dn2o, pa, dt&
2376 & , dtb, h2oexp, h2oexpb, conexp, conexpb, co2exp, co2expb, n2oexp, &
2377 & n2oexpb)
2378  IMPLICIT NONE
2379  INTEGER :: np, k
2380 !---- input parameters -----
2381  REAL*8 :: dh2o(0:np), dcont(0:np), dn2o(0:np)
2382  REAL*8 :: dh2ob(0:np), dcontb(0:np)
2383  REAL*8 :: dco2(0:np), pa(0:np), dt(0:np)
2384  REAL*8 :: dtb(0:np)
2385 !---- output parameters -----
2386  REAL*8 :: h2oexp(0:np, 5), conexp(0:np), co2exp(0:np, 6), n2oexp(0:np&
2387 & , 2)
2388  REAL*8 :: h2oexpb(0:np, 5), conexpb(0:np), co2expb(0:np, 6), n2oexpb(0&
2389 & :np, 2)
2390 !---- temporary arrays -----
2391  REAL*8 :: xx, xx1, xx2, xx3
2392  REAL*8 :: xxb, xx1b, xx2b, xx3b
2393  INTRINSIC exp
2394  REAL*8 :: tempb2
2395  REAL*8 :: tempb1
2396  REAL*8 :: tempb0
2397  REAL*8 :: tempb
2398 !**********************************************************************
2399  DO k=0,np
2400 !-----Compute scaled h2o-line amount for Band 10 (Eq. 4.4 and Table 3).
2401  CALL pushreal8(xx)
2402  xx = dh2o(k)*(pa(k)/500.0)*(1.+(0.0149+6.20e-5*dt(k))*dt(k))
2403 !-----six exponentials by powers of 8
2404  h2oexp(k, 1) = exp(-(xx*0.10624))
2405  xx = h2oexp(k, 1)*h2oexp(k, 1)
2406  CALL pushreal8(xx)
2407  xx = xx*xx
2408  CALL pushreal8(h2oexp(k, 2))
2409  h2oexp(k, 2) = xx*xx
2410  CALL pushreal8(xx)
2411  xx = h2oexp(k, 2)*h2oexp(k, 2)
2412  CALL pushreal8(xx)
2413  xx = xx*xx
2414  CALL pushreal8(h2oexp(k, 3))
2415  h2oexp(k, 3) = xx*xx
2416  CALL pushreal8(xx)
2417  xx = h2oexp(k, 3)*h2oexp(k, 3)
2418  CALL pushreal8(xx)
2419  xx = xx*xx
2420  CALL pushreal8(h2oexp(k, 4))
2421  h2oexp(k, 4) = xx*xx
2422  CALL pushreal8(xx)
2423  xx = h2oexp(k, 4)*h2oexp(k, 4)
2424  CALL pushreal8(xx)
2425  xx = xx*xx
2426 !-----one exponential of h2o continuum for sub-band 3a (Table 9).
2427 !-----Scaled co2 amount for the Band 10 (Eq. 4.4, Tables 3 and 6).
2428  CALL pushreal8(xx)
2429  xx = dco2(k)*(pa(k)/300.0)**0.5*(1.+(0.0179+1.02e-4*dt(k))*dt(k))
2430 !-----six exponentials by powers of 8
2431  co2exp(k, 1) = exp(-(xx*2.656e-5))
2432  xx = co2exp(k, 1)*co2exp(k, 1)
2433  CALL pushreal8(xx)
2434  xx = xx*xx
2435  CALL pushreal8(co2exp(k, 2))
2436  co2exp(k, 2) = xx*xx
2437  CALL pushreal8(xx)
2438  xx = co2exp(k, 2)*co2exp(k, 2)
2439  CALL pushreal8(xx)
2440  xx = xx*xx
2441  CALL pushreal8(co2exp(k, 3))
2442  co2exp(k, 3) = xx*xx
2443  CALL pushreal8(xx)
2444  xx = co2exp(k, 3)*co2exp(k, 3)
2445  CALL pushreal8(xx)
2446  xx = xx*xx
2447  CALL pushreal8(co2exp(k, 4))
2448  co2exp(k, 4) = xx*xx
2449  CALL pushreal8(xx)
2450  xx = co2exp(k, 4)*co2exp(k, 4)
2451  CALL pushreal8(xx)
2452  xx = xx*xx
2453  CALL pushreal8(co2exp(k, 5))
2454  co2exp(k, 5) = xx*xx
2455  CALL pushreal8(xx)
2456  xx = co2exp(k, 5)*co2exp(k, 5)
2457  CALL pushreal8(xx)
2458  xx = xx*xx
2459 !-----Compute the scaled n2o amount for Band 10 (Table 5).
2460  CALL pushreal8(xx)
2461  xx = dn2o(k)*(1.+(1.4476e-3+3.6656e-6*dt(k))*dt(k))
2462 !-----Two exponentials by powers of 58
2463  n2oexp(k, 1) = exp(-(xx*0.25238))
2464  xx = n2oexp(k, 1)*n2oexp(k, 1)
2465  CALL pushreal8(xx1)
2466  xx1 = xx*xx
2467  CALL pushreal8(xx1)
2468  xx1 = xx1*xx1
2469  END DO
2470  DO k=np,0,-1
2471  xx2 = xx1*xx1
2472  xx3 = xx2*xx2
2473  tempb = xx2*xx3*n2oexpb(k, 2)
2474  tempb0 = xx*xx1*n2oexpb(k, 2)
2475  xxb = xx1*tempb
2476  xx3b = xx2*tempb0
2477  xx2b = 2*xx2*xx3b + xx3*tempb0
2478  xx1b = 2*xx1*xx2b + xx*tempb
2479  n2oexpb(k, 2) = 0.0_8
2480  CALL popreal8(xx1)
2481  xx1b = 2*xx1*xx1b
2482  CALL popreal8(xx1)
2483  xxb = xxb + 2*xx*xx1b
2484  n2oexpb(k, 1) = n2oexpb(k, 1) + 2*n2oexp(k, 1)*xxb
2485  xx = dn2o(k)*(1.+(1.4476e-3+3.6656e-6*dt(k))*dt(k))
2486  xxb = -(exp(-(0.25238*xx))*0.25238*n2oexpb(k, 1))
2487  n2oexpb(k, 1) = 0.0_8
2488  CALL popreal8(xx)
2489  dtb(k) = dtb(k) + (dn2o(k)*(3.6656e-6*dt(k)+1.4476e-3)+dt(k)*dn2o(k)&
2490 & *3.6656e-6)*xxb
2491  xxb = 2*xx*co2expb(k, 6)
2492  co2expb(k, 6) = 0.0_8
2493  CALL popreal8(xx)
2494  xxb = 2*xx*xxb
2495  CALL popreal8(xx)
2496  co2expb(k, 5) = co2expb(k, 5) + 2*co2exp(k, 5)*xxb
2497  CALL popreal8(co2exp(k, 5))
2498  xxb = 2*xx*co2expb(k, 5)
2499  co2expb(k, 5) = 0.0_8
2500  CALL popreal8(xx)
2501  xxb = 2*xx*xxb
2502  CALL popreal8(xx)
2503  co2expb(k, 4) = co2expb(k, 4) + 2*co2exp(k, 4)*xxb
2504  CALL popreal8(co2exp(k, 4))
2505  xxb = 2*xx*co2expb(k, 4)
2506  co2expb(k, 4) = 0.0_8
2507  CALL popreal8(xx)
2508  xxb = 2*xx*xxb
2509  CALL popreal8(xx)
2510  co2expb(k, 3) = co2expb(k, 3) + 2*co2exp(k, 3)*xxb
2511  CALL popreal8(co2exp(k, 3))
2512  xxb = 2*xx*co2expb(k, 3)
2513  co2expb(k, 3) = 0.0_8
2514  CALL popreal8(xx)
2515  xxb = 2*xx*xxb
2516  CALL popreal8(xx)
2517  co2expb(k, 2) = co2expb(k, 2) + 2*co2exp(k, 2)*xxb
2518  CALL popreal8(co2exp(k, 2))
2519  xxb = 2*xx*co2expb(k, 2)
2520  co2expb(k, 2) = 0.0_8
2521  CALL popreal8(xx)
2522  xxb = 2*xx*xxb
2523  co2expb(k, 1) = co2expb(k, 1) + 2*co2exp(k, 1)*xxb
2524  xx = dco2(k)*(pa(k)/300.0)**0.5*(1.+(0.0179+1.02e-4*dt(k))*dt(k))
2525  xxb = -(exp(-(2.656e-5*xx))*2.656e-5*co2expb(k, 1))
2526  co2expb(k, 1) = 0.0_8
2527  CALL popreal8(xx)
2528  tempb1 = (pa(k)/300.0)**0.5*dco2(k)*xxb
2529  dcontb(k) = dcontb(k) - exp(-(109.0*dcont(k)))*109.0*conexpb(k)
2530  conexpb(k) = 0.0_8
2531  xxb = 2*xx*h2oexpb(k, 5)
2532  h2oexpb(k, 5) = 0.0_8
2533  CALL popreal8(xx)
2534  xxb = 2*xx*xxb
2535  CALL popreal8(xx)
2536  h2oexpb(k, 4) = h2oexpb(k, 4) + 2*h2oexp(k, 4)*xxb
2537  CALL popreal8(h2oexp(k, 4))
2538  xxb = 2*xx*h2oexpb(k, 4)
2539  h2oexpb(k, 4) = 0.0_8
2540  CALL popreal8(xx)
2541  xxb = 2*xx*xxb
2542  CALL popreal8(xx)
2543  h2oexpb(k, 3) = h2oexpb(k, 3) + 2*h2oexp(k, 3)*xxb
2544  CALL popreal8(h2oexp(k, 3))
2545  xxb = 2*xx*h2oexpb(k, 3)
2546  h2oexpb(k, 3) = 0.0_8
2547  CALL popreal8(xx)
2548  xxb = 2*xx*xxb
2549  CALL popreal8(xx)
2550  h2oexpb(k, 2) = h2oexpb(k, 2) + 2*h2oexp(k, 2)*xxb
2551  CALL popreal8(h2oexp(k, 2))
2552  xxb = 2*xx*h2oexpb(k, 2)
2553  h2oexpb(k, 2) = 0.0_8
2554  CALL popreal8(xx)
2555  xxb = 2*xx*xxb
2556  h2oexpb(k, 1) = h2oexpb(k, 1) + 2*h2oexp(k, 1)*xxb
2557  xx = dh2o(k)*(pa(k)/500.0)*(1.+(0.0149+6.20e-5*dt(k))*dt(k))
2558  xxb = -(exp(-(0.10624*xx))*0.10624*h2oexpb(k, 1))
2559  h2oexpb(k, 1) = 0.0_8
2560  CALL popreal8(xx)
2561  tempb2 = pa(k)*dh2o(k)*xxb/500.0
2562  dtb(k) = dtb(k) + (6.20e-5*dt(k)+dt(k)*6.20e-5+0.0149)*tempb2 + (&
2563 & 1.02e-4*dt(k)+dt(k)*1.02e-4+0.0179)*tempb1
2564  dh2ob(k) = dh2ob(k) + ((6.20e-5*dt(k)+0.0149)*dt(k)+1.)*pa(k)*xxb/&
2565 & 500.0
2566  END DO
2567 END SUBROUTINE b10exps_b
2568 
2569 ! Differentiation of tablup in reverse (adjoint) mode:
2570 ! gradient of useful results: s1 s2 s3 dt dw tran
2571 ! with respect to varying inputs: s1 s2 s3 dt dw tran
2572 !**********************************************************************
2573 SUBROUTINE tablup_b(nx1, nh1, dw, dwb, p, dt, dtb, s1, s1b, s2, s2b, s3&
2574 & , s3b, w1, p1, dwe, dpe, coef1, coef2, coef3, tran, tranb)
2575  IMPLICIT NONE
2576  INTEGER :: nx1, nh1
2577 !---- input parameters -----
2578  REAL*8 :: w1, p1, dwe, dpe
2579  REAL*8 :: dw, p, dt
2580  REAL*8 :: dwb, dtb
2581  REAL*8 :: coef1(nx1, nh1), coef2(nx1, nh1), coef3(nx1, nh1)
2582 !---- update parameter -----
2583  REAL*8 :: s1, s2, s3, tran
2584  REAL*8 :: s1b, s2b, s3b, tranb
2585 !---- temporary variables -----
2586  REAL*8 :: we, pe, fw, fp, pa, pb, pc, ax, ba, bb, t1, ca, cb, t2
2587  REAL*8 :: web, peb, fwb, fpb, pab, pbb, pcb, axb, bab, bbb, t1b, cab, &
2588 & cbb, t2b
2589  REAL*8 :: x1, x2, x3, xx, x1c
2590  REAL*8 :: x1b, x2b, x3b, xxb, x1cb
2591  INTEGER :: iw, ip
2592  INTRINSIC log10
2593  INTRINSIC real
2594  INTRINSIC min
2595  INTRINSIC int
2596  INTRINSIC max
2597  INTEGER :: branch
2598  REAL*8 :: tempb0
2599  REAL*8 :: tempb
2600  REAL*8 :: y2
2601  REAL*8 :: y1
2602 !-----Compute effective pressure (x2) and temperature (x3) following
2603 ! Eqs. (8.28) and (8.29)
2604  s1 = s1 + dw
2605  s2 = s2 + p*dw
2606  s3 = s3 + dt*dw
2607  x1 = s1
2608  x1c = 1.0/s1
2609  x2 = s2*x1c
2610  x3 = s3*x1c
2611 !-----normalize we and pe
2612 ! we=(log10(x1)-w1)/dwe
2613 ! pe=(log10(x2)-p1)/dpe
2614  we = (log10(x1)-w1)*dwe
2615  pe = (log10(x2)-p1)*dpe
2616  y1 = REAL(nh1 - 1)
2617  IF (we .GT. y1) THEN
2618  we = y1
2619  CALL pushcontrol1b(0)
2620  ELSE
2621  CALL pushcontrol1b(1)
2622  we = we
2623  END IF
2624  y2 = REAL(nx1 - 1)
2625  IF (pe .GT. y2) THEN
2626  pe = y2
2627  CALL pushcontrol1b(0)
2628  ELSE
2629  CALL pushcontrol1b(1)
2630  pe = pe
2631  END IF
2632 !-----assign iw and ip and compute the distance of we and pe
2633 ! from iw and ip.
2634  iw = int(we + 1.0)
2635  IF (iw .GT. nh1 - 1) THEN
2636  iw = nh1 - 1
2637  ELSE
2638  iw = iw
2639  END IF
2640  IF (iw .LT. 2) THEN
2641  iw = 2
2642  ELSE
2643  iw = iw
2644  END IF
2645  fw = we - REAL(iw - 1)
2646  ip = int(pe + 1.0)
2647  IF (ip .GT. nx1 - 1) THEN
2648  ip = nx1 - 1
2649  ELSE
2650  ip = ip
2651  END IF
2652  IF (ip .LT. 1) THEN
2653  ip = 1
2654  ELSE
2655  ip = ip
2656  END IF
2657  fp = pe - REAL(ip - 1)
2658 !-----linear interpolation in pressure
2659  pa = coef1(ip, iw-1) + (coef1(ip+1, iw-1)-coef1(ip, iw-1))*fp
2660  pb = coef1(ip, iw) + (coef1(ip+1, iw)-coef1(ip, iw))*fp
2661  pc = coef1(ip, iw+1) + (coef1(ip+1, iw+1)-coef1(ip, iw+1))*fp
2662 !-----quadratic interpolation in absorber amount for coef1
2663  ax = ((pc+pa)*fw+(pc-pa))*fw*0.5 + pb*(1.-fw*fw)
2664 !-----linear interpolation in absorber amount for coef2 and coef3
2665  ba = coef2(ip, iw) + (coef2(ip+1, iw)-coef2(ip, iw))*fp
2666  bb = coef2(ip, iw+1) + (coef2(ip+1, iw+1)-coef2(ip, iw+1))*fp
2667  t1 = ba + (bb-ba)*fw
2668  ca = coef3(ip, iw) + (coef3(ip+1, iw)-coef3(ip, iw))*fp
2669  cb = coef3(ip, iw+1) + (coef3(ip+1, iw+1)-coef3(ip, iw+1))*fp
2670  t2 = ca + (cb-ca)*fw
2671 !-----update the total transmittance between levels k1 and k2
2672  xx = ax + (t1+t2*x3)*x3
2673  IF (xx .GT. 0.9999999) THEN
2674  xx = 0.9999999
2675  CALL pushcontrol1b(0)
2676  ELSE
2677  CALL pushcontrol1b(1)
2678  xx = xx
2679  END IF
2680  IF (xx .LT. 0.0000001) THEN
2681  xx = 0.0000001
2682  CALL pushcontrol1b(0)
2683  ELSE
2684  CALL pushcontrol1b(1)
2685  xx = xx
2686  END IF
2687  xxb = tran*tranb
2688  tranb = xx*tranb
2689  CALL popcontrol1b(branch)
2690  IF (branch .EQ. 0) xxb = 0.0_8
2691  CALL popcontrol1b(branch)
2692  IF (branch .EQ. 0) xxb = 0.0_8
2693  tempb = x3*xxb
2694  axb = xxb
2695  t1b = tempb
2696  t2b = x3*tempb
2697  x3b = (t1+t2*x3)*xxb + t2*tempb
2698  cab = (1.0_8-fw)*t2b
2699  cbb = fw*t2b
2700  bab = (1.0_8-fw)*t1b
2701  bbb = fw*t1b
2702  tempb0 = 0.5*fw*axb
2703  fwb = (bb-ba)*t1b + (0.5*((pc+pa)*fw+pc-pa)-pb*2*fw)*axb + (pc+pa)*&
2704 & tempb0 + (cb-ca)*t2b
2705  pcb = (fw+1.0_8)*tempb0
2706  pab = (fw-1.0)*tempb0
2707  pbb = (1.-fw**2)*axb
2708  fpb = (coef3(ip+1, iw)-coef3(ip, iw))*cab + (coef2(ip+1, iw)-coef2(ip&
2709 & , iw))*bab + (coef1(ip+1, iw)-coef1(ip, iw))*pbb + (coef1(ip+1, iw-1&
2710 & )-coef1(ip, iw-1))*pab + (coef1(ip+1, iw+1)-coef1(ip, iw+1))*pcb + (&
2711 & coef2(ip+1, iw+1)-coef2(ip, iw+1))*bbb + (coef3(ip+1, iw+1)-coef3(ip&
2712 & , iw+1))*cbb
2713  peb = fpb
2714  web = fwb
2715  CALL popcontrol1b(branch)
2716  IF (branch .EQ. 0) peb = 0.0_8
2717  CALL popcontrol1b(branch)
2718  IF (branch .EQ. 0) web = 0.0_8
2719  x2b = dpe*peb/(x2*log(10.0))
2720  x1b = dwe*web/(x1*log(10.0))
2721  s3b = s3b + x1c*x3b
2722  x1cb = s2*x2b + s3*x3b
2723  s2b = s2b + x1c*x2b
2724  s1b = s1b + x1b - x1cb/s1**2
2725  dtb = dtb + dw*s3b
2726  dwb = dwb + p*s2b + s1b + dt*s3b
2727 END SUBROUTINE tablup_b
2728 
2729 ! Differentiation of h2okdis in reverse (adjoint) mode:
2730 ! gradient of useful results: tran tcon h2oexp th2o conexp
2731 ! with respect to varying inputs: tcon h2oexp th2o conexp
2732 !**********************************************************************
2733 SUBROUTINE h2okdis_b(ib, np, k, fkw, gkw, ne, h2oexp, h2oexpb, conexp, &
2734 & conexpb, th2o, th2ob, tcon, tconb, tran, tranb)
2735  IMPLICIT NONE
2736 !---- input parameters ------
2737  INTEGER :: ib, ne, np, k
2738  REAL*8 :: h2oexp(0:np, 6), conexp(0:np, 3)
2739  REAL*8 :: h2oexpb(0:np, 6), conexpb(0:np, 3)
2740  REAL*8 :: fkw(6, 9), gkw(6, 3)
2741 !---- updated parameters -----
2742  REAL*8 :: th2o(6), tcon(3), tran
2743  REAL*8 :: th2ob(6), tconb(3), tranb
2744 !---- temporary arrays -----
2745  REAL*8 :: trnth2o
2746  REAL*8 :: trnth2ob
2747  REAL*8 :: tempb2
2748  REAL*8 :: tempb1
2749  REAL*8 :: tempb0
2750  REAL*8 :: tempb
2751 !-----tco2 are the six exp factors between levels k1 and k2
2752 ! tran is the updated total transmittance between levels k1 and k2
2753 !-----th2o is the 6 exp factors between levels k1 and k2 due to
2754 ! h2o line absorption.
2755 !-----tcon is the 3 exp factors between levels k1 and k2 due to
2756 ! h2o continuum absorption.
2757 !-----trnth2o is the total transmittance between levels k1 and k2 due
2758 ! to both line and continuum absorption.
2759 !-----Compute th2o following Eq. (8.23).
2760  CALL pushreal8(th2o(1))
2761  th2o(1) = th2o(1)*h2oexp(k, 1)
2762  CALL pushreal8(th2o(2))
2763  th2o(2) = th2o(2)*h2oexp(k, 2)
2764  CALL pushreal8(th2o(3))
2765  th2o(3) = th2o(3)*h2oexp(k, 3)
2766  CALL pushreal8(th2o(4))
2767  th2o(4) = th2o(4)*h2oexp(k, 4)
2768  CALL pushreal8(th2o(5))
2769  th2o(5) = th2o(5)*h2oexp(k, 5)
2770  CALL pushreal8(th2o(6))
2771  th2o(6) = th2o(6)*h2oexp(k, 6)
2772  IF (ne .EQ. 0) THEN
2773  trnth2ob = tran*tranb
2774  th2ob(1) = th2ob(1) + fkw(1, ib)*trnth2ob
2775  th2ob(2) = th2ob(2) + fkw(2, ib)*trnth2ob
2776  th2ob(3) = th2ob(3) + fkw(3, ib)*trnth2ob
2777  th2ob(4) = th2ob(4) + fkw(4, ib)*trnth2ob
2778  th2ob(5) = th2ob(5) + fkw(5, ib)*trnth2ob
2779  th2ob(6) = th2ob(6) + fkw(6, ib)*trnth2ob
2780  ELSE IF (ne .EQ. 1) THEN
2781 !-----Compute trnh2o following Eqs. (8.25) and (4.27).
2782  CALL pushreal8(tcon(1))
2783  tcon(1) = tcon(1)*conexp(k, 1)
2784  trnth2ob = tran*tranb
2785  tempb = tcon(1)*trnth2ob
2786  th2ob(1) = th2ob(1) + fkw(1, ib)*tempb
2787  th2ob(2) = th2ob(2) + fkw(2, ib)*tempb
2788  th2ob(3) = th2ob(3) + fkw(3, ib)*tempb
2789  th2ob(4) = th2ob(4) + fkw(4, ib)*tempb
2790  th2ob(5) = th2ob(5) + fkw(5, ib)*tempb
2791  th2ob(6) = th2ob(6) + fkw(6, ib)*tempb
2792  tconb(1) = tconb(1) + (fkw(1, ib)*th2o(1)+fkw(2, ib)*th2o(2)+fkw(3, &
2793 & ib)*th2o(3)+fkw(4, ib)*th2o(4)+fkw(5, ib)*th2o(5)+fkw(6, ib)*th2o(&
2794 & 6))*trnth2ob
2795  CALL popreal8(tcon(1))
2796  conexpb(k, 1) = conexpb(k, 1) + tcon(1)*tconb(1)
2797  tconb(1) = conexp(k, 1)*tconb(1)
2798  ELSE
2799 !-----For band 3. This band is divided into 3 subbands.
2800  CALL pushreal8(tcon(1))
2801  tcon(1) = tcon(1)*conexp(k, 1)
2802  CALL pushreal8(tcon(2))
2803  tcon(2) = tcon(2)*conexp(k, 2)
2804  CALL pushreal8(tcon(3))
2805  tcon(3) = tcon(3)*conexp(k, 3)
2806 !-----Compute trnh2o following Eqs. (4.29) and (8.25).
2807  trnth2ob = tran*tranb
2808  tempb0 = tcon(1)*trnth2ob
2809  tempb1 = tcon(2)*trnth2ob
2810  tempb2 = tcon(3)*trnth2ob
2811  th2ob(1) = th2ob(1) + gkw(1, 3)*tempb2 + gkw(1, 2)*tempb1 + gkw(1, 1&
2812 & )*tempb0
2813  th2ob(2) = th2ob(2) + gkw(2, 3)*tempb2 + gkw(2, 2)*tempb1 + gkw(2, 1&
2814 & )*tempb0
2815  th2ob(3) = th2ob(3) + gkw(3, 3)*tempb2 + gkw(3, 2)*tempb1 + gkw(3, 1&
2816 & )*tempb0
2817  th2ob(4) = th2ob(4) + gkw(4, 3)*tempb2 + gkw(4, 2)*tempb1 + gkw(4, 1&
2818 & )*tempb0
2819  th2ob(5) = th2ob(5) + gkw(5, 3)*tempb2 + gkw(5, 2)*tempb1 + gkw(5, 1&
2820 & )*tempb0
2821  th2ob(6) = th2ob(6) + gkw(6, 3)*tempb2 + gkw(6, 2)*tempb1 + gkw(6, 1&
2822 & )*tempb0
2823  tconb(1) = tconb(1) + (gkw(1, 1)*th2o(1)+gkw(2, 1)*th2o(2)+gkw(3, 1)&
2824 & *th2o(3)+gkw(4, 1)*th2o(4)+gkw(5, 1)*th2o(5)+gkw(6, 1)*th2o(6))*&
2825 & trnth2ob
2826  tconb(2) = tconb(2) + (gkw(1, 2)*th2o(1)+gkw(2, 2)*th2o(2)+gkw(3, 2)&
2827 & *th2o(3)+gkw(4, 2)*th2o(4)+gkw(5, 2)*th2o(5)+gkw(6, 2)*th2o(6))*&
2828 & trnth2ob
2829  tconb(3) = tconb(3) + (gkw(1, 3)*th2o(1)+gkw(2, 3)*th2o(2)+gkw(3, 3)&
2830 & *th2o(3)+gkw(4, 3)*th2o(4)+gkw(5, 3)*th2o(5)+gkw(6, 3)*th2o(6))*&
2831 & trnth2ob
2832  CALL popreal8(tcon(3))
2833  conexpb(k, 3) = conexpb(k, 3) + tcon(3)*tconb(3)
2834  tconb(3) = conexp(k, 3)*tconb(3)
2835  CALL popreal8(tcon(2))
2836  conexpb(k, 2) = conexpb(k, 2) + tcon(2)*tconb(2)
2837  tconb(2) = conexp(k, 2)*tconb(2)
2838  CALL popreal8(tcon(1))
2839  conexpb(k, 1) = conexpb(k, 1) + tcon(1)*tconb(1)
2840  tconb(1) = conexp(k, 1)*tconb(1)
2841  END IF
2842  CALL popreal8(th2o(6))
2843  h2oexpb(k, 6) = h2oexpb(k, 6) + th2o(6)*th2ob(6)
2844  th2ob(6) = h2oexp(k, 6)*th2ob(6)
2845  CALL popreal8(th2o(5))
2846  h2oexpb(k, 5) = h2oexpb(k, 5) + th2o(5)*th2ob(5)
2847  th2ob(5) = h2oexp(k, 5)*th2ob(5)
2848  CALL popreal8(th2o(4))
2849  h2oexpb(k, 4) = h2oexpb(k, 4) + th2o(4)*th2ob(4)
2850  th2ob(4) = h2oexp(k, 4)*th2ob(4)
2851  CALL popreal8(th2o(3))
2852  h2oexpb(k, 3) = h2oexpb(k, 3) + th2o(3)*th2ob(3)
2853  th2ob(3) = h2oexp(k, 3)*th2ob(3)
2854  CALL popreal8(th2o(2))
2855  h2oexpb(k, 2) = h2oexpb(k, 2) + th2o(2)*th2ob(2)
2856  th2ob(2) = h2oexp(k, 2)*th2ob(2)
2857  CALL popreal8(th2o(1))
2858  h2oexpb(k, 1) = h2oexpb(k, 1) + th2o(1)*th2ob(1)
2859  th2ob(1) = h2oexp(k, 1)*th2ob(1)
2860 END SUBROUTINE h2okdis_b
2861 
2862 ! Differentiation of n2okdis in reverse (adjoint) mode:
2863 ! gradient of useful results: tran tn2o n2oexp
2864 ! with respect to varying inputs: tran tn2o n2oexp
2865 !**********************************************************************
2866 SUBROUTINE n2okdis_b(ib, np, k, n2oexp, n2oexpb, tn2o, tn2ob, tran, &
2867 & tranb)
2868  IMPLICIT NONE
2869  INTEGER :: ib, np, k
2870 !---- input parameters -----
2871  REAL*8 :: n2oexp(0:np, 4)
2872  REAL*8 :: n2oexpb(0:np, 4)
2873 !---- updated parameters -----
2874  REAL*8 :: tn2o(4), tran
2875  REAL*8 :: tn2ob(4), tranb
2876 !---- temporary arrays -----
2877  REAL*8 :: xc
2878  REAL*8 :: xcb
2879  INTEGER :: branch
2880 !-----tn2o is computed from Eq. (8.23).
2881 ! xc is the total n2o transmittance computed from (8.25)
2882 ! The k-distribution functions are given in Table 5.
2883 !-----band 6
2884  IF (ib .EQ. 6) THEN
2885  CALL pushreal8(tn2o(1))
2886  tn2o(1) = tn2o(1)*n2oexp(k, 1)
2887  xc = 0.940414*tn2o(1)
2888  CALL pushreal8(tn2o(2))
2889  tn2o(2) = tn2o(2)*n2oexp(k, 2)
2890  xc = xc + 0.059586*tn2o(2)
2891 !-----band 7
2892  CALL pushcontrol1b(0)
2893  ELSE
2894  CALL pushreal8(tn2o(1))
2895  tn2o(1) = tn2o(1)*n2oexp(k, 1)
2896  xc = 0.561961*tn2o(1)
2897  CALL pushreal8(tn2o(2))
2898  tn2o(2) = tn2o(2)*n2oexp(k, 2)
2899  xc = xc + 0.138707*tn2o(2)
2900  CALL pushreal8(tn2o(3))
2901  tn2o(3) = tn2o(3)*n2oexp(k, 3)
2902  xc = xc + 0.240670*tn2o(3)
2903  CALL pushreal8(tn2o(4))
2904  tn2o(4) = tn2o(4)*n2oexp(k, 4)
2905  xc = xc + 0.058662*tn2o(4)
2906  CALL pushcontrol1b(1)
2907  END IF
2908  xcb = tran*tranb
2909  tranb = xc*tranb
2910  CALL popcontrol1b(branch)
2911  IF (branch .EQ. 0) THEN
2912  tn2ob(2) = tn2ob(2) + 0.059586*xcb
2913  CALL popreal8(tn2o(2))
2914  n2oexpb(k, 2) = n2oexpb(k, 2) + tn2o(2)*tn2ob(2)
2915  tn2ob(2) = n2oexp(k, 2)*tn2ob(2)
2916  tn2ob(1) = tn2ob(1) + 0.940414*xcb
2917  CALL popreal8(tn2o(1))
2918  n2oexpb(k, 1) = n2oexpb(k, 1) + tn2o(1)*tn2ob(1)
2919  tn2ob(1) = n2oexp(k, 1)*tn2ob(1)
2920  ELSE
2921  tn2ob(4) = tn2ob(4) + 0.058662*xcb
2922  CALL popreal8(tn2o(4))
2923  n2oexpb(k, 4) = n2oexpb(k, 4) + tn2o(4)*tn2ob(4)
2924  tn2ob(4) = n2oexp(k, 4)*tn2ob(4)
2925  tn2ob(3) = tn2ob(3) + 0.240670*xcb
2926  CALL popreal8(tn2o(3))
2927  n2oexpb(k, 3) = n2oexpb(k, 3) + tn2o(3)*tn2ob(3)
2928  tn2ob(3) = n2oexp(k, 3)*tn2ob(3)
2929  tn2ob(2) = tn2ob(2) + 0.138707*xcb
2930  CALL popreal8(tn2o(2))
2931  n2oexpb(k, 2) = n2oexpb(k, 2) + tn2o(2)*tn2ob(2)
2932  tn2ob(2) = n2oexp(k, 2)*tn2ob(2)
2933  tn2ob(1) = tn2ob(1) + 0.561961*xcb
2934  CALL popreal8(tn2o(1))
2935  n2oexpb(k, 1) = n2oexpb(k, 1) + tn2o(1)*tn2ob(1)
2936  tn2ob(1) = n2oexp(k, 1)*tn2ob(1)
2937  END IF
2938 END SUBROUTINE n2okdis_b
2939 
2940 ! Differentiation of ch4kdis in reverse (adjoint) mode:
2941 ! gradient of useful results: tran ch4exp tch4
2942 ! with respect to varying inputs: tran ch4exp tch4
2943 !**********************************************************************
2944 SUBROUTINE ch4kdis_b(ib, np, k, ch4exp, ch4expb, tch4, tch4b, tran, &
2945 & tranb)
2946  IMPLICIT NONE
2947  INTEGER :: ib, np, k
2948 !---- input parameters -----
2949  REAL*8 :: ch4exp(0:np, 4)
2950  REAL*8 :: ch4expb(0:np, 4)
2951 !---- updated parameters -----
2952  REAL*8 :: tch4(4), tran
2953  REAL*8 :: tch4b(4), tranb
2954 !---- temporary arrays -----
2955  REAL*8 :: xc
2956  REAL*8 :: xcb
2957  INTEGER :: branch
2958 !-----tch4 is computed from Eq. (8.23).
2959 ! xc is the total ch4 transmittance computed from (8.25)
2960 ! The k-distribution functions are given in Table 5.
2961 !-----band 6
2962  IF (ib .EQ. 6) THEN
2963  CALL pushreal8(tch4(1))
2964  tch4(1) = tch4(1)*ch4exp(k, 1)
2965  xc = tch4(1)
2966 !-----band 7
2967  CALL pushcontrol1b(0)
2968  ELSE
2969  CALL pushreal8(tch4(1))
2970  tch4(1) = tch4(1)*ch4exp(k, 1)
2971  xc = 0.610650*tch4(1)
2972  CALL pushreal8(tch4(2))
2973  tch4(2) = tch4(2)*ch4exp(k, 2)
2974  xc = xc + 0.280212*tch4(2)
2975  CALL pushreal8(tch4(3))
2976  tch4(3) = tch4(3)*ch4exp(k, 3)
2977  xc = xc + 0.107349*tch4(3)
2978  CALL pushreal8(tch4(4))
2979  tch4(4) = tch4(4)*ch4exp(k, 4)
2980  xc = xc + 0.001789*tch4(4)
2981  CALL pushcontrol1b(1)
2982  END IF
2983  xcb = tran*tranb
2984  tranb = xc*tranb
2985  CALL popcontrol1b(branch)
2986  IF (branch .EQ. 0) THEN
2987  tch4b(1) = tch4b(1) + xcb
2988  CALL popreal8(tch4(1))
2989  ch4expb(k, 1) = ch4expb(k, 1) + tch4(1)*tch4b(1)
2990  tch4b(1) = ch4exp(k, 1)*tch4b(1)
2991  ELSE
2992  tch4b(4) = tch4b(4) + 0.001789*xcb
2993  CALL popreal8(tch4(4))
2994  ch4expb(k, 4) = ch4expb(k, 4) + tch4(4)*tch4b(4)
2995  tch4b(4) = ch4exp(k, 4)*tch4b(4)
2996  tch4b(3) = tch4b(3) + 0.107349*xcb
2997  CALL popreal8(tch4(3))
2998  ch4expb(k, 3) = ch4expb(k, 3) + tch4(3)*tch4b(3)
2999  tch4b(3) = ch4exp(k, 3)*tch4b(3)
3000  tch4b(2) = tch4b(2) + 0.280212*xcb
3001  CALL popreal8(tch4(2))
3002  ch4expb(k, 2) = ch4expb(k, 2) + tch4(2)*tch4b(2)
3003  tch4b(2) = ch4exp(k, 2)*tch4b(2)
3004  tch4b(1) = tch4b(1) + 0.610650*xcb
3005  CALL popreal8(tch4(1))
3006  ch4expb(k, 1) = ch4expb(k, 1) + tch4(1)*tch4b(1)
3007  tch4b(1) = ch4exp(k, 1)*tch4b(1)
3008  END IF
3009 END SUBROUTINE ch4kdis_b
3010 
3011 ! Differentiation of comkdis in reverse (adjoint) mode:
3012 ! gradient of useful results: tran tcom comexp
3013 ! with respect to varying inputs: tran tcom comexp
3014 !**********************************************************************
3015 SUBROUTINE comkdis_b(ib, np, k, comexp, comexpb, tcom, tcomb, tran, &
3016 & tranb)
3017  IMPLICIT NONE
3018  INTEGER :: ib, np, k
3019 !---- input parameters -----
3020  REAL*8 :: comexp(0:np, 6)
3021  REAL*8 :: comexpb(0:np, 6)
3022 !---- updated parameters -----
3023  REAL*8 :: tcom(6), tran
3024  REAL*8 :: tcomb(6), tranb
3025 !---- temporary arrays -----
3026  REAL*8 :: xc
3027  REAL*8 :: xcb
3028  INTEGER :: branch
3029 !-----tcom is computed from Eq. (8.23).
3030 ! xc is the total co2 transmittance computed from (8.25)
3031 ! The k-distribution functions are given in Table 6.
3032 !-----band 4
3033  IF (ib .EQ. 4) THEN
3034  CALL pushreal8(tcom(1))
3035  tcom(1) = tcom(1)*comexp(k, 1)
3036  xc = 0.12159*tcom(1)
3037  CALL pushreal8(tcom(2))
3038  tcom(2) = tcom(2)*comexp(k, 2)
3039  xc = xc + 0.24359*tcom(2)
3040  CALL pushreal8(tcom(3))
3041  tcom(3) = tcom(3)*comexp(k, 3)
3042  xc = xc + 0.24981*tcom(3)
3043  CALL pushreal8(tcom(4))
3044  tcom(4) = tcom(4)*comexp(k, 4)
3045  xc = xc + 0.26427*tcom(4)
3046  CALL pushreal8(tcom(5))
3047  tcom(5) = tcom(5)*comexp(k, 5)
3048  xc = xc + 0.07807*tcom(5)
3049  CALL pushreal8(tcom(6))
3050  tcom(6) = tcom(6)*comexp(k, 6)
3051  xc = xc + 0.04267*tcom(6)
3052 !-----band 5
3053  CALL pushcontrol1b(0)
3054  ELSE
3055  CALL pushreal8(tcom(1))
3056  tcom(1) = tcom(1)*comexp(k, 1)
3057  xc = 0.06869*tcom(1)
3058  CALL pushreal8(tcom(2))
3059  tcom(2) = tcom(2)*comexp(k, 2)
3060  xc = xc + 0.14795*tcom(2)
3061  CALL pushreal8(tcom(3))
3062  tcom(3) = tcom(3)*comexp(k, 3)
3063  xc = xc + 0.19512*tcom(3)
3064  CALL pushreal8(tcom(4))
3065  tcom(4) = tcom(4)*comexp(k, 4)
3066  xc = xc + 0.33446*tcom(4)
3067  CALL pushreal8(tcom(5))
3068  tcom(5) = tcom(5)*comexp(k, 5)
3069  xc = xc + 0.17199*tcom(5)
3070  CALL pushreal8(tcom(6))
3071  tcom(6) = tcom(6)*comexp(k, 6)
3072  xc = xc + 0.08179*tcom(6)
3073  CALL pushcontrol1b(1)
3074  END IF
3075  xcb = tran*tranb
3076  tranb = xc*tranb
3077  CALL popcontrol1b(branch)
3078  IF (branch .EQ. 0) THEN
3079  tcomb(6) = tcomb(6) + 0.04267*xcb
3080  CALL popreal8(tcom(6))
3081  comexpb(k, 6) = comexpb(k, 6) + tcom(6)*tcomb(6)
3082  tcomb(6) = comexp(k, 6)*tcomb(6)
3083  tcomb(5) = tcomb(5) + 0.07807*xcb
3084  CALL popreal8(tcom(5))
3085  comexpb(k, 5) = comexpb(k, 5) + tcom(5)*tcomb(5)
3086  tcomb(5) = comexp(k, 5)*tcomb(5)
3087  tcomb(4) = tcomb(4) + 0.26427*xcb
3088  CALL popreal8(tcom(4))
3089  comexpb(k, 4) = comexpb(k, 4) + tcom(4)*tcomb(4)
3090  tcomb(4) = comexp(k, 4)*tcomb(4)
3091  tcomb(3) = tcomb(3) + 0.24981*xcb
3092  CALL popreal8(tcom(3))
3093  comexpb(k, 3) = comexpb(k, 3) + tcom(3)*tcomb(3)
3094  tcomb(3) = comexp(k, 3)*tcomb(3)
3095  tcomb(2) = tcomb(2) + 0.24359*xcb
3096  CALL popreal8(tcom(2))
3097  comexpb(k, 2) = comexpb(k, 2) + tcom(2)*tcomb(2)
3098  tcomb(2) = comexp(k, 2)*tcomb(2)
3099  tcomb(1) = tcomb(1) + 0.12159*xcb
3100  CALL popreal8(tcom(1))
3101  comexpb(k, 1) = comexpb(k, 1) + tcom(1)*tcomb(1)
3102  tcomb(1) = comexp(k, 1)*tcomb(1)
3103  ELSE
3104  tcomb(6) = tcomb(6) + 0.08179*xcb
3105  CALL popreal8(tcom(6))
3106  comexpb(k, 6) = comexpb(k, 6) + tcom(6)*tcomb(6)
3107  tcomb(6) = comexp(k, 6)*tcomb(6)
3108  tcomb(5) = tcomb(5) + 0.17199*xcb
3109  CALL popreal8(tcom(5))
3110  comexpb(k, 5) = comexpb(k, 5) + tcom(5)*tcomb(5)
3111  tcomb(5) = comexp(k, 5)*tcomb(5)
3112  tcomb(4) = tcomb(4) + 0.33446*xcb
3113  CALL popreal8(tcom(4))
3114  comexpb(k, 4) = comexpb(k, 4) + tcom(4)*tcomb(4)
3115  tcomb(4) = comexp(k, 4)*tcomb(4)
3116  tcomb(3) = tcomb(3) + 0.19512*xcb
3117  CALL popreal8(tcom(3))
3118  comexpb(k, 3) = comexpb(k, 3) + tcom(3)*tcomb(3)
3119  tcomb(3) = comexp(k, 3)*tcomb(3)
3120  tcomb(2) = tcomb(2) + 0.14795*xcb
3121  CALL popreal8(tcom(2))
3122  comexpb(k, 2) = comexpb(k, 2) + tcom(2)*tcomb(2)
3123  tcomb(2) = comexp(k, 2)*tcomb(2)
3124  tcomb(1) = tcomb(1) + 0.06869*xcb
3125  CALL popreal8(tcom(1))
3126  comexpb(k, 1) = comexpb(k, 1) + tcom(1)*tcomb(1)
3127  tcomb(1) = comexp(k, 1)*tcomb(1)
3128  END IF
3129 END SUBROUTINE comkdis_b
3130 
3131 ! Differentiation of cfckdis in reverse (adjoint) mode:
3132 ! gradient of useful results: tran tcfc cfcexp
3133 ! with respect to varying inputs: tran tcfc cfcexp
3134 !**********************************************************************
3135 SUBROUTINE cfckdis_b(np, k, cfcexp, cfcexpb, tcfc, tcfcb, tran, tranb)
3136  IMPLICIT NONE
3137 !---- input parameters -----
3138  INTEGER :: k, np
3139  REAL*8 :: cfcexp(0:np)
3140  REAL*8 :: cfcexpb(0:np)
3141 !---- updated parameters -----
3142  REAL*8 :: tcfc, tran
3143  REAL*8 :: tcfcb, tranb
3144 !-----tcfc is the exp factors between levels k1 and k2.
3145  CALL pushreal8(tcfc)
3146  tcfc = tcfc*cfcexp(k)
3147  tcfcb = tcfcb + tran*tranb
3148  tranb = tcfc*tranb
3149  CALL popreal8(tcfc)
3150  cfcexpb(k) = cfcexpb(k) + tcfc*tcfcb
3151  tcfcb = cfcexp(k)*tcfcb
3152 END SUBROUTINE cfckdis_b
3153 
3154 ! Differentiation of b10kdis in reverse (adjoint) mode:
3155 ! gradient of useful results: tran tn2o co2exp tcon h2oexp
3156 ! n2oexp th2o tco2 conexp
3157 ! with respect to varying inputs: tn2o co2exp tcon h2oexp n2oexp
3158 ! th2o tco2 conexp
3159 !**********************************************************************
3160 SUBROUTINE b10kdis_b(np, k, h2oexp, h2oexpb, conexp, conexpb, co2exp, &
3161 & co2expb, n2oexp, n2oexpb, th2o, th2ob, tcon, tconb, tco2, tco2b, tn2o&
3162 & , tn2ob, tran, tranb)
3163  IMPLICIT NONE
3164  INTEGER :: np, k
3165 !---- input parameters -----
3166  REAL*8 :: h2oexp(0:np, 5), conexp(0:np), co2exp(0:np, 6), n2oexp(0:np&
3167 & , 2)
3168  REAL*8 :: h2oexpb(0:np, 5), conexpb(0:np), co2expb(0:np, 6), n2oexpb(0&
3169 & :np, 2)
3170 !---- updated parameters -----
3171  REAL*8 :: th2o(6), tcon(3), tco2(6), tn2o(4), tran
3172  REAL*8 :: th2ob(6), tconb(3), tco2b(6), tn2ob(4), tranb
3173 !---- temporary arrays -----
3174  REAL*8 :: xx
3175  REAL*8 :: xxb
3176 !-----For h2o line. The k-distribution functions are given in Table 4.
3177  CALL pushreal8(th2o(1))
3178  th2o(1) = th2o(1)*h2oexp(k, 1)
3179  xx = 0.3153*th2o(1)
3180  CALL pushreal8(th2o(2))
3181  th2o(2) = th2o(2)*h2oexp(k, 2)
3182  xx = xx + 0.4604*th2o(2)
3183  CALL pushreal8(th2o(3))
3184  th2o(3) = th2o(3)*h2oexp(k, 3)
3185  xx = xx + 0.1326*th2o(3)
3186  CALL pushreal8(th2o(4))
3187  th2o(4) = th2o(4)*h2oexp(k, 4)
3188  xx = xx + 0.0798*th2o(4)
3189  CALL pushreal8(th2o(5))
3190  th2o(5) = th2o(5)*h2oexp(k, 5)
3191  xx = xx + 0.0119*th2o(5)
3192  tran = xx
3193 !-----For h2o continuum. Note that conexp(k,3) is for subband 3a.
3194  CALL pushreal8(tcon(1))
3195  tcon(1) = tcon(1)*conexp(k)
3196  CALL pushreal8(tran)
3197  tran = tran*tcon(1)
3198 !-----For co2 (Table 6)
3199  CALL pushreal8(tco2(1))
3200  tco2(1) = tco2(1)*co2exp(k, 1)
3201  xx = 0.2673*tco2(1)
3202  CALL pushreal8(tco2(2))
3203  tco2(2) = tco2(2)*co2exp(k, 2)
3204  xx = xx + 0.2201*tco2(2)
3205  CALL pushreal8(tco2(3))
3206  tco2(3) = tco2(3)*co2exp(k, 3)
3207  xx = xx + 0.2106*tco2(3)
3208  CALL pushreal8(tco2(4))
3209  tco2(4) = tco2(4)*co2exp(k, 4)
3210  xx = xx + 0.2409*tco2(4)
3211  CALL pushreal8(tco2(5))
3212  tco2(5) = tco2(5)*co2exp(k, 5)
3213  xx = xx + 0.0196*tco2(5)
3214  CALL pushreal8(tco2(6))
3215  tco2(6) = tco2(6)*co2exp(k, 6)
3216  xx = xx + 0.0415*tco2(6)
3217  CALL pushreal8(tran)
3218  tran = tran*xx
3219 !-----For n2o (Table 5)
3220  CALL pushreal8(tn2o(1))
3221  tn2o(1) = tn2o(1)*n2oexp(k, 1)
3222  CALL pushreal8(xx)
3223  xx = 0.970831*tn2o(1)
3224  CALL pushreal8(tn2o(2))
3225  tn2o(2) = tn2o(2)*n2oexp(k, 2)
3226  xx = xx + 0.029169*tn2o(2)
3227  xxb = tran*tranb
3228  tranb = (xx-1.0)*tranb
3229  tn2ob(2) = tn2ob(2) + 0.029169*xxb
3230  CALL popreal8(tn2o(2))
3231  n2oexpb(k, 2) = n2oexpb(k, 2) + tn2o(2)*tn2ob(2)
3232  tn2ob(2) = n2oexp(k, 2)*tn2ob(2)
3233  CALL popreal8(xx)
3234  tn2ob(1) = tn2ob(1) + 0.970831*xxb
3235  CALL popreal8(tn2o(1))
3236  n2oexpb(k, 1) = n2oexpb(k, 1) + tn2o(1)*tn2ob(1)
3237  tn2ob(1) = n2oexp(k, 1)*tn2ob(1)
3238  CALL popreal8(tran)
3239  xxb = tran*tranb
3240  tranb = xx*tranb
3241  tco2b(6) = tco2b(6) + 0.0415*xxb
3242  CALL popreal8(tco2(6))
3243  co2expb(k, 6) = co2expb(k, 6) + tco2(6)*tco2b(6)
3244  tco2b(6) = co2exp(k, 6)*tco2b(6)
3245  tco2b(5) = tco2b(5) + 0.0196*xxb
3246  CALL popreal8(tco2(5))
3247  co2expb(k, 5) = co2expb(k, 5) + tco2(5)*tco2b(5)
3248  tco2b(5) = co2exp(k, 5)*tco2b(5)
3249  tco2b(4) = tco2b(4) + 0.2409*xxb
3250  CALL popreal8(tco2(4))
3251  co2expb(k, 4) = co2expb(k, 4) + tco2(4)*tco2b(4)
3252  tco2b(4) = co2exp(k, 4)*tco2b(4)
3253  tco2b(3) = tco2b(3) + 0.2106*xxb
3254  CALL popreal8(tco2(3))
3255  co2expb(k, 3) = co2expb(k, 3) + tco2(3)*tco2b(3)
3256  tco2b(3) = co2exp(k, 3)*tco2b(3)
3257  tco2b(2) = tco2b(2) + 0.2201*xxb
3258  CALL popreal8(tco2(2))
3259  co2expb(k, 2) = co2expb(k, 2) + tco2(2)*tco2b(2)
3260  tco2b(2) = co2exp(k, 2)*tco2b(2)
3261  tco2b(1) = tco2b(1) + 0.2673*xxb
3262  CALL popreal8(tco2(1))
3263  co2expb(k, 1) = co2expb(k, 1) + tco2(1)*tco2b(1)
3264  tco2b(1) = co2exp(k, 1)*tco2b(1)
3265  CALL popreal8(tran)
3266  tconb(1) = tconb(1) + tran*tranb
3267  tranb = tcon(1)*tranb
3268  CALL popreal8(tcon(1))
3269  conexpb(k) = conexpb(k) + tcon(1)*tconb(1)
3270  tconb(1) = conexp(k)*tconb(1)
3271  xxb = tranb
3272  th2ob(5) = th2ob(5) + 0.0119*xxb
3273  CALL popreal8(th2o(5))
3274  h2oexpb(k, 5) = h2oexpb(k, 5) + th2o(5)*th2ob(5)
3275  th2ob(5) = h2oexp(k, 5)*th2ob(5)
3276  th2ob(4) = th2ob(4) + 0.0798*xxb
3277  CALL popreal8(th2o(4))
3278  h2oexpb(k, 4) = h2oexpb(k, 4) + th2o(4)*th2ob(4)
3279  th2ob(4) = h2oexp(k, 4)*th2ob(4)
3280  th2ob(3) = th2ob(3) + 0.1326*xxb
3281  CALL popreal8(th2o(3))
3282  h2oexpb(k, 3) = h2oexpb(k, 3) + th2o(3)*th2ob(3)
3283  th2ob(3) = h2oexp(k, 3)*th2ob(3)
3284  th2ob(2) = th2ob(2) + 0.4604*xxb
3285  CALL popreal8(th2o(2))
3286  h2oexpb(k, 2) = h2oexpb(k, 2) + th2o(2)*th2ob(2)
3287  th2ob(2) = h2oexp(k, 2)*th2ob(2)
3288  th2ob(1) = th2ob(1) + 0.3153*xxb
3289  CALL popreal8(th2o(1))
3290  h2oexpb(k, 1) = h2oexpb(k, 1) + th2o(1)*th2ob(1)
3291  th2ob(1) = h2oexp(k, 1)*th2ob(1)
3292 END SUBROUTINE b10kdis_b
3293 
3294 ! Differentiation of cldovlp in reverse (adjoint) mode:
3295 ! gradient of useful results: cldhi ett cldlw enn cldmd
3296 ! with respect to varying inputs: cldhi ett cldlw enn cldmd
3297 !mjs
3298 !***********************************************************************
3299 SUBROUTINE cldovlp_b(np, k1, k2, ict, icb, icx, ncld, enn, ennb, ett, &
3300 & ettb, cldhi, cldhib, cldmd, cldmdb, cldlw, cldlwb)
3301  IMPLICIT NONE
3302  INTEGER, INTENT(IN) :: np, k1, k2, ict, icb, icx(0:np), ncld(3)
3303  REAL*8, INTENT(IN) :: enn(0:np), ett(0:np)
3304  REAL*8 :: ennb(0:np), ettb(0:np)
3305  REAL*8, INTENT(INOUT) :: cldhi, cldmd, cldlw
3306  REAL*8 :: cldhib, cldmdb, cldlwb
3307  INTEGER :: j, k, km, kx
3308  INTEGER :: branch
3309  km = k2 - 1
3310  IF (km .LT. ict) THEN
3311 ! do high clouds
3312  kx = ncld(1)
3313  IF (kx .EQ. 1 .OR. cldhi .EQ. 0.) THEN
3314  ennb(km) = ennb(km) + cldhib
3315  ELSE
3316 !if ( (kx.lt.1 .or. kx.gt.1) .and. abs(cldhi) .gt. 0.) then
3317  CALL pushreal8(cldhi)
3318  cldhi = 0.0
3319  IF (kx .NE. 0) THEN
3320  DO k=ict-kx,ict-1
3321  j = icx(k)
3322  IF (j .GE. k1 .AND. j .LE. km) THEN
3323  CALL pushreal8(cldhi)
3324  cldhi = enn(j) + ett(j)*cldhi
3325  CALL pushcontrol1b(1)
3326  ELSE
3327  CALL pushcontrol1b(0)
3328  END IF
3329  END DO
3330  DO k=ict-1,ict-kx,-1
3331  CALL popcontrol1b(branch)
3332  IF (branch .NE. 0) THEN
3333  j = icx(k)
3334  CALL popreal8(cldhi)
3335  ennb(j) = ennb(j) + cldhib
3336  ettb(j) = ettb(j) + cldhi*cldhib
3337  cldhib = ett(j)*cldhib
3338  END IF
3339  END DO
3340  END IF
3341  CALL popreal8(cldhi)
3342  END IF
3343  cldhib = 0.0_8
3344  ELSE IF (km .GE. ict .AND. km .LT. icb) THEN
3345 ! do middle clouds
3346  kx = ncld(2)
3347  IF (kx .EQ. 1 .OR. cldmd .EQ. 0.) THEN
3348  ennb(km) = ennb(km) + cldmdb
3349  ELSE
3350 !if ( (kx.lt.1 .or. kx.gt.1) .and. abs(cldhi) .gt. 0.) then
3351  CALL pushreal8(cldmd)
3352  cldmd = 0.0
3353  IF (kx .NE. 0) THEN
3354  DO k=icb-kx,icb-1
3355  j = icx(k)
3356  IF (j .GE. k1 .AND. j .LE. km) THEN
3357  CALL pushreal8(cldmd)
3358  cldmd = enn(j) + ett(j)*cldmd
3359  CALL pushcontrol1b(1)
3360  ELSE
3361  CALL pushcontrol1b(0)
3362  END IF
3363  END DO
3364  DO k=icb-1,icb-kx,-1
3365  CALL popcontrol1b(branch)
3366  IF (branch .NE. 0) THEN
3367  j = icx(k)
3368  CALL popreal8(cldmd)
3369  ennb(j) = ennb(j) + cldmdb
3370  ettb(j) = ettb(j) + cldmd*cldmdb
3371  cldmdb = ett(j)*cldmdb
3372  END IF
3373  END DO
3374  END IF
3375  CALL popreal8(cldmd)
3376  END IF
3377  cldmdb = 0.0_8
3378  ELSE
3379 ! do low clouds
3380  kx = ncld(3)
3381  IF (kx .EQ. 1 .OR. cldlw .EQ. 0.) THEN
3382  ennb(km) = ennb(km) + cldlwb
3383  ELSE
3384 !if ( (kx.lt.1 .or. kx.gt.1) .and. abs(cldhi) .gt. 0.) then
3385  CALL pushreal8(cldlw)
3386  cldlw = 0.0
3387  IF (kx .NE. 0) THEN
3388  DO k=np+1-kx,np
3389  j = icx(k)
3390  IF (j .GE. k1 .AND. j .LE. km) THEN
3391  CALL pushreal8(cldlw)
3392  cldlw = enn(j) + ett(j)*cldlw
3393  CALL pushcontrol1b(1)
3394  ELSE
3395  CALL pushcontrol1b(0)
3396  END IF
3397  END DO
3398  DO k=np,np+1-kx,-1
3399  CALL popcontrol1b(branch)
3400  IF (branch .NE. 0) THEN
3401  j = icx(k)
3402  CALL popreal8(cldlw)
3403  ennb(j) = ennb(j) + cldlwb
3404  ettb(j) = ettb(j) + cldlw*cldlwb
3405  cldlwb = ett(j)*cldlwb
3406  END IF
3407  END DO
3408  END IF
3409  CALL popreal8(cldlw)
3410  END IF
3411  cldlwb = 0.0_8
3412  END IF
3413 END SUBROUTINE cldovlp_b
3414 
3415 ! Differentiation of getirtau1 in reverse (adjoint) mode:
3416 ! gradient of useful results: hydromets tcldlyr fcld enn
3417 ! reff
3418 ! with respect to varying inputs: hydromets tcldlyr fcld enn
3419 ! reff
3420 SUBROUTINE getirtau1_b(ib, nlevs, dp, fcld, fcldb, reff, reffb, &
3421 & hydromets, hydrometsb, tcldlyr, tcldlyrb, enn, ennb, aib_ir1, awb_ir1&
3422 & , aiw_ir1, aww_ir1, aig_ir1, awg_ir1, cons_grav)
3423  IMPLICIT NONE
3424 ! !INPUT PARAMETERS:
3425 ! Band number
3426  INTEGER, INTENT(IN) :: ib
3427 ! Number of levels
3428  INTEGER, INTENT(IN) :: nlevs
3429 ! Delta pressure in Pa
3430  REAL*8, INTENT(IN) :: dp(nlevs)
3431 ! Cloud fraction (used sometimes)
3432  REAL*8, INTENT(IN) :: fcld(nlevs)
3433  REAL*8 :: fcldb(nlevs)
3434 ! Effective radius (microns)
3435  REAL*8, INTENT(IN) :: reff(nlevs, 4)
3436  REAL*8 :: reffb(nlevs, 4)
3437 ! Hydrometeors (kg/kg)
3438  REAL*8, INTENT(IN) :: hydromets(nlevs, 4)
3439  REAL*8 :: hydrometsb(nlevs, 4)
3440  REAL*8, INTENT(IN) :: aib_ir1(3, 10), awb_ir1(4, 10), aiw_ir1(4, 10)
3441  REAL*8, INTENT(IN) :: aww_ir1(4, 10), aig_ir1(4, 10), awg_ir1(4, 10)
3442  REAL*8, INTENT(IN) :: cons_grav
3443 ! !OUTPUT PARAMETERS:
3444 ! Flux transmissivity?
3445  REAL*8 :: tcldlyr(0:nlevs)
3446  REAL*8 :: tcldlyrb(0:nlevs)
3447 ! Flux transmissivity of a cloud layer?
3448  REAL*8 :: enn(0:nlevs)
3449  REAL*8 :: ennb(0:nlevs)
3450 ! !DESCRIPTION:
3451 ! Compute in-cloud or grid mean optical depths for infrared wavelengths
3452 ! Slots for reff, hydrometeors and tauall are as follows:
3453 ! 1 Cloud Ice
3454 ! 2 Cloud Liquid
3455 ! 3 Falling Liquid (Rain)
3456 ! 4 Falling Ice (Snow)
3457 !
3458 ! In the below calculations, the constants used in the tau calculation are in
3459 ! m$^2$ g$^{-1}$ and m$^2$ g$^{-1}$ $\mu$m. Thus, we must convert the kg contained in the
3460 ! pressure (Pa = kg m$^{-1}$ s$^{-2}$) to grams.
3461 !
3462 ! !REVISION HISTORY:
3463 ! 2011.11.18 MAT moved to Radiation_Shared and revised arg list, units
3464 !
3465 !EOP
3466 !------------------------------------------------------------------------------
3467 !BOC
3468  INTEGER :: k
3469  REAL*8 :: taucld1, taucld2, taucld3, taucld4
3470  REAL*8 :: taucld1b, taucld2b, taucld3b, taucld4b
3471  REAL*8 :: g1, g2, g3, g4, gg
3472  REAL*8 :: g1b, g2b, g3b, g4b, ggb
3473  REAL*8 :: w1, w2, w3, w4, ww
3474  REAL*8 :: w1b, w2b, w3b, w4b, wwb
3475  REAL*8 :: ff, tauc
3476  REAL*8 :: ffb, taucb
3477  REAL*8 :: reff_snow
3478  REAL*8 :: reff_snowb
3479  INTRINSIC min
3480  INTRINSIC abs
3481  INTRINSIC max
3482  INTRINSIC exp
3483 !-----Compute cloud optical thickness following Eqs. (6.4a,b) and (6.7)
3484 ! Rain optical thickness is set to 0.00307 /(gm/m**2).
3485 ! It is for a specific drop size distribution provided by Q. Fu.
3486  INTEGER :: branch
3487  REAL*8 :: temp3
3488  REAL*8 :: temp2
3489  REAL*8 :: temp1
3490  REAL*8 :: temp0
3491  REAL*8 :: tempb9
3492  REAL*8 :: tempb8
3493  REAL*8 :: tempb7
3494  REAL*8 :: tempb6
3495  REAL*8 :: tempb5
3496  REAL*8 :: tempb4
3497  REAL*8 :: tempb3
3498  REAL*8 :: tempb2
3499  REAL*8 :: tempb1
3500  REAL*8 :: tempb0
3501  REAL*8 :: tempb12
3502  REAL*8 :: tempb11
3503  REAL*8 :: tempb10
3504  REAL*8 :: temp16
3505  REAL*8 :: temp15
3506  REAL*8 :: temp14
3507  REAL*8 :: temp13
3508  REAL*8 :: temp12
3509  REAL*8 :: temp11
3510  REAL*8 :: temp10
3511  REAL*8 :: max1b
3512  REAL*8 :: tempb
3513  REAL*8 :: abs0
3514  REAL*8 :: temp
3515  REAL*8 :: max1
3516  REAL*8 :: temp9
3517  REAL*8 :: temp8
3518  REAL*8 :: temp7
3519  REAL*8 :: temp6
3520  REAL*8 :: temp5
3521  REAL*8 :: temp4
3522  DO k=1,nlevs
3523  IF (reff(k, 1) .LE. 0.0) THEN
3524  CALL pushreal8(taucld1)
3525  taucld1 = 0.0
3526  CALL pushcontrol1b(1)
3527  ELSE
3528  CALL pushreal8(taucld1)
3529  taucld1 = dp(k)*1.0e3/cons_grav*hydromets(k, 1)*(aib_ir1(1, ib)+&
3530 & aib_ir1(2, ib)/reff(k, 1)**aib_ir1(3, ib))
3531  CALL pushcontrol1b(0)
3532  END IF
3533  CALL pushreal8(taucld2)
3534  taucld2 = dp(k)*1.0e3/cons_grav*hydromets(k, 2)*(awb_ir1(1, ib)+(&
3535 & awb_ir1(2, ib)+(awb_ir1(3, ib)+awb_ir1(4, ib)*reff(k, 2))*reff(k, &
3536 & 2))*reff(k, 2))
3537  taucld3 = 0.00307*(dp(k)*1.0e3/cons_grav*hydromets(k, 3))
3538  IF (reff(k, 4) .GT. 112.0) THEN
3539  CALL pushreal8(reff_snow)
3540  reff_snow = 112.0
3541  CALL pushcontrol1b(0)
3542  ELSE
3543  CALL pushreal8(reff_snow)
3544  reff_snow = reff(k, 4)
3545  CALL pushcontrol1b(1)
3546  END IF
3547  IF (reff_snow .LE. 0.0) THEN
3548  CALL pushreal8(taucld4)
3549  taucld4 = 0.0
3550  CALL pushcontrol1b(0)
3551  ELSE
3552  CALL pushreal8(taucld4)
3553  taucld4 = dp(k)*1.0e3/cons_grav*hydromets(k, 4)*(aib_ir1(1, ib)+&
3554 & aib_ir1(2, ib)/reff_snow**aib_ir1(3, ib))
3555  CALL pushcontrol1b(1)
3556  END IF
3557 !-----Compute cloud single-scattering albedo and asymmetry factor for
3558 ! a mixture of ice particles and liquid drops following
3559 ! Eqs. (6.5), (6.6), (6.15) and (6.16).
3560 ! Single-scattering albedo and asymmetry factor of rain are set
3561 ! to 0.54 and 0.95, respectively, based on the information provided
3562 ! by Prof. Qiang Fu.
3563  CALL pushreal8(tauc)
3564  tauc = taucld1 + taucld2 + taucld3 + taucld4
3565  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
3566  CALL pushreal8(w1)
3567  w1 = taucld1*(aiw_ir1(1, ib)+(aiw_ir1(2, ib)+(aiw_ir1(3, ib)+&
3568 & aiw_ir1(4, ib)*reff(k, 1))*reff(k, 1))*reff(k, 1))
3569  CALL pushreal8(w2)
3570  w2 = taucld2*(aww_ir1(1, ib)+(aww_ir1(2, ib)+(aww_ir1(3, ib)+&
3571 & aww_ir1(4, ib)*reff(k, 2))*reff(k, 2))*reff(k, 2))
3572  w3 = taucld3*0.54
3573  CALL pushreal8(w4)
3574  w4 = taucld4*(aiw_ir1(1, ib)+(aiw_ir1(2, ib)+(aiw_ir1(3, ib)+&
3575 & aiw_ir1(4, ib)*reff_snow)*reff_snow)*reff_snow)
3576  ww = (w1+w2+w3+w4)/tauc
3577  CALL pushreal8(g1)
3578  g1 = w1*(aig_ir1(1, ib)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(4&
3579 & , ib)*reff(k, 1))*reff(k, 1))*reff(k, 1))
3580  CALL pushreal8(g2)
3581  g2 = w2*(awg_ir1(1, ib)+(awg_ir1(2, ib)+(awg_ir1(3, ib)+awg_ir1(4&
3582 & , ib)*reff(k, 2))*reff(k, 2))*reff(k, 2))
3583  g3 = w3*0.95
3584  CALL pushreal8(g4)
3585  g4 = w4*(aig_ir1(1, ib)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(4&
3586 & , ib)*reff_snow)*reff_snow)*reff_snow)
3587  IF (w1 + w2 + w3 + w4 .GE. 0.) THEN
3588  abs0 = w1 + w2 + w3 + w4
3589  ELSE
3590  abs0 = -(w1+w2+w3+w4)
3591  END IF
3592  IF (abs0 .GT. 0.0) THEN
3593  CALL pushreal8(gg)
3594  gg = (g1+g2+g3+g4)/(w1+w2+w3+w4)
3595  CALL pushcontrol1b(1)
3596  ELSE
3597  CALL pushreal8(gg)
3598  gg = 0.5
3599  CALL pushcontrol1b(0)
3600  END IF
3601 !-----Parameterization of LW scattering following Eqs. (6.11)
3602 ! and (6.12).
3603  ff = 0.5 + (0.3739+(0.0076+0.1185*gg)*gg)*gg
3604  IF (1. - ww*ff .LT. 0.0) THEN
3605  CALL pushreal8(max1)
3606  max1 = 0.0
3607  CALL pushcontrol1b(0)
3608  ELSE
3609  CALL pushreal8(max1)
3610  max1 = 1. - ww*ff
3611  CALL pushcontrol1b(1)
3612  END IF
3613 !ALT: temporary protection against negative cloud optical thickness
3614  CALL pushreal8(tauc)
3615  tauc = max1*tauc
3616 !-----compute cloud diffuse transmittance. It is approximated by using
3617 ! a diffusivity factor of 1.66.
3618  CALL pushreal8(tcldlyr(k))
3619  tcldlyr(k) = exp(-(1.66*tauc))
3620 ! N in the documentation (6.13)
3621  CALL pushcontrol1b(1)
3622  ELSE
3623  CALL pushcontrol1b(0)
3624  END IF
3625  END DO
3626  DO k=nlevs,1,-1
3627  CALL popcontrol1b(branch)
3628  IF (branch .EQ. 0) THEN
3629  ennb(k) = 0.0_8
3630  tcldlyrb(k) = 0.0_8
3631  reff_snowb = 0.0_8
3632  taucld1b = 0.0_8
3633  taucld2b = 0.0_8
3634  taucld3b = 0.0_8
3635  taucld4b = 0.0_8
3636  taucb = 0.0_8
3637  ELSE
3638  fcldb(k) = fcldb(k) + (1.0-tcldlyr(k))*ennb(k)
3639  tcldlyrb(k) = tcldlyrb(k) - fcld(k)*ennb(k)
3640  ennb(k) = 0.0_8
3641  CALL popreal8(tcldlyr(k))
3642  taucb = -(exp(-(1.66*tauc))*1.66*tcldlyrb(k))
3643  tcldlyrb(k) = 0.0_8
3644  CALL popreal8(tauc)
3645  max1b = tauc*taucb
3646  taucb = max1*taucb
3647  taucld3 = 0.00307*(dp(k)*1.0e3/cons_grav*hydromets(k, 3))
3648  w3 = taucld3*0.54
3649  ww = (w1+w2+w3+w4)/tauc
3650  ff = 0.5 + (0.3739+(0.0076+0.1185*gg)*gg)*gg
3651  CALL popcontrol1b(branch)
3652  IF (branch .EQ. 0) THEN
3653  CALL popreal8(max1)
3654  wwb = 0.0_8
3655  ffb = 0.0_8
3656  ELSE
3657  CALL popreal8(max1)
3658  wwb = -(ff*max1b)
3659  ffb = -(ww*max1b)
3660  END IF
3661  ggb = ((0.1185*gg+0.0076)*gg+gg*(0.1185*gg+0.0076)+gg**2*0.1185+&
3662 & 0.3739)*ffb
3663  g3 = w3*0.95
3664  CALL popcontrol1b(branch)
3665  IF (branch .EQ. 0) THEN
3666  CALL popreal8(gg)
3667  w1b = 0.0_8
3668  w2b = 0.0_8
3669  w3b = 0.0_8
3670  w4b = 0.0_8
3671  g1b = 0.0_8
3672  g2b = 0.0_8
3673  g3b = 0.0_8
3674  g4b = 0.0_8
3675  ELSE
3676  CALL popreal8(gg)
3677  tempb11 = ggb/(w1+w2+w3+w4)
3678  tempb12 = -((g1+g2+g3+g4)*tempb11/(w1+w2+w3+w4))
3679  g1b = tempb11
3680  g2b = tempb11
3681  g3b = tempb11
3682  g4b = tempb11
3683  w1b = tempb12
3684  w2b = tempb12
3685  w3b = tempb12
3686  w4b = tempb12
3687  END IF
3688  tempb5 = wwb/tauc
3689  CALL popreal8(g4)
3690  temp16 = aig_ir1(3, ib) + aig_ir1(4, ib)*reff_snow
3691  temp15 = aig_ir1(2, ib) + temp16*reff_snow
3692  tempb4 = w4*g4b
3693  w4b = w4b + tempb5 + (aig_ir1(1, ib)+temp15*reff_snow)*g4b
3694  w3b = w3b + tempb5 + 0.95*g3b
3695  CALL popreal8(g2)
3696  temp14 = awg_ir1(3, ib) + awg_ir1(4, ib)*reff(k, 2)
3697  temp13 = awg_ir1(2, ib) + temp14*reff(k, 2)
3698  tempb7 = w2*reff(k, 2)*g2b
3699  w2b = w2b + tempb5 + (awg_ir1(1, ib)+temp13*reff(k, 2))*g2b
3700  reffb(k, 2) = reffb(k, 2) + w2*temp13*g2b + (temp14+reff(k, 2)*&
3701 & awg_ir1(4, ib))*tempb7
3702  CALL popreal8(g1)
3703  temp12 = aig_ir1(3, ib) + aig_ir1(4, ib)*reff(k, 1)
3704  temp11 = aig_ir1(2, ib) + temp12*reff(k, 1)
3705  tempb8 = w1*reff(k, 1)*g1b
3706  w1b = w1b + tempb5 + (aig_ir1(1, ib)+temp11*reff(k, 1))*g1b
3707  reffb(k, 1) = reffb(k, 1) + w1*temp11*g1b + (temp12+reff(k, 1)*&
3708 & aig_ir1(4, ib))*tempb8
3709  taucb = taucb - (w1+w2+w3+w4)*tempb5/tauc
3710  CALL popreal8(w4)
3711  temp10 = aiw_ir1(3, ib) + aiw_ir1(4, ib)*reff_snow
3712  temp9 = aiw_ir1(2, ib) + temp10*reff_snow
3713  tempb6 = taucld4*w4b
3714  reff_snowb = (temp9+reff_snow*temp10+reff_snow**2*aiw_ir1(4, ib))*&
3715 & tempb6 + (temp15+reff_snow*temp16+reff_snow**2*aig_ir1(4, ib))*&
3716 & tempb4
3717  taucld4b = (aiw_ir1(1, ib)+temp9*reff_snow)*w4b
3718  taucld3b = 0.54*w3b
3719  CALL popreal8(w2)
3720  temp8 = aww_ir1(3, ib) + aww_ir1(4, ib)*reff(k, 2)
3721  temp7 = aww_ir1(2, ib) + temp8*reff(k, 2)
3722  tempb9 = taucld2*reff(k, 2)*w2b
3723  taucld2b = (aww_ir1(1, ib)+temp7*reff(k, 2))*w2b
3724  reffb(k, 2) = reffb(k, 2) + taucld2*temp7*w2b + (temp8+reff(k, 2)*&
3725 & aww_ir1(4, ib))*tempb9
3726  CALL popreal8(w1)
3727  temp6 = aiw_ir1(3, ib) + aiw_ir1(4, ib)*reff(k, 1)
3728  temp5 = aiw_ir1(2, ib) + temp6*reff(k, 1)
3729  tempb10 = taucld1*reff(k, 1)*w1b
3730  taucld1b = (aiw_ir1(1, ib)+temp5*reff(k, 1))*w1b
3731  reffb(k, 1) = reffb(k, 1) + taucld1*temp5*w1b + (temp6+reff(k, 1)*&
3732 & aiw_ir1(4, ib))*tempb10
3733  END IF
3734  CALL popreal8(tauc)
3735  taucld1b = taucld1b + taucb
3736  taucld2b = taucld2b + taucb
3737  taucld3b = taucld3b + taucb
3738  taucld4b = taucld4b + taucb
3739  CALL popcontrol1b(branch)
3740  IF (branch .EQ. 0) THEN
3741  CALL popreal8(taucld4)
3742  ELSE
3743  CALL popreal8(taucld4)
3744  temp4 = reff_snow**aib_ir1(3, ib)
3745  temp3 = aib_ir1(2, ib)/temp4
3746  tempb3 = dp(k)*1.0e3*taucld4b
3747  hydrometsb(k, 4) = hydrometsb(k, 4) + (aib_ir1(1, ib)+temp3)*&
3748 & tempb3/cons_grav
3749  IF (.NOT.(reff_snow .LE. 0.0 .AND. (aib_ir1(3, ib) .EQ. 0.0 .OR. &
3750 & aib_ir1(3, ib) .NE. int(aib_ir1(3, ib))))) reff_snowb = &
3751 & reff_snowb - temp3*hydromets(k, 4)*aib_ir1(3, ib)*reff_snow**(&
3752 & aib_ir1(3, ib)-1)*tempb3/(temp4*cons_grav)
3753  END IF
3754  CALL popcontrol1b(branch)
3755  IF (branch .EQ. 0) THEN
3756  CALL popreal8(reff_snow)
3757  ELSE
3758  CALL popreal8(reff_snow)
3759  reffb(k, 4) = reffb(k, 4) + reff_snowb
3760  END IF
3761  hydrometsb(k, 3) = hydrometsb(k, 3) + dp(k)*1.0e3*0.00307*taucld3b/&
3762 & cons_grav
3763  CALL popreal8(taucld2)
3764  temp2 = awb_ir1(3, ib) + awb_ir1(4, ib)*reff(k, 2)
3765  temp1 = awb_ir1(2, ib) + temp2*reff(k, 2)
3766  tempb0 = dp(k)*1.0e3*taucld2b
3767  tempb1 = hydromets(k, 2)*tempb0/cons_grav
3768  tempb2 = reff(k, 2)*tempb1
3769  hydrometsb(k, 2) = hydrometsb(k, 2) + (awb_ir1(1, ib)+temp1*reff(k, &
3770 & 2))*tempb0/cons_grav
3771  reffb(k, 2) = reffb(k, 2) + temp1*tempb1 + (temp2+reff(k, 2)*awb_ir1&
3772 & (4, ib))*tempb2
3773  CALL popcontrol1b(branch)
3774  IF (branch .EQ. 0) THEN
3775  CALL popreal8(taucld1)
3776  temp0 = reff(k, 1)**aib_ir1(3, ib)
3777  temp = aib_ir1(2, ib)/temp0
3778  tempb = dp(k)*1.0e3*taucld1b
3779  hydrometsb(k, 1) = hydrometsb(k, 1) + (aib_ir1(1, ib)+temp)*tempb/&
3780 & cons_grav
3781  IF (.NOT.(reff(k, 1) .LE. 0.0 .AND. (aib_ir1(3, ib) .EQ. 0.0 .OR. &
3782 & aib_ir1(3, ib) .NE. int(aib_ir1(3, ib))))) reffb(k, 1) = reffb&
3783 & (k, 1) - temp*hydromets(k, 1)*aib_ir1(3, ib)*reff(k, 1)**(&
3784 & aib_ir1(3, ib)-1)*tempb/(temp0*cons_grav)
3785  ELSE
3786  CALL popreal8(taucld1)
3787  END IF
3788  END DO
3789  ennb(0) = 0.0_8
3790  tcldlyrb(0) = 0.0_8
3791 END SUBROUTINE getirtau1_b
3792 
3793 end module irrad_ad
irradmod::b10exps
subroutine, public b10exps(np, dh2o, dcont, dco2, dn2o, pa, dt, h2oexp, conexp, co2exp, n2oexp)
Definition: irrad.F90:1246
irradmod::h2okdis
subroutine, public h2okdis(ib, np, k, fkw, gkw, ne, h2oexp, conexp, th2o, tcon, tran)
Definition: irrad.F90:1484
popinteger4
subroutine popinteger4(x)
Definition: adBuffer.f:541
irrad_ad::b10kdis_b
subroutine b10kdis_b(np, k, h2oexp, h2oexpb, conexp, conexpb, co2exp, co2expb, n2oexp, n2oexpb, th2o, th2ob, tcon, tconb, tco2, tco2b, tn2o, tn2ob, tran, tranb)
Definition: irrad_ad.F90:3163
popcontrol2b
subroutine popcontrol2b(cc)
Definition: adBuffer.f:146
irradmod::ch4exps
subroutine, public ch4exps(ib, np, dch4, pa, dt, ch4exp)
Definition: irrad.F90:1063
irradmod::cfckdis
subroutine, public cfckdis(np, k, cfcexp, tcfc, tran)
Definition: irrad.F90:1827
irrad_ad
Definition: irrad_ad.F90:1
irradmod::ch4kdis
subroutine, public ch4kdis(ib, np, k, ch4exp, tch4, tran)
Definition: irrad.F90:1682
popinteger4array
void popinteger4array(int *x, int n)
Definition: adStack.c:335
popreal8array
void popreal8array(double *x, int n)
Definition: adStack.c:375
irradmod::cfcexps
subroutine, public cfcexps(ib, np, a1, b1, fk1, a2, b2, fk2, dcfc, dt, cfcexp)
Definition: irrad.F90:1188
type_kinds::ip
integer, parameter, public ip
Definition: Type_Kinds.f90:97
irradmod::n2oexps
subroutine, public n2oexps(ib, np, dn2o, pa, dt, n2oexp)
Definition: irrad.F90:991
irradmod::n2okdis
subroutine, public n2okdis(ib, np, k, n2oexp, tn2o, tran)
Definition: irrad.F90:1612
irradmod::comexps
subroutine, public comexps(ib, np, dcom, dt, comexp)
Definition: irrad.F90:1132
pushcontrol1b
subroutine pushcontrol1b(cc)
Definition: adBuffer.f:115
irradmod::cldovlp
subroutine, public cldovlp(np, k1, k2, ict, icb, icx, ncld, enn, ett, cldhi, cldmd, cldlw)
Definition: irrad.F90:1973
irradmod::planck
subroutine, public planck(ibn, cb, t, xlayer)
Definition: irrad.F90:819
irrad_ad::tablup_b
subroutine tablup_b(nx1, nh1, dw, dwb, p, dt, dtb, s1, s1b, s2, s2b, s3, s3b, w1, p1, dwe, dpe, coef1, coef2, coef3, tran, tranb)
Definition: irrad_ad.F90:2575
irradmod::conexps
subroutine, public conexps(ib, np, dcont, xke, conexp)
Definition: irrad.F90:940
irrad_ad::planck_b
subroutine planck_b(ibn, cb, t, tb, xlayer, xlayerb)
Definition: irrad_ad.F90:1869
pushinteger4array
void pushinteger4array(int *x, int n)
Definition: adStack.c:332
pushcontrol2b
subroutine pushcontrol2b(cc)
Definition: adBuffer.f:140
irrad_ad::n2oexps_b
subroutine n2oexps_b(ib, np, dn2o, pa, dt, dtb, n2oexp, n2oexpb)
Definition: irrad_ad.F90:2061
pushreal8array
void pushreal8array(double *x, int n)
Definition: adStack.c:372
popreal8
subroutine popreal8(x)
Definition: adBuffer.f:820
irrad_ad::cfcexps_b
subroutine cfcexps_b(ib, np, a1, b1, fk1, a2, b2, fk2, dcfc, dt, dtb, cfcexp, cfcexpb)
Definition: irrad_ad.F90:2329
irradmod::h2oexps
subroutine, public h2oexps(ib, np, dh2o, pa, dt, xkw, aw, bw, pm, mw, h2oexp)
Definition: irrad.F90:855
irrad_ad::conexps_b
subroutine conexps_b(ib, np, dcont, dcontb, xke, conexp, conexpb)
Definition: irrad_ad.F90:2015
irradmod::tablup
subroutine, public tablup(nx1, nh1, dw, p, dt, s1, s2, s3, w1, p1, dwe, dpe, coef1, coef2, coef3, tran)
Definition: irrad.F90:1364
irradmod::sfcflux
subroutine, public sfcflux(ibn, cb, dcb, ns, fs, tg, eg, tv, ev, rv, bs, dbs, rflxs)
Definition: irrad.F90:2067
irrad_ad::b10exps_b
subroutine b10exps_b(np, dh2o, dh2ob, dcont, dcontb, dco2, dn2o, pa, dt, dtb, h2oexp, h2oexpb, conexp, conexpb, co2exp, co2expb, n2oexp, n2oexpb)
Definition: irrad_ad.F90:2378
irrad_ad::comexps_b
subroutine comexps_b(ib, np, dcom, dt, dtb, comexp, comexpb)
Definition: irrad_ad.F90:2252
irrad_ad::cldovlp_b
subroutine cldovlp_b(np, k1, k2, ict, icb, icx, ncld, enn, ennb, ett, ettb, cldhi, cldhib, cldmd, cldmdb, cldlw, cldlwb)
Definition: irrad_ad.F90:3301
irradmod::getirtau1
subroutine, public getirtau1(ib, nlevs, dp, fcld, reff, hydromets, tcldlyr, enn, aib_ir1, awb_ir1, aiw_ir1, aww_ir1, aig_ir1, awg_ir1, CONS_GRAV)
Definition: irrad.F90:2243
irradmod::b10kdis
subroutine, public b10kdis(np, k, h2oexp, conexp, co2exp, n2oexp, th2o, tcon, tco2, tn2o, tran)
Definition: irrad.F90:1864
irrad_ad::comkdis_b
subroutine comkdis_b(ib, np, k, comexp, comexpb, tcom, tcomb, tran, tranb)
Definition: irrad_ad.F90:3017
pushreal8
subroutine pushreal8(x)
Definition: adBuffer.f:763
irrad_ad::getirtau1_b
subroutine getirtau1_b(ib, nlevs, dp, fcld, fcldb, reff, reffb, hydromets, hydrometsb, tcldlyr, tcldlyrb, enn, ennb, aib_ir1, awb_ir1, aiw_ir1, aww_ir1, aig_ir1, awg_ir1, cons_grav)
Definition: irrad_ad.F90:3423
popcontrol1b
subroutine popcontrol1b(cc)
Definition: adBuffer.f:120
irrad_ad::h2oexps_b
subroutine h2oexps_b(ib, np, dh2o, dh2ob, pa, dt, dtb, xkw, aw, bw, pm, mw, h2oexp, h2oexpb)
Definition: irrad_ad.F90:1898
max
#define max(a, b)
Definition: mosaic_util.h:33
irrad_ad::ch4exps_b
subroutine ch4exps_b(ib, np, dch4, pa, dt, dtb, ch4exp, ch4expb)
Definition: irrad_ad.F90:2162
irrad_ad::n2okdis_b
subroutine n2okdis_b(ib, np, k, n2oexp, n2oexpb, tn2o, tn2ob, tran, tranb)
Definition: irrad_ad.F90:2868
popcontrol4b
subroutine popcontrol4b(cc)
Definition: adBuffer.f:207
irradmod::comkdis
subroutine, public comkdis(ib, np, k, comexp, tcom, tran)
Definition: irrad.F90:1749
irrad_ad::h2okdis_b
subroutine h2okdis_b(ib, np, k, fkw, gkw, ne, h2oexp, h2oexpb, conexp, conexpb, th2o, th2ob, tcon, tconb, tran, tranb)
Definition: irrad_ad.F90:2735
irrad_ad::cfckdis_b
subroutine cfckdis_b(np, k, cfcexp, cfcexpb, tcfc, tcfcb, tran, tranb)
Definition: irrad_ad.F90:3136
irrad_ad::irrad_b
subroutine, public irrad_b(np, ple_dev, ta_dev, ta_devb, wa_dev, wa_devb, oa_dev, oa_devb, tb_dev, tb_devb, co2, trace, n2o_dev, ch4_dev, cfc11_dev, cfc12_dev, cfc22_dev, cwc_dev, cwc_devb, fcld_dev, fcld_devb, ict, icb, reff_dev, reff_devb, ns, fs_dev, tg_dev, eg_dev, tv_dev, ev_dev, rv_dev, na, nb, taua_dev, taua_devb, ssaa_dev, ssaa_devb, asya_dev, asya_devb, flxu_dev, flxu_devb, flxd_dev, flxd_devb, dfdts_dev, dfdts_devb, aib_ir, awb_ir, aiw_ir, aww_ir, aig_ir, awg_ir, xkw, xke, mw, aw, bw, pm, fkw, gkw, cb, dcb, w11, w12, w13, p11, p12, p13, dwe, dpe, c1, c2, c3, oo1, oo2, oo3, h11, h12, h13, h21, h22, h23, h81, h82, h83)
Definition: irrad_ad.F90:22
pushcontrol4b
subroutine pushcontrol4b(cc)
Definition: adBuffer.f:199
irrad_ad::ch4kdis_b
subroutine ch4kdis_b(ib, np, k, ch4exp, ch4expb, tch4, tch4b, tran, tranb)
Definition: irrad_ad.F90:2946
min
#define min(a, b)
Definition: mosaic_util.h:32
irradmod::mkicx
subroutine, public mkicx(np, ict, icb, enn, icx, ncld)
Definition: irrad.F90:2185
pushinteger4
subroutine pushinteger4(x)
Definition: adBuffer.f:484