FV3 Bundle
irrad_tl.F90
Go to the documentation of this file.
1 module irrad_tl
2 
3 use irradmod, only: sfcflux, mkicx
4 
5 IMPLICIT NONE
6 
7 PRIVATE
8 PUBLIC :: irrad_d
9 
10 contains
11 
12 SUBROUTINE irrad_d(np, ple_dev, ta_dev, ta_devd, wa_dev, wa_devd, oa_dev&
13 & , oa_devd, tb_dev, tb_devd, co2, trace, n2o_dev, ch4_dev, cfc11_dev, &
14 & cfc12_dev, cfc22_dev, cwc_dev, cwc_devd, fcld_dev, fcld_devd, ict, icb&
15 & , reff_dev, reff_devd, ns, fs_dev, tg_dev, eg_dev, tv_dev, ev_dev, &
16 & rv_dev, na, nb, taua_dev, taua_devd, ssaa_dev, ssaa_devd, asya_dev, &
17 & asya_devd, flxu_dev, flxu_devd, flxd_dev, flxd_devd, dfdts_dev, &
18 & dfdts_devd, 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, INTENT(IN) :: tb_devd
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), INTENT(IN) :: ta_devd, wa_devd, oa_devd, &
46 & fcld_devd
47  REAL*8, DIMENSION(np), INTENT(IN) :: n2o_dev, ch4_dev, cfc11_dev, &
48 & cfc12_dev, cfc22_dev
49  REAL*8, DIMENSION(np + 1), INTENT(IN) :: ple_dev
50 !Rank 3 (surface types)
51  REAL*8, DIMENSION(ns), INTENT(IN) :: fs_dev, tg_dev, tv_dev
52  REAL*8, DIMENSION(ns, 10), INTENT(IN) :: eg_dev, ev_dev, rv_dev
53 !Rank 3 (diagnostic cloud parts)
54  REAL*8, DIMENSION(np, 4), INTENT(IN) :: cwc_dev, reff_dev
55  REAL*8, DIMENSION(np, 4), INTENT(IN) :: cwc_devd, reff_devd
56 !Rank 3 (aerosols)
57  REAL*8, DIMENSION(np, nb), INTENT(INOUT) :: taua_dev, ssaa_dev, &
58 & asya_dev
59  REAL*8, DIMENSION(np, nb), INTENT(INOUT) :: taua_devd, ssaa_devd, &
60 & asya_devd
61 !----- OUPUTS -----
62  REAL*8, DIMENSION(np + 1), INTENT(OUT) :: flxu_dev
63  REAL*8, DIMENSION(np+1), INTENT(OUT) :: flxu_devd
64  REAL*8, DIMENSION(np + 1), INTENT(OUT) :: flxd_dev
65  REAL*8, DIMENSION(np+1), INTENT(OUT) :: flxd_devd
66  REAL*8, DIMENSION(np + 1), INTENT(OUT) :: dfdts_dev
67  REAL*8, DIMENSION(np+1), INTENT(OUT) :: dfdts_devd
68 !----- LOCALS -----
69  REAL*8, PARAMETER :: cons_grav=9.80665
70  INTEGER, PARAMETER :: nx1=26
71  INTEGER, PARAMETER :: no1=21
72  INTEGER, PARAMETER :: nc1=30
73  INTEGER, PARAMETER :: nh1=31
74 !Temporary arrays
75  REAL*8 :: pa(0:np), dt(0:np)
76  REAL*8 :: pad(0:np), dtd(0:np)
77  REAL*8 :: x1, x2, x3
78  REAL*8 :: x1d, x2d, x3d
79  REAL*8 :: dh2o(0:np), dcont(0:np), dco2(0:np), do3(0:np)
80  REAL*8 :: dh2od(0:np), dcontd(0:np), dco2d(0:np), do3d(0:np)
81  REAL*8 :: dn2o(0:np), dch4(0:np)
82  REAL*8 :: dn2od(0:np), dch4d(0:np)
83  REAL*8 :: df11(0:np), df12(0:np), df22(0:np)
84  REAL*8 :: df11d(0:np), df12d(0:np), df22d(0:np)
85  REAL*8 :: th2o(6), tcon(3), tco2(6)
86  REAL*8 :: th2od(6), tcond(3), tco2d(6)
87  REAL*8 :: tn2o(4), tch4(4), tcom(6)
88  REAL*8 :: tn2od(4), tch4d(4), tcomd(6)
89  REAL*8 :: tf11, tf12, tf22
90  REAL*8 :: tf11d, tf12d, tf22d
91  REAL*8 :: blayer(0:np+1), blevel(0:np+1)
92  REAL*8 :: blayerd(0:np+1), bleveld(0:np+1)
93  REAL*8 :: bd(0:np+1), bu(0:np+1)
94  REAL*8 :: bdd(0:np+1), bud(0:np+1)
95  REAL*8 :: bs, dbs, rflxs
96  REAL*8 :: dp(0:np)
97  REAL*8 :: dpd(0:np)
98  REAL*8 :: trant, tranal
99  REAL*8 :: trantd, tranald
100  REAL*8 :: transfc(0:np+1)
101  REAL*8 :: transfcd(0:np+1)
102  REAL*8 :: flxu(0:np+1), flxd(0:np+1)
103  REAL*8 :: flxud(0:np+1), flxdd(0:np+1)
104  REAL*8 :: taerlyr(0:np)
105  REAL*8 :: taerlyrd(0:np)
106 !OVERCAST
107  INTEGER :: ncld(3)
108  INTEGER :: icx(0:np)
109 !OVERCAST
110  INTEGER :: idx, rc
111  INTEGER :: k, l, ip, iw, ibn, ik, iq, isb, k1, k2, ne
112  REAL*8 :: enn(0:np)
113  REAL*8 :: ennd(0:np)
114  REAL*8 :: cldhi, cldmd, cldlw, tcldlyr(0:np), fclr
115  REAL*8 :: cldhid, cldmdd, cldlwd, tcldlyrd(0:np), fclrd
116  REAL*8 :: x, xx, yy, p1, a1, b1, fk1, a2, b2, fk2
117  REAL*8 :: xxd, yyd
118  REAL*8 :: w1, ff
119  REAL*8 :: ffd
120  LOGICAL :: oznbnd, co2bnd, h2otable, conbnd, n2obnd
121  LOGICAL :: ch4bnd, combnd, f11bnd, f12bnd, f22bnd, b10bnd
122  LOGICAL :: do_aerosol
123 !Temp arrays and variables for consolidation of tables
124  INTEGER, PARAMETER :: max_num_tables=17
125  REAL*8 :: exptbl(0:np, max_num_tables)
126  REAL*8 :: exptbld(0:np, max_num_tables)
127  TYPE band_table
128  INTEGER :: start
129  INTEGER :: end
130  END TYPE band_table
131  TYPE(band_table) :: h2oexp
132  TYPE(band_table) :: conexp
133  TYPE(band_table) :: co2exp
134  TYPE(band_table) :: n2oexp
135  TYPE(band_table) :: ch4exp
136  TYPE(band_table) :: comexp
137  TYPE(band_table) :: f11exp
138  TYPE(band_table) :: f12exp
139  TYPE(band_table) :: f22exp
140 !Variables for new getirtau routine
141  REAL*8 :: dp_pa(np)
142  REAL*8 :: dp_pad(np)
143  REAL*8 :: fcld_col(np)
144  REAL*8 :: fcld_cold(np)
145  REAL*8 :: reff_col(np, 4)
146  REAL*8 :: reff_cold(np, 4)
147  REAL*8 :: cwc_col(np, 4)
148  REAL*8 :: cwc_cold(np, 4)
149  REAL*8 :: h2oexp_tmp(0:np, 5), conexp_tmp(0:np), co2exp_tmp(0:np, 6), &
150 & n2oexp_tmp(0:np, 2)
151  REAL*8 :: h2oexp_tmpd(0:np, 5), co2exp_tmpd(0:np, 6), n2oexp_tmpd(0:np&
152 & , 2)
153  INTRINSIC max
154  INTRINSIC exp
155  INTRINSIC min
156  INTRINSIC log
157  dh2od = 0.0_8
158  dcontd = 0.0_8
159  dtd = 0.0_8
160  reff_cold = 0.0_8
161  fcld_cold = 0.0_8
162  cwc_cold = 0.0_8
163  do3d = 0.0_8
164 !BEGIN CALCULATIONS ...
165 !compute layer pressure (pa) and layer temperature minus 250K (dt)
166  DO k=1,np
167  pad(k) = 0.0_8
168  pa(k) = 0.5*(ple_dev(k+1)+ple_dev(k))*0.01
169  dpd(k) = 0.0_8
170  dp(k) = (ple_dev(k+1)-ple_dev(k))*0.01
171 ! dp in Pascals for getirtau
172  dp_pad(k) = 0.0_8
173  dp_pa(k) = ple_dev(k+1) - ple_dev(k)
174  dtd(k) = ta_devd(k)
175  dt(k) = ta_dev(k) - 250.0
176 !compute layer absorber amount
177 !dh2o : water vapor amount (g/cm**2)
178 !dcont: scaled water vapor amount for continuum absorption
179 ! (g/cm**2)
180 !dco2 : co2 amount (cm-atm)stp
181 !do3 : o3 amount (cm-atm)stp
182 !dn2o : n2o amount (cm-atm)stp
183 !dch4 : ch4 amount (cm-atm)stp
184 !df11 : cfc11 amount (cm-atm)stp
185 !df12 : cfc12 amount (cm-atm)stp
186 !df22 : cfc22 amount (cm-atm)stp
187 !the factor 1.02 is equal to 1000/980
188 !factors 789 and 476 are for unit conversion
189 !the factor 0.001618 is equal to 1.02/(.622*1013.25)
190 !the factor 6.081 is equal to 1800/296
191  dh2od(k) = 1.02*dp(k)*wa_devd(k)
192  dh2o(k) = 1.02*wa_dev(k)*dp(k)
193  do3d(k) = 476.*dp(k)*oa_devd(k)
194  do3(k) = 476.*oa_dev(k)*dp(k)
195  dco2d(k) = 0.0_8
196  dco2(k) = 789.*co2*dp(k)
197  dch4d(k) = 0.0_8
198  dch4(k) = 789.*ch4_dev(k)*dp(k)
199  dn2od(k) = 0.0_8
200  dn2o(k) = 789.*n2o_dev(k)*dp(k)
201  df11d(k) = 0.0_8
202  df11(k) = 789.*cfc11_dev(k)*dp(k)
203  df12d(k) = 0.0_8
204  df12(k) = 789.*cfc12_dev(k)*dp(k)
205  df22d(k) = 0.0_8
206  df22(k) = 789.*cfc22_dev(k)*dp(k)
207  IF (dh2o(k) .LT. 1.e-10) THEN
208  dh2od(k) = 0.0_8
209  dh2o(k) = 1.e-10
210  ELSE
211  dh2o(k) = dh2o(k)
212  END IF
213  IF (do3(k) .LT. 1.e-6) THEN
214  do3d(k) = 0.0_8
215  do3(k) = 1.e-6
216  ELSE
217  do3(k) = do3(k)
218  END IF
219  IF (dco2(k) .LT. 1.e-4) THEN
220  dco2d(k) = 0.0_8
221  dco2(k) = 1.e-4
222  ELSE
223  dco2d(k) = 0.0_8
224  dco2(k) = dco2(k)
225  END IF
226 !Compute scaled water vapor amount for h2o continuum absorption
227 !following eq. (4.21).
228  xxd = pa(k)*0.001618*dp(k)*(wa_devd(k)*wa_dev(k)+wa_dev(k)*wa_devd(k&
229 & ))
230  xx = pa(k)*0.001618*wa_dev(k)*wa_dev(k)*dp(k)
231  dcontd(k) = xxd*exp(1800./ta_dev(k)-6.081) - xx*1800.*ta_devd(k)*exp&
232 & (1800./ta_dev(k)-6.081)/ta_dev(k)**2
233  dcont(k) = xx*exp(1800./ta_dev(k)-6.081)
234 !Fill the reff, cwc, and fcld for the column
235  fcld_cold(k) = fcld_devd(k)
236  fcld_col(k) = fcld_dev(k)
237  DO l=1,4
238  reff_cold(k, l) = reff_devd(k, l)
239  reff_col(k, l) = reff_dev(k, l)
240  cwc_cold(k, l) = cwc_devd(k, l)
241  cwc_col(k, l) = cwc_dev(k, l)
242  END DO
243  END DO
244  IF (ple_dev(1)*0.01 .LT. 0.005) THEN
245  dpd(0) = 0.0_8
246  dp(0) = 0.005
247  ELSE
248  dpd(0) = 0.0_8
249  dp(0) = ple_dev(1)*0.01
250  END IF
251  pad(0) = 0.0_8
252  pa(0) = 0.5*dp(0)
253  dtd(0) = ta_devd(1)
254  dt(0) = ta_dev(1) - 250.0
255  dh2od(0) = 1.02*dp(0)*wa_devd(1)
256  dh2o(0) = 1.02*wa_dev(1)*dp(0)
257  do3d(0) = 476.*dp(0)*oa_devd(1)
258  do3(0) = 476.*oa_dev(1)*dp(0)
259  dco2d(0) = 0.0_8
260  dco2(0) = 789.*co2*dp(0)
261  dch4d(0) = 0.0_8
262  dch4(0) = 789.*ch4_dev(1)*dp(0)
263  dn2od(0) = 0.0_8
264  dn2o(0) = 789.*n2o_dev(1)*dp(0)
265  df11d(0) = 0.0_8
266  df11(0) = 789.*cfc11_dev(1)*dp(0)
267  df12d(0) = 0.0_8
268  df12(0) = 789.*cfc12_dev(1)*dp(0)
269  df22d(0) = 0.0_8
270  df22(0) = 789.*cfc22_dev(1)*dp(0)
271  IF (dh2o(0) .LT. 1.e-10) THEN
272  dh2od(0) = 0.0_8
273  dh2o(0) = 1.e-10
274  ELSE
275  dh2o(0) = dh2o(0)
276  END IF
277  IF (do3(0) .LT. 1.e-6) THEN
278  do3d(0) = 0.0_8
279  do3(0) = 1.e-6
280  ELSE
281  do3(0) = do3(0)
282  END IF
283  IF (dco2(0) .LT. 1.e-4) THEN
284  dco2d(0) = 0.0_8
285  dco2(0) = 1.e-4
286  ELSE
287  dco2d(0) = 0.0_8
288  dco2(0) = dco2(0)
289  END IF
290  xxd = pa(0)*0.001618*dp(0)*(wa_devd(1)*wa_dev(1)+wa_dev(1)*wa_devd(1))
291  xx = pa(0)*0.001618*wa_dev(1)*wa_dev(1)*dp(0)
292  dcontd(0) = xxd*exp(1800./ta_dev(1)-6.081) - xx*1800.*ta_devd(1)*exp(&
293 & 1800./ta_dev(1)-6.081)/ta_dev(1)**2
294  dcont(0) = xx*exp(1800./ta_dev(1)-6.081)
295 !The surface (np+1) is treated as a layer filled with black clouds.
296 !transfc is the transmittance between the surface and a pressure
297 !level.
298  transfcd(np+1) = 0.0_8
299  transfc(np+1) = 1.0
300 !Initialize fluxes
301  DO k=1,np+1
302  flxu_devd(k) = 0.0_8
303  flxu_dev(k) = 0.0
304  flxd_devd(k) = 0.0_8
305  flxd_dev(k) = 0.0
306  dfdts_devd(k) = 0.0_8
307  dfdts_dev(k) = 0.0
308  END DO
309  n2oexp_tmpd = 0.0_8
310  blayerd = 0.0_8
311  tn2od = 0.0_8
312  transfcd = 0.0_8
313  bdd = 0.0_8
314  h2oexp_tmpd = 0.0_8
315  tcomd = 0.0_8
316  tf11d = 0.0_8
317  tf12d = 0.0_8
318  trantd = 0.0_8
319  bud = 0.0_8
320  th2od = 0.0_8
321  tcldlyrd = 0.0_8
322  tch4d = 0.0_8
323  tf22d = 0.0_8
324  ennd = 0.0_8
325  fclrd = 0.0_8
326  tco2d = 0.0_8
327  co2exp_tmpd = 0.0_8
328  taerlyrd = 0.0_8
329  bleveld = 0.0_8
330 !Integration over spectral bands
331  DO ibn=1,10
332  IF (ibn .EQ. 10 .AND. (.NOT.trace)) THEN
333  GOTO 100
334  ELSE
335 !if h2otable, compute h2o (line) transmittance using table look-up.
336 !if conbnd, compute h2o (continuum) transmittance in bands 2-7.
337 !if co2bnd, compute co2 transmittance in band 3.
338 !if oznbnd, compute o3 transmittance in band 5.
339 !if n2obnd, compute n2o transmittance in bands 6 and 7.
340 !if ch4bnd, compute ch4 transmittance in bands 6 and 7.
341 !if combnd, compute co2-minor transmittance in bands 4 and 5.
342 !if f11bnd, compute cfc11 transmittance in bands 4 and 5.
343 !if f12bnd, compute cfc12 transmittance in bands 4 and 6.
344 !if f22bnd, compute cfc22 transmittance in bands 4 and 6.
345 !if b10bnd, compute flux reduction due to n2o in band 10.
346  h2otable = (ibn .EQ. 1 .OR. ibn .EQ. 2) .OR. ibn .EQ. 8
347  conbnd = ibn .GE. 2 .AND. ibn .LE. 7
348  co2bnd = ibn .EQ. 3
349  oznbnd = ibn .EQ. 5
350  n2obnd = ibn .EQ. 6 .OR. ibn .EQ. 7
351  ch4bnd = ibn .EQ. 6 .OR. ibn .EQ. 7
352  combnd = ibn .EQ. 4 .OR. ibn .EQ. 5
353  f11bnd = ibn .EQ. 4 .OR. ibn .EQ. 5
354  f12bnd = ibn .EQ. 4 .OR. ibn .EQ. 6
355  f22bnd = ibn .EQ. 4 .OR. ibn .EQ. 6
356  b10bnd = ibn .EQ. 10
357  do_aerosol = na .GT. 0
358  exptbl = 0.0
359 !Control packing of the new exponential tables by band
360  SELECT CASE (ibn)
361  CASE (2)
362  conexp%start = 1
363  conexp%end = 1
364  CASE (3)
365  h2oexp%start = 1
366  h2oexp%end = 6
367  conexp%start = 7
368  conexp%end = 9
369  CASE (4)
370  h2oexp%start = 1
371  h2oexp%end = 6
372  conexp%start = 7
373  conexp%end = 7
374  comexp%start = 8
375  comexp%end = 13
376  f11exp%start = 14
377  f11exp%end = 14
378  f12exp%start = 15
379  f12exp%end = 15
380  f22exp%start = 16
381  f22exp%end = 16
382  CASE (5)
383  h2oexp%start = 1
384  h2oexp%end = 6
385  conexp%start = 7
386  conexp%end = 7
387  comexp%start = 8
388  comexp%end = 13
389  f11exp%start = 14
390  f11exp%end = 14
391  CASE (6)
392  h2oexp%start = 1
393  h2oexp%end = 6
394  conexp%start = 7
395  conexp%end = 7
396  n2oexp%start = 8
397  n2oexp%end = 11
398  ch4exp%start = 12
399  ch4exp%end = 15
400  f12exp%start = 16
401  f12exp%end = 16
402  f22exp%start = 17
403  f22exp%end = 17
404  CASE (7)
405  h2oexp%start = 1
406  h2oexp%end = 6
407  conexp%start = 7
408  conexp%end = 7
409  n2oexp%start = 8
410  n2oexp%end = 11
411  ch4exp%start = 12
412  ch4exp%end = 15
413  CASE (9)
414  h2oexp%start = 1
415  h2oexp%end = 6
416  CASE (10)
417  h2oexp%start = 1
418  h2oexp%end = 5
419  conexp%start = 6
420  conexp%end = 6
421  co2exp%start = 7
422  co2exp%end = 12
423  n2oexp%start = 13
424  n2oexp%end = 14
425  END SELECT
426 !blayer is the spectrally integrated planck flux of the mean layer
427 !temperature derived from eq. (3.11)
428 !The fitting for the planck flux is valid for the range 160-345 K.
429  DO k=1,np
430  CALL planck_d(ibn, cb, ta_dev(k), ta_devd(k), blayer(k), blayerd&
431 & (k))
432  END DO
433 !Index "0" is the layer above the top of the atmosphere.
434  blayerd(0) = blayerd(1)
435  blayer(0) = blayer(1)
436  bleveld(0) = blayerd(1)
437  blevel(0) = blayer(1)
438 !Surface emission and reflectivity. See Section 9.
439 !bs and dbs include the effect of surface emissivity.
440  CALL sfcflux(ibn, cb, dcb, ns, fs_dev, tg_dev, eg_dev, tv_dev, &
441 & ev_dev, rv_dev, bs, dbs, rflxs)
442  blayerd(np+1) = 0.0_8
443  blayer(np+1) = bs
444 !interpolate Planck function at model levels (linear in p)
445  DO k=2,np
446  bleveld(k) = (dp(k)*blayerd(k-1)+dp(k-1)*blayerd(k))/(dp(k-1)+dp&
447 & (k))
448  blevel(k) = (blayer(k-1)*dp(k)+blayer(k)*dp(k-1))/(dp(k-1)+dp(k)&
449 & )
450  END DO
451 !Extrapolate blevel(1) from blayer(2) and blayer(1)
452  bleveld(1) = blayerd(1) + dp(1)*(blayerd(1)-blayerd(2))/(dp(1)+dp(&
453 & 2))
454  blevel(1) = blayer(1) + (blayer(1)-blayer(2))*dp(1)/(dp(1)+dp(2))
455  bleveld(0) = bleveld(1)
456  blevel(0) = blevel(1)
457 !If the surface air temperature tb is known, compute blevel(np+1)
458  CALL planck_d(ibn, cb, tb_dev, tb_devd, blevel(np+1), bleveld(np+1&
459 & ))
460 !Compute cloud optical thickness following Eqs. (6.4a,b) and (6.7)
461 ! NOTE: dp_pa is only dims(1:np) as the 0'th level isn't needed in getirtau.
462 ! Plus, the pressures in getirtau *MUST* be in Pascals.
463 ! Slots for reff, hydrometeors and tauall are as follows:
464 ! 1 Cloud Ice
465 ! 2 Cloud Liquid
466 ! 3 Falling Liquid (Rain)
467 ! 4 Falling Ice (Snow)
468  CALL getirtau1_d(ibn, np, dp_pa, fcld_col, fcld_cold, reff_col, &
469 & reff_cold, cwc_col, cwc_cold, tcldlyr, tcldlyrd, enn, &
470 & ennd, aib_ir, awb_ir, aiw_ir, aww_ir, aig_ir, awg_ir, &
471 & cons_grav)
472  DO k=0,np
473  icx(k) = k
474  END DO
475  CALL mkicx(np, ict, icb, enn, icx, ncld)
476 !Compute optical thickness, single-scattering albedo and asymmetry
477 !factor for a mixture of "na" aerosol types. Eqs. (7.1)-(7.3)
478  IF (do_aerosol) THEN
479  taerlyrd(0) = 0.0_8
480  taerlyr(0) = 1.0
481  DO k=1,np
482 !-----taerlyr is the aerosol diffuse transmittance
483  taerlyrd(k) = 0.0_8
484  taerlyr(k) = 1.0
485  IF (taua_dev(k, ibn) .GT. 0.001) THEN
486  IF (ssaa_dev(k, ibn) .GT. 0.001) THEN
487  asya_devd(k, ibn) = (asya_devd(k, ibn)*ssaa_dev(k, ibn)-&
488 & asya_dev(k, ibn)*ssaa_devd(k, ibn))/ssaa_dev(k, ibn)**2
489  asya_dev(k, ibn) = asya_dev(k, ibn)/ssaa_dev(k, ibn)
490  ssaa_devd(k, ibn) = (ssaa_devd(k, ibn)*taua_dev(k, ibn)-&
491 & ssaa_dev(k, ibn)*taua_devd(k, ibn))/taua_dev(k, ibn)**2
492  ssaa_dev(k, ibn) = ssaa_dev(k, ibn)/taua_dev(k, ibn)
493 !Parameterization of aerosol scattering following
494  ffd = (0.1185*asya_devd(k, ibn)*asya_dev(k, ibn)+(0.0076+&
495 & 0.1185*asya_dev(k, ibn))*asya_devd(k, ibn))*asya_dev(k, &
496 & ibn) + (.3739+(0.0076+0.1185*asya_dev(k, ibn))*asya_dev(&
497 & k, ibn))*asya_devd(k, ibn)
498  ff = .5 + (.3739+(0.0076+0.1185*asya_dev(k, ibn))*asya_dev&
499 & (k, ibn))*asya_dev(k, ibn)
500  taua_devd(k, ibn) = taua_devd(k, ibn)*(1.-ssaa_dev(k, ibn)&
501 & *ff) + taua_dev(k, ibn)*(-(ssaa_devd(k, ibn)*ff)-&
502 & ssaa_dev(k, ibn)*ffd)
503  taua_dev(k, ibn) = taua_dev(k, ibn)*(1.-ssaa_dev(k, ibn)*&
504 & ff)
505  END IF
506  taerlyrd(k) = -(1.66*taua_devd(k, ibn)*exp(-(1.66*taua_dev(k&
507 & , ibn))))
508  taerlyr(k) = exp(-(1.66*taua_dev(k, ibn)))
509  END IF
510  END DO
511  END IF
512 !Compute the exponential terms (Eq. 8.21) at each layer due to
513 !water vapor line absorption when k-distribution is used
514  IF (.NOT.h2otable .AND. (.NOT.b10bnd)) THEN
515  CALL h2oexps_d(ibn, np, dh2o, dh2od, pa, dt, dtd, xkw, aw, bw, &
516 & pm, mw, exptbl(:, h2oexp%start:h2oexp%end), exptbld(:, &
517 & h2oexp%start:h2oexp%end))
518  ELSE
519  exptbld = 0.0_8
520  END IF
521 !compute the exponential terms (Eq. 4.24) at each layer due to
522 !water vapor continuum absorption.
523  ne = 0
524  IF (conbnd) THEN
525  ne = 1
526  IF (ibn .EQ. 3) ne = 3
527  CALL conexps_d(ibn, np, dcont, dcontd, xke, exptbl(:, conexp%&
528 & start:conexp%end), exptbld(:, conexp%start:conexp%end))
529  END IF
530  IF (trace) THEN
531 !compute the exponential terms at each layer due to n2o absorption
532  IF (n2obnd) CALL n2oexps_d(ibn, np, dn2o, pa, dt, dtd, exptbl(:&
533 & , n2oexp%start:n2oexp%end), exptbld(:, &
534 & n2oexp%start:n2oexp%end))
535 !compute the exponential terms at each layer due to ch4 absorption
536  IF (ch4bnd) CALL ch4exps_d(ibn, np, dch4, pa, dt, dtd, exptbl(:&
537 & , ch4exp%start:ch4exp%end), exptbld(:, &
538 & ch4exp%start:ch4exp%end))
539 !Compute the exponential terms due to co2 minor absorption
540  IF (combnd) CALL comexps_d(ibn, np, dco2, dt, dtd, exptbl(:, &
541 & comexp%start:comexp%end), exptbld(:, comexp&
542 & %start:comexp%end))
543 !Compute the exponential terms due to cfc11 absorption.
544 !The values of the parameters are given in Table 7.
545  IF (f11bnd) THEN
546  a1 = 1.26610e-3
547  b1 = 3.55940e-6
548  fk1 = 1.89736e+1
549  a2 = 8.19370e-4
550  b2 = 4.67810e-6
551  fk2 = 1.01487e+1
552  CALL cfcexps_d(ibn, np, a1, b1, fk1, a2, b2, fk2, df11, dt, &
553 & dtd, exptbl(:, f11exp%start:f11exp%end), exptbld(:, &
554 & f11exp%start:f11exp%end))
555  END IF
556 !Compute the exponential terms due to cfc12 absorption.
557  IF (f12bnd) THEN
558  a1 = 8.77370e-4
559  b1 = -5.88440e-6
560  fk1 = 1.58104e+1
561  a2 = 8.62000e-4
562  b2 = -4.22500e-6
563  fk2 = 3.70107e+1
564  CALL cfcexps_d(ibn, np, a1, b1, fk1, a2, b2, fk2, df12, dt, &
565 & dtd, exptbl(:, f12exp%start:f12exp%end), exptbld(:, &
566 & f12exp%start:f12exp%end))
567  END IF
568 !Compute the exponential terms due to cfc22 absorption.
569  IF (f22bnd) THEN
570  a1 = 9.65130e-4
571  b1 = 1.31280e-5
572  fk1 = 6.18536e+0
573  a2 = -3.00010e-5
574  b2 = 5.25010e-7
575  fk2 = 3.27912e+1
576  CALL cfcexps_d(ibn, np, a1, b1, fk1, a2, b2, fk2, df22, dt, &
577 & dtd, exptbl(:, f22exp%start:f22exp%end), exptbld(:, &
578 & f22exp%start:f22exp%end))
579  END IF
580 !Compute the exponential terms at each layer in band 10 due to
581 !h2o line and continuum, co2, and n2o absorption
582  IF (b10bnd) THEN
583  CALL b10exps_d(np, dh2o, dh2od, dcont, dcontd, dco2, dn2o, pa&
584 & , dt, dtd, h2oexp_tmp, h2oexp_tmpd, exptbl(:, conexp%&
585 & start:conexp%end), exptbld(:, conexp%start:conexp%end&
586 & ), co2exp_tmp, co2exp_tmpd, n2oexp_tmp, n2oexp_tmpd)
587  exptbld(:, h2oexp%start:h2oexp%end) = h2oexp_tmpd
588  exptbl(:, h2oexp%start:h2oexp%end) = h2oexp_tmp
589  exptbld(:, co2exp%start:co2exp%end) = co2exp_tmpd
590  exptbl(:, co2exp%start:co2exp%end) = co2exp_tmp
591  exptbld(:, n2oexp%start:n2oexp%end) = n2oexp_tmpd
592  exptbl(:, n2oexp%start:n2oexp%end) = n2oexp_tmp
593  END IF
594  END IF
595 !blayer(np+1) includes the effect of surface emissivity.
596 ! ALT: this was undefined, check with Max if 0.0 is good value
597  bud(0) = 0.0_8
598  bu(0) = 0.0
599  bdd(0) = blayerd(1)
600  bd(0) = blayer(1)
601  bud(np+1) = blayerd(np+1)
602  bu(np+1) = blayer(np+1)
603 !do-loop 1500 is for computing upward (bu) and downward (bd)
604 !Here, trant is the transmittance of the layer k2-1.
605  DO k2=1,np+1
606 !for h2o line transmission
607  IF (.NOT.h2otable) THEN
608  th2o = 1.0
609  th2od = 0.0_8
610  END IF
611 !for h2o continuum transmission
612  tcon = 1.0
613  x1 = 0.0
614  x2 = 0.0
615  x3 = 0.0
616  trant = 1.0
617  IF (h2otable) THEN
618 !Compute water vapor transmittance using table look-up.
619  IF (ibn .EQ. 1) THEN
620  trantd = 0.0_8
621  x3d = 0.0_8
622  x2d = 0.0_8
623  x1d = 0.0_8
624  CALL tablup_d(nx1, nh1, dh2o(k2-1), dh2od(k2-1), pa(k2-1), &
625 & dt(k2-1), dtd(k2-1), x1, x1d, x2, x2d, x3, x3d, w11&
626 & , p11, dwe, dpe, h11, h12, h13, trant, trantd)
627  ELSE
628  trantd = 0.0_8
629  x1d = 0.0_8
630  x2d = 0.0_8
631  x3d = 0.0_8
632  END IF
633  IF (ibn .EQ. 2) CALL tablup_d(nx1, nh1, dh2o(k2-1), dh2od(k2-1&
634 & ), pa(k2-1), dt(k2-1), dtd(k2-1), x1, &
635 & x1d, x2, x2d, x3, x3d, w11, p11, dwe, &
636 & dpe, h21, h22, h23, trant, trantd)
637  IF (ibn .EQ. 8) CALL tablup_d(nx1, nh1, dh2o(k2-1), dh2od(k2-1&
638 & ), pa(k2-1), dt(k2-1), dtd(k2-1), x1, &
639 & x1d, x2, x2d, x3, x3d, w11, p11, dwe, &
640 & dpe, h81, h82, h83, trant, trantd)
641 !for water vapor continuum absorption
642  IF (conbnd) THEN
643 ! Only the first exp
644  tcond = 0.0_8
645  tcond(1) = tcon(1)*exptbld(k2-1, conexp%start)
646  tcon(1) = tcon(1)*exptbl(k2-1, conexp%start)
647  trantd = trantd*tcon(1) + trant*tcond(1)
648  trant = trant*tcon(1)
649  END IF
650  ELSE IF (.NOT.b10bnd) THEN
651 !compute water vapor transmittance using k-distribution
652  tcond = 0.0_8
653  CALL h2okdis_d(ibn, np, k2 - 1, fkw, gkw, ne, exptbl(:, h2oexp&
654 & %start:h2oexp%end), exptbld(:, h2oexp%start:h2oexp%&
655 & end), exptbl(:, conexp%start:conexp%end), exptbld(:, &
656 & conexp%start:conexp%end), th2o, th2od, tcon, tcond, &
657 & trant, trantd)
658  x1d = 0.0_8
659  x2d = 0.0_8
660  x3d = 0.0_8
661  ELSE
662  trantd = 0.0_8
663  x1d = 0.0_8
664  x2d = 0.0_8
665  x3d = 0.0_8
666  END IF
667  IF (co2bnd) THEN
668 !Compute co2 transmittance using table look-up method
669  dco2d = 0.0_8
670  CALL tablup_d(nx1, nc1, dco2(k2-1), dco2d(k2-1), pa(k2-1), dt(&
671 & k2-1), dtd(k2-1), x1, x1d, x2, x2d, x3, x3d, w12, p12&
672 & , dwe, dpe, c1, c2, c3, trant, trantd)
673  END IF
674 !Always use table look-up to compute o3 transmittance.
675  IF (oznbnd) CALL tablup_d(nx1, no1, do3(k2-1), do3d(k2-1), pa(k2&
676 & -1), dt(k2-1), dtd(k2-1), x1, x1d, x2, x2d, &
677 & x3, x3d, w13, p13, dwe, dpe, oo1, oo2, oo3, &
678 & trant, trantd)
679 !include aerosol effect
680  IF (do_aerosol) THEN
681  trantd = trantd*taerlyr(k2-1) + trant*taerlyrd(k2-1)
682  trant = trant*taerlyr(k2-1)
683  END IF
684 !Compute upward (bu) and downward (bd) emission of the layer k2-1
685  xxd = (1.-enn(k2-1))*trantd - ennd(k2-1)*trant
686  xx = (1.-enn(k2-1))*trant
687  IF (0.9999 .GT. xx) THEN
688  yyd = xxd
689  yy = xx
690  ELSE
691  yy = 0.9999
692  yyd = 0.0_8
693  END IF
694  IF (0.00001 .LT. yy) THEN
695  yy = yy
696  ELSE
697  yy = 0.00001
698  yyd = 0.0_8
699  END IF
700  xxd = ((bleveld(k2-1)-bleveld(k2))*log(yy)-(blevel(k2-1)-blevel(&
701 & k2))*yyd/yy)/log(yy)**2
702  xx = (blevel(k2-1)-blevel(k2))/log(yy)
703  bdd(k2-1) = ((bleveld(k2)-bleveld(k2-1)*yy-blevel(k2-1)*yyd)*(&
704 & 1.0-yy)+(blevel(k2)-blevel(k2-1)*yy)*yyd)/(1.0-yy)**2 - xxd
705  bd(k2-1) = (blevel(k2)-blevel(k2-1)*yy)/(1.0-yy) - xx
706  bud(k2-1) = bleveld(k2-1) + bleveld(k2) - bdd(k2-1)
707  bu(k2-1) = blevel(k2-1) + blevel(k2) - bd(k2-1)
708  END DO
709 !Initialize fluxes
710  flxu = 0.0
711  flxd = 0.0
712  flxud = 0.0_8
713  flxdd = 0.0_8
714 !Compute upward and downward fluxes for each spectral band, ibn.
715  DO k1=0,np
716 !initialization
717  cldlw = 0.0
718  cldmd = 0.0
719  cldhi = 0.0
720  tranal = 1.0
721 !for h2o line transmission
722  IF (.NOT.h2otable) THEN
723  th2o = 1.0
724  th2od = 0.0_8
725  END IF
726 !for h2o continuum transmission
727  tcon = 1.0
728  IF (trace) THEN
729 !Add trace gases contribution
730  IF (n2obnd) THEN
731 !n2o
732  tn2o = 1.0
733  tn2od = 0.0_8
734  END IF
735  IF (ch4bnd) THEN
736 !ch4
737  tch4 = 1.0
738  tch4d = 0.0_8
739  END IF
740  IF (combnd) THEN
741 !co2-minor
742  tcom = 1.0
743  tcomd = 0.0_8
744  END IF
745  IF (f11bnd) THEN
746 !cfc-11
747  tf11 = 1.0
748  tf11d = 0.0_8
749  END IF
750  IF (f12bnd) THEN
751 !cfc-12
752  tf12 = 1.0
753  tf12d = 0.0_8
754  END IF
755  IF (f22bnd) THEN
756 !cfc-22
757  tf22 = 1.0
758  tf22d = 0.0_8
759  END IF
760  IF (b10bnd) THEN
761 !
762  th2o = 1.0
763  tco2 = 1.0
764  tcond(1) = 0.0_8
765  tcon(1) = 1.0
766  tn2o = 1.0
767  tn2od = 0.0_8
768  th2od = 0.0_8
769  tco2d = 0.0_8
770  END IF
771  END IF
772  x1 = 0.0
773  x2 = 0.0
774  x3 = 0.0
775  cldhid = 0.0_8
776  tcond = 0.0_8
777  tranald = 0.0_8
778  x1d = 0.0_8
779  x2d = 0.0_8
780  x3d = 0.0_8
781  cldlwd = 0.0_8
782  cldmdd = 0.0_8
783  DO k2=k1+1,np+1
784  trant = 1.0
785  fclr = 1.0
786  IF (h2otable) THEN
787 !Compute water vapor transmittance using table look-up.
788  IF (ibn .EQ. 1) THEN
789  trantd = 0.0_8
790  CALL tablup_d(nx1, nh1, dh2o(k2-1), dh2od(k2-1), pa(k2-1)&
791 & , dt(k2-1), dtd(k2-1), x1, x1d, x2, x2d, x3, x3d, &
792 & w11, p11, dwe, dpe, h11, h12, h13, trant, trantd)
793  ELSE
794  trantd = 0.0_8
795  END IF
796  IF (ibn .EQ. 2) CALL tablup_d(nx1, nh1, dh2o(k2-1), dh2od(k2&
797 & -1), pa(k2-1), dt(k2-1), dtd(k2-1), &
798 & x1, x1d, x2, x2d, x3, x3d, w11, p11&
799 & , dwe, dpe, h21, h22, h23, trant, &
800 & trantd)
801  IF (ibn .EQ. 8) CALL tablup_d(nx1, nh1, dh2o(k2-1), dh2od(k2&
802 & -1), pa(k2-1), dt(k2-1), dtd(k2-1), &
803 & x1, x1d, x2, x2d, x3, x3d, w11, p11&
804 & , dwe, dpe, h81, h82, h83, trant, &
805 & trantd)
806  IF (conbnd) THEN
807 ! Only the first exp
808  tcond(1) = tcond(1)*exptbl(k2-1, conexp%start) + tcon(1)*&
809 & exptbld(k2-1, conexp%start)
810  tcon(1) = tcon(1)*exptbl(k2-1, conexp%start)
811  trantd = trantd*tcon(1) + trant*tcond(1)
812  trant = trant*tcon(1)
813  END IF
814  ELSE IF (.NOT.b10bnd) THEN
815 !Compute water vapor transmittance using k-distribution
816  CALL h2okdis_d(ibn, np, k2 - 1, fkw, gkw, ne, exptbl(:, &
817 & h2oexp%start:h2oexp%end), exptbld(:, h2oexp%start:&
818 & h2oexp%end), exptbl(:, conexp%start:conexp%end), &
819 & exptbld(:, conexp%start:conexp%end), th2o, th2od, &
820 & tcon, tcond, trant, trantd)
821  ELSE
822  trantd = 0.0_8
823  END IF
824  IF (co2bnd) THEN
825 !Compute co2 transmittance using table look-up method.
826  dco2d = 0.0_8
827  CALL tablup_d(nx1, nc1, dco2(k2-1), dco2d(k2-1), pa(k2-1), &
828 & dt(k2-1), dtd(k2-1), x1, x1d, x2, x2d, x3, x3d, w12&
829 & , p12, dwe, dpe, c1, c2, c3, trant, trantd)
830  END IF
831  IF (oznbnd) CALL tablup_d(nx1, no1, do3(k2-1), do3d(k2-1), pa(&
832 & k2-1), dt(k2-1), dtd(k2-1), x1, x1d, x2, &
833 & x2d, x3, x3d, w13, p13, dwe, dpe, oo1, oo2&
834 & , oo3, trant, trantd)
835 !Always use table look-up to compute o3 transmittanc
836  IF (trace) THEN
837 !Add trace gas effects
838  IF (n2obnd) CALL n2okdis_d(ibn, np, k2 - 1, exptbl(:, n2oexp&
839 & %start:n2oexp%end), exptbld(:, n2oexp%&
840 & start:n2oexp%end), tn2o, tn2od, trant, &
841 & trantd)
842 !n2o
843  IF (ch4bnd) CALL ch4kdis_d(ibn, np, k2 - 1, exptbl(:, ch4exp&
844 & %start:ch4exp%end), exptbld(:, ch4exp%&
845 & start:ch4exp%end), tch4, tch4d, trant, &
846 & trantd)
847 !ch4
848  IF (combnd) CALL comkdis_d(ibn, np, k2 - 1, exptbl(:, comexp&
849 & %start:comexp%end), exptbld(:, comexp%&
850 & start:comexp%end), tcom, tcomd, trant, &
851 & trantd)
852 !co2-minor
853  IF (f11bnd) CALL cfckdis_d(np, k2 - 1, exptbl(:, f11exp%&
854 & start:f11exp%end), exptbld(:, f11exp%&
855 & start:f11exp%end), tf11, tf11d, trant, &
856 & trantd)
857 !cfc11
858  IF (f12bnd) CALL cfckdis_d(np, k2 - 1, exptbl(:, f12exp%&
859 & start:f12exp%end), exptbld(:, f12exp%&
860 & start:f12exp%end), tf12, tf12d, trant, &
861 & trantd)
862 !cfc12
863  IF (f22bnd) CALL cfckdis_d(np, k2 - 1, exptbl(:, f22exp%&
864 & start:f22exp%end), exptbld(:, f22exp%&
865 & start:f22exp%end), tf22, tf22d, trant, &
866 & trantd)
867 !CFC22
868  IF (b10bnd) CALL b10kdis_d(np, k2 - 1, exptbl(:, h2oexp%&
869 & start:h2oexp%end), exptbld(:, h2oexp%&
870 & start:h2oexp%end), exptbl(:, conexp%&
871 & start:conexp%end), exptbld(:, conexp%&
872 & start:conexp%end), exptbl(:, co2exp%&
873 & start:co2exp%end), exptbld(:, co2exp%&
874 & start:co2exp%end), exptbl(:, n2oexp%&
875 & start:n2oexp%end), exptbld(:, n2oexp%&
876 & start:n2oexp%end), th2o, th2od, tcon, &
877 & tcond, tco2, tco2d, tn2o, tn2od, trant&
878 & , trantd)
879  END IF
880  IF (do_aerosol) THEN
881  tranald = tranald*taerlyr(k2-1) + tranal*taerlyrd(k2-1)
882  tranal = tranal*taerlyr(k2-1)
883  trantd = trantd*tranal + trant*tranald
884  trant = trant*tranal
885  END IF
886  IF (enn(k2-1) .GE. 0.001) CALL cldovlp_d(np, k1, k2, ict, icb&
887 & , icx, ncld, enn, ennd, &
888 & tcldlyr, tcldlyrd, cldhi, &
889 & cldhid, cldmd, cldmdd, &
890 & cldlw, cldlwd)
891  fclrd = (-(cldhid*(1.0-cldmd))-(1.0-cldhi)*cldmdd)*(1.0-cldlw)&
892 & - (1.0-cldhi)*(1.0-cldmd)*cldlwd
893  fclr = (1.0-cldhi)*(1.0-cldmd)*(1.0-cldlw)
894  IF (k2 .EQ. k1 + 1 .AND. ibn .NE. 10) THEN
895  flxud(k1) = flxud(k1) - bud(k1)
896  flxu(k1) = flxu(k1) - bu(k1)
897  flxdd(k2) = flxdd(k2) + bdd(k1)
898  flxd(k2) = flxd(k2) + bd(k1)
899  END IF
900  xxd = trantd*(bu(k2-1)-bu(k2)) + trant*(bud(k2-1)-bud(k2))
901  xx = trant*(bu(k2-1)-bu(k2))
902  flxud(k1) = flxud(k1) + xxd*fclr + xx*fclrd
903  flxu(k1) = flxu(k1) + xx*fclr
904  IF (k1 .EQ. 0) THEN
905 !mjs bd(-1) is not defined
906  xxd = -(trantd*bd(k1)+trant*bdd(k1))
907  xx = -(trant*bd(k1))
908  ELSE
909  xxd = trantd*(bd(k1-1)-bd(k1)) + trant*(bdd(k1-1)-bdd(k1))
910  xx = trant*(bd(k1-1)-bd(k1))
911  END IF
912  flxdd(k2) = flxdd(k2) + xxd*fclr + xx*fclrd
913  flxd(k2) = flxd(k2) + xx*fclr
914  END DO
915  transfcd(k1) = trantd*fclr + trant*fclrd
916  transfc(k1) = trant*fclr
917  IF (k1 .GT. 0) THEN
918  dfdts_devd(k1) = dfdts_devd(k1) - dbs*transfcd(k1)
919  dfdts_dev(k1) = dfdts_dev(k1) - dbs*transfc(k1)
920  END IF
921  END DO
922  IF (.NOT.b10bnd) THEN
923 !surface emission
924  flxud(np+1) = -blayerd(np+1)
925  flxu(np+1) = -blayer(np+1)
926  dfdts_dev(np+1) = dfdts_dev(np+1) - dbs
927  DO k=1,np+1
928  flxud(k) = flxud(k) - rflxs*(flxdd(np+1)*transfc(k)+flxd(np+1)&
929 & *transfcd(k))
930  flxu(k) = flxu(k) - flxd(np+1)*transfc(k)*rflxs
931  END DO
932  END IF
933 !Summation of fluxes over spectral bands
934  DO k=1,np+1
935  flxu_devd(k) = flxu_devd(k) + flxud(k)
936  flxu_dev(k) = flxu_dev(k) + flxu(k)
937  flxd_devd(k) = flxd_devd(k) + flxdd(k)
938  flxd_dev(k) = flxd_dev(k) + flxd(k)
939  END DO
940  END IF
941  END DO
942  GOTO 110
943  100 RETURN
944  110 CONTINUE
945 END SUBROUTINE irrad_d
946 
947 ! Differentiation of planck in forward (tangent) mode:
948 ! variations of useful results: xlayer
949 ! with respect to varying inputs: t
950 !***********************************************************************
951 SUBROUTINE planck_d(ibn, cb, t, td, xlayer, xlayerd)
952  IMPLICIT NONE
953 ! spectral band index
954  INTEGER, INTENT(IN) :: ibn
955 ! Planck table coefficients
956  REAL*8, INTENT(IN) :: cb(6, 10)
957 ! temperature (K)
958  REAL*8, INTENT(IN) :: t
959  REAL*8, INTENT(IN) :: td
960 ! planck flux (w/m2)
961  REAL*8, INTENT(OUT) :: xlayer
962  REAL*8, INTENT(OUT) :: xlayerd
963  xlayerd = td*(t*(t*(t*(t*cb(6, ibn)+cb(5, ibn))+cb(4, ibn))+cb(3, ibn)&
964 & )+cb(2, ibn)) + t*(td*(t*(t*(t*cb(6, ibn)+cb(5, ibn))+cb(4, ibn))+cb&
965 & (3, ibn))+t*(td*(t*(t*cb(6, ibn)+cb(5, ibn))+cb(4, ibn))+t*(td*(t*cb&
966 & (6, ibn)+cb(5, ibn))+t*cb(6, ibn)*td)))
967  xlayer = t*(t*(t*(t*(t*cb(6, ibn)+cb(5, ibn))+cb(4, ibn))+cb(3, ibn))+&
968 & cb(2, ibn)) + cb(1, ibn)
969 END SUBROUTINE planck_d
970 
971 ! Differentiation of h2oexps in forward (tangent) mode:
972 ! variations of useful results: h2oexp
973 ! with respect to varying inputs: dh2o dt
974 !**********************************************************************
975 SUBROUTINE h2oexps_d(ib, np, dh2o, dh2od, pa, dt, dtd, xkw, aw, bw, pm, &
976 & mw, h2oexp, h2oexpd)
977  IMPLICIT NONE
978  INTEGER :: ib, np, ik, k
979 !---- input parameters ------
980  REAL*8 :: dh2o(0:np), pa(0:np), dt(0:np)
981  REAL*8 :: dh2od(0:np), dtd(0:np)
982 !---- output parameters -----
983  REAL*8 :: h2oexp(0:np, 6)
984  REAL*8 :: h2oexpd(0:np, 6)
985 !---- static data -----
986  INTEGER :: mw(9)
987  REAL*8 :: xkw(9), aw(9), bw(9), pm(9)
988 !---- temporary arrays -----
989  REAL*8 :: xh
990  REAL*8 :: xhd
991  INTRINSIC exp
992  REAL*8 :: pwx1
993  REAL*8 :: pwr1
994  h2oexpd = 0.0_8
995 !**********************************************************************
996 ! note that the 3 sub-bands in band 3 use the same set of xkw, aw,
997 ! and bw, therefore, h2oexp for these sub-bands are identical.
998 !**********************************************************************
999  DO k=0,np
1000 !-----xh is the scaled water vapor amount for line absorption
1001 ! computed from Eq. (4.4).
1002  pwx1 = pa(k)/500.
1003  pwr1 = pwx1**pm(ib)
1004  xhd = pwr1*(dh2od(k)*(1.+(aw(ib)+bw(ib)*dt(k))*dt(k))+dh2o(k)*(bw(ib&
1005 & )*dtd(k)*dt(k)+(aw(ib)+bw(ib)*dt(k))*dtd(k)))
1006  xh = dh2o(k)*pwr1*(1.+(aw(ib)+bw(ib)*dt(k))*dt(k))
1007 !-----h2oexp is the water vapor transmittance of the layer k
1008 ! due to line absorption
1009  h2oexpd(k, 1) = -(xkw(ib)*xhd*exp(-(xh*xkw(ib))))
1010  h2oexp(k, 1) = exp(-(xh*xkw(ib)))
1011 !-----compute transmittances from Eq. (8.22)
1012  DO ik=2,6
1013  IF (mw(ib) .EQ. 6) THEN
1014  xhd = h2oexpd(k, ik-1)*h2oexp(k, ik-1) + h2oexp(k, ik-1)*h2oexpd&
1015 & (k, ik-1)
1016  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)
1017  h2oexpd(k, ik) = (xhd*xh+xh*xhd)*xh + xh**2*xhd
1018  h2oexp(k, ik) = xh*xh*xh
1019  ELSE IF (mw(ib) .EQ. 8) THEN
1020  xhd = h2oexpd(k, ik-1)*h2oexp(k, ik-1) + h2oexp(k, ik-1)*h2oexpd&
1021 & (k, ik-1)
1022  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)
1023  xhd = xhd*xh + xh*xhd
1024  xh = xh*xh
1025  h2oexpd(k, ik) = xhd*xh + xh*xhd
1026  h2oexp(k, ik) = xh*xh
1027  ELSE IF (mw(ib) .EQ. 9) THEN
1028  xhd = (h2oexpd(k, ik-1)*h2oexp(k, ik-1)+h2oexp(k, ik-1)*h2oexpd(&
1029 & k, ik-1))*h2oexp(k, ik-1) + h2oexp(k, ik-1)**2*h2oexpd(k, ik-1&
1030 & )
1031  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)*h2oexp(k, ik-1)
1032  h2oexpd(k, ik) = (xhd*xh+xh*xhd)*xh + xh**2*xhd
1033  h2oexp(k, ik) = xh*xh*xh
1034  ELSE
1035  xhd = h2oexpd(k, ik-1)*h2oexp(k, ik-1) + h2oexp(k, ik-1)*h2oexpd&
1036 & (k, ik-1)
1037  xh = h2oexp(k, ik-1)*h2oexp(k, ik-1)
1038  xhd = xhd*xh + xh*xhd
1039  xh = xh*xh
1040  xhd = xhd*xh + xh*xhd
1041  xh = xh*xh
1042  h2oexpd(k, ik) = xhd*xh + xh*xhd
1043  h2oexp(k, ik) = xh*xh
1044  END IF
1045  END DO
1046  END DO
1047 END SUBROUTINE h2oexps_d
1048 
1049 ! Differentiation of conexps in forward (tangent) mode:
1050 ! variations of useful results: conexp
1051 ! with respect to varying inputs: dcont conexp
1052 !**********************************************************************
1053 SUBROUTINE conexps_d(ib, np, dcont, dcontd, xke, conexp, conexpd)
1054  IMPLICIT NONE
1055  INTEGER :: ib, np, k
1056 !---- input parameters ------
1057  REAL*8 :: dcont(0:np)
1058  REAL*8 :: dcontd(0:np)
1059 !---- updated parameters -----
1060  REAL*8 :: conexp(0:np, 3)
1061  REAL*8 :: conexpd(0:np, 3)
1062 !---- static data -----
1063  REAL*8 :: xke(9)
1064  INTRINSIC exp
1065 !****************************************************************
1066  DO k=0,np
1067  conexpd(k, 1) = -(xke(ib)*dcontd(k)*exp(-(dcont(k)*xke(ib))))
1068  conexp(k, 1) = exp(-(dcont(k)*xke(ib)))
1069 !-----The absorption coefficients for sub-bands 3b and 3a are, respectively,
1070 ! two and four times the absorption coefficient for sub-band 3c (Table 9).
1071 ! Note that conexp(3) is for sub-band 3a.
1072  IF (ib .EQ. 3) THEN
1073  conexpd(k, 2) = conexpd(k, 1)*conexp(k, 1) + conexp(k, 1)*conexpd(&
1074 & k, 1)
1075  conexp(k, 2) = conexp(k, 1)*conexp(k, 1)
1076  conexpd(k, 3) = conexpd(k, 2)*conexp(k, 2) + conexp(k, 2)*conexpd(&
1077 & k, 2)
1078  conexp(k, 3) = conexp(k, 2)*conexp(k, 2)
1079  END IF
1080  END DO
1081 END SUBROUTINE conexps_d
1082 
1083 ! Differentiation of n2oexps in forward (tangent) mode:
1084 ! variations of useful results: n2oexp
1085 ! with respect to varying inputs: dt n2oexp
1086 !**********************************************************************
1087 SUBROUTINE n2oexps_d(ib, np, dn2o, pa, dt, dtd, n2oexp, n2oexpd)
1088  IMPLICIT NONE
1089  INTEGER :: ib, k, np
1090 !---- input parameters -----
1091  REAL*8 :: dn2o(0:np), pa(0:np), dt(0:np)
1092  REAL*8 :: dtd(0:np)
1093 !---- output parameters -----
1094  REAL*8 :: n2oexp(0:np, 4)
1095  REAL*8 :: n2oexpd(0:np, 4)
1096 !---- temporary arrays -----
1097  REAL*8 :: xc, xc1, xc2
1098  REAL*8 :: xcd, xc1d, xc2d
1099  INTRINSIC exp
1100 !-----Scaling and absorption data are given in Table 5.
1101 ! Transmittances are computed using Eqs. (8.21) and (8.22).
1102  DO k=0,np
1103 !-----four exponential by powers of 21 for band 6.
1104  IF (ib .EQ. 6) THEN
1105  xcd = dn2o(k)*(4.3750e-6*dtd(k)*dt(k)+(1.9297e-3+4.3750e-6*dt(k))*&
1106 & dtd(k))
1107  xc = dn2o(k)*(1.+(1.9297e-3+4.3750e-6*dt(k))*dt(k))
1108  n2oexpd(k, 1) = -(6.31582e-2*xcd*exp(-(xc*6.31582e-2)))
1109  n2oexp(k, 1) = exp(-(xc*6.31582e-2))
1110  xcd = (n2oexpd(k, 1)*n2oexp(k, 1)+n2oexp(k, 1)*n2oexpd(k, 1))*&
1111 & n2oexp(k, 1) + n2oexp(k, 1)**2*n2oexpd(k, 1)
1112  xc = n2oexp(k, 1)*n2oexp(k, 1)*n2oexp(k, 1)
1113  xc1d = xcd*xc + xc*xcd
1114  xc1 = xc*xc
1115  xc2d = xc1d*xc1 + xc1*xc1d
1116  xc2 = xc1*xc1
1117  n2oexpd(k, 2) = (xcd*xc1+xc*xc1d)*xc2 + xc*xc1*xc2d
1118  n2oexp(k, 2) = xc*xc1*xc2
1119 !-----four exponential by powers of 8 for band 7
1120  ELSE
1121  xcd = dn2o(k)*(pa(k)/500.0)**0.48*(7.4838e-6*dtd(k)*dt(k)+(&
1122 & 1.3804e-3+7.4838e-6*dt(k))*dtd(k))
1123  xc = dn2o(k)*(pa(k)/500.0)**0.48*(1.+(1.3804e-3+7.4838e-6*dt(k))*&
1124 & dt(k))
1125  n2oexpd(k, 1) = -(5.35779e-2*xcd*exp(-(xc*5.35779e-2)))
1126  n2oexp(k, 1) = exp(-(xc*5.35779e-2))
1127  xcd = n2oexpd(k, 1)*n2oexp(k, 1) + n2oexp(k, 1)*n2oexpd(k, 1)
1128  xc = n2oexp(k, 1)*n2oexp(k, 1)
1129  xcd = xcd*xc + xc*xcd
1130  xc = xc*xc
1131  n2oexpd(k, 2) = xcd*xc + xc*xcd
1132  n2oexp(k, 2) = xc*xc
1133  xcd = n2oexpd(k, 2)*n2oexp(k, 2) + n2oexp(k, 2)*n2oexpd(k, 2)
1134  xc = n2oexp(k, 2)*n2oexp(k, 2)
1135  xcd = xcd*xc + xc*xcd
1136  xc = xc*xc
1137  n2oexpd(k, 3) = xcd*xc + xc*xcd
1138  n2oexp(k, 3) = xc*xc
1139  xcd = n2oexpd(k, 3)*n2oexp(k, 3) + n2oexp(k, 3)*n2oexpd(k, 3)
1140  xc = n2oexp(k, 3)*n2oexp(k, 3)
1141  xcd = xcd*xc + xc*xcd
1142  xc = xc*xc
1143  n2oexpd(k, 4) = xcd*xc + xc*xcd
1144  n2oexp(k, 4) = xc*xc
1145  END IF
1146  END DO
1147 END SUBROUTINE n2oexps_d
1148 
1149 ! Differentiation of ch4exps in forward (tangent) mode:
1150 ! variations of useful results: ch4exp
1151 ! with respect to varying inputs: dt ch4exp
1152 !**********************************************************************
1153 SUBROUTINE ch4exps_d(ib, np, dch4, pa, dt, dtd, ch4exp, ch4expd)
1154  IMPLICIT NONE
1155  INTEGER :: ib, np, k
1156 !---- input parameters -----
1157  REAL*8 :: dch4(0:np), pa(0:np), dt(0:np)
1158  REAL*8 :: dtd(0:np)
1159 !---- output parameters -----
1160  REAL*8 :: ch4exp(0:np, 4)
1161  REAL*8 :: ch4expd(0:np, 4)
1162 !---- temporary arrays -----
1163  REAL*8 :: xc
1164  REAL*8 :: xcd
1165  INTRINSIC exp
1166 !***** Scaling and absorption data are given in Table 5 *****
1167  DO k=0,np
1168 !-----four exponentials for band 6
1169  IF (ib .EQ. 6) THEN
1170  xcd = dch4(k)*(1.5826e-4*dtd(k)*dt(k)+(1.7007e-2+1.5826e-4*dt(k))*&
1171 & dtd(k))
1172  xc = dch4(k)*(1.+(1.7007e-2+1.5826e-4*dt(k))*dt(k))
1173  ch4expd(k, 1) = -(5.80708e-3*xcd*exp(-(xc*5.80708e-3)))
1174  ch4exp(k, 1) = exp(-(xc*5.80708e-3))
1175 !-----four exponentials by powers of 12 for band 7
1176  ELSE
1177  xcd = dch4(k)*(pa(k)/500.0)**0.65*((5.9590e-4-2.2931e-6*dt(k))*dtd&
1178 & (k)-2.2931e-6*dtd(k)*dt(k))
1179  xc = dch4(k)*(pa(k)/500.0)**0.65*(1.+(5.9590e-4-2.2931e-6*dt(k))*&
1180 & dt(k))
1181  ch4expd(k, 1) = -(6.29247e-2*xcd*exp(-(xc*6.29247e-2)))
1182  ch4exp(k, 1) = exp(-(xc*6.29247e-2))
1183  xcd = (ch4expd(k, 1)*ch4exp(k, 1)+ch4exp(k, 1)*ch4expd(k, 1))*&
1184 & ch4exp(k, 1) + ch4exp(k, 1)**2*ch4expd(k, 1)
1185  xc = ch4exp(k, 1)*ch4exp(k, 1)*ch4exp(k, 1)
1186  xcd = xcd*xc + xc*xcd
1187  xc = xc*xc
1188  ch4expd(k, 2) = xcd*xc + xc*xcd
1189  ch4exp(k, 2) = xc*xc
1190  xcd = (ch4expd(k, 2)*ch4exp(k, 2)+ch4exp(k, 2)*ch4expd(k, 2))*&
1191 & ch4exp(k, 2) + ch4exp(k, 2)**2*ch4expd(k, 2)
1192  xc = ch4exp(k, 2)*ch4exp(k, 2)*ch4exp(k, 2)
1193  xcd = xcd*xc + xc*xcd
1194  xc = xc*xc
1195  ch4expd(k, 3) = xcd*xc + xc*xcd
1196  ch4exp(k, 3) = xc*xc
1197  xcd = (ch4expd(k, 3)*ch4exp(k, 3)+ch4exp(k, 3)*ch4expd(k, 3))*&
1198 & ch4exp(k, 3) + ch4exp(k, 3)**2*ch4expd(k, 3)
1199  xc = ch4exp(k, 3)*ch4exp(k, 3)*ch4exp(k, 3)
1200  xcd = xcd*xc + xc*xcd
1201  xc = xc*xc
1202  ch4expd(k, 4) = xcd*xc + xc*xcd
1203  ch4exp(k, 4) = xc*xc
1204  END IF
1205  END DO
1206 END SUBROUTINE ch4exps_d
1207 
1208 ! Differentiation of comexps in forward (tangent) mode:
1209 ! variations of useful results: comexp
1210 ! with respect to varying inputs: dt comexp
1211 !**********************************************************************
1212 SUBROUTINE comexps_d(ib, np, dcom, dt, dtd, comexp, comexpd)
1213  IMPLICIT NONE
1214  INTEGER :: ib, ik, np, k
1215 !---- input parameters -----
1216  REAL*8 :: dcom(0:np), dt(0:np)
1217  REAL*8 :: dtd(0:np)
1218 !---- output parameters -----
1219  REAL*8 :: comexp(0:np, 6)
1220  REAL*8 :: comexpd(0:np, 6)
1221 !---- temporary arrays -----
1222  REAL*8 :: xc
1223  REAL*8 :: xcd
1224  INTRINSIC exp
1225  xcd = 0.0_8
1226 !***** Scaling and absorpton data are given in Table 6 *****
1227  DO k=0,np
1228  IF (ib .EQ. 4) THEN
1229  xcd = dcom(k)*(4.0447e-4*dtd(k)*dt(k)+(3.5775e-2+4.0447e-4*dt(k))*&
1230 & dtd(k))
1231  xc = dcom(k)*(1.+(3.5775e-2+4.0447e-4*dt(k))*dt(k))
1232  END IF
1233  IF (ib .EQ. 5) THEN
1234  xcd = dcom(k)*(3.7401e-4*dtd(k)*dt(k)+(3.4268e-2+3.7401e-4*dt(k))*&
1235 & dtd(k))
1236  xc = dcom(k)*(1.+(3.4268e-2+3.7401e-4*dt(k))*dt(k))
1237  END IF
1238  comexpd(k, 1) = -(1.922e-7*xcd*exp(-(xc*1.922e-7)))
1239  comexp(k, 1) = exp(-(xc*1.922e-7))
1240  DO ik=2,6
1241  xcd = comexpd(k, ik-1)*comexp(k, ik-1) + comexp(k, ik-1)*comexpd(k&
1242 & , ik-1)
1243  xc = comexp(k, ik-1)*comexp(k, ik-1)
1244  xcd = xcd*xc + xc*xcd
1245  xc = xc*xc
1246  comexpd(k, ik) = xcd*comexp(k, ik-1) + xc*comexpd(k, ik-1)
1247  comexp(k, ik) = xc*comexp(k, ik-1)
1248  END DO
1249  END DO
1250 END SUBROUTINE comexps_d
1251 
1252 ! Differentiation of cfcexps in forward (tangent) mode:
1253 ! variations of useful results: cfcexp
1254 ! with respect to varying inputs: dt cfcexp
1255 !**********************************************************************
1256 SUBROUTINE cfcexps_d(ib, np, a1, b1, fk1, a2, b2, fk2, dcfc, dt, dtd, &
1257 & cfcexp, cfcexpd)
1258  IMPLICIT NONE
1259  INTEGER :: ib, np, k
1260 !---- input parameters -----
1261  REAL*8 :: dcfc(0:np), dt(0:np)
1262  REAL*8 :: dtd(0:np)
1263 !---- output parameters -----
1264  REAL*8 :: cfcexp(0:np)
1265  REAL*8 :: cfcexpd(0:np)
1266 !---- static data -----
1267  REAL*8 :: a1, b1, fk1, a2, b2, fk2
1268 !---- temporary arrays -----
1269  REAL*8 :: xf
1270  REAL*8 :: xfd
1271  INTRINSIC exp
1272 !**********************************************************************
1273  DO k=0,np
1274 !-----compute the scaled cfc amount (xf) and exponential (cfcexp)
1275  IF (ib .EQ. 4) THEN
1276  xfd = dcfc(k)*(b1*dtd(k)*dt(k)+(a1+b1*dt(k))*dtd(k))
1277  xf = dcfc(k)*(1.+(a1+b1*dt(k))*dt(k))
1278  cfcexpd(k) = -(fk1*xfd*exp(-(xf*fk1)))
1279  cfcexp(k) = exp(-(xf*fk1))
1280  ELSE
1281  xfd = dcfc(k)*(b2*dtd(k)*dt(k)+(a2+b2*dt(k))*dtd(k))
1282  xf = dcfc(k)*(1.+(a2+b2*dt(k))*dt(k))
1283  cfcexpd(k) = -(fk2*xfd*exp(-(xf*fk2)))
1284  cfcexp(k) = exp(-(xf*fk2))
1285  END IF
1286  END DO
1287 END SUBROUTINE cfcexps_d
1288 
1289 ! Differentiation of b10exps in forward (tangent) mode:
1290 ! variations of useful results: co2exp h2oexp n2oexp conexp
1291 ! with respect to varying inputs: dh2o dcont dt co2exp h2oexp
1292 ! n2oexp conexp
1293 !**********************************************************************
1294 SUBROUTINE b10exps_d(np, dh2o, dh2od, dcont, dcontd, dco2, dn2o, pa, dt&
1295 & , dtd, h2oexp, h2oexpd, conexp, conexpd, co2exp, co2expd, n2oexp, &
1296 & n2oexpd)
1297  IMPLICIT NONE
1298  INTEGER :: np, k
1299 !---- input parameters -----
1300  REAL*8 :: dh2o(0:np), dcont(0:np), dn2o(0:np)
1301  REAL*8 :: dh2od(0:np), dcontd(0:np)
1302  REAL*8 :: dco2(0:np), pa(0:np), dt(0:np)
1303  REAL*8 :: dtd(0:np)
1304 !---- output parameters -----
1305  REAL*8 :: h2oexp(0:np, 5), conexp(0:np), co2exp(0:np, 6), n2oexp(0:np&
1306 & , 2)
1307  REAL*8 :: h2oexpd(0:np, 5), conexpd(0:np), co2expd(0:np, 6), n2oexpd(0&
1308 & :np, 2)
1309 !---- temporary arrays -----
1310  REAL*8 :: xx, xx1, xx2, xx3
1311  REAL*8 :: xxd, xx1d, xx2d, xx3d
1312  INTRINSIC exp
1313 !**********************************************************************
1314  DO k=0,np
1315 !-----Compute scaled h2o-line amount for Band 10 (Eq. 4.4 and Table 3).
1316  xxd = pa(k)*(dh2od(k)*(1.+(0.0149+6.20e-5*dt(k))*dt(k))+dh2o(k)*(&
1317 & 6.20e-5*dtd(k)*dt(k)+(0.0149+6.20e-5*dt(k))*dtd(k)))/500.0
1318  xx = dh2o(k)*(pa(k)/500.0)*(1.+(0.0149+6.20e-5*dt(k))*dt(k))
1319 !-----six exponentials by powers of 8
1320  h2oexpd(k, 1) = -(0.10624*xxd*exp(-(xx*0.10624)))
1321  h2oexp(k, 1) = exp(-(xx*0.10624))
1322  xxd = h2oexpd(k, 1)*h2oexp(k, 1) + h2oexp(k, 1)*h2oexpd(k, 1)
1323  xx = h2oexp(k, 1)*h2oexp(k, 1)
1324  xxd = xxd*xx + xx*xxd
1325  xx = xx*xx
1326  h2oexpd(k, 2) = xxd*xx + xx*xxd
1327  h2oexp(k, 2) = xx*xx
1328  xxd = h2oexpd(k, 2)*h2oexp(k, 2) + h2oexp(k, 2)*h2oexpd(k, 2)
1329  xx = h2oexp(k, 2)*h2oexp(k, 2)
1330  xxd = xxd*xx + xx*xxd
1331  xx = xx*xx
1332  h2oexpd(k, 3) = xxd*xx + xx*xxd
1333  h2oexp(k, 3) = xx*xx
1334  xxd = h2oexpd(k, 3)*h2oexp(k, 3) + h2oexp(k, 3)*h2oexpd(k, 3)
1335  xx = h2oexp(k, 3)*h2oexp(k, 3)
1336  xxd = xxd*xx + xx*xxd
1337  xx = xx*xx
1338  h2oexpd(k, 4) = xxd*xx + xx*xxd
1339  h2oexp(k, 4) = xx*xx
1340  xxd = h2oexpd(k, 4)*h2oexp(k, 4) + h2oexp(k, 4)*h2oexpd(k, 4)
1341  xx = h2oexp(k, 4)*h2oexp(k, 4)
1342  xxd = xxd*xx + xx*xxd
1343  xx = xx*xx
1344  h2oexpd(k, 5) = xxd*xx + xx*xxd
1345  h2oexp(k, 5) = xx*xx
1346 !-----one exponential of h2o continuum for sub-band 3a (Table 9).
1347  conexpd(k) = -(109.0*dcontd(k)*exp(-(dcont(k)*109.0)))
1348  conexp(k) = exp(-(dcont(k)*109.0))
1349 !-----Scaled co2 amount for the Band 10 (Eq. 4.4, Tables 3 and 6).
1350  xxd = dco2(k)*(pa(k)/300.0)**0.5*(1.02e-4*dtd(k)*dt(k)+(0.0179+&
1351 & 1.02e-4*dt(k))*dtd(k))
1352  xx = dco2(k)*(pa(k)/300.0)**0.5*(1.+(0.0179+1.02e-4*dt(k))*dt(k))
1353 !-----six exponentials by powers of 8
1354  co2expd(k, 1) = -(2.656e-5*xxd*exp(-(xx*2.656e-5)))
1355  co2exp(k, 1) = exp(-(xx*2.656e-5))
1356  xxd = co2expd(k, 1)*co2exp(k, 1) + co2exp(k, 1)*co2expd(k, 1)
1357  xx = co2exp(k, 1)*co2exp(k, 1)
1358  xxd = xxd*xx + xx*xxd
1359  xx = xx*xx
1360  co2expd(k, 2) = xxd*xx + xx*xxd
1361  co2exp(k, 2) = xx*xx
1362  xxd = co2expd(k, 2)*co2exp(k, 2) + co2exp(k, 2)*co2expd(k, 2)
1363  xx = co2exp(k, 2)*co2exp(k, 2)
1364  xxd = xxd*xx + xx*xxd
1365  xx = xx*xx
1366  co2expd(k, 3) = xxd*xx + xx*xxd
1367  co2exp(k, 3) = xx*xx
1368  xxd = co2expd(k, 3)*co2exp(k, 3) + co2exp(k, 3)*co2expd(k, 3)
1369  xx = co2exp(k, 3)*co2exp(k, 3)
1370  xxd = xxd*xx + xx*xxd
1371  xx = xx*xx
1372  co2expd(k, 4) = xxd*xx + xx*xxd
1373  co2exp(k, 4) = xx*xx
1374  xxd = co2expd(k, 4)*co2exp(k, 4) + co2exp(k, 4)*co2expd(k, 4)
1375  xx = co2exp(k, 4)*co2exp(k, 4)
1376  xxd = xxd*xx + xx*xxd
1377  xx = xx*xx
1378  co2expd(k, 5) = xxd*xx + xx*xxd
1379  co2exp(k, 5) = xx*xx
1380  xxd = co2expd(k, 5)*co2exp(k, 5) + co2exp(k, 5)*co2expd(k, 5)
1381  xx = co2exp(k, 5)*co2exp(k, 5)
1382  xxd = xxd*xx + xx*xxd
1383  xx = xx*xx
1384  co2expd(k, 6) = xxd*xx + xx*xxd
1385  co2exp(k, 6) = xx*xx
1386 !-----Compute the scaled n2o amount for Band 10 (Table 5).
1387  xxd = dn2o(k)*(3.6656e-6*dtd(k)*dt(k)+(1.4476e-3+3.6656e-6*dt(k))*&
1388 & dtd(k))
1389  xx = dn2o(k)*(1.+(1.4476e-3+3.6656e-6*dt(k))*dt(k))
1390 !-----Two exponentials by powers of 58
1391  n2oexpd(k, 1) = -(0.25238*xxd*exp(-(xx*0.25238)))
1392  n2oexp(k, 1) = exp(-(xx*0.25238))
1393  xxd = n2oexpd(k, 1)*n2oexp(k, 1) + n2oexp(k, 1)*n2oexpd(k, 1)
1394  xx = n2oexp(k, 1)*n2oexp(k, 1)
1395  xx1d = xxd*xx + xx*xxd
1396  xx1 = xx*xx
1397  xx1d = xx1d*xx1 + xx1*xx1d
1398  xx1 = xx1*xx1
1399  xx2d = xx1d*xx1 + xx1*xx1d
1400  xx2 = xx1*xx1
1401  xx3d = xx2d*xx2 + xx2*xx2d
1402  xx3 = xx2*xx2
1403  n2oexpd(k, 2) = (xxd*xx1+xx*xx1d)*xx2*xx3 + xx*xx1*(xx2d*xx3+xx2*&
1404 & xx3d)
1405  n2oexp(k, 2) = xx*xx1*xx2*xx3
1406  END DO
1407 END SUBROUTINE b10exps_d
1408 
1409 ! Differentiation of tablup in forward (tangent) mode:
1410 ! variations of useful results: s1 s2 s3 tran
1411 ! with respect to varying inputs: s1 s2 s3 dt dw tran
1412 !**********************************************************************
1413 SUBROUTINE tablup_d(nx1, nh1, dw, dwd, p, dt, dtd, s1, s1d, s2, s2d, s3&
1414 & , s3d, w1, p1, dwe, dpe, coef1, coef2, coef3, tran, trand)
1415  IMPLICIT NONE
1416  INTEGER :: nx1, nh1
1417 !---- input parameters -----
1418  REAL*8 :: w1, p1, dwe, dpe
1419  REAL*8 :: dw, p, dt
1420  REAL*8 :: dwd, dtd
1421  REAL*8 :: coef1(nx1, nh1), coef2(nx1, nh1), coef3(nx1, nh1)
1422 !---- update parameter -----
1423  REAL*8 :: s1, s2, s3, tran
1424  REAL*8 :: s1d, s2d, s3d, trand
1425 !---- temporary variables -----
1426  REAL*8 :: we, pe, fw, fp, pa, pb, pc, ax, ba, bb, t1, ca, cb, t2
1427  REAL*8 :: wed, ped, fwd, fpd, pad, pbd, pcd, axd, bad, bbd, t1d, cad, &
1428 & cbd, t2d
1429  REAL*8 :: x1, x2, x3, xx, x1c
1430  REAL*8 :: x1d, x2d, x3d, xxd, x1cd
1431  INTEGER :: iw, ip
1432  INTRINSIC log10
1433  INTRINSIC real
1434  INTRINSIC min
1435  INTRINSIC int
1436  INTRINSIC max
1437  REAL*8 :: y2
1438  REAL*8 :: y1
1439 !-----Compute effective pressure (x2) and temperature (x3) following
1440 ! Eqs. (8.28) and (8.29)
1441  s1d = s1d + dwd
1442  s1 = s1 + dw
1443  s2d = s2d + p*dwd
1444  s2 = s2 + p*dw
1445  s3d = s3d + dtd*dw + dt*dwd
1446  s3 = s3 + dt*dw
1447  x1d = s1d
1448  x1 = s1
1449  x1cd = -(s1d/s1**2)
1450  x1c = 1.0/s1
1451  x2d = s2d*x1c + s2*x1cd
1452  x2 = s2*x1c
1453  x3d = s3d*x1c + s3*x1cd
1454  x3 = s3*x1c
1455 !-----normalize we and pe
1456 ! we=(log10(x1)-w1)/dwe
1457 ! pe=(log10(x2)-p1)/dpe
1458  wed = dwe*x1d/(x1*log(10.0))
1459  we = (log10(x1)-w1)*dwe
1460  ped = dpe*x2d/(x2*log(10.0))
1461  pe = (log10(x2)-p1)*dpe
1462  y1 = REAL(nh1 - 1)
1463  IF (we .GT. y1) THEN
1464  we = y1
1465  wed = 0.0_8
1466  ELSE
1467  we = we
1468  END IF
1469  y2 = REAL(nx1 - 1)
1470  IF (pe .GT. y2) THEN
1471  pe = y2
1472  ped = 0.0_8
1473  ELSE
1474  pe = pe
1475  END IF
1476 !-----assign iw and ip and compute the distance of we and pe
1477 ! from iw and ip.
1478  iw = int(we + 1.0)
1479  IF (iw .GT. nh1 - 1) THEN
1480  iw = nh1 - 1
1481  ELSE
1482  iw = iw
1483  END IF
1484  IF (iw .LT. 2) THEN
1485  iw = 2
1486  ELSE
1487  iw = iw
1488  END IF
1489  fwd = wed
1490  fw = we - REAL(iw - 1)
1491  ip = int(pe + 1.0)
1492  IF (ip .GT. nx1 - 1) THEN
1493  ip = nx1 - 1
1494  ELSE
1495  ip = ip
1496  END IF
1497  IF (ip .LT. 1) THEN
1498  ip = 1
1499  ELSE
1500  ip = ip
1501  END IF
1502  fpd = ped
1503  fp = pe - REAL(ip - 1)
1504 !-----linear interpolation in pressure
1505  pad = (coef1(ip+1, iw-1)-coef1(ip, iw-1))*fpd
1506  pa = coef1(ip, iw-1) + (coef1(ip+1, iw-1)-coef1(ip, iw-1))*fp
1507  pbd = (coef1(ip+1, iw)-coef1(ip, iw))*fpd
1508  pb = coef1(ip, iw) + (coef1(ip+1, iw)-coef1(ip, iw))*fp
1509  pcd = (coef1(ip+1, iw+1)-coef1(ip, iw+1))*fpd
1510  pc = coef1(ip, iw+1) + (coef1(ip+1, iw+1)-coef1(ip, iw+1))*fp
1511 !-----quadratic interpolation in absorber amount for coef1
1512  axd = 0.5*(((pcd+pad)*fw+(pc+pa)*fwd+pcd-pad)*fw+((pc+pa)*fw+(pc-pa))*&
1513 & fwd) + pbd*(1.-fw*fw) + pb*(-(fwd*fw)-fw*fwd)
1514  ax = ((pc+pa)*fw+(pc-pa))*fw*0.5 + pb*(1.-fw*fw)
1515 !-----linear interpolation in absorber amount for coef2 and coef3
1516  bad = (coef2(ip+1, iw)-coef2(ip, iw))*fpd
1517  ba = coef2(ip, iw) + (coef2(ip+1, iw)-coef2(ip, iw))*fp
1518  bbd = (coef2(ip+1, iw+1)-coef2(ip, iw+1))*fpd
1519  bb = coef2(ip, iw+1) + (coef2(ip+1, iw+1)-coef2(ip, iw+1))*fp
1520  t1d = bad + (bbd-bad)*fw + (bb-ba)*fwd
1521  t1 = ba + (bb-ba)*fw
1522  cad = (coef3(ip+1, iw)-coef3(ip, iw))*fpd
1523  ca = coef3(ip, iw) + (coef3(ip+1, iw)-coef3(ip, iw))*fp
1524  cbd = (coef3(ip+1, iw+1)-coef3(ip, iw+1))*fpd
1525  cb = coef3(ip, iw+1) + (coef3(ip+1, iw+1)-coef3(ip, iw+1))*fp
1526  t2d = cad + (cbd-cad)*fw + (cb-ca)*fwd
1527  t2 = ca + (cb-ca)*fw
1528 !-----update the total transmittance between levels k1 and k2
1529  xxd = axd + (t1d+t2d*x3+t2*x3d)*x3 + (t1+t2*x3)*x3d
1530  xx = ax + (t1+t2*x3)*x3
1531  IF (xx .GT. 0.9999999) THEN
1532  xx = 0.9999999
1533  xxd = 0.0_8
1534  ELSE
1535  xx = xx
1536  END IF
1537  IF (xx .LT. 0.0000001) THEN
1538  xx = 0.0000001
1539  xxd = 0.0_8
1540  ELSE
1541  xx = xx
1542  END IF
1543  trand = trand*xx + tran*xxd
1544  tran = tran*xx
1545 END SUBROUTINE tablup_d
1546 
1547 ! Differentiation of h2okdis in forward (tangent) mode:
1548 ! variations of useful results: tran tcon th2o
1549 ! with respect to varying inputs: tcon h2oexp th2o conexp
1550 !**********************************************************************
1551 SUBROUTINE h2okdis_d(ib, np, k, fkw, gkw, ne, h2oexp, h2oexpd, conexp, &
1552 & conexpd, th2o, th2od, tcon, tcond, tran, trand)
1553  IMPLICIT NONE
1554 !---- input parameters ------
1555  INTEGER :: ib, ne, np, k
1556  REAL*8 :: h2oexp(0:np, 6), conexp(0:np, 3)
1557  REAL*8 :: h2oexpd(0:np, 6), conexpd(0:np, 3)
1558  REAL*8 :: fkw(6, 9), gkw(6, 3)
1559 !---- updated parameters -----
1560  REAL*8 :: th2o(6), tcon(3), tran
1561  REAL*8 :: th2od(6), tcond(3), trand
1562 !---- temporary arrays -----
1563  REAL*8 :: trnth2o
1564  REAL*8 :: trnth2od
1565 !-----tco2 are the six exp factors between levels k1 and k2
1566 ! tran is the updated total transmittance between levels k1 and k2
1567 !-----th2o is the 6 exp factors between levels k1 and k2 due to
1568 ! h2o line absorption.
1569 !-----tcon is the 3 exp factors between levels k1 and k2 due to
1570 ! h2o continuum absorption.
1571 !-----trnth2o is the total transmittance between levels k1 and k2 due
1572 ! to both line and continuum absorption.
1573 !-----Compute th2o following Eq. (8.23).
1574  th2od(1) = th2od(1)*h2oexp(k, 1) + th2o(1)*h2oexpd(k, 1)
1575  th2o(1) = th2o(1)*h2oexp(k, 1)
1576  th2od(2) = th2od(2)*h2oexp(k, 2) + th2o(2)*h2oexpd(k, 2)
1577  th2o(2) = th2o(2)*h2oexp(k, 2)
1578  th2od(3) = th2od(3)*h2oexp(k, 3) + th2o(3)*h2oexpd(k, 3)
1579  th2o(3) = th2o(3)*h2oexp(k, 3)
1580  th2od(4) = th2od(4)*h2oexp(k, 4) + th2o(4)*h2oexpd(k, 4)
1581  th2o(4) = th2o(4)*h2oexp(k, 4)
1582  th2od(5) = th2od(5)*h2oexp(k, 5) + th2o(5)*h2oexpd(k, 5)
1583  th2o(5) = th2o(5)*h2oexp(k, 5)
1584  th2od(6) = th2od(6)*h2oexp(k, 6) + th2o(6)*h2oexpd(k, 6)
1585  th2o(6) = th2o(6)*h2oexp(k, 6)
1586  IF (ne .EQ. 0) THEN
1587 !-----Compute trnh2o following Eq. (8.25). fkw is given in Table 4.
1588  trnth2od = fkw(1, ib)*th2od(1) + fkw(2, ib)*th2od(2) + fkw(3, ib)*&
1589 & th2od(3) + fkw(4, ib)*th2od(4) + fkw(5, ib)*th2od(5) + fkw(6, ib)*&
1590 & th2od(6)
1591  trnth2o = fkw(1, ib)*th2o(1) + fkw(2, ib)*th2o(2) + fkw(3, ib)*th2o(&
1592 & 3) + fkw(4, ib)*th2o(4) + fkw(5, ib)*th2o(5) + fkw(6, ib)*th2o(6)
1593  trand = tran*trnth2od
1594  tran = tran*trnth2o
1595  ELSE IF (ne .EQ. 1) THEN
1596 !-----Compute trnh2o following Eqs. (8.25) and (4.27).
1597  tcond(1) = tcond(1)*conexp(k, 1) + tcon(1)*conexpd(k, 1)
1598  tcon(1) = tcon(1)*conexp(k, 1)
1599  trnth2od = (fkw(1, ib)*th2od(1)+fkw(2, ib)*th2od(2)+fkw(3, ib)*th2od&
1600 & (3)+fkw(4, ib)*th2od(4)+fkw(5, ib)*th2od(5)+fkw(6, ib)*th2od(6))*&
1601 & tcon(1) + (fkw(1, ib)*th2o(1)+fkw(2, ib)*th2o(2)+fkw(3, ib)*th2o(3&
1602 & )+fkw(4, ib)*th2o(4)+fkw(5, ib)*th2o(5)+fkw(6, ib)*th2o(6))*tcond(&
1603 & 1)
1604  trnth2o = (fkw(1, ib)*th2o(1)+fkw(2, ib)*th2o(2)+fkw(3, ib)*th2o(3)+&
1605 & fkw(4, ib)*th2o(4)+fkw(5, ib)*th2o(5)+fkw(6, ib)*th2o(6))*tcon(1)
1606  trand = tran*trnth2od
1607  tran = tran*trnth2o
1608  ELSE
1609 !-----For band 3. This band is divided into 3 subbands.
1610  tcond(1) = tcond(1)*conexp(k, 1) + tcon(1)*conexpd(k, 1)
1611  tcon(1) = tcon(1)*conexp(k, 1)
1612  tcond(2) = tcond(2)*conexp(k, 2) + tcon(2)*conexpd(k, 2)
1613  tcon(2) = tcon(2)*conexp(k, 2)
1614  tcond(3) = tcond(3)*conexp(k, 3) + tcon(3)*conexpd(k, 3)
1615  tcon(3) = tcon(3)*conexp(k, 3)
1616 !-----Compute trnh2o following Eqs. (4.29) and (8.25).
1617  trnth2od = (gkw(1, 1)*th2od(1)+gkw(2, 1)*th2od(2)+gkw(3, 1)*th2od(3)&
1618 & +gkw(4, 1)*th2od(4)+gkw(5, 1)*th2od(5)+gkw(6, 1)*th2od(6))*tcon(1)&
1619 & + (gkw(1, 1)*th2o(1)+gkw(2, 1)*th2o(2)+gkw(3, 1)*th2o(3)+gkw(4, 1)&
1620 & *th2o(4)+gkw(5, 1)*th2o(5)+gkw(6, 1)*th2o(6))*tcond(1) + (gkw(1, 2&
1621 & )*th2od(1)+gkw(2, 2)*th2od(2)+gkw(3, 2)*th2od(3)+gkw(4, 2)*th2od(4&
1622 & )+gkw(5, 2)*th2od(5)+gkw(6, 2)*th2od(6))*tcon(2) + (gkw(1, 2)*th2o&
1623 & (1)+gkw(2, 2)*th2o(2)+gkw(3, 2)*th2o(3)+gkw(4, 2)*th2o(4)+gkw(5, 2&
1624 & )*th2o(5)+gkw(6, 2)*th2o(6))*tcond(2) + (gkw(1, 3)*th2od(1)+gkw(2&
1625 & , 3)*th2od(2)+gkw(3, 3)*th2od(3)+gkw(4, 3)*th2od(4)+gkw(5, 3)*&
1626 & th2od(5)+gkw(6, 3)*th2od(6))*tcon(3) + (gkw(1, 3)*th2o(1)+gkw(2, 3&
1627 & )*th2o(2)+gkw(3, 3)*th2o(3)+gkw(4, 3)*th2o(4)+gkw(5, 3)*th2o(5)+&
1628 & gkw(6, 3)*th2o(6))*tcond(3)
1629  trnth2o = (gkw(1, 1)*th2o(1)+gkw(2, 1)*th2o(2)+gkw(3, 1)*th2o(3)+gkw&
1630 & (4, 1)*th2o(4)+gkw(5, 1)*th2o(5)+gkw(6, 1)*th2o(6))*tcon(1) + (gkw&
1631 & (1, 2)*th2o(1)+gkw(2, 2)*th2o(2)+gkw(3, 2)*th2o(3)+gkw(4, 2)*th2o(&
1632 & 4)+gkw(5, 2)*th2o(5)+gkw(6, 2)*th2o(6))*tcon(2) + (gkw(1, 3)*th2o(&
1633 & 1)+gkw(2, 3)*th2o(2)+gkw(3, 3)*th2o(3)+gkw(4, 3)*th2o(4)+gkw(5, 3)&
1634 & *th2o(5)+gkw(6, 3)*th2o(6))*tcon(3)
1635  trand = tran*trnth2od
1636  tran = tran*trnth2o
1637  END IF
1638 END SUBROUTINE h2okdis_d
1639 
1640 ! Differentiation of n2okdis in forward (tangent) mode:
1641 ! variations of useful results: tran tn2o
1642 ! with respect to varying inputs: tran tn2o n2oexp
1643 !**********************************************************************
1644 SUBROUTINE n2okdis_d(ib, np, k, n2oexp, n2oexpd, tn2o, tn2od, tran, &
1645 & trand)
1646  IMPLICIT NONE
1647  INTEGER :: ib, np, k
1648 !---- input parameters -----
1649  REAL*8 :: n2oexp(0:np, 4)
1650  REAL*8 :: n2oexpd(0:np, 4)
1651 !---- updated parameters -----
1652  REAL*8 :: tn2o(4), tran
1653  REAL*8 :: tn2od(4), trand
1654 !---- temporary arrays -----
1655  REAL*8 :: xc
1656  REAL*8 :: xcd
1657 !-----tn2o is computed from Eq. (8.23).
1658 ! xc is the total n2o transmittance computed from (8.25)
1659 ! The k-distribution functions are given in Table 5.
1660 !-----band 6
1661  IF (ib .EQ. 6) THEN
1662  tn2od(1) = tn2od(1)*n2oexp(k, 1) + tn2o(1)*n2oexpd(k, 1)
1663  tn2o(1) = tn2o(1)*n2oexp(k, 1)
1664  xcd = 0.940414*tn2od(1)
1665  xc = 0.940414*tn2o(1)
1666  tn2od(2) = tn2od(2)*n2oexp(k, 2) + tn2o(2)*n2oexpd(k, 2)
1667  tn2o(2) = tn2o(2)*n2oexp(k, 2)
1668  xcd = xcd + 0.059586*tn2od(2)
1669  xc = xc + 0.059586*tn2o(2)
1670 !-----band 7
1671  ELSE
1672  tn2od(1) = tn2od(1)*n2oexp(k, 1) + tn2o(1)*n2oexpd(k, 1)
1673  tn2o(1) = tn2o(1)*n2oexp(k, 1)
1674  xcd = 0.561961*tn2od(1)
1675  xc = 0.561961*tn2o(1)
1676  tn2od(2) = tn2od(2)*n2oexp(k, 2) + tn2o(2)*n2oexpd(k, 2)
1677  tn2o(2) = tn2o(2)*n2oexp(k, 2)
1678  xcd = xcd + 0.138707*tn2od(2)
1679  xc = xc + 0.138707*tn2o(2)
1680  tn2od(3) = tn2od(3)*n2oexp(k, 3) + tn2o(3)*n2oexpd(k, 3)
1681  tn2o(3) = tn2o(3)*n2oexp(k, 3)
1682  xcd = xcd + 0.240670*tn2od(3)
1683  xc = xc + 0.240670*tn2o(3)
1684  tn2od(4) = tn2od(4)*n2oexp(k, 4) + tn2o(4)*n2oexpd(k, 4)
1685  tn2o(4) = tn2o(4)*n2oexp(k, 4)
1686  xcd = xcd + 0.058662*tn2od(4)
1687  xc = xc + 0.058662*tn2o(4)
1688  END IF
1689  trand = trand*xc + tran*xcd
1690  tran = tran*xc
1691 END SUBROUTINE n2okdis_d
1692 
1693 ! Differentiation of ch4kdis in forward (tangent) mode:
1694 ! variations of useful results: tran tch4
1695 ! with respect to varying inputs: tran ch4exp tch4
1696 !**********************************************************************
1697 SUBROUTINE ch4kdis_d(ib, np, k, ch4exp, ch4expd, tch4, tch4d, tran, &
1698 & trand)
1699  IMPLICIT NONE
1700  INTEGER :: ib, np, k
1701 !---- input parameters -----
1702  REAL*8 :: ch4exp(0:np, 4)
1703  REAL*8 :: ch4expd(0:np, 4)
1704 !---- updated parameters -----
1705  REAL*8 :: tch4(4), tran
1706  REAL*8 :: tch4d(4), trand
1707 !---- temporary arrays -----
1708  REAL*8 :: xc
1709  REAL*8 :: xcd
1710 !-----tch4 is computed from Eq. (8.23).
1711 ! xc is the total ch4 transmittance computed from (8.25)
1712 ! The k-distribution functions are given in Table 5.
1713 !-----band 6
1714  IF (ib .EQ. 6) THEN
1715  tch4d(1) = tch4d(1)*ch4exp(k, 1) + tch4(1)*ch4expd(k, 1)
1716  tch4(1) = tch4(1)*ch4exp(k, 1)
1717  xcd = tch4d(1)
1718  xc = tch4(1)
1719 !-----band 7
1720  ELSE
1721  tch4d(1) = tch4d(1)*ch4exp(k, 1) + tch4(1)*ch4expd(k, 1)
1722  tch4(1) = tch4(1)*ch4exp(k, 1)
1723  xcd = 0.610650*tch4d(1)
1724  xc = 0.610650*tch4(1)
1725  tch4d(2) = tch4d(2)*ch4exp(k, 2) + tch4(2)*ch4expd(k, 2)
1726  tch4(2) = tch4(2)*ch4exp(k, 2)
1727  xcd = xcd + 0.280212*tch4d(2)
1728  xc = xc + 0.280212*tch4(2)
1729  tch4d(3) = tch4d(3)*ch4exp(k, 3) + tch4(3)*ch4expd(k, 3)
1730  tch4(3) = tch4(3)*ch4exp(k, 3)
1731  xcd = xcd + 0.107349*tch4d(3)
1732  xc = xc + 0.107349*tch4(3)
1733  tch4d(4) = tch4d(4)*ch4exp(k, 4) + tch4(4)*ch4expd(k, 4)
1734  tch4(4) = tch4(4)*ch4exp(k, 4)
1735  xcd = xcd + 0.001789*tch4d(4)
1736  xc = xc + 0.001789*tch4(4)
1737  END IF
1738  trand = trand*xc + tran*xcd
1739  tran = tran*xc
1740 END SUBROUTINE ch4kdis_d
1741 
1742 ! Differentiation of comkdis in forward (tangent) mode:
1743 ! variations of useful results: tran tcom
1744 ! with respect to varying inputs: tran tcom comexp
1745 !**********************************************************************
1746 SUBROUTINE comkdis_d(ib, np, k, comexp, comexpd, tcom, tcomd, tran, &
1747 & trand)
1748  IMPLICIT NONE
1749  INTEGER :: ib, np, k
1750 !---- input parameters -----
1751  REAL*8 :: comexp(0:np, 6)
1752  REAL*8 :: comexpd(0:np, 6)
1753 !---- updated parameters -----
1754  REAL*8 :: tcom(6), tran
1755  REAL*8 :: tcomd(6), trand
1756 !---- temporary arrays -----
1757  REAL*8 :: xc
1758  REAL*8 :: xcd
1759 !-----tcom is computed from Eq. (8.23).
1760 ! xc is the total co2 transmittance computed from (8.25)
1761 ! The k-distribution functions are given in Table 6.
1762 !-----band 4
1763  IF (ib .EQ. 4) THEN
1764  tcomd(1) = tcomd(1)*comexp(k, 1) + tcom(1)*comexpd(k, 1)
1765  tcom(1) = tcom(1)*comexp(k, 1)
1766  xcd = 0.12159*tcomd(1)
1767  xc = 0.12159*tcom(1)
1768  tcomd(2) = tcomd(2)*comexp(k, 2) + tcom(2)*comexpd(k, 2)
1769  tcom(2) = tcom(2)*comexp(k, 2)
1770  xcd = xcd + 0.24359*tcomd(2)
1771  xc = xc + 0.24359*tcom(2)
1772  tcomd(3) = tcomd(3)*comexp(k, 3) + tcom(3)*comexpd(k, 3)
1773  tcom(3) = tcom(3)*comexp(k, 3)
1774  xcd = xcd + 0.24981*tcomd(3)
1775  xc = xc + 0.24981*tcom(3)
1776  tcomd(4) = tcomd(4)*comexp(k, 4) + tcom(4)*comexpd(k, 4)
1777  tcom(4) = tcom(4)*comexp(k, 4)
1778  xcd = xcd + 0.26427*tcomd(4)
1779  xc = xc + 0.26427*tcom(4)
1780  tcomd(5) = tcomd(5)*comexp(k, 5) + tcom(5)*comexpd(k, 5)
1781  tcom(5) = tcom(5)*comexp(k, 5)
1782  xcd = xcd + 0.07807*tcomd(5)
1783  xc = xc + 0.07807*tcom(5)
1784  tcomd(6) = tcomd(6)*comexp(k, 6) + tcom(6)*comexpd(k, 6)
1785  tcom(6) = tcom(6)*comexp(k, 6)
1786  xcd = xcd + 0.04267*tcomd(6)
1787  xc = xc + 0.04267*tcom(6)
1788 !-----band 5
1789  ELSE
1790  tcomd(1) = tcomd(1)*comexp(k, 1) + tcom(1)*comexpd(k, 1)
1791  tcom(1) = tcom(1)*comexp(k, 1)
1792  xcd = 0.06869*tcomd(1)
1793  xc = 0.06869*tcom(1)
1794  tcomd(2) = tcomd(2)*comexp(k, 2) + tcom(2)*comexpd(k, 2)
1795  tcom(2) = tcom(2)*comexp(k, 2)
1796  xcd = xcd + 0.14795*tcomd(2)
1797  xc = xc + 0.14795*tcom(2)
1798  tcomd(3) = tcomd(3)*comexp(k, 3) + tcom(3)*comexpd(k, 3)
1799  tcom(3) = tcom(3)*comexp(k, 3)
1800  xcd = xcd + 0.19512*tcomd(3)
1801  xc = xc + 0.19512*tcom(3)
1802  tcomd(4) = tcomd(4)*comexp(k, 4) + tcom(4)*comexpd(k, 4)
1803  tcom(4) = tcom(4)*comexp(k, 4)
1804  xcd = xcd + 0.33446*tcomd(4)
1805  xc = xc + 0.33446*tcom(4)
1806  tcomd(5) = tcomd(5)*comexp(k, 5) + tcom(5)*comexpd(k, 5)
1807  tcom(5) = tcom(5)*comexp(k, 5)
1808  xcd = xcd + 0.17199*tcomd(5)
1809  xc = xc + 0.17199*tcom(5)
1810  tcomd(6) = tcomd(6)*comexp(k, 6) + tcom(6)*comexpd(k, 6)
1811  tcom(6) = tcom(6)*comexp(k, 6)
1812  xcd = xcd + 0.08179*tcomd(6)
1813  xc = xc + 0.08179*tcom(6)
1814  END IF
1815  trand = trand*xc + tran*xcd
1816  tran = tran*xc
1817 END SUBROUTINE comkdis_d
1818 
1819 ! Differentiation of cfckdis in forward (tangent) mode:
1820 ! variations of useful results: tran tcfc
1821 ! with respect to varying inputs: tran tcfc cfcexp
1822 !**********************************************************************
1823 SUBROUTINE cfckdis_d(np, k, cfcexp, cfcexpd, tcfc, tcfcd, tran, trand)
1824  IMPLICIT NONE
1825 !---- input parameters -----
1826  INTEGER :: k, np
1827  REAL*8 :: cfcexp(0:np)
1828  REAL*8 :: cfcexpd(0:np)
1829 !---- updated parameters -----
1830  REAL*8 :: tcfc, tran
1831  REAL*8 :: tcfcd, trand
1832 !-----tcfc is the exp factors between levels k1 and k2.
1833  tcfcd = tcfcd*cfcexp(k) + tcfc*cfcexpd(k)
1834  tcfc = tcfc*cfcexp(k)
1835  trand = trand*tcfc + tran*tcfcd
1836  tran = tran*tcfc
1837 END SUBROUTINE cfckdis_d
1838 
1839 ! Differentiation of b10kdis in forward (tangent) mode:
1840 ! variations of useful results: tran tn2o tcon th2o tco2
1841 ! with respect to varying inputs: tn2o co2exp tcon h2oexp n2oexp
1842 ! th2o tco2 conexp
1843 !**********************************************************************
1844 SUBROUTINE b10kdis_d(np, k, h2oexp, h2oexpd, conexp, conexpd, co2exp, &
1845 & co2expd, n2oexp, n2oexpd, th2o, th2od, tcon, tcond, tco2, tco2d, tn2o&
1846 & , tn2od, tran, trand)
1847  IMPLICIT NONE
1848  INTEGER :: np, k
1849 !---- input parameters -----
1850  REAL*8 :: h2oexp(0:np, 5), conexp(0:np), co2exp(0:np, 6), n2oexp(0:np&
1851 & , 2)
1852  REAL*8 :: h2oexpd(0:np, 5), conexpd(0:np), co2expd(0:np, 6), n2oexpd(0&
1853 & :np, 2)
1854 !---- updated parameters -----
1855  REAL*8 :: th2o(6), tcon(3), tco2(6), tn2o(4), tran
1856  REAL*8 :: th2od(6), tcond(3), tco2d(6), tn2od(4), trand
1857 !---- temporary arrays -----
1858  REAL*8 :: xx
1859  REAL*8 :: xxd
1860 !-----For h2o line. The k-distribution functions are given in Table 4.
1861  th2od(1) = th2od(1)*h2oexp(k, 1) + th2o(1)*h2oexpd(k, 1)
1862  th2o(1) = th2o(1)*h2oexp(k, 1)
1863  xxd = 0.3153*th2od(1)
1864  xx = 0.3153*th2o(1)
1865  th2od(2) = th2od(2)*h2oexp(k, 2) + th2o(2)*h2oexpd(k, 2)
1866  th2o(2) = th2o(2)*h2oexp(k, 2)
1867  xxd = xxd + 0.4604*th2od(2)
1868  xx = xx + 0.4604*th2o(2)
1869  th2od(3) = th2od(3)*h2oexp(k, 3) + th2o(3)*h2oexpd(k, 3)
1870  th2o(3) = th2o(3)*h2oexp(k, 3)
1871  xxd = xxd + 0.1326*th2od(3)
1872  xx = xx + 0.1326*th2o(3)
1873  th2od(4) = th2od(4)*h2oexp(k, 4) + th2o(4)*h2oexpd(k, 4)
1874  th2o(4) = th2o(4)*h2oexp(k, 4)
1875  xxd = xxd + 0.0798*th2od(4)
1876  xx = xx + 0.0798*th2o(4)
1877  th2od(5) = th2od(5)*h2oexp(k, 5) + th2o(5)*h2oexpd(k, 5)
1878  th2o(5) = th2o(5)*h2oexp(k, 5)
1879  xxd = xxd + 0.0119*th2od(5)
1880  xx = xx + 0.0119*th2o(5)
1881  trand = xxd
1882  tran = xx
1883 !-----For h2o continuum. Note that conexp(k,3) is for subband 3a.
1884  tcond(1) = tcond(1)*conexp(k) + tcon(1)*conexpd(k)
1885  tcon(1) = tcon(1)*conexp(k)
1886  trand = trand*tcon(1) + tran*tcond(1)
1887  tran = tran*tcon(1)
1888 !-----For co2 (Table 6)
1889  tco2d(1) = tco2d(1)*co2exp(k, 1) + tco2(1)*co2expd(k, 1)
1890  tco2(1) = tco2(1)*co2exp(k, 1)
1891  xxd = 0.2673*tco2d(1)
1892  xx = 0.2673*tco2(1)
1893  tco2d(2) = tco2d(2)*co2exp(k, 2) + tco2(2)*co2expd(k, 2)
1894  tco2(2) = tco2(2)*co2exp(k, 2)
1895  xxd = xxd + 0.2201*tco2d(2)
1896  xx = xx + 0.2201*tco2(2)
1897  tco2d(3) = tco2d(3)*co2exp(k, 3) + tco2(3)*co2expd(k, 3)
1898  tco2(3) = tco2(3)*co2exp(k, 3)
1899  xxd = xxd + 0.2106*tco2d(3)
1900  xx = xx + 0.2106*tco2(3)
1901  tco2d(4) = tco2d(4)*co2exp(k, 4) + tco2(4)*co2expd(k, 4)
1902  tco2(4) = tco2(4)*co2exp(k, 4)
1903  xxd = xxd + 0.2409*tco2d(4)
1904  xx = xx + 0.2409*tco2(4)
1905  tco2d(5) = tco2d(5)*co2exp(k, 5) + tco2(5)*co2expd(k, 5)
1906  tco2(5) = tco2(5)*co2exp(k, 5)
1907  xxd = xxd + 0.0196*tco2d(5)
1908  xx = xx + 0.0196*tco2(5)
1909  tco2d(6) = tco2d(6)*co2exp(k, 6) + tco2(6)*co2expd(k, 6)
1910  tco2(6) = tco2(6)*co2exp(k, 6)
1911  xxd = xxd + 0.0415*tco2d(6)
1912  xx = xx + 0.0415*tco2(6)
1913  trand = trand*xx + tran*xxd
1914  tran = tran*xx
1915 !-----For n2o (Table 5)
1916  tn2od(1) = tn2od(1)*n2oexp(k, 1) + tn2o(1)*n2oexpd(k, 1)
1917  tn2o(1) = tn2o(1)*n2oexp(k, 1)
1918  xxd = 0.970831*tn2od(1)
1919  xx = 0.970831*tn2o(1)
1920  tn2od(2) = tn2od(2)*n2oexp(k, 2) + tn2o(2)*n2oexpd(k, 2)
1921  tn2o(2) = tn2o(2)*n2oexp(k, 2)
1922  xxd = xxd + 0.029169*tn2od(2)
1923  xx = xx + 0.029169*tn2o(2)
1924  trand = trand*(xx-1.0) + tran*xxd
1925  tran = tran*(xx-1.0)
1926 END SUBROUTINE b10kdis_d
1927 
1928 ! Differentiation of cldovlp in forward (tangent) mode:
1929 ! variations of useful results: cldhi cldlw cldmd
1930 ! with respect to varying inputs: cldhi ett cldlw enn cldmd
1931 !mjs
1932 !***********************************************************************
1933 SUBROUTINE cldovlp_d(np, k1, k2, ict, icb, icx, ncld, enn, ennd, ett, &
1934 & ettd, cldhi, cldhid, cldmd, cldmdd, cldlw, cldlwd)
1935  IMPLICIT NONE
1936  INTEGER, INTENT(IN) :: np, k1, k2, ict, icb, icx(0:np), ncld(3)
1937  REAL*8, INTENT(IN) :: enn(0:np), ett(0:np)
1938  REAL*8, INTENT(IN) :: ennd(0:np), ettd(0:np)
1939  REAL*8, INTENT(INOUT) :: cldhi, cldmd, cldlw
1940  REAL*8, INTENT(INOUT) :: cldhid, cldmdd, cldlwd
1941  INTEGER :: j, k, km, kx
1942  km = k2 - 1
1943  IF (km .LT. ict) THEN
1944 ! do high clouds
1945  kx = ncld(1)
1946  IF (kx .EQ. 1 .OR. cldhi .EQ. 0.) THEN
1947  cldhid = ennd(km)
1948  cldhi = enn(km)
1949  ELSE
1950 !if ( (kx.lt.1 .or. kx.gt.1) .and. abs(cldhi) .gt. 0.) then
1951  cldhi = 0.0
1952  IF (kx .NE. 0) THEN
1953  cldhid = 0.0_8
1954  DO k=ict-kx,ict-1
1955  j = icx(k)
1956  IF (j .GE. k1 .AND. j .LE. km) THEN
1957  cldhid = ennd(j) + ettd(j)*cldhi + ett(j)*cldhid
1958  cldhi = enn(j) + ett(j)*cldhi
1959  END IF
1960  END DO
1961  ELSE
1962  cldhid = 0.0_8
1963  END IF
1964  END IF
1965  ELSE IF (km .GE. ict .AND. km .LT. icb) THEN
1966 ! do middle clouds
1967  kx = ncld(2)
1968  IF (kx .EQ. 1 .OR. cldmd .EQ. 0.) THEN
1969  cldmdd = ennd(km)
1970  cldmd = enn(km)
1971  ELSE
1972 !if ( (kx.lt.1 .or. kx.gt.1) .and. abs(cldhi) .gt. 0.) then
1973  cldmd = 0.0
1974  IF (kx .NE. 0) THEN
1975  cldmdd = 0.0_8
1976  DO k=icb-kx,icb-1
1977  j = icx(k)
1978  IF (j .GE. k1 .AND. j .LE. km) THEN
1979  cldmdd = ennd(j) + ettd(j)*cldmd + ett(j)*cldmdd
1980  cldmd = enn(j) + ett(j)*cldmd
1981  END IF
1982  END DO
1983  ELSE
1984  cldmdd = 0.0_8
1985  END IF
1986  END IF
1987  ELSE
1988 ! do low clouds
1989  kx = ncld(3)
1990  IF (kx .EQ. 1 .OR. cldlw .EQ. 0.) THEN
1991  cldlwd = ennd(km)
1992  cldlw = enn(km)
1993  ELSE
1994 !if ( (kx.lt.1 .or. kx.gt.1) .and. abs(cldhi) .gt. 0.) then
1995  cldlw = 0.0
1996  IF (kx .NE. 0) THEN
1997  cldlwd = 0.0_8
1998  DO k=np+1-kx,np
1999  j = icx(k)
2000  IF (j .GE. k1 .AND. j .LE. km) THEN
2001  cldlwd = ennd(j) + ettd(j)*cldlw + ett(j)*cldlwd
2002  cldlw = enn(j) + ett(j)*cldlw
2003  END IF
2004  END DO
2005  ELSE
2006  cldlwd = 0.0_8
2007  END IF
2008  END IF
2009  END IF
2010 END SUBROUTINE cldovlp_d
2011 
2012 ! Differentiation of getirtau1 in forward (tangent) mode:
2013 ! variations of useful results: tcldlyr enn
2014 ! with respect to varying inputs: hydromets tcldlyr fcld enn
2015 ! reff
2016 SUBROUTINE getirtau1_d(ib, nlevs, dp, fcld, fcldd, reff, reffd, &
2017 & hydromets, hydrometsd, tcldlyr, tcldlyrd, enn, ennd, aib_ir1, awb_ir1&
2018 & , aiw_ir1, aww_ir1, aig_ir1, awg_ir1, cons_grav)
2019  IMPLICIT NONE
2020 ! !INPUT PARAMETERS:
2021 ! Band number
2022  INTEGER, INTENT(IN) :: ib
2023 ! Number of levels
2024  INTEGER, INTENT(IN) :: nlevs
2025 ! Delta pressure in Pa
2026  REAL*8, INTENT(IN) :: dp(nlevs)
2027 ! Cloud fraction (used sometimes)
2028  REAL*8, INTENT(IN) :: fcld(nlevs)
2029  REAL*8, INTENT(IN) :: fcldd(nlevs)
2030 ! Effective radius (microns)
2031  REAL*8, INTENT(IN) :: reff(nlevs, 4)
2032  REAL*8, INTENT(IN) :: reffd(nlevs, 4)
2033 ! Hydrometeors (kg/kg)
2034  REAL*8, INTENT(IN) :: hydromets(nlevs, 4)
2035  REAL*8, INTENT(IN) :: hydrometsd(nlevs, 4)
2036  REAL*8, INTENT(IN) :: aib_ir1(3, 10), awb_ir1(4, 10), aiw_ir1(4, 10)
2037  REAL*8, INTENT(IN) :: aww_ir1(4, 10), aig_ir1(4, 10), awg_ir1(4, 10)
2038  REAL*8, INTENT(IN) :: cons_grav
2039 ! !OUTPUT PARAMETERS:
2040 ! Flux transmissivity?
2041  REAL*8, INTENT(OUT) :: tcldlyr(0:nlevs)
2042  REAL*8, INTENT(OUT) :: tcldlyrd(0:nlevs)
2043 ! Flux transmissivity of a cloud layer?
2044  REAL*8, INTENT(OUT) :: enn(0:nlevs)
2045  REAL*8, INTENT(OUT) :: ennd(0:nlevs)
2046 ! !DESCRIPTION:
2047 ! Compute in-cloud or grid mean optical depths for infrared wavelengths
2048 ! Slots for reff, hydrometeors and tauall are as follows:
2049 ! 1 Cloud Ice
2050 ! 2 Cloud Liquid
2051 ! 3 Falling Liquid (Rain)
2052 ! 4 Falling Ice (Snow)
2053 !
2054 ! In the below calculations, the constants used in the tau calculation are in
2055 ! m$^2$ g$^{-1}$ and m$^2$ g$^{-1}$ $\mu$m. Thus, we must convert the kg contained in the
2056 ! pressure (Pa = kg m$^{-1}$ s$^{-2}$) to grams.
2057 !
2058 ! !REVISION HISTORY:
2059 ! 2011.11.18 MAT moved to Radiation_Shared and revised arg list, units
2060 !
2061 !EOP
2062 !------------------------------------------------------------------------------
2063 !BOC
2064  INTEGER :: k
2065  REAL*8 :: taucld1, taucld2, taucld3, taucld4
2066  REAL*8 :: taucld1d, taucld2d, taucld3d, taucld4d
2067  REAL*8 :: g1, g2, g3, g4, gg
2068  REAL*8 :: g1d, g2d, g3d, g4d, ggd
2069  REAL*8 :: w1, w2, w3, w4, ww
2070  REAL*8 :: w1d, w2d, w3d, w4d, wwd
2071  REAL*8 :: ff, tauc
2072  REAL*8 :: ffd, taucd
2073  REAL*8 :: reff_snow
2074  REAL*8 :: reff_snowd
2075  INTRINSIC min
2076  INTRINSIC abs
2077  INTRINSIC max
2078  INTRINSIC exp
2079  REAL*8 :: pwr1
2080  REAL*8 :: pwr1d
2081  REAL*8 :: max1d
2082  REAL*8 :: abs0
2083  REAL*8 :: max1
2084 !-----Compute cloud optical thickness following Eqs. (6.4a,b) and (6.7)
2085 ! Rain optical thickness is set to 0.00307 /(gm/m**2).
2086 ! It is for a specific drop size distribution provided by Q. Fu.
2087  tcldlyrd(0) = 0.0_8
2088  tcldlyr(0) = 1.0
2089  ennd(0) = 0.0_8
2090  enn(0) = 0.0
2091  DO k=1,nlevs
2092  IF (reff(k, 1) .LE. 0.0) THEN
2093  taucld1 = 0.0
2094  taucld1d = 0.0_8
2095  ELSE
2096  IF (reff(k, 1) .GT. 0.0 .OR. (reff(k, 1) .LT. 0.0 .AND. aib_ir1(3&
2097 & , ib) .EQ. int(aib_ir1(3, ib)))) THEN
2098  pwr1d = aib_ir1(3, ib)*reff(k, 1)**(aib_ir1(3, ib)-1)*reffd(k, 1&
2099 & )
2100  ELSE IF (reff(k, 1) .EQ. 0.0 .AND. aib_ir1(3, ib) .EQ. 1.0) THEN
2101  pwr1d = reffd(k, 1)
2102  ELSE
2103  pwr1d = 0.0
2104  END IF
2105  pwr1 = reff(k, 1)**aib_ir1(3, ib)
2106  taucld1d = dp(k)*1.0e3*(hydrometsd(k, 1)*(aib_ir1(1, ib)+aib_ir1(2&
2107 & , ib)/pwr1)-hydromets(k, 1)*aib_ir1(2, ib)*pwr1d/pwr1**2)/&
2108 & cons_grav
2109  taucld1 = dp(k)*1.0e3/cons_grav*hydromets(k, 1)*(aib_ir1(1, ib)+&
2110 & aib_ir1(2, ib)/pwr1)
2111  END IF
2112  taucld2d = dp(k)*1.0e3*(hydrometsd(k, 2)*(awb_ir1(1, ib)+(awb_ir1(2&
2113 & , ib)+(awb_ir1(3, ib)+awb_ir1(4, ib)*reff(k, 2))*reff(k, 2))*reff(&
2114 & k, 2))+hydromets(k, 2)*((awb_ir1(4, ib)*reffd(k, 2)*reff(k, 2)+(&
2115 & awb_ir1(3, ib)+awb_ir1(4, ib)*reff(k, 2))*reffd(k, 2))*reff(k, 2)+&
2116 & (awb_ir1(2, ib)+(awb_ir1(3, ib)+awb_ir1(4, ib)*reff(k, 2))*reff(k&
2117 & , 2))*reffd(k, 2)))/cons_grav
2118  taucld2 = dp(k)*1.0e3/cons_grav*hydromets(k, 2)*(awb_ir1(1, ib)+(&
2119 & awb_ir1(2, ib)+(awb_ir1(3, ib)+awb_ir1(4, ib)*reff(k, 2))*reff(k, &
2120 & 2))*reff(k, 2))
2121  taucld3d = 0.00307*dp(k)*1.0e3*hydrometsd(k, 3)/cons_grav
2122  taucld3 = 0.00307*(dp(k)*1.0e3/cons_grav*hydromets(k, 3))
2123  IF (reff(k, 4) .GT. 112.0) THEN
2124  reff_snow = 112.0
2125  reff_snowd = 0.0_8
2126  ELSE
2127  reff_snowd = reffd(k, 4)
2128  reff_snow = reff(k, 4)
2129  END IF
2130  IF (reff_snow .LE. 0.0) THEN
2131  taucld4 = 0.0
2132  taucld4d = 0.0_8
2133  ELSE
2134  IF (reff_snow .GT. 0.0 .OR. (reff_snow .LT. 0.0 .AND. aib_ir1(3, &
2135 & ib) .EQ. int(aib_ir1(3, ib)))) THEN
2136  pwr1d = aib_ir1(3, ib)*reff_snow**(aib_ir1(3, ib)-1)*reff_snowd
2137  ELSE IF (reff_snow .EQ. 0.0 .AND. aib_ir1(3, ib) .EQ. 1.0) THEN
2138  pwr1d = reff_snowd
2139  ELSE
2140  pwr1d = 0.0
2141  END IF
2142  pwr1 = reff_snow**aib_ir1(3, ib)
2143  taucld4d = dp(k)*1.0e3*(hydrometsd(k, 4)*(aib_ir1(1, ib)+aib_ir1(2&
2144 & , ib)/pwr1)-hydromets(k, 4)*aib_ir1(2, ib)*pwr1d/pwr1**2)/&
2145 & cons_grav
2146  taucld4 = dp(k)*1.0e3/cons_grav*hydromets(k, 4)*(aib_ir1(1, ib)+&
2147 & aib_ir1(2, ib)/pwr1)
2148  END IF
2149 !-----Compute cloud single-scattering albedo and asymmetry factor for
2150 ! a mixture of ice particles and liquid drops following
2151 ! Eqs. (6.5), (6.6), (6.15) and (6.16).
2152 ! Single-scattering albedo and asymmetry factor of rain are set
2153 ! to 0.54 and 0.95, respectively, based on the information provided
2154 ! by Prof. Qiang Fu.
2155  taucd = taucld1d + taucld2d + taucld3d + taucld4d
2156  tauc = taucld1 + taucld2 + taucld3 + taucld4
2157  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
2158  w1d = taucld1d*(aiw_ir1(1, ib)+(aiw_ir1(2, ib)+(aiw_ir1(3, ib)+&
2159 & aiw_ir1(4, ib)*reff(k, 1))*reff(k, 1))*reff(k, 1)) + taucld1*((&
2160 & aiw_ir1(4, ib)*reffd(k, 1)*reff(k, 1)+(aiw_ir1(3, ib)+aiw_ir1(4&
2161 & , ib)*reff(k, 1))*reffd(k, 1))*reff(k, 1)+(aiw_ir1(2, ib)+(&
2162 & aiw_ir1(3, ib)+aiw_ir1(4, ib)*reff(k, 1))*reff(k, 1))*reffd(k, 1&
2163 & ))
2164  w1 = taucld1*(aiw_ir1(1, ib)+(aiw_ir1(2, ib)+(aiw_ir1(3, ib)+&
2165 & aiw_ir1(4, ib)*reff(k, 1))*reff(k, 1))*reff(k, 1))
2166  w2d = taucld2d*(aww_ir1(1, ib)+(aww_ir1(2, ib)+(aww_ir1(3, ib)+&
2167 & aww_ir1(4, ib)*reff(k, 2))*reff(k, 2))*reff(k, 2)) + taucld2*((&
2168 & aww_ir1(4, ib)*reffd(k, 2)*reff(k, 2)+(aww_ir1(3, ib)+aww_ir1(4&
2169 & , ib)*reff(k, 2))*reffd(k, 2))*reff(k, 2)+(aww_ir1(2, ib)+(&
2170 & aww_ir1(3, ib)+aww_ir1(4, ib)*reff(k, 2))*reff(k, 2))*reffd(k, 2&
2171 & ))
2172  w2 = taucld2*(aww_ir1(1, ib)+(aww_ir1(2, ib)+(aww_ir1(3, ib)+&
2173 & aww_ir1(4, ib)*reff(k, 2))*reff(k, 2))*reff(k, 2))
2174  w3d = 0.54*taucld3d
2175  w3 = taucld3*0.54
2176  w4d = taucld4d*(aiw_ir1(1, ib)+(aiw_ir1(2, ib)+(aiw_ir1(3, ib)+&
2177 & aiw_ir1(4, ib)*reff_snow)*reff_snow)*reff_snow) + taucld4*((&
2178 & aiw_ir1(4, ib)*reff_snowd*reff_snow+(aiw_ir1(3, ib)+aiw_ir1(4, &
2179 & ib)*reff_snow)*reff_snowd)*reff_snow+(aiw_ir1(2, ib)+(aiw_ir1(3&
2180 & , ib)+aiw_ir1(4, ib)*reff_snow)*reff_snow)*reff_snowd)
2181  w4 = taucld4*(aiw_ir1(1, ib)+(aiw_ir1(2, ib)+(aiw_ir1(3, ib)+&
2182 & aiw_ir1(4, ib)*reff_snow)*reff_snow)*reff_snow)
2183  wwd = ((w1d+w2d+w3d+w4d)*tauc-(w1+w2+w3+w4)*taucd)/tauc**2
2184  ww = (w1+w2+w3+w4)/tauc
2185  g1d = w1d*(aig_ir1(1, ib)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(&
2186 & 4, ib)*reff(k, 1))*reff(k, 1))*reff(k, 1)) + w1*((aig_ir1(4, ib)&
2187 & *reffd(k, 1)*reff(k, 1)+(aig_ir1(3, ib)+aig_ir1(4, ib)*reff(k, 1&
2188 & ))*reffd(k, 1))*reff(k, 1)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+&
2189 & aig_ir1(4, ib)*reff(k, 1))*reff(k, 1))*reffd(k, 1))
2190  g1 = w1*(aig_ir1(1, ib)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(4&
2191 & , ib)*reff(k, 1))*reff(k, 1))*reff(k, 1))
2192  g2d = w2d*(awg_ir1(1, ib)+(awg_ir1(2, ib)+(awg_ir1(3, ib)+awg_ir1(&
2193 & 4, ib)*reff(k, 2))*reff(k, 2))*reff(k, 2)) + w2*((awg_ir1(4, ib)&
2194 & *reffd(k, 2)*reff(k, 2)+(awg_ir1(3, ib)+awg_ir1(4, ib)*reff(k, 2&
2195 & ))*reffd(k, 2))*reff(k, 2)+(awg_ir1(2, ib)+(awg_ir1(3, ib)+&
2196 & awg_ir1(4, ib)*reff(k, 2))*reff(k, 2))*reffd(k, 2))
2197  g2 = w2*(awg_ir1(1, ib)+(awg_ir1(2, ib)+(awg_ir1(3, ib)+awg_ir1(4&
2198 & , ib)*reff(k, 2))*reff(k, 2))*reff(k, 2))
2199  g3d = 0.95*w3d
2200  g3 = w3*0.95
2201  g4d = w4d*(aig_ir1(1, ib)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(&
2202 & 4, ib)*reff_snow)*reff_snow)*reff_snow) + w4*((aig_ir1(4, ib)*&
2203 & reff_snowd*reff_snow+(aig_ir1(3, ib)+aig_ir1(4, ib)*reff_snow)*&
2204 & reff_snowd)*reff_snow+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(4&
2205 & , ib)*reff_snow)*reff_snow)*reff_snowd)
2206  g4 = w4*(aig_ir1(1, ib)+(aig_ir1(2, ib)+(aig_ir1(3, ib)+aig_ir1(4&
2207 & , ib)*reff_snow)*reff_snow)*reff_snow)
2208  IF (w1 + w2 + w3 + w4 .GE. 0.) THEN
2209  abs0 = w1 + w2 + w3 + w4
2210  ELSE
2211  abs0 = -(w1+w2+w3+w4)
2212  END IF
2213  IF (abs0 .GT. 0.0) THEN
2214  ggd = ((g1d+g2d+g3d+g4d)*(w1+w2+w3+w4)-(g1+g2+g3+g4)*(w1d+w2d+&
2215 & w3d+w4d))/(w1+w2+w3+w4)**2
2216  gg = (g1+g2+g3+g4)/(w1+w2+w3+w4)
2217  ELSE
2218  gg = 0.5
2219  ggd = 0.0_8
2220  END IF
2221 !-----Parameterization of LW scattering following Eqs. (6.11)
2222 ! and (6.12).
2223  ffd = (0.1185*ggd*gg+(0.0076+0.1185*gg)*ggd)*gg + (0.3739+(0.0076+&
2224 & 0.1185*gg)*gg)*ggd
2225  ff = 0.5 + (0.3739+(0.0076+0.1185*gg)*gg)*gg
2226  IF (1. - ww*ff .LT. 0.0) THEN
2227  max1 = 0.0
2228  max1d = 0.0_8
2229  ELSE
2230  max1d = -(wwd*ff) - ww*ffd
2231  max1 = 1. - ww*ff
2232  END IF
2233 !ALT: temporary protection against negative cloud optical thickness
2234  taucd = max1d*tauc + max1*taucd
2235  tauc = max1*tauc
2236 !-----compute cloud diffuse transmittance. It is approximated by using
2237 ! a diffusivity factor of 1.66.
2238  tcldlyrd(k) = -(1.66*taucd*exp(-(1.66*tauc)))
2239  tcldlyr(k) = exp(-(1.66*tauc))
2240 ! N in the documentation (6.13)
2241  ennd(k) = fcldd(k)*(1.0-tcldlyr(k)) - fcld(k)*tcldlyrd(k)
2242  enn(k) = fcld(k)*(1.0-tcldlyr(k))
2243  ELSE
2244  tcldlyrd(k) = 0.0_8
2245  tcldlyr(k) = 1.0
2246  ennd(k) = 0.0_8
2247  enn(k) = 0.0
2248  END IF
2249  END DO
2250 END SUBROUTINE getirtau1_d
2251 
2252 end module irrad_tl
irrad_tl::n2okdis_d
subroutine n2okdis_d(ib, np, k, n2oexp, n2oexpd, tn2o, tn2od, tran, trand)
Definition: irrad_tl.F90:1646
irrad_tl::h2oexps_d
subroutine h2oexps_d(ib, np, dh2o, dh2od, pa, dt, dtd, xkw, aw, bw, pm, mw, h2oexp, h2oexpd)
Definition: irrad_tl.F90:977
irrad_tl::planck_d
subroutine planck_d(ibn, cb, t, td, xlayer, xlayerd)
Definition: irrad_tl.F90:952
irrad_tl::tablup_d
subroutine tablup_d(nx1, nh1, dw, dwd, p, dt, dtd, s1, s1d, s2, s2d, s3, s3d, w1, p1, dwe, dpe, coef1, coef2, coef3, tran, trand)
Definition: irrad_tl.F90:1415
type_kinds::ip
integer, parameter, public ip
Definition: Type_Kinds.f90:97
irrad_tl::getirtau1_d
subroutine getirtau1_d(ib, nlevs, dp, fcld, fcldd, reff, reffd, hydromets, hydrometsd, tcldlyr, tcldlyrd, enn, ennd, aib_ir1, awb_ir1, aiw_ir1, aww_ir1, aig_ir1, awg_ir1, cons_grav)
Definition: irrad_tl.F90:2019
irrad_tl::irrad_d
subroutine, public irrad_d(np, ple_dev, ta_dev, ta_devd, wa_dev, wa_devd, oa_dev, oa_devd, tb_dev, tb_devd, co2, trace, n2o_dev, ch4_dev, cfc11_dev, cfc12_dev, cfc22_dev, cwc_dev, cwc_devd, fcld_dev, fcld_devd, ict, icb, reff_dev, reff_devd, ns, fs_dev, tg_dev, eg_dev, tv_dev, ev_dev, rv_dev, na, nb, taua_dev, taua_devd, ssaa_dev, ssaa_devd, asya_dev, asya_devd, flxu_dev, flxu_devd, flxd_dev, flxd_devd, dfdts_dev, dfdts_devd, 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_tl.F90:22
irrad_tl::n2oexps_d
subroutine n2oexps_d(ib, np, dn2o, pa, dt, dtd, n2oexp, n2oexpd)
Definition: irrad_tl.F90:1088
irrad_tl::comkdis_d
subroutine comkdis_d(ib, np, k, comexp, comexpd, tcom, tcomd, tran, trand)
Definition: irrad_tl.F90:1748
irrad_tl::ch4kdis_d
subroutine ch4kdis_d(ib, np, k, ch4exp, ch4expd, tch4, tch4d, tran, trand)
Definition: irrad_tl.F90:1699
irrad_tl::cldovlp_d
subroutine cldovlp_d(np, k1, k2, ict, icb, icx, ncld, enn, ennd, ett, ettd, cldhi, cldhid, cldmd, cldmdd, cldlw, cldlwd)
Definition: irrad_tl.F90:1935
irradmod::sfcflux
subroutine, public sfcflux(ibn, cb, dcb, ns, fs, tg, eg, tv, ev, rv, bs, dbs, rflxs)
Definition: irrad.F90:2067
irrad_tl::cfckdis_d
subroutine cfckdis_d(np, k, cfcexp, cfcexpd, tcfc, tcfcd, tran, trand)
Definition: irrad_tl.F90:1824
irrad_tl::h2okdis_d
subroutine h2okdis_d(ib, np, k, fkw, gkw, ne, h2oexp, h2oexpd, conexp, conexpd, th2o, th2od, tcon, tcond, tran, trand)
Definition: irrad_tl.F90:1553
irrad_tl::b10kdis_d
subroutine b10kdis_d(np, k, h2oexp, h2oexpd, conexp, conexpd, co2exp, co2expd, n2oexp, n2oexpd, th2o, th2od, tcon, tcond, tco2, tco2d, tn2o, tn2od, tran, trand)
Definition: irrad_tl.F90:1847
irrad_tl
Definition: irrad_tl.F90:1
irrad_tl::cfcexps_d
subroutine cfcexps_d(ib, np, a1, b1, fk1, a2, b2, fk2, dcfc, dt, dtd, cfcexp, cfcexpd)
Definition: irrad_tl.F90:1258
max
#define max(a, b)
Definition: mosaic_util.h:33
irrad_tl::conexps_d
subroutine conexps_d(ib, np, dcont, dcontd, xke, conexp, conexpd)
Definition: irrad_tl.F90:1054
irrad_tl::comexps_d
subroutine comexps_d(ib, np, dcom, dt, dtd, comexp, comexpd)
Definition: irrad_tl.F90:1213
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
irrad_tl::b10exps_d
subroutine b10exps_d(np, dh2o, dh2od, dcont, dcontd, dco2, dn2o, pa, dt, dtd, h2oexp, h2oexpd, conexp, conexpd, co2exp, co2expd, n2oexp, n2oexpd)
Definition: irrad_tl.F90:1297
irrad_tl::ch4exps_d
subroutine ch4exps_d(ib, np, dch4, pa, dt, dtd, ch4exp, ch4expd)
Definition: irrad_tl.F90:1154