FV3 Bundle
sorad_ad.F90
Go to the documentation of this file.
1 module sorad_ad
2 
3 USE soradmod
4 
5 IMPLICIT NONE
6 
7 PRIVATE
8 PUBLIC :: sorad_b
9 
10 contains
11 
12 SUBROUTINE sorad_b(m, np, nb, cosz_dev, pl_dev, ta_dev, ta_devb, wa_dev&
13 & , wa_devb, oa_dev, oa_devb, co2, cwc_dev, cwc_devb, fcld_dev, &
14 & fcld_devb, ict, icb, reff_dev, reff_devb, hk_uv, hk_ir, taua_dev, &
15 & taua_devb, ssaa_dev, ssaa_devb, asya_dev, asya_devb, rsuvbm_dev, &
16 & rsuvdf_dev, rsirbm_dev, rsirdf_dev, flx_dev, flx_devb, cons_grav, &
17 & wk_uv, zk_uv, ry_uv, xk_ir, ry_ir, cah, coa, aig_uv, awg_uv, arg_uv, &
18 & aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, &
19 & ara_nir, aig_nir, awg_nir, arg_nir, caib, caif)
20  IMPLICIT NONE
21 !end do RUN_LOOP
22 ! Parameters
23 ! ----------
24  INTEGER, PARAMETER :: nu=43
25  INTEGER, PARAMETER :: nw=37
26  INTEGER, PARAMETER :: nx=62
27  INTEGER, PARAMETER :: ny=101
28  INTEGER, PARAMETER :: nband_uv=5
29  INTEGER, PARAMETER :: nk_ir=10
30  INTEGER, PARAMETER :: nband_ir=3
31  INTEGER, PARAMETER :: nband=nband_uv+nband_ir
32  REAL*8, PARAMETER :: dsm=0.602
33 !-----input values
34 !-----input parameters
35  INTEGER :: m, np, ict, icb, nb
36  REAL*8 :: cosz_dev(m), pl_dev(m, np+1), ta_dev(m, np), wa_dev(m, np), &
37 & oa_dev(m, np), co2
38  REAL*8 :: ta_devb(m, np), wa_devb(m, np), oa_devb(m, np)
39  REAL*8 :: cwc_dev(m, np, 4), fcld_dev(m, np), reff_dev(m, np, 4), &
40 & hk_uv(5), hk_ir(3, 10)
41  REAL*8 :: cwc_devb(m, np, 4), fcld_devb(m, np), reff_devb(m, np, 4)
42  REAL*8 :: rsuvbm_dev, rsuvdf_dev, rsirbm_dev, rsirdf_dev
43  REAL*8 :: taua_dev(m, np, nb)
44  REAL*8 :: taua_devb(m, np, nb)
45  REAL*8 :: ssaa_dev(m, np, nb)
46  REAL*8 :: ssaa_devb(m, np, nb)
47  REAL*8 :: asya_dev(m, np, nb)
48  REAL*8 :: asya_devb(m, np, nb)
49  LOGICAL :: overcast
50 ! Constants
51  REAL*8, INTENT(IN) :: wk_uv(5), zk_uv(5), ry_uv(5)
52  REAL*8, INTENT(IN) :: xk_ir(10), ry_ir(3)
53  REAL*8, INTENT(IN) :: cah(43, 37), coa(62, 101)
54  REAL*8, INTENT(IN) :: aig_uv(3), awg_uv(3), arg_uv(3)
55  REAL*8, INTENT(IN) :: aib_uv, awb_uv(2), arb_uv(2)
56  REAL*8, INTENT(IN) :: aib_nir, awb_nir(3, 2), arb_nir(3, 2)
57  REAL*8, INTENT(IN) :: aia_nir(3, 3), awa_nir(3, 3), ara_nir(3, 3)
58  REAL*8, INTENT(IN) :: aig_nir(3, 3), awg_nir(3, 3), arg_nir(3, 3)
59  REAL*8, INTENT(IN) :: caib(11, 9, 11), caif(9, 11)
60  REAL*8, INTENT(IN) :: cons_grav
61 !-----output parameters
62  REAL*8 :: flx_dev(m, np+1), flc_dev(m, np+1)
63  REAL*8 :: flx_devb(m, np+1)
64  REAL*8 :: flxu_dev(m, np+1), flcu_dev(m, np+1)
65  REAL*8 :: fdiruv_dev(m), fdifuv_dev(m)
66  REAL*8 :: fdirpar_dev(m), fdifpar_dev(m)
67  REAL*8 :: fdirir_dev(m), fdifir_dev(m)
68  REAL*8 :: flx_sfc_band_dev(m, nband)
69 !-----temporary arrays
70  INTEGER :: i, j, k, l, in, ntop
71  REAL*8 :: dp(np), wh(np), oh(np)
72  REAL*8 :: whb(np), ohb(np)
73  REAL*8 :: scal(np)
74  REAL*8 :: swh(np+1), so2(np+1), df(0:np+1)
75  REAL*8 :: swhb(np+1), dfb(0:np+1)
76  REAL*8 :: scal0, wvtoa, o3toa, pa
77  REAL*8 :: wvtoab, o3toab
78  REAL*8 :: snt, cnt, x, xx4, xtoa
79  REAL*8 :: xx4b
80  REAL*8 :: dp_pa(np)
81 !-----parameters for co2 transmission tables
82  REAL*8 :: w1, dw, u1, du
83  INTEGER :: ib, rc
84  REAL*8 :: tauclb(np), tauclf(np), asycl(np)
85  REAL*8 :: tauclbb(np), tauclfb(np), asyclb(np)
86  REAL*8 :: taubeam(np, 4), taudiff(np, 4)
87  REAL*8 :: taubeamb(np, 4), taudiffb(np, 4)
88  REAL*8 :: fcld_col(np)
89  REAL*8 :: fcld_colb(np)
90  REAL*8 :: cwc_col(np, 4)
91  REAL*8 :: cwc_colb(np, 4)
92  REAL*8 :: reff_col(np, 4)
93  REAL*8 :: reff_colb(np, 4)
94  REAL*8 :: taurs, tauoz, tauwv
95  REAL*8 :: tauozb, tauwvb
96  REAL*8 :: tausto, ssatau, asysto
97  REAL*8 :: taustob, ssataub, asystob
98  REAL*8 :: tautob, ssatob, asytob
99  REAL*8 :: tautobb, ssatobb, asytobb
100  REAL*8 :: tautof, ssatof, asytof
101  REAL*8 :: tautofb, ssatofb, asytofb
102  REAL*8 :: rr(0:np+1, 2), tt(0:np+1, 2), td(0:np+1, 2)
103  REAL*8 :: rrb(0:np+1, 2), ttb(0:np+1, 2), tdb(0:np+1, 2)
104  REAL*8 :: rs(0:np+1, 2), ts(0:np+1, 2)
105  REAL*8 :: rsb(0:np+1, 2), tsb(0:np+1, 2)
106  REAL*8 :: fall(np+1), fclr(np+1), fsdir, fsdif
107  REAL*8 :: fallb(np+1)
108  REAL*8 :: fupa(np+1), fupc(np+1)
109  REAL*8 :: cc1, cc2, cc3
110  REAL*8 :: cc1b, cc2b, cc3b
111  REAL*8 :: rrt, ttt, tdt, rst, tst
112  REAL*8 :: rrtb, tttb, tdtb, rstb, tstb
113  INTEGER :: iv, ik
114  REAL*8 :: ssacl(np)
115  REAL*8 :: ssaclb(np)
116  INTEGER :: im
117  INTEGER :: ic, iw
118  REAL*8 :: ulog, wlog, dc, dd, x0, x1, x2, y0, y1, y2, du2, dw2
119  REAL*8 :: wlogb, ddb, x2b, y2b
120  INTEGER :: ih
121 !if (overcast == true) then
122 !real(8) :: rra(0:np+1),rxa(0:np+1)
123 !real(8) :: ttaold,tdaold,rsaold
124 !real(8) :: ttanew,tdanew,rsanew
125 !else
126  REAL*8 :: rra(0:np+1, 2, 2), tta(0:np, 2, 2)
127  REAL*8 :: rrab(0:np+1, 2, 2), ttab(0:np, 2, 2)
128  REAL*8 :: tda(0:np, 2, 2)
129  REAL*8 :: tdab(0:np, 2, 2)
130  REAL*8 :: rsa(0:np, 2, 2), rxa(0:np+1, 2, 2)
131  REAL*8 :: rsab(0:np, 2, 2), rxab(0:np+1, 2, 2)
132 !endif
133  REAL*8 :: flxdn
134  REAL*8 :: flxdnb
135  REAL*8 :: fdndir, fdndif, fupdif
136  REAL*8 :: fdndirb, fdndifb, fupdifb
137  REAL*8 :: denm, yy
138  REAL*8 :: denmb, yyb
139  INTEGER :: is
140  REAL*8 :: ch, cm, ct
141  REAL*8 :: chb, cmb, ctb
142  INTEGER :: foundtop
143  REAL*8 :: dftop
144  REAL*8 :: dftopb
145 !-----Variables for aerosols
146  INTEGER :: ii, jj, irhp1, an
147  REAL*8 :: dum
148  REAL*8 :: dumb
149  INTRINSIC max
150  INTRINSIC exp
151  INTRINSIC min
152  INTRINSIC sqrt
153  INTRINSIC real
154  INTRINSIC log10
155  INTRINSIC int
156  INTRINSIC abs
157  INTRINSIC epsilon
158  REAL :: result1
159  INTEGER :: branch
160  REAL*8 :: temp3
161  REAL*8 :: temp29
162  REAL*8 :: tempb52
163  REAL*8 :: temp2
164  REAL*8 :: temp28
165  REAL*8 :: tempb51
166  REAL*8 :: temp1
167  REAL*8 :: temp27
168  REAL*8 :: tempb50
169  REAL*8 :: temp0
170  REAL*8 :: temp26
171  REAL*8 :: temp25
172  REAL*8 :: temp24
173  REAL*8 :: temp23
174  REAL*8 :: temp22
175  REAL*8 :: temp21
176  REAL*8 :: temp20
177  REAL*8 :: tempb9
178  REAL*8 :: tempb8
179  REAL*8 :: tempb7
180  REAL*8 :: tempb6
181  REAL*8 :: tempb5
182  REAL*8 :: tempb4
183  REAL*8 :: tempb19
184  REAL*8 :: tempb3
185  REAL*8 :: tempb18
186  REAL*8 :: tempb2
187  REAL*8 :: tempb17
188  REAL*8 :: tempb1
189  REAL*8 :: tempb16
190  REAL*8 :: tempb0
191  REAL*8 :: tempb15
192  REAL*8 :: tempb14
193  REAL*8 :: tempb13
194  REAL*8 :: tempb12
195  REAL*8 :: tempb49
196  REAL*8 :: x6
197  REAL*8 :: tempb11
198  REAL*8 :: tempb48
199  REAL*8 :: x5
200  REAL*8 :: tempb10
201  REAL*8 :: tempb47
202  REAL*8 :: x4
203  REAL*8 :: tempb46
204  REAL*8 :: x3
205  REAL*8 :: tempb45
206  REAL*8 :: tempb44
207  REAL*8 :: tempb43
208  REAL*8 :: temp19
209  REAL*8 :: tempb42
210  REAL*8 :: temp18
211  REAL*8 :: tempb41
212  REAL*8 :: temp17
213  REAL*8 :: tempb40
214  REAL*8 :: temp16
215  REAL*8 :: temp15
216  REAL*8 :: temp14
217  REAL*8 :: temp13
218  REAL*8 :: temp12
219  REAL*8 :: temp11
220  REAL*8 :: temp10
221  REAL*8 :: tempb
222  REAL*8 :: tempb39
223  REAL*8 :: tempb38
224  REAL*8 :: tempb37
225  REAL*8 :: tempb36
226  REAL*8 :: tempb35
227  REAL*8 :: tempb34
228  REAL*8 :: tempb33
229  REAL*8 :: tempb32
230  REAL*8 :: tempb31
231  REAL*8 :: tempb30
232  REAL*8 :: x4b
233  REAL*8 :: temp38
234  REAL*8 :: temp37
235  REAL*8 :: temp36
236  REAL*8 :: temp35
237  REAL*8 :: temp34
238  REAL*8 :: temp33
239  REAL*8 :: temp32
240  REAL*8 :: temp31
241  REAL*8 :: temp30
242  REAL*8 :: abs0
243  REAL*8 :: tempb29
244  REAL*8 :: tempb28
245  REAL*8 :: tempb27
246  REAL*8 :: tempb26
247  REAL*8 :: tempb25
248  REAL*8 :: temp
249  REAL*8 :: tempb24
250  REAL*8 :: tempb23
251  REAL*8 :: tempb22
252  REAL*8 :: temp9
253  REAL*8 :: tempb21
254  REAL*8 :: temp8
255  REAL*8 :: tempb20
256  REAL*8 :: temp7
257  REAL*8 :: tempb56
258  REAL*8 :: temp6
259  REAL*8 :: tempb55
260  REAL*8 :: temp5
261  REAL*8 :: tempb54
262  REAL*8 :: temp4
263  REAL*8 :: tempb53
264  i = 1
265 !RUN_LOOP: do i=1,m
266  ntop = 0
267 !-----Beginning of sorad code
268 !-----wvtoa and o3toa are the water vapor and o3 amounts of the region
269 ! above the pl(1) level.
270 ! snt is the secant of the solar zenith angle
271  snt = 1.0/cosz_dev(i)
272  IF (pl_dev(i, 1) .LT. 1.e-3) THEN
273  xtoa = 1.e-3
274  ELSE
275  xtoa = pl_dev(i, 1)
276  END IF
277  scal0 = xtoa*(0.5*xtoa/300.)**.8
278  o3toa = 1.02*oa_dev(i, 1)*xtoa*466.7 + 1.0e-8
279  wvtoa = 1.02*wa_dev(i, 1)*scal0*(1.0+0.00135*(ta_dev(i, 1)-240.)) + &
280 & 1.0e-9
281  swh(1) = wvtoa
282  DO k=1,np
283 !-----compute layer thickness. indices for the surface level and
284 ! surface layer are np+1 and np, respectively.
285  dp(k) = pl_dev(i, k+1) - pl_dev(i, k)
286 ! dp in pascals
287  dp_pa(k) = dp(k)*100.
288 !-----compute scaled water vapor amount following Eqs. (3.3) and (3.5)
289 ! unit is g/cm**2
290 !
291  pa = 0.5*(pl_dev(i, k)+pl_dev(i, k+1))
292  scal(k) = dp(k)*(pa/300.)**.8
293  wh(k) = 1.02*wa_dev(i, k)*scal(k)*(1.+0.00135*(ta_dev(i, k)-240.)) +&
294 & 1.e-9
295  swh(k+1) = swh(k) + wh(k)
296 !-----compute ozone amount, unit is (cm-atm)stp
297 ! the number 466.7 is the unit conversion factor
298 ! from g/cm**2 to (cm-atm)stp
299  oh(k) = 1.02*oa_dev(i, k)*dp(k)*466.7 + 1.e-8
300 !-----Fill the reff, cwc, and fcld for the column
301  fcld_col(k) = fcld_dev(i, k)
302  DO l=1,4
303  reff_col(k, l) = reff_dev(i, k, l)
304  cwc_col(k, l) = cwc_dev(i, k, l)
305  END DO
306  END DO
307 !-----Initialize temporary arrays to zero to avoid UMR
308  rr = 0.0
309  tt = 0.0
310  td = 0.0
311  rs = 0.0
312  ts = 0.0
313  rra = 0.0
314  rxa = 0.0
315 !if( OVERCAST == .false. ) then
316  tta = 0.0
317  tda = 0.0
318  rsa = 0.0
319 !endif
320 !-----initialize fluxes for all-sky (flx), clear-sky (flc), and
321 ! flux reduction (df)
322 !
323  DO k=1,np+1
324  flx_dev(i, k) = 0.
325  END DO
326 !-----Begin inline of SOLUV
327 !-----compute solar uv and par fluxes
328 !-----initialize fdiruv, fdifuv, surface reflectances and transmittances.
329 ! the reflectance and transmittance of the clear and cloudy portions
330 ! of a layer are denoted by 1 and 2, respectively.
331 ! cc is the maximum cloud cover in each of the high, middle, and low
332 ! cloud groups.
333 ! 1/dsm=1/cos(53) = 1.66
334  rr(np+1, 1) = rsuvbm_dev
335  rr(np+1, 2) = rsuvbm_dev
336  rs(np+1, 1) = rsuvdf_dev
337  rs(np+1, 2) = rsuvdf_dev
338  td(np+1, 1) = 0.0
339  td(np+1, 2) = 0.0
340  tt(np+1, 1) = 0.0
341  tt(np+1, 2) = 0.0
342  ts(np+1, 1) = 0.0
343  ts(np+1, 2) = 0.0
344  rr(0, 1) = 0.0
345  rr(0, 2) = 0.0
346  rs(0, 1) = 0.0
347  rs(0, 2) = 0.0
348 ! td(0,1)=1.0
349 ! td(0,2)=1.0
350  tt(0, 1) = 1.0
351  tt(0, 2) = 1.0
352  ts(0, 1) = 1.0
353  ts(0, 2) = 1.0
354  cc1 = 0.0
355  cc2 = 0.0
356  cc3 = 0.0
357 !-----options for scaling cloud optical thickness
358 !if ( OVERCAST == .true. ) then
359 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10) and
360 ! compute cloud asymmetry factor
361 ! Note: the cloud optical properties are assumed to be independent
362 ! of spectral bands in the UV and PAR regions.
363 ! for a mixture of ice and liquid particles.
364 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
365 !
366 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
367 ! respectively.
368 ! call getvistau1(np,cosz_dev(i),dp_pa,fcld_col,reff_col,cwc_col,0,0,&
369 ! taubeam,taudiff,asycl, &
370 ! aig_uv, awg_uv, arg_uv, &
371 ! aib_uv, awb_uv, arb_uv, &
372 ! aib_nir, awb_nir, arb_nir, &
373 ! aia_nir, awa_nir, ara_nir, &
374 ! aig_nir, awg_nir, arg_nir, &
375 ! caib, caif, &
376 ! CONS_GRAV )
377 !else
378 !-----Compute scaled cloud optical thickness. Eqs. (4.6) and (4.10) and
379 ! compute cloud asymmetry factor
380 ! Note: the cloud optical properties are assumed to be independent
381 ! of spectral bands in the UV and PAR regions.
382 ! for a mixture of ice and liquid particles.
383 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
384 !
385 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
386 ! respectively.
387  CALL getvistau1(np, cosz_dev(i), dp_pa, fcld_col, reff_col, cwc_col, &
388 & ict, icb, taubeam, taudiff, asycl, aig_uv, awg_uv, arg_uv, &
389 & aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, &
390 & awa_nir, ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, &
391 & cons_grav)
392 !-----clouds within each of the high, middle, and low clouds are
393 ! assumed to be maximally overlapped, and the cloud cover (cc)
394 ! for a group (high, middle, or low) is the maximum cloud cover
395 ! of all the layers within a group
396 ! The cc1,2,3 are still needed in the flux calculations below
397 !MAT---DO NOT FUSE THIS LOOP
398 !MAT---Loop must run to completion so that cc[1,2,3] are correct.
399  DO k=1,np
400  IF (k .LT. ict) THEN
401  IF (cc1 .LT. fcld_dev(i, k)) THEN
402  cc1 = fcld_dev(i, k)
403  CALL pushcontrol3b(4)
404  ELSE
405  CALL pushcontrol3b(5)
406  cc1 = cc1
407  END IF
408  ELSE IF (k .LT. icb) THEN
409  IF (cc2 .LT. fcld_dev(i, k)) THEN
410  cc2 = fcld_dev(i, k)
411  CALL pushcontrol3b(2)
412  ELSE
413  CALL pushcontrol3b(3)
414  cc2 = cc2
415  END IF
416  ELSE IF (cc3 .LT. fcld_dev(i, k)) THEN
417  cc3 = fcld_dev(i, k)
418  CALL pushcontrol3b(0)
419  ELSE
420  CALL pushcontrol3b(1)
421  cc3 = cc3
422  END IF
423  END DO
424 !MAT---DO NOT FUSE THIS LOOP
425 !endif !overcast
426  DO k=1,np
427  tauclb(k) = taubeam(k, 1) + taubeam(k, 2) + taubeam(k, 3) + taubeam(&
428 & k, 4)
429  tauclf(k) = taudiff(k, 1) + taudiff(k, 2) + taudiff(k, 3) + taudiff(&
430 & k, 4)
431  END DO
432 !-----integration over spectral bands
433 !-----Compute optical thickness, single-scattering albedo and asymmetry
434 ! factor for a mixture of "na" aerosol types. [Eqs. (4.16)-(4.18)]
435  DO ib=1,nband_uv
436 !-----compute direct beam transmittances of the layer above pl(1)
437  CALL pushreal8(td(0, 1))
438  td(0, 1) = exp(-((wvtoa*wk_uv(ib)+o3toa*zk_uv(ib))/cosz_dev(i)))
439  CALL pushreal8(td(0, 2))
440  td(0, 2) = td(0, 1)
441  DO k=1,np
442 !-----compute clear-sky optical thickness, single scattering albedo,
443 ! and asymmetry factor (Eqs. 6.2-6.4)
444  taurs = ry_uv(ib)*dp(k)
445  tauoz = zk_uv(ib)*oh(k)
446  tauwv = wk_uv(ib)*wh(k)
447  tausto = taurs + tauoz + tauwv + taua_dev(i, k, ib) + 1.0e-7
448  ssatau = ssaa_dev(i, k, ib) + taurs
449  asysto = asya_dev(i, k, ib)
450  CALL pushreal8(tautob)
451  tautob = tausto
452  asytob = asysto/ssatau
453  CALL pushreal8(ssatob)
454  ssatob = ssatau/tautob + 1.0e-8
455  IF (ssatob .GT. 0.999999) THEN
456  ssatob = 0.999999
457  CALL pushcontrol1b(0)
458  ELSE
459  CALL pushcontrol1b(1)
460  ssatob = ssatob
461  END IF
462 !-----for direct incident radiation
463  CALL deledd(tautob, ssatob, asytob, cosz_dev(i), rrt, ttt, tdt)
464 !-----diffuse incident radiation is approximated by beam radiation with
465 ! an incident angle of 53 degrees, Eqs. (6.5) and (6.6)
466  CALL deledd(tautob, ssatob, asytob, dsm, rst, tst, dum)
467  CALL pushreal8(rr(k, 1))
468  rr(k, 1) = rrt
469  CALL pushreal8(tt(k, 1))
470  tt(k, 1) = ttt
471  CALL pushreal8(td(k, 1))
472  td(k, 1) = tdt
473  CALL pushreal8(rs(k, 1))
474  rs(k, 1) = rst
475  CALL pushreal8(ts(k, 1))
476  ts(k, 1) = tst
477 !-----compute reflectance and transmittance of the cloudy portion
478 ! of a layer
479 !-----for direct incident radiation
480 ! The effective layer optical properties. Eqs. (6.2)-(6.4)
481  tautob = tausto + tauclb(k)
482  CALL pushreal8(ssatob)
483  ssatob = (ssatau+tauclb(k))/tautob + 1.0e-8
484  IF (ssatob .GT. 0.999999) THEN
485  ssatob = 0.999999
486  CALL pushcontrol1b(0)
487  ELSE
488  CALL pushcontrol1b(1)
489  ssatob = ssatob
490  END IF
491  asytob = (asysto+asycl(k)*tauclb(k))/(ssatob*tautob)
492 !-----for diffuse incident radiation
493  CALL pushreal8(tautof)
494  tautof = tausto + tauclf(k)
495  CALL pushreal8(ssatof)
496  ssatof = (ssatau+tauclf(k))/tautof + 1.0e-8
497  IF (ssatof .GT. 0.999999) THEN
498  ssatof = 0.999999
499  CALL pushcontrol1b(0)
500  ELSE
501  ssatof = ssatof
502  CALL pushcontrol1b(1)
503  END IF
504  asytof = (asysto+asycl(k)*tauclf(k))/(ssatof*tautof)
505 !-----for direct incident radiation
506 ! note that the cloud optical thickness is scaled differently
507 ! for direct and diffuse insolation, Eqs. (7.3) and (7.4).
508  CALL deledd(tautob, ssatob, asytob, cosz_dev(i), rrt, ttt, tdt)
509 !-----diffuse incident radiation is approximated by beam radiation
510 ! with an incident angle of 53 degrees, Eqs. (6.5) and (6.6)
511  CALL deledd(tautof, ssatof, asytof, dsm, rst, tst, dum)
512  CALL pushreal8(rr(k, 2))
513  rr(k, 2) = rrt
514  CALL pushreal8(tt(k, 2))
515  tt(k, 2) = ttt
516  CALL pushreal8(td(k, 2))
517  td(k, 2) = tdt
518  CALL pushreal8(rs(k, 2))
519  rs(k, 2) = rst
520  CALL pushreal8(ts(k, 2))
521  ts(k, 2) = tst
522  END DO
523 !-----flux calculations
524 ! initialize clear-sky flux (fclr), all-sky flux (fall),
525 ! and surface downward fluxes (fsdir and fsdif)
526  DO k=1,np+1
527  fall(k) = 0.0
528  END DO
529 !if ( OVERCAST == .true. ) then
530 !-----Inline CLDFLXY
531 !-----for clear- and all-sky flux calculations when fractional
532 ! cloud cover is either 0 or 1.
533 !-----layers are added one at a time, going up from the surface.
534 ! rra is the composite reflectance illuminated by beam radiation
535 ! rxa is the composite reflectance illuminated from above
536 ! by diffuse radiation
537 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
538 !-----compute transmittances and reflectances for a composite of
539 ! layers. layers are added one at a time, going down from the top.
540 ! tda is the composite direct transmittance illuminated by
541 ! beam radiation
542 ! tta is the composite total transmittance illuminated by
543 ! beam radiation
544 ! rsa is the composite reflectance illuminated from below
545 ! by diffuse radiation
546 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
547 !-----compute fluxes following Eq. (6.15) for fupdif and
548 ! Eq. (6.16) for (fdndir+fdndif)
549 ! fdndir is the direct downward flux
550 ! fdndif is the diffuse downward flux
551 ! fupdif is the diffuse upward flux
552 !-----ih=1 for clear sky; ih=2 for cloudy sky.
553 !-----First set is ih = 1
554 ! rra(np+1)=rr(np+1,1)
555 ! rxa(np+1)=rs(np+1,1)
556 !
557 ! do k=np,0,-1
558 ! denm=ts(k,1)/(1.-rs(k,1)*rxa(k+1))
559 ! rra(k)=rr(k,1)+(td(k,1)*rra(k+1)+(tt(k,1)-td(k,1))*rxa(k+1))*denm
560 ! rxa(k)=rs(k,1)+ts(k,1)*rxa(k+1)*denm
561 ! end do
562 !
563 ! do k=1,np+1
564 ! if (k <= np) then
565 ! if (k == 1) then
566 ! tdaold = td(0,1)
567 ! ttaold = tt(0,1)
568 ! rsaold = rs(0,1)
569 !
570 ! tdanew = 0.0
571 ! ttanew = 0.0
572 ! rsanew = 0.0
573 ! end if
574 ! denm=ts(k,1)/(1.-rsaold*rs(k,1))
575 ! tdanew=tdaold*td(k,1)
576 ! ttanew=tdaold*tt(k,1)+(tdaold*rsaold*rr(k,1)+ttaold-tdaold)*denm
577 ! rsanew=rs(k,1)+ts(k,1)*rsaold*denm
578 ! end if
579 !
580 ! denm=1./(1.-rsaold*rxa(k))
581 ! fdndir=tdaold
582 ! xx4=tdaold*rra(k)
583 ! yy=ttaold-tdaold
584 ! fdndif=(xx4*rsaold+yy)*denm
585 ! fupdif=(xx4+yy*rxa(k))*denm
586 ! flxdn=fdndir+fdndif-fupdif
587 ! fupc(k)=fupdif
588 ! fclr(k)=flxdn
589 !
590 ! tdaold = tdanew
591 ! ttaold = ttanew
592 ! rsaold = rsanew
593 !
594 ! tdanew = 0.0
595 ! ttanew = 0.0
596 ! rsanew = 0.0
597 ! end do
598 !
599 !!-----Second set is ih = 2
600 !
601 ! rra(np+1)=rr(np+1,2)
602 ! rxa(np+1)=rs(np+1,2)
603 !
604 ! do k=np,0,-1
605 ! denm=ts(k,2)/(1.-rs(k,2)*rxa(k+1))
606 ! rra(k)=rr(k,2)+(td(k,2)*rra(k+1)+(tt(k,2)-td(k,2))*rxa(k+1))*denm
607 ! rxa(k)=rs(k,2)+ts(k,2)*rxa(k+1)*denm
608 ! end do
609 !
610 ! do k=1,np+1
611 ! if (k <= np) then
612 ! if (k == 1) then
613 ! tdaold = td(0,2)
614 ! ttaold = tt(0,2)
615 ! rsaold = rs(0,2)
616 ! tdanew = 0.0
617 ! ttanew = 0.0
618 ! rsanew = 0.0
619 ! end if
620 ! denm=ts(k,2)/(1.-rsaold*rs(k,2))
621 ! tdanew=tdaold*td(k,2)
622 ! ttanew=tdaold*tt(k,2)+(tdaold*rsaold*rr(k,2)+ttaold-tdaold)*denm
623 ! rsanew=rs(k,2)+ts(k,2)*rsaold*denm
624 ! end if
625 !
626 ! denm=1./(1.-rsaold*rxa(k))
627 ! fdndir=tdaold
628 ! xx4=tdaold*rra(k)
629 ! yy=ttaold-tdaold
630 ! fdndif=(xx4*rsaold+yy)*denm
631 ! fupdif=(xx4+yy*rxa(k))*denm
632 ! flxdn=fdndir+fdndif-fupdif
633 !
634 ! fupa(k)=fupdif
635 ! fall(k)=flxdn
636 !
637 ! tdaold = tdanew
638 ! ttaold = ttanew
639 ! rsaold = rsanew
640 !
641 ! tdanew = 0.0
642 ! ttanew = 0.0
643 ! rsanew = 0.0
644 ! end do
645 !
646 ! fsdir=fdndir
647 ! fsdif=fdndif
648 !
649 !!-----End CLDFLXY inline
650 !
651 !else
652 !-----for clear- and all-sky flux calculations when fractional
653 ! cloud cover is allowed to be between 0 and 1.
654 ! the all-sky flux, fall is the summation inside the brackets
655 ! of Eq. (7.11)
656 !-----Inline CLDFLX
657 !-----compute transmittances and reflectances for a composite of
658 ! layers. layers are added one at a time, going down from the top.
659 ! tda is the composite direct transmittance illuminated
660 ! by beam radiation
661 ! tta is the composite total transmittance illuminated by
662 ! beam radiation
663 ! rsa is the composite reflectance illuminated from below
664 ! by diffuse radiation
665 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
666 !-----To save memory space, tda, tta, and rsa are pre-computed
667 ! for k<icb. The dimension of these parameters is (m,np,2,2).
668 ! It would have been (m,np,2,2,2) if these parameters were
669 ! computed for all k's.
670 !-----for high clouds
671 ! ih=1 for clear-sky condition, ih=2 for cloudy-sky condition
672  DO ih=1,2
673  CALL pushreal8(tda(0, ih, 1))
674  tda(0, ih, 1) = td(0, ih)
675  CALL pushreal8(tta(0, ih, 1))
676  tta(0, ih, 1) = tt(0, ih)
677  CALL pushreal8(rsa(0, ih, 1))
678  rsa(0, ih, 1) = rs(0, ih)
679  CALL pushreal8(tda(0, ih, 2))
680  tda(0, ih, 2) = td(0, ih)
681  CALL pushreal8(tta(0, ih, 2))
682  tta(0, ih, 2) = tt(0, ih)
683  CALL pushreal8(rsa(0, ih, 2))
684  rsa(0, ih, 2) = rs(0, ih)
685  DO k=1,ict-1
686  CALL pushreal8(denm)
687  denm = ts(k, ih)/(1.-rsa(k-1, ih, 1)*rs(k, ih))
688  CALL pushreal8(tda(k, ih, 1))
689  tda(k, ih, 1) = tda(k-1, ih, 1)*td(k, ih)
690  CALL pushreal8(tta(k, ih, 1))
691  tta(k, ih, 1) = tda(k-1, ih, 1)*tt(k, ih) + (tda(k-1, ih, 1)*rsa&
692 & (k-1, ih, 1)*rr(k, ih)+tta(k-1, ih, 1)-tda(k-1, ih, 1))*denm
693  CALL pushreal8(rsa(k, ih, 1))
694  rsa(k, ih, 1) = rs(k, ih) + ts(k, ih)*rsa(k-1, ih, 1)*denm
695  CALL pushreal8(tda(k, ih, 2))
696  tda(k, ih, 2) = tda(k, ih, 1)
697  CALL pushreal8(tta(k, ih, 2))
698  tta(k, ih, 2) = tta(k, ih, 1)
699  CALL pushreal8(rsa(k, ih, 2))
700  rsa(k, ih, 2) = rsa(k, ih, 1)
701  END DO
702 ! k loop
703 !-----for middle clouds
704 ! im=1 for clear-sky condition, im=2 for cloudy-sky condition
705  DO k=ict,icb-1
706  DO im=1,2
707  CALL pushreal8(denm)
708  denm = ts(k, im)/(1.-rsa(k-1, ih, im)*rs(k, im))
709  CALL pushreal8(tda(k, ih, im))
710  tda(k, ih, im) = tda(k-1, ih, im)*td(k, im)
711  CALL pushreal8(tta(k, ih, im))
712  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, im) + (tda(k-1, ih, im&
713 & )*rsa(k-1, ih, im)*rr(k, im)+tta(k-1, ih, im)-tda(k-1, ih, &
714 & im))*denm
715  CALL pushreal8(rsa(k, ih, im))
716  rsa(k, ih, im) = rs(k, im) + ts(k, im)*rsa(k-1, ih, im)*denm
717  END DO
718  END DO
719  END DO
720 ! im loop
721 ! k loop
722 ! ih loop
723 !-----layers are added one at a time, going up from the surface.
724 ! rra is the composite reflectance illuminated by beam radiation
725 ! rxa is the composite reflectance illuminated from above
726 ! by diffuse radiation
727 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
728 !-----To save memory space, rra and rxa are pre-computed for k>=icb.
729 ! the dimension of these parameters is (m,np,2,2). It would have
730 ! been (m,np,2,2,2) if these parameters were computed for all k's.
731 !-----for the low clouds
732 ! is=1 for clear-sky condition, is=2 for cloudy-sky condition
733  DO is=1,2
734  CALL pushreal8(rra(np+1, 1, is))
735  rra(np+1, 1, is) = rr(np+1, is)
736  CALL pushreal8(rxa(np+1, 1, is))
737  rxa(np+1, 1, is) = rs(np+1, is)
738  CALL pushreal8(rra(np+1, 2, is))
739  rra(np+1, 2, is) = rr(np+1, is)
740  CALL pushreal8(rxa(np+1, 2, is))
741  rxa(np+1, 2, is) = rs(np+1, is)
742  DO k=np,icb,-1
743  CALL pushreal8(denm)
744  denm = ts(k, is)/(1.-rs(k, is)*rxa(k+1, 1, is))
745  CALL pushreal8(rra(k, 1, is))
746  rra(k, 1, is) = rr(k, is) + (td(k, is)*rra(k+1, 1, is)+(tt(k, is&
747 & )-td(k, is))*rxa(k+1, 1, is))*denm
748  CALL pushreal8(rxa(k, 1, is))
749  rxa(k, 1, is) = rs(k, is) + ts(k, is)*rxa(k+1, 1, is)*denm
750  CALL pushreal8(rra(k, 2, is))
751  rra(k, 2, is) = rra(k, 1, is)
752  CALL pushreal8(rxa(k, 2, is))
753  rxa(k, 2, is) = rxa(k, 1, is)
754  END DO
755 ! k loop
756 !-----for middle clouds
757  DO k=icb-1,ict,-1
758  DO im=1,2
759  CALL pushreal8(denm)
760  denm = ts(k, im)/(1.-rs(k, im)*rxa(k+1, im, is))
761  CALL pushreal8(rra(k, im, is))
762  rra(k, im, is) = rr(k, im) + (td(k, im)*rra(k+1, im, is)+(tt(k&
763 & , im)-td(k, im))*rxa(k+1, im, is))*denm
764  CALL pushreal8(rxa(k, im, is))
765  rxa(k, im, is) = rs(k, im) + ts(k, im)*rxa(k+1, im, is)*denm
766  END DO
767  END DO
768  END DO
769 ! im loop
770 ! k loop
771 ! is loop
772 !-----integration over eight sky situations.
773 ! ih, im, is denote high, middle and low cloud groups.
774  DO ih=1,2
775 !-----clear portion
776  IF (ih .EQ. 1) THEN
777  CALL pushreal8(ch)
778  ch = 1.0 - cc1
779 !-----cloudy portion
780  CALL pushcontrol1b(1)
781  ELSE
782  CALL pushreal8(ch)
783  ch = cc1
784  CALL pushcontrol1b(0)
785  END IF
786  DO im=1,2
787 !-----clear portion
788  IF (im .EQ. 1) THEN
789  CALL pushreal8(cm)
790  cm = ch*(1.0-cc2)
791 !-----cloudy portion
792  CALL pushcontrol1b(1)
793  ELSE
794  CALL pushreal8(cm)
795  cm = ch*cc2
796  CALL pushcontrol1b(0)
797  END IF
798  DO is=1,2
799 !-----clear portion
800  IF (is .EQ. 1) THEN
801  CALL pushreal8(ct)
802  ct = cm*(1.0-cc3)
803 !-----cloudy portion
804  CALL pushcontrol1b(1)
805  ELSE
806  CALL pushreal8(ct)
807  ct = cm*cc3
808  CALL pushcontrol1b(0)
809  END IF
810 !-----add one layer at a time, going down.
811  DO k=icb,np
812  CALL pushreal8(denm)
813  denm = ts(k, is)/(1.-rsa(k-1, ih, im)*rs(k, is))
814  CALL pushreal8(tda(k, ih, im))
815  tda(k, ih, im) = tda(k-1, ih, im)*td(k, is)
816  CALL pushreal8(tta(k, ih, im))
817  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, is) + (tda(k-1, ih, &
818 & im)*rr(k, is)*rsa(k-1, ih, im)+tta(k-1, ih, im)-tda(k-1, &
819 & ih, im))*denm
820  CALL pushreal8(rsa(k, ih, im))
821  rsa(k, ih, im) = rs(k, is) + ts(k, is)*rsa(k-1, ih, im)*denm
822  END DO
823 ! k loop
824 !-----add one layer at a time, going up.
825  DO k=ict-1,0,-1
826  CALL pushreal8(denm)
827  denm = ts(k, ih)/(1.-rs(k, ih)*rxa(k+1, im, is))
828  CALL pushreal8(rra(k, im, is))
829  rra(k, im, is) = rr(k, ih) + (td(k, ih)*rra(k+1, im, is)+(tt&
830 & (k, ih)-td(k, ih))*rxa(k+1, im, is))*denm
831  CALL pushreal8(rxa(k, im, is))
832  rxa(k, im, is) = rs(k, ih) + ts(k, ih)*rxa(k+1, im, is)*denm
833  END DO
834 ! k loop
835 !-----compute fluxes following Eq. (6.15) for fupdif and
836 ! Eq. (6.16) for (fdndir+fdndif)
837 ! fdndir is the direct downward flux
838 ! fdndif is the diffuse downward flux
839 ! fupdif is the diffuse upward flux
840  DO k=1,np+1
841  CALL pushreal8(denm)
842  denm = 1./(1.-rsa(k-1, ih, im)*rxa(k, im, is))
843  fdndir = tda(k-1, ih, im)
844  xx4 = tda(k-1, ih, im)*rra(k, im, is)
845  yy = tta(k-1, ih, im) - tda(k-1, ih, im)
846  fdndif = (xx4*rsa(k-1, ih, im)+yy)*denm
847  fupdif = (xx4+yy*rxa(k, im, is))*denm
848  flxdn = fdndir + fdndif - fupdif
849 !-----summation of fluxes over all sky situations;
850 ! the term in the brackets of Eq. (7.11)
851  fall(k) = fall(k) + flxdn*ct
852  END DO
853  END DO
854  END DO
855  END DO
856 ! is loop
857 ! im loop
858 ! ih loop
859 !-----End CLDFLX inline
860 !endif !overcast
861 !-----flux integration, Eq. (6.1)
862  DO k=1,np+1
863  flx_dev(i, k) = flx_dev(i, k) + fall(k)*hk_uv(ib)
864  END DO
865  END DO
866 !-----Inline SOLIR
867 !-----compute and update solar ir fluxes
868  CALL pushreal8(rr(np+1, 1))
869  rr(np+1, 1) = rsirbm_dev
870  CALL pushreal8(rr(np+1, 2))
871  rr(np+1, 2) = rsirbm_dev
872  CALL pushreal8(rs(np+1, 1))
873  rs(np+1, 1) = rsirdf_dev
874  CALL pushreal8(rs(np+1, 2))
875  rs(np+1, 2) = rsirdf_dev
876  CALL pushreal8(td(np+1, 1))
877  td(np+1, 1) = 0.0
878  CALL pushreal8(td(np+1, 2))
879  td(np+1, 2) = 0.0
880  CALL pushreal8(tt(np+1, 1))
881  tt(np+1, 1) = 0.0
882  CALL pushreal8(tt(np+1, 2))
883  tt(np+1, 2) = 0.0
884  CALL pushreal8(ts(np+1, 1))
885  ts(np+1, 1) = 0.0
886  CALL pushreal8(ts(np+1, 2))
887  ts(np+1, 2) = 0.0
888  CALL pushreal8(rr(0, 1))
889  rr(0, 1) = 0.0
890  CALL pushreal8(rr(0, 2))
891  rr(0, 2) = 0.0
892  CALL pushreal8(rs(0, 1))
893  rs(0, 1) = 0.0
894  CALL pushreal8(rs(0, 2))
895  rs(0, 2) = 0.0
896 ! td(0,1)=1.0
897 ! td(0,2)=1.0
898  CALL pushreal8(tt(0, 1))
899  tt(0, 1) = 1.0
900  CALL pushreal8(tt(0, 2))
901  tt(0, 2) = 1.0
902  CALL pushreal8(ts(0, 1))
903  ts(0, 1) = 1.0
904  CALL pushreal8(ts(0, 2))
905  ts(0, 2) = 1.0
906  cc1 = 0.0
907  CALL pushreal8(cc2)
908  cc2 = 0.0
909  CALL pushreal8(cc3)
910  cc3 = 0.0
911 !-----integration over spectral bands
912 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10)
913 ! The indices 1, 2, 3 are for ice, water, rain particles,
914 ! respectively.
915  DO ib=1,nband_ir
916  CALL pushinteger4(iv)
917  iv = ib + 5
918 !-----options for scaling cloud optical thickness
919 !if ( OVERCAST == .true. ) then
920 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10) and
921 ! compute cloud single scattering albedo and asymmetry factor
922 ! for a mixture of ice and liquid particles.
923 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
924 !
925 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
926 ! respectively.
927 ! call getnirtau1(ib,np,cosz_dev(i),dp_pa,fcld_col,reff_col,cwc_col,0,0,&
928 ! taubeam,taudiff,asycl,ssacl, &
929 ! aig_uv, awg_uv, arg_uv, &
930 ! aib_uv, awb_uv, arb_uv, &
931 ! aib_nir, awb_nir, arb_nir, &
932 ! aia_nir, awa_nir, ara_nir, &
933 ! aig_nir, awg_nir, arg_nir, &
934 ! caib, caif, &
935 ! CONS_GRAV )
936 !else
937 !-----Compute scaled cloud optical thickness. Eqs. (4.6) and (4.10) and
938 ! compute cloud single scattering albedo and asymmetry factor
939 ! for a mixture of ice and liquid particles.
940 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
941 !
942 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
943 ! respectively.
944  CALL pushreal8array(ssacl, np)
945  CALL pushreal8array(asycl, np)
946  CALL getnirtau1(ib, np, cosz_dev(i), dp_pa, fcld_col, reff_col, &
947 & cwc_col, ict, icb, taubeam, taudiff, asycl, ssacl, aig_uv&
948 & , awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir&
949 & , arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, &
950 & arg_nir, caib, caif, cons_grav)
951 !-----clouds within each of the high, middle, and low clouds are
952 ! assumed to be maximally overlapped, and the cloud cover (cc)
953 ! for a group (high, middle, or low) is the maximum cloud cover
954 ! of all the layers within a group
955 !MAT--DO NOT FUSE THIS LOOP
956 !MAT Loop must run to completion so that cc[1,2,3] are correct.
957  DO k=1,np
958  IF (k .LT. ict) THEN
959  IF (cc1 .LT. fcld_dev(i, k)) THEN
960  cc1 = fcld_dev(i, k)
961  CALL pushcontrol3b(4)
962  ELSE
963  CALL pushcontrol3b(5)
964  cc1 = cc1
965  END IF
966  ELSE IF (k .LT. icb) THEN
967  IF (cc2 .LT. fcld_dev(i, k)) THEN
968  CALL pushreal8(cc2)
969  cc2 = fcld_dev(i, k)
970  CALL pushcontrol3b(2)
971  ELSE
972  CALL pushreal8(cc2)
973  cc2 = cc2
974  CALL pushcontrol3b(3)
975  END IF
976  ELSE IF (cc3 .LT. fcld_dev(i, k)) THEN
977  CALL pushreal8(cc3)
978  cc3 = fcld_dev(i, k)
979  CALL pushcontrol3b(0)
980  ELSE
981  CALL pushreal8(cc3)
982  cc3 = cc3
983  CALL pushcontrol3b(1)
984  END IF
985  END DO
986 !MAT--DO NOT FUSE THIS LOOP
987 !endif !overcast
988  DO k=1,np
989  CALL pushreal8(tauclb(k))
990  tauclb(k) = taubeam(k, 1) + taubeam(k, 2) + taubeam(k, 3) + &
991 & taubeam(k, 4)
992  CALL pushreal8(tauclf(k))
993  tauclf(k) = taudiff(k, 1) + taudiff(k, 2) + taudiff(k, 3) + &
994 & taudiff(k, 4)
995  END DO
996 !-----integration over the k-distribution function
997  DO ik=1,nk_ir
998 !-----compute direct beam transmittances of the layer above pl(1)
999  CALL pushreal8(td(0, 1))
1000  td(0, 1) = exp(-(wvtoa*xk_ir(ik)/cosz_dev(i)))
1001  CALL pushreal8(td(0, 2))
1002  td(0, 2) = td(0, 1)
1003  DO k=1,np
1004  taurs = ry_ir(ib)*dp(k)
1005  tauwv = xk_ir(ik)*wh(k)
1006 !-----compute clear-sky optical thickness, single scattering albedo,
1007 ! and asymmetry factor. Eqs.(6.2)-(6.4)
1008  tausto = taurs + tauwv + taua_dev(i, k, iv) + 1.0e-7
1009  ssatau = ssaa_dev(i, k, iv) + taurs + 1.0e-8
1010  asysto = asya_dev(i, k, iv)
1011  CALL pushreal8(tautob)
1012  tautob = tausto
1013  asytob = asysto/ssatau
1014  CALL pushreal8(ssatob)
1015  ssatob = ssatau/tautob + 1.0e-8
1016  IF (ssatob .GT. 0.999999) THEN
1017  ssatob = 0.999999
1018  CALL pushcontrol1b(0)
1019  ELSE
1020  CALL pushcontrol1b(1)
1021  ssatob = ssatob
1022  END IF
1023 !-----Compute reflectance and transmittance of the clear portion
1024 ! of a layer
1025 !-----for direct incident radiation
1026  CALL deledd(tautob, ssatob, asytob, cosz_dev(i), rrt, ttt, tdt)
1027 !-----diffuse incident radiation is approximated by beam radiation with
1028 ! an incident angle of 53 degrees, Eqs. (6.5) and (6.6)
1029  CALL deledd(tautob, ssatob, asytob, dsm, rst, tst, dum)
1030  CALL pushreal8(rr(k, 1))
1031  rr(k, 1) = rrt
1032  CALL pushreal8(tt(k, 1))
1033  tt(k, 1) = ttt
1034  CALL pushreal8(td(k, 1))
1035  td(k, 1) = tdt
1036  CALL pushreal8(rs(k, 1))
1037  rs(k, 1) = rst
1038  CALL pushreal8(ts(k, 1))
1039  ts(k, 1) = tst
1040 !-----compute reflectance and transmittance of the cloudy portion
1041 ! of a layer
1042 !-----for direct incident radiation. Eqs.(6.2)-(6.4)
1043  tautob = tausto + tauclb(k)
1044  CALL pushreal8(ssatob)
1045  ssatob = (ssatau+ssacl(k)*tauclb(k))/tautob + 1.0e-8
1046  IF (ssatob .GT. 0.999999) THEN
1047  ssatob = 0.999999
1048  CALL pushcontrol1b(0)
1049  ELSE
1050  CALL pushcontrol1b(1)
1051  ssatob = ssatob
1052  END IF
1053  asytob = (asysto+asycl(k)*ssacl(k)*tauclb(k))/(ssatob*tautob)
1054 !-----for diffuse incident radiation
1055  CALL pushreal8(tautof)
1056  tautof = tausto + tauclf(k)
1057  CALL pushreal8(ssatof)
1058  ssatof = (ssatau+ssacl(k)*tauclf(k))/tautof + 1.0e-8
1059  IF (ssatof .GT. 0.999999) THEN
1060  ssatof = 0.999999
1061  CALL pushcontrol1b(0)
1062  ELSE
1063  ssatof = ssatof
1064  CALL pushcontrol1b(1)
1065  END IF
1066  asytof = (asysto+asycl(k)*ssacl(k)*tauclf(k))/(ssatof*tautof)
1067 !-----for direct incident radiation
1068  CALL deledd(tautob, ssatob, asytob, cosz_dev(i), rrt, ttt, tdt)
1069 !-----diffuse incident radiation is approximated by beam radiation with
1070 ! an incident angle of 53 degrees, Eqs.(6.5) and (6.6)
1071  CALL deledd(tautof, ssatof, asytof, dsm, rst, tst, dum)
1072  CALL pushreal8(rr(k, 2))
1073  rr(k, 2) = rrt
1074  CALL pushreal8(tt(k, 2))
1075  tt(k, 2) = ttt
1076  CALL pushreal8(td(k, 2))
1077  td(k, 2) = tdt
1078  CALL pushreal8(rs(k, 2))
1079  rs(k, 2) = rst
1080  CALL pushreal8(ts(k, 2))
1081  ts(k, 2) = tst
1082  END DO
1083 !-----FLUX CALCULATIONS
1084 ! initialize clear-sky flux (fclr), all-sky flux (fall),
1085 ! and surface downward fluxes (fsdir and fsdif)
1086  DO k=1,np+1
1087  fall(k) = 0.0
1088  END DO
1089 !-----for clear- and all-sky flux calculations when fractional
1090 ! cloud cover is either 0 or 1.
1091 !if ( OVERCAST == .true. ) then
1092 !-----Inline CLDFLXY
1093 !-----layers are added one at a time, going up from the surface.
1094 ! rra is the composite reflectance illuminated by beam radiation
1095 ! rxa is the composite reflectance illuminated from above
1096 ! by diffuse radiation
1097 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
1098 !-----compute transmittances and reflectances for a composite of
1099 ! layers. layers are added one at a time, going down from the top.
1100 ! tda is the composite direct transmittance illuminated by
1101 ! beam radiation
1102 ! tta is the composite total transmittance illuminated by
1103 ! beam radiation
1104 ! rsa is the composite reflectance illuminated from below
1105 ! by diffuse radiation
1106 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
1107 !-----compute fluxes following Eq. (6.15) for fupdif and
1108 ! Eq. (6.16) for (fdndir+fdndif)
1109 ! fdndir is the direct downward flux
1110 ! fdndif is the diffuse downward flux
1111 ! fupdif is the diffuse upward flux
1112 !-----ih=1 for clear sky; ih=2 for cloudy sky.
1113 !-----First set is ih = 1
1114 ! rra(np+1)=rr(np+1,1)
1115 ! rxa(np+1)=rs(np+1,1)
1116 !
1117 ! do k=np,0,-1
1118 ! denm=ts(k,1)/(1.-rs(k,1)*rxa(k+1))
1119 ! rra(k)=rr(k,1)+(td(k,1)*rra(k+1)+(tt(k,1)-td(k,1))*rxa(k+1))*denm
1120 ! rxa(k)=rs(k,1)+ts(k,1)*rxa(k+1)*denm
1121 ! end do
1122 !
1123 ! do k=1,np+1
1124 ! if (k <= np) then
1125 ! if (k == 1) then
1126 ! tdaold = td(0,1)
1127 ! ttaold = tt(0,1)
1128 ! rsaold = rs(0,1)
1129 !
1130 ! tdanew = 0.0
1131 ! ttanew = 0.0
1132 ! rsanew = 0.0
1133 ! end if
1134 ! denm=ts(k,1)/(1.-rsaold*rs(k,1))
1135 ! tdanew=tdaold*td(k,1)
1136 ! ttanew=tdaold*tt(k,1)+(tdaold*rsaold*rr(k,1)+ttaold-tdaold)*denm
1137 ! rsanew=rs(k,1)+ts(k,1)*rsaold*denm
1138 ! end if
1139 !
1140 ! denm=1./(1.-rsaold*rxa(k))
1141 ! fdndir=tdaold
1142 ! xx4=tdaold*rra(k)
1143 ! yy=ttaold-tdaold
1144 ! fdndif=(xx4*rsaold+yy)*denm
1145 ! fupdif=(xx4+yy*rxa(k))*denm
1146 ! flxdn=fdndir+fdndif-fupdif
1147 !
1148 ! fupc(k)=fupdif
1149 ! fclr(k)=flxdn
1150 !
1151 ! tdaold = tdanew
1152 ! ttaold = ttanew
1153 ! rsaold = rsanew
1154 !
1155 ! tdanew = 0.0
1156 ! ttanew = 0.0
1157 ! rsanew = 0.0
1158 ! end do
1159 !
1160 !!-----Second set is ih = 2
1161 !
1162 ! rra(np+1)=rr(np+1,2)
1163 ! rxa(np+1)=rs(np+1,2)
1164 !
1165 ! do k=np,0,-1
1166 ! denm=ts(k,2)/(1.-rs(k,2)*rxa(k+1))
1167 ! rra(k)=rr(k,2)+(td(k,2)*rra(k+1)+(tt(k,2)-td(k,2))*rxa(k+1))*denm
1168 ! rxa(k)=rs(k,2)+ts(k,2)*rxa(k+1)*denm
1169 ! end do
1170 !
1171 ! do k=1,np+1
1172 ! if (k <= np) then
1173 ! if (k == 1) then
1174 ! tdaold = td(0,2)
1175 ! ttaold = tt(0,2)
1176 ! rsaold = rs(0,2)
1177 !
1178 ! tdanew = 0.0
1179 ! ttanew = 0.0
1180 ! rsanew = 0.0
1181 ! end if
1182 ! denm=ts(k,2)/(1.-rsaold*rs(k,2))
1183 ! tdanew=tdaold*td(k,2)
1184 ! ttanew=tdaold*tt(k,2)+(tdaold*rsaold*rr(k,2)+ttaold-tdaold)*denm
1185 ! rsanew=rs(k,2)+ts(k,2)*rsaold*denm
1186 ! end if
1187 !
1188 ! denm=1./(1.-rsaold*rxa(k))
1189 ! fdndir=tdaold
1190 ! xx4=tdaold*rra(k)
1191 ! yy=ttaold-tdaold
1192 ! fdndif=(xx4*rsaold+yy)*denm
1193 ! fupdif=(xx4+yy*rxa(k))*denm
1194 ! flxdn=fdndir+fdndif-fupdif
1195 !
1196 ! fupa(k)=fupdif
1197 ! fall(k)=flxdn
1198 !
1199 ! tdaold = tdanew
1200 ! ttaold = ttanew
1201 ! rsaold = rsanew
1202 !
1203 ! tdanew = 0.0
1204 ! ttanew = 0.0
1205 ! rsanew = 0.0
1206 ! end do
1207 !
1208 ! fsdir=fdndir
1209 ! fsdif=fdndif
1210 !
1211 !!-----End CLDFLXY inline
1212 !
1213 !else
1214 !-----for clear- and all-sky flux calculations when fractional
1215 ! cloud cover is allowed to be between 0 and 1.
1216 ! the all-sky flux, fall is the summation inside the brackets
1217 ! of Eq. (7.11)
1218 !-----Inline CLDFLX
1219 !-----compute transmittances and reflectances for a composite of
1220 ! layers. layers are added one at a time, going down from the top.
1221 ! tda is the composite direct transmittance illuminated
1222 ! by beam radiation
1223 ! tta is the composite total transmittance illuminated by
1224 ! beam radiation
1225 ! rsa is the composite reflectance illuminated from below
1226 ! by diffuse radiation
1227 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
1228 !-----To save memory space, tda, tta, and rsa are pre-computed
1229 ! for k<icb. The dimension of these parameters is (m,np,2,2).
1230 ! It would have been (m,np,2,2,2) if these parameters were
1231 ! computed for all k's.
1232 !-----for high clouds
1233 ! ih=1 for clear-sky condition, ih=2 for cloudy-sky condition
1234  DO ih=1,2
1235  CALL pushreal8(tda(0, ih, 1))
1236  tda(0, ih, 1) = td(0, ih)
1237  CALL pushreal8(tta(0, ih, 1))
1238  tta(0, ih, 1) = tt(0, ih)
1239  CALL pushreal8(rsa(0, ih, 1))
1240  rsa(0, ih, 1) = rs(0, ih)
1241  CALL pushreal8(tda(0, ih, 2))
1242  tda(0, ih, 2) = td(0, ih)
1243  CALL pushreal8(tta(0, ih, 2))
1244  tta(0, ih, 2) = tt(0, ih)
1245  CALL pushreal8(rsa(0, ih, 2))
1246  rsa(0, ih, 2) = rs(0, ih)
1247  DO k=1,ict-1
1248  CALL pushreal8(denm)
1249  denm = ts(k, ih)/(1.-rsa(k-1, ih, 1)*rs(k, ih))
1250  CALL pushreal8(tda(k, ih, 1))
1251  tda(k, ih, 1) = tda(k-1, ih, 1)*td(k, ih)
1252  CALL pushreal8(tta(k, ih, 1))
1253  tta(k, ih, 1) = tda(k-1, ih, 1)*tt(k, ih) + (tda(k-1, ih, 1)*&
1254 & rsa(k-1, ih, 1)*rr(k, ih)+tta(k-1, ih, 1)-tda(k-1, ih, 1))*&
1255 & denm
1256  CALL pushreal8(rsa(k, ih, 1))
1257  rsa(k, ih, 1) = rs(k, ih) + ts(k, ih)*rsa(k-1, ih, 1)*denm
1258  CALL pushreal8(tda(k, ih, 2))
1259  tda(k, ih, 2) = tda(k, ih, 1)
1260  CALL pushreal8(tta(k, ih, 2))
1261  tta(k, ih, 2) = tta(k, ih, 1)
1262  CALL pushreal8(rsa(k, ih, 2))
1263  rsa(k, ih, 2) = rsa(k, ih, 1)
1264  END DO
1265 ! k loop
1266 !-----for middle clouds
1267 ! im=1 for clear-sky condition, im=2 for cloudy-sky condition
1268  DO k=ict,icb-1
1269  DO im=1,2
1270  CALL pushreal8(denm)
1271  denm = ts(k, im)/(1.-rsa(k-1, ih, im)*rs(k, im))
1272  CALL pushreal8(tda(k, ih, im))
1273  tda(k, ih, im) = tda(k-1, ih, im)*td(k, im)
1274  CALL pushreal8(tta(k, ih, im))
1275  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, im) + (tda(k-1, ih, &
1276 & im)*rsa(k-1, ih, im)*rr(k, im)+tta(k-1, ih, im)-tda(k-1, &
1277 & ih, im))*denm
1278  CALL pushreal8(rsa(k, ih, im))
1279  rsa(k, ih, im) = rs(k, im) + ts(k, im)*rsa(k-1, ih, im)*denm
1280  END DO
1281  END DO
1282  END DO
1283 ! im loop
1284 ! k loop
1285 ! ih loop
1286 !-----layers are added one at a time, going up from the surface.
1287 ! rra is the composite reflectance illuminated by beam radiation
1288 ! rxa is the composite reflectance illuminated from above
1289 ! by diffuse radiation
1290 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
1291 !-----To save memory space, rra and rxa are pre-computed for k>=icb.
1292 ! the dimension of these parameters is (m,np,2,2). It would have
1293 ! been (m,np,2,2,2) if these parameters were computed for all k's.
1294 !-----for the low clouds
1295 ! is=1 for clear-sky condition, is=2 for cloudy-sky condition
1296  DO is=1,2
1297  CALL pushreal8(rra(np+1, 1, is))
1298  rra(np+1, 1, is) = rr(np+1, is)
1299  CALL pushreal8(rxa(np+1, 1, is))
1300  rxa(np+1, 1, is) = rs(np+1, is)
1301  CALL pushreal8(rra(np+1, 2, is))
1302  rra(np+1, 2, is) = rr(np+1, is)
1303  CALL pushreal8(rxa(np+1, 2, is))
1304  rxa(np+1, 2, is) = rs(np+1, is)
1305  DO k=np,icb,-1
1306  CALL pushreal8(denm)
1307  denm = ts(k, is)/(1.-rs(k, is)*rxa(k+1, 1, is))
1308  CALL pushreal8(rra(k, 1, is))
1309  rra(k, 1, is) = rr(k, is) + (td(k, is)*rra(k+1, 1, is)+(tt(k, &
1310 & is)-td(k, is))*rxa(k+1, 1, is))*denm
1311  CALL pushreal8(rxa(k, 1, is))
1312  rxa(k, 1, is) = rs(k, is) + ts(k, is)*rxa(k+1, 1, is)*denm
1313  CALL pushreal8(rra(k, 2, is))
1314  rra(k, 2, is) = rra(k, 1, is)
1315  CALL pushreal8(rxa(k, 2, is))
1316  rxa(k, 2, is) = rxa(k, 1, is)
1317  END DO
1318 ! k loop
1319 !-----for middle clouds
1320  DO k=icb-1,ict,-1
1321  DO im=1,2
1322  CALL pushreal8(denm)
1323  denm = ts(k, im)/(1.-rs(k, im)*rxa(k+1, im, is))
1324  CALL pushreal8(rra(k, im, is))
1325  rra(k, im, is) = rr(k, im) + (td(k, im)*rra(k+1, im, is)+(tt&
1326 & (k, im)-td(k, im))*rxa(k+1, im, is))*denm
1327  CALL pushreal8(rxa(k, im, is))
1328  rxa(k, im, is) = rs(k, im) + ts(k, im)*rxa(k+1, im, is)*denm
1329  END DO
1330  END DO
1331  END DO
1332 ! im loop
1333 ! k loop
1334 ! is loop
1335 !-----integration over eight sky situations.
1336 ! ih, im, is denote high, middle and low cloud groups.
1337  DO ih=1,2
1338 !-----clear portion
1339  IF (ih .EQ. 1) THEN
1340  CALL pushreal8(ch)
1341  ch = 1.0 - cc1
1342 !-----cloudy portion
1343  CALL pushcontrol1b(1)
1344  ELSE
1345  CALL pushreal8(ch)
1346  ch = cc1
1347  CALL pushcontrol1b(0)
1348  END IF
1349  DO im=1,2
1350 !-----clear portion
1351  IF (im .EQ. 1) THEN
1352  CALL pushreal8(cm)
1353  cm = ch*(1.0-cc2)
1354 !-----cloudy portion
1355  CALL pushcontrol1b(1)
1356  ELSE
1357  CALL pushreal8(cm)
1358  cm = ch*cc2
1359  CALL pushcontrol1b(0)
1360  END IF
1361  DO is=1,2
1362 !-----clear portion
1363  IF (is .EQ. 1) THEN
1364  CALL pushreal8(ct)
1365  ct = cm*(1.0-cc3)
1366 !-----cloudy portion
1367  CALL pushcontrol1b(1)
1368  ELSE
1369  CALL pushreal8(ct)
1370  ct = cm*cc3
1371  CALL pushcontrol1b(0)
1372  END IF
1373 !-----add one layer at a time, going down.
1374  DO k=icb,np
1375  CALL pushreal8(denm)
1376  denm = ts(k, is)/(1.-rsa(k-1, ih, im)*rs(k, is))
1377  CALL pushreal8(tda(k, ih, im))
1378  tda(k, ih, im) = tda(k-1, ih, im)*td(k, is)
1379  CALL pushreal8(tta(k, ih, im))
1380  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, is) + (tda(k-1, ih&
1381 & , im)*rr(k, is)*rsa(k-1, ih, im)+tta(k-1, ih, im)-tda(k-&
1382 & 1, ih, im))*denm
1383  CALL pushreal8(rsa(k, ih, im))
1384  rsa(k, ih, im) = rs(k, is) + ts(k, is)*rsa(k-1, ih, im)*&
1385 & denm
1386  END DO
1387 ! k loop
1388 !-----add one layer at a time, going up.
1389  DO k=ict-1,0,-1
1390  CALL pushreal8(denm)
1391  denm = ts(k, ih)/(1.-rs(k, ih)*rxa(k+1, im, is))
1392  CALL pushreal8(rra(k, im, is))
1393  rra(k, im, is) = rr(k, ih) + (td(k, ih)*rra(k+1, im, is)+(&
1394 & tt(k, ih)-td(k, ih))*rxa(k+1, im, is))*denm
1395  CALL pushreal8(rxa(k, im, is))
1396  rxa(k, im, is) = rs(k, ih) + ts(k, ih)*rxa(k+1, im, is)*&
1397 & denm
1398  END DO
1399 ! k loop
1400 !-----compute fluxes following Eq. (6.15) for fupdif and
1401 ! Eq. (6.16) for (fdndir+fdndif)
1402 ! fdndir is the direct downward flux
1403 ! fdndif is the diffuse downward flux
1404 ! fupdif is the diffuse upward flux
1405  DO k=1,np+1
1406  CALL pushreal8(denm)
1407  denm = 1./(1.-rsa(k-1, ih, im)*rxa(k, im, is))
1408  fdndir = tda(k-1, ih, im)
1409  xx4 = tda(k-1, ih, im)*rra(k, im, is)
1410  yy = tta(k-1, ih, im) - tda(k-1, ih, im)
1411  fdndif = (xx4*rsa(k-1, ih, im)+yy)*denm
1412  fupdif = (xx4+yy*rxa(k, im, is))*denm
1413  flxdn = fdndir + fdndif - fupdif
1414 !-----summation of fluxes over all sky situations;
1415 ! the term in the brackets of Eq. (7.11)
1416  fall(k) = fall(k) + flxdn*ct
1417  END DO
1418  END DO
1419  END DO
1420  END DO
1421 ! is loop
1422 ! im loop
1423 ! ih loop
1424 !-----End CLDFLX inline
1425 !endif !overcast
1426 !-----flux integration following Eq. (6.1)
1427  DO k=1,np+1
1428  flx_dev(i, k) = flx_dev(i, k) + fall(k)*hk_ir(ib, ik)
1429  END DO
1430  END DO
1431  END DO
1432 ! ik loop
1433 !-----compute pressure-scaled o2 amount following Eq. (3.5) with f=1.
1434 ! unit is (cm-atm)stp. 165.22 = (1000/980)*23.14%*(22400/32)
1435 ! compute flux reduction due to oxygen following Eq. (3.18). 0.0633 is the
1436 ! fraction of insolation contained in the oxygen bands
1437  df(0) = 0.0
1438  cnt = 165.22*snt
1439  so2(1) = scal0*cnt
1440 ! LLT increased parameter 145 to 155 to enhance effect
1441  df(1) = 0.0633*(1.-exp(-(0.000155*sqrt(so2(1)))))
1442  DO k=1,np
1443  so2(k+1) = so2(k) + scal(k)*cnt
1444 ! LLT increased parameter 145 to 155 to enhance effect
1445  df(k+1) = 0.0633*(1.0-exp(-(0.000155*sqrt(so2(k+1)))))
1446  END DO
1447 !-----for solar heating due to co2 scaling follows Eq(3.5) with f=1.
1448 ! unit is (cm-atm)stp. 789 = (1000/980)*(44/28.97)*(22400/44)
1449  so2(1) = 789.*co2*scal0
1450  DO k=1,np
1451  so2(k+1) = so2(k) + 789.*co2*scal(k)
1452  END DO
1453 !-----The updated flux reduction for co2 absorption in Band 7 where absorption due to
1454 ! water vapor and co2 are both moderate. df is given by the second term on the
1455 ! right-hand-side of Eq. (3.24) divided by So. so2 and swh are the co2 and
1456 ! water vapor amounts integrated from the top of the atmosphere
1457  u1 = -3.0
1458  du = 0.15
1459  w1 = -4.0
1460  dw = 0.15
1461 !-----Inline RFLX
1462  x0 = u1 + REAL(nu)*du
1463  y0 = w1 + REAL(nw)*dw
1464  x1 = u1 - 0.5*du
1465  y1 = w1 - 0.5*dw
1466  DO k=1,np+1
1467  x3 = log10(so2(k)*snt)
1468  IF (x3 .GT. x0) THEN
1469  ulog = x0
1470  ELSE
1471  ulog = x3
1472  END IF
1473  x4 = log10(swh(k)*snt)
1474  IF (x4 .GT. y0) THEN
1475  wlog = y0
1476  CALL pushcontrol1b(0)
1477  ELSE
1478  wlog = x4
1479  CALL pushcontrol1b(1)
1480  END IF
1481  CALL pushinteger4(ic)
1482  ic = int((ulog-x1)/du + 1.)
1483  CALL pushinteger4(iw)
1484  iw = int((wlog-y1)/dw + 1.)
1485  IF (ic .LT. 2) ic = 2
1486  IF (iw .LT. 2) iw = 2
1487  IF (ic .GT. nu) ic = nu
1488  IF (iw .GT. nw) iw = nw
1489  dc = ulog - REAL(ic-2)*du - u1
1490  dd = wlog - REAL(iw-2)*dw - w1
1491  x2 = cah(ic-1, iw-1) + (cah(ic-1, iw)-cah(ic-1, iw-1))/dw*dd
1492  y2 = x2 + (cah(ic, iw-1)-cah(ic-1, iw-1))/du*dc
1493  IF (y2 .LT. 0.0) THEN
1494  y2 = 0.0
1495  CALL pushcontrol1b(0)
1496  ELSE
1497  CALL pushcontrol1b(1)
1498  y2 = y2
1499  END IF
1500 ! LLT increase CO2 effect to help reduce cold tropopause bias
1501  df(k) = df(k) + 1.5*y2
1502  END DO
1503 !-----df is the updated flux reduction for co2 absorption
1504 ! in Band 8 where the co2 absorption has a large impact
1505 ! on the heating of middle atmosphere. From the table
1506 ! given by Eq. (3.19)
1507  u1 = 0.000250
1508  du = 0.000050
1509  w1 = -2.0
1510  dw = 0.05
1511 !-----Inline RFLX
1512  x0 = u1 + REAL(nx)*du
1513  y0 = w1 + REAL(ny)*dw
1514  x1 = u1 - 0.5*du
1515  y1 = w1 - 0.5*dw
1516  DO k=1,np+1
1517  IF (co2*snt .GT. x0) THEN
1518  ulog = x0
1519  ELSE
1520  ulog = co2*snt
1521  END IF
1522  x5 = log10(pl_dev(i, k))
1523  IF (x5 .GT. y0) THEN
1524  wlog = y0
1525  ELSE
1526  wlog = x5
1527  END IF
1528  CALL pushinteger4(ic)
1529  ic = int((ulog-x1)/du + 1.)
1530  CALL pushinteger4(iw)
1531  iw = int((wlog-y1)/dw + 1.)
1532  IF (ic .LT. 2) ic = 2
1533  IF (iw .LT. 2) iw = 2
1534  IF (ic .GT. nx) ic = nx
1535  IF (iw .GT. ny) iw = ny
1536  dc = ulog - REAL(ic-2)*du - u1
1537  dd = wlog - REAL(iw-2)*dw - w1
1538  x2 = coa(ic-1, iw-1) + (coa(ic-1, iw)-coa(ic-1, iw-1))/dw*dd
1539  y2 = x2 + (coa(ic, iw-1)-coa(ic-1, iw-1))/du*dc
1540  IF (y2 .LT. 0.0) THEN
1541  y2 = 0.0
1542  ELSE
1543  y2 = y2
1544  END IF
1545 ! LLT increase CO2 effect to help reduce cold tropopause bias
1546  df(k) = df(k) + 1.5*y2
1547  END DO
1548 !-----adjust the o2-co2 reduction below cloud top following Eq. (6.18)
1549  foundtop = 0
1550  DO k=1,np
1551  IF (fcld_dev(i, k) .GT. 0.02 .AND. foundtop .EQ. 0) THEN
1552  foundtop = 1
1553  ntop = k
1554  END IF
1555  END DO
1556  IF (foundtop .EQ. 0) ntop = np + 1
1557  dftop = df(ntop)
1558  DO k=1,np+1
1559  IF (k .GT. ntop) THEN
1560  xx4 = flx_dev(i, k)/flx_dev(i, ntop)
1561  CALL pushreal8(df(k))
1562  df(k) = dftop + xx4*(df(k)-dftop)
1563  CALL pushcontrol1b(1)
1564  ELSE
1565  CALL pushcontrol1b(0)
1566  END IF
1567  END DO
1568 !-----update the net fluxes
1569  DO k=1,np+1
1570  IF (df(k) .GT. flx_dev(i, k) - 1.0e-8) THEN
1571  CALL pushcontrol1b(0)
1572  ELSE
1573  CALL pushcontrol1b(1)
1574  END IF
1575  END DO
1576  dfb = 0.0_8
1577  DO k=np+1,1,-1
1578  dfb(k) = dfb(k) - flx_devb(i, k)
1579  CALL popcontrol1b(branch)
1580  IF (branch .EQ. 0) THEN
1581  flx_devb(i, k) = flx_devb(i, k) + dfb(k)
1582  dfb(k) = 0.0_8
1583  END IF
1584  END DO
1585  dftopb = 0.0_8
1586  DO k=np+1,1,-1
1587  CALL popcontrol1b(branch)
1588  IF (branch .NE. 0) THEN
1589  xx4 = flx_dev(i, k)/flx_dev(i, ntop)
1590  CALL popreal8(df(k))
1591  dftopb = dftopb + (1.0_8-xx4)*dfb(k)
1592  xx4b = (df(k)-dftop)*dfb(k)
1593  dfb(k) = xx4*dfb(k)
1594  tempb56 = xx4b/flx_dev(i, ntop)
1595  flx_devb(i, k) = flx_devb(i, k) + tempb56
1596  flx_devb(i, ntop) = flx_devb(i, ntop) - flx_dev(i, k)*tempb56/&
1597 & flx_dev(i, ntop)
1598  END IF
1599  END DO
1600  dfb(ntop) = dfb(ntop) + dftopb
1601  DO k=np+1,1,-1
1602  CALL popinteger4(iw)
1603  CALL popinteger4(ic)
1604  END DO
1605  dw = 0.15
1606  swhb = 0.0_8
1607  DO k=np+1,1,-1
1608  y2b = 1.5*dfb(k)
1609  CALL popcontrol1b(branch)
1610  IF (branch .EQ. 0) y2b = 0.0_8
1611  x2b = y2b
1612  ddb = (cah(ic-1, iw)-cah(ic-1, iw-1))*x2b/dw
1613  wlogb = ddb
1614  CALL popinteger4(iw)
1615  CALL popinteger4(ic)
1616  CALL popcontrol1b(branch)
1617  IF (branch .EQ. 0) THEN
1618  x4b = 0.0_8
1619  ELSE
1620  x4b = wlogb
1621  END IF
1622  swhb(k) = swhb(k) + x4b/(swh(k)*log(10.0))
1623  END DO
1624  tsb = 0.0_8
1625  ttb = 0.0_8
1626  rsab = 0.0_8
1627  reff_colb = 0.0_8
1628  rrb = 0.0_8
1629  rsb = 0.0_8
1630  wvtoab = 0.0_8
1631  fallb = 0.0_8
1632  cc1b = 0.0_8
1633  cc2b = 0.0_8
1634  cc3b = 0.0_8
1635  asyclb = 0.0_8
1636  ttab = 0.0_8
1637  rxab = 0.0_8
1638  tdab = 0.0_8
1639  ssaclb = 0.0_8
1640  fcld_colb = 0.0_8
1641  cwc_colb = 0.0_8
1642  rrab = 0.0_8
1643  tauclbb = 0.0_8
1644  tauclfb = 0.0_8
1645  whb = 0.0_8
1646  tdb = 0.0_8
1647  DO ib=nband_ir,1,-1
1648  DO ik=nk_ir,1,-1
1649  DO k=np+1,1,-1
1650  fallb(k) = fallb(k) + hk_ir(ib, ik)*flx_devb(i, k)
1651  END DO
1652  DO ih=2,1,-1
1653  chb = 0.0_8
1654  DO im=2,1,-1
1655  cmb = 0.0_8
1656  DO is=2,1,-1
1657  ctb = 0.0_8
1658  DO k=np+1,1,-1
1659  yy = tta(k-1, ih, im) - tda(k-1, ih, im)
1660  denm = 1./(1.-rsa(k-1, ih, im)*rxa(k, im, is))
1661  xx4 = tda(k-1, ih, im)*rra(k, im, is)
1662  fupdif = (xx4+yy*rxa(k, im, is))*denm
1663  fdndif = (xx4*rsa(k-1, ih, im)+yy)*denm
1664  fdndir = tda(k-1, ih, im)
1665  flxdn = fdndir + fdndif - fupdif
1666  flxdnb = ct*fallb(k)
1667  ctb = ctb + flxdn*fallb(k)
1668  fdndirb = flxdnb
1669  fdndifb = flxdnb
1670  fupdifb = -flxdnb
1671  tempb53 = denm*fupdifb
1672  denmb = (xx4*rsa(k-1, ih, im)+yy)*fdndifb + (xx4+yy*rxa(k&
1673 & , im, is))*fupdifb
1674  tempb54 = denm*fdndifb
1675  xx4b = rsa(k-1, ih, im)*tempb54 + tempb53
1676  yyb = tempb54 + rxa(k, im, is)*tempb53
1677  ttab(k-1, ih, im) = ttab(k-1, ih, im) + yyb
1678  tdab(k-1, ih, im) = tdab(k-1, ih, im) + rra(k, im, is)*&
1679 & xx4b + fdndirb - yyb
1680  rrab(k, im, is) = rrab(k, im, is) + tda(k-1, ih, im)*xx4b
1681  CALL popreal8(denm)
1682  temp38 = -(rsa(k-1, ih, im)*rxa(k, im, is)) + 1.
1683  tempb55 = -(denmb/temp38**2)
1684  rxab(k, im, is) = rxab(k, im, is) + yy*tempb53 - rsa(k-1, &
1685 & ih, im)*tempb55
1686  rsab(k-1, ih, im) = rsab(k-1, ih, im) + xx4*tempb54 - rxa(&
1687 & k, im, is)*tempb55
1688  END DO
1689  DO k=0,ict-1,1
1690  CALL popreal8(rxa(k, im, is))
1691  tempb50 = rxa(k+1, im, is)*rxab(k, im, is)
1692  CALL popreal8(rra(k, im, is))
1693  tempb52 = denm*rrab(k, im, is)
1694  temp37 = rxa(k+1, im, is)
1695  temp36 = tt(k, ih) - td(k, ih)
1696  rrb(k, ih) = rrb(k, ih) + rrab(k, im, is)
1697  tdb(k, ih) = tdb(k, ih) + (rra(k+1, im, is)-temp37)*&
1698 & tempb52
1699  rrab(k+1, im, is) = rrab(k+1, im, is) + td(k, ih)*tempb52
1700  denmb = (td(k, ih)*rra(k+1, im, is)+temp36*temp37)*rrab(k&
1701 & , im, is) + ts(k, ih)*tempb50
1702  ttb(k, ih) = ttb(k, ih) + temp37*tempb52
1703  rrab(k, im, is) = 0.0_8
1704  temp35 = -(rs(k, ih)*rxa(k+1, im, is)) + 1.
1705  tsb(k, ih) = tsb(k, ih) + denmb/temp35 + denm*tempb50
1706  tempb51 = -(ts(k, ih)*denmb/temp35**2)
1707  rsb(k, ih) = rsb(k, ih) + rxab(k, im, is) - rxa(k+1, im, &
1708 & is)*tempb51
1709  rxab(k+1, im, is) = rxab(k+1, im, is) + ts(k, ih)*denm*&
1710 & rxab(k, im, is)
1711  rxab(k, im, is) = 0.0_8
1712  rxab(k+1, im, is) = rxab(k+1, im, is) + temp36*tempb52 - &
1713 & rs(k, ih)*tempb51
1714  CALL popreal8(denm)
1715  END DO
1716  DO k=np,icb,-1
1717  CALL popreal8(rsa(k, ih, im))
1718  tempb47 = rsa(k-1, ih, im)*rsab(k, ih, im)
1719  CALL popreal8(tta(k, ih, im))
1720  tempb49 = denm*ttab(k, ih, im)
1721  temp34 = rsa(k-1, ih, im)
1722  temp33 = tda(k-1, ih, im)*rr(k, is)
1723  tdab(k-1, ih, im) = tdab(k-1, ih, im) + td(k, is)*tdab(k, &
1724 & ih, im) + (temp34*rr(k, is)-1.0)*tempb49 + tt(k, is)*&
1725 & ttab(k, ih, im)
1726  ttb(k, is) = ttb(k, is) + tda(k-1, ih, im)*ttab(k, ih, im)
1727  rrb(k, is) = rrb(k, is) + temp34*tda(k-1, ih, im)*tempb49
1728  ttab(k-1, ih, im) = ttab(k-1, ih, im) + tempb49
1729  denmb = (temp33*temp34+tta(k-1, ih, im)-tda(k-1, ih, im))*&
1730 & ttab(k, ih, im) + ts(k, is)*tempb47
1731  ttab(k, ih, im) = 0.0_8
1732  CALL popreal8(tda(k, ih, im))
1733  tdb(k, is) = tdb(k, is) + tda(k-1, ih, im)*tdab(k, ih, im)
1734  tdab(k, ih, im) = 0.0_8
1735  temp32 = -(rsa(k-1, ih, im)*rs(k, is)) + 1.
1736  tsb(k, is) = tsb(k, is) + denmb/temp32 + denm*tempb47
1737  tempb48 = -(ts(k, is)*denmb/temp32**2)
1738  rsb(k, is) = rsb(k, is) + rsab(k, ih, im) - rsa(k-1, ih, &
1739 & im)*tempb48
1740  rsab(k-1, ih, im) = rsab(k-1, ih, im) + ts(k, is)*denm*&
1741 & rsab(k, ih, im)
1742  rsab(k, ih, im) = 0.0_8
1743  rsab(k-1, ih, im) = rsab(k-1, ih, im) + temp33*tempb49 - &
1744 & rs(k, is)*tempb48
1745  CALL popreal8(denm)
1746  END DO
1747  CALL popcontrol1b(branch)
1748  IF (branch .EQ. 0) THEN
1749  CALL popreal8(ct)
1750  cmb = cmb + cc3*ctb
1751  cc3b = cc3b + cm*ctb
1752  ELSE
1753  CALL popreal8(ct)
1754  cmb = cmb + (1.0-cc3)*ctb
1755  cc3b = cc3b - cm*ctb
1756  END IF
1757  END DO
1758  CALL popcontrol1b(branch)
1759  IF (branch .EQ. 0) THEN
1760  CALL popreal8(cm)
1761  chb = chb + cc2*cmb
1762  cc2b = cc2b + ch*cmb
1763  ELSE
1764  CALL popreal8(cm)
1765  chb = chb + (1.0-cc2)*cmb
1766  cc2b = cc2b - ch*cmb
1767  END IF
1768  END DO
1769  CALL popcontrol1b(branch)
1770  IF (branch .EQ. 0) THEN
1771  CALL popreal8(ch)
1772  cc1b = cc1b + chb
1773  ELSE
1774  CALL popreal8(ch)
1775  cc1b = cc1b - chb
1776  END IF
1777  END DO
1778  DO is=2,1,-1
1779  DO k=ict,icb-1,1
1780  DO im=2,1,-1
1781  CALL popreal8(rxa(k, im, is))
1782  tempb44 = rxa(k+1, im, is)*rxab(k, im, is)
1783  CALL popreal8(rra(k, im, is))
1784  tempb46 = denm*rrab(k, im, is)
1785  temp31 = rxa(k+1, im, is)
1786  temp30 = tt(k, im) - td(k, im)
1787  rrb(k, im) = rrb(k, im) + rrab(k, im, is)
1788  tdb(k, im) = tdb(k, im) + (rra(k+1, im, is)-temp31)*tempb46
1789  rrab(k+1, im, is) = rrab(k+1, im, is) + td(k, im)*tempb46
1790  denmb = (td(k, im)*rra(k+1, im, is)+temp30*temp31)*rrab(k, &
1791 & im, is) + ts(k, im)*tempb44
1792  ttb(k, im) = ttb(k, im) + temp31*tempb46
1793  rrab(k, im, is) = 0.0_8
1794  temp29 = -(rs(k, im)*rxa(k+1, im, is)) + 1.
1795  tsb(k, im) = tsb(k, im) + denmb/temp29 + denm*tempb44
1796  tempb45 = -(ts(k, im)*denmb/temp29**2)
1797  rsb(k, im) = rsb(k, im) + rxab(k, im, is) - rxa(k+1, im, is)&
1798 & *tempb45
1799  rxab(k+1, im, is) = rxab(k+1, im, is) + ts(k, im)*denm*rxab(&
1800 & k, im, is)
1801  rxab(k, im, is) = 0.0_8
1802  rxab(k+1, im, is) = rxab(k+1, im, is) + temp30*tempb46 - rs(&
1803 & k, im)*tempb45
1804  CALL popreal8(denm)
1805  END DO
1806  END DO
1807  DO k=icb,np,1
1808  CALL popreal8(rxa(k, 2, is))
1809  rxab(k, 1, is) = rxab(k, 1, is) + rxab(k, 2, is)
1810  rxab(k, 2, is) = 0.0_8
1811  CALL popreal8(rra(k, 2, is))
1812  rrab(k, 1, is) = rrab(k, 1, is) + rrab(k, 2, is)
1813  rrab(k, 2, is) = 0.0_8
1814  CALL popreal8(rxa(k, 1, is))
1815  tempb41 = rxa(k+1, 1, is)*rxab(k, 1, is)
1816  CALL popreal8(rra(k, 1, is))
1817  tempb43 = denm*rrab(k, 1, is)
1818  temp28 = rxa(k+1, 1, is)
1819  temp27 = tt(k, is) - td(k, is)
1820  rrb(k, is) = rrb(k, is) + rrab(k, 1, is)
1821  tdb(k, is) = tdb(k, is) + (rra(k+1, 1, is)-temp28)*tempb43
1822  rrab(k+1, 1, is) = rrab(k+1, 1, is) + td(k, is)*tempb43
1823  denmb = (td(k, is)*rra(k+1, 1, is)+temp27*temp28)*rrab(k, 1, &
1824 & is) + ts(k, is)*tempb41
1825  ttb(k, is) = ttb(k, is) + temp28*tempb43
1826  rrab(k, 1, is) = 0.0_8
1827  temp26 = -(rs(k, is)*rxa(k+1, 1, is)) + 1.
1828  tsb(k, is) = tsb(k, is) + denmb/temp26 + denm*tempb41
1829  tempb42 = -(ts(k, is)*denmb/temp26**2)
1830  rsb(k, is) = rsb(k, is) + rxab(k, 1, is) - rxa(k+1, 1, is)*&
1831 & tempb42
1832  rxab(k+1, 1, is) = rxab(k+1, 1, is) + ts(k, is)*denm*rxab(k, 1&
1833 & , is)
1834  rxab(k, 1, is) = 0.0_8
1835  rxab(k+1, 1, is) = rxab(k+1, 1, is) + temp27*tempb43 - rs(k, &
1836 & is)*tempb42
1837  CALL popreal8(denm)
1838  END DO
1839  CALL popreal8(rxa(np+1, 2, is))
1840  rsb(np+1, is) = rsb(np+1, is) + rxab(np+1, 2, is)
1841  rxab(np+1, 2, is) = 0.0_8
1842  CALL popreal8(rra(np+1, 2, is))
1843  rrb(np+1, is) = rrb(np+1, is) + rrab(np+1, 2, is)
1844  rrab(np+1, 2, is) = 0.0_8
1845  CALL popreal8(rxa(np+1, 1, is))
1846  rsb(np+1, is) = rsb(np+1, is) + rxab(np+1, 1, is)
1847  rxab(np+1, 1, is) = 0.0_8
1848  CALL popreal8(rra(np+1, 1, is))
1849  rrb(np+1, is) = rrb(np+1, is) + rrab(np+1, 1, is)
1850  rrab(np+1, 1, is) = 0.0_8
1851  END DO
1852  DO ih=2,1,-1
1853  DO k=icb-1,ict,-1
1854  DO im=2,1,-1
1855  CALL popreal8(rsa(k, ih, im))
1856  tempb38 = rsa(k-1, ih, im)*rsab(k, ih, im)
1857  CALL popreal8(tta(k, ih, im))
1858  tempb40 = denm*ttab(k, ih, im)
1859  temp25 = rsa(k-1, ih, im)
1860  temp24 = tda(k-1, ih, im)*rr(k, im)
1861  tdab(k-1, ih, im) = tdab(k-1, ih, im) + td(k, im)*tdab(k, ih&
1862 & , im) + (temp25*rr(k, im)-1.0)*tempb40 + tt(k, im)*ttab(k&
1863 & , ih, im)
1864  ttb(k, im) = ttb(k, im) + tda(k-1, ih, im)*ttab(k, ih, im)
1865  rrb(k, im) = rrb(k, im) + temp25*tda(k-1, ih, im)*tempb40
1866  ttab(k-1, ih, im) = ttab(k-1, ih, im) + tempb40
1867  denmb = (temp24*temp25+tta(k-1, ih, im)-tda(k-1, ih, im))*&
1868 & ttab(k, ih, im) + ts(k, im)*tempb38
1869  ttab(k, ih, im) = 0.0_8
1870  CALL popreal8(tda(k, ih, im))
1871  tdb(k, im) = tdb(k, im) + tda(k-1, ih, im)*tdab(k, ih, im)
1872  tdab(k, ih, im) = 0.0_8
1873  temp23 = -(rsa(k-1, ih, im)*rs(k, im)) + 1.
1874  tsb(k, im) = tsb(k, im) + denmb/temp23 + denm*tempb38
1875  tempb39 = -(ts(k, im)*denmb/temp23**2)
1876  rsb(k, im) = rsb(k, im) + rsab(k, ih, im) - rsa(k-1, ih, im)&
1877 & *tempb39
1878  rsab(k-1, ih, im) = rsab(k-1, ih, im) + ts(k, im)*denm*rsab(&
1879 & k, ih, im)
1880  rsab(k, ih, im) = 0.0_8
1881  rsab(k-1, ih, im) = rsab(k-1, ih, im) + temp24*tempb40 - rs(&
1882 & k, im)*tempb39
1883  CALL popreal8(denm)
1884  END DO
1885  END DO
1886  DO k=ict-1,1,-1
1887  CALL popreal8(rsa(k, ih, 2))
1888  rsab(k, ih, 1) = rsab(k, ih, 1) + rsab(k, ih, 2)
1889  rsab(k, ih, 2) = 0.0_8
1890  CALL popreal8(tta(k, ih, 2))
1891  ttab(k, ih, 1) = ttab(k, ih, 1) + ttab(k, ih, 2)
1892  ttab(k, ih, 2) = 0.0_8
1893  CALL popreal8(tda(k, ih, 2))
1894  tdab(k, ih, 1) = tdab(k, ih, 1) + tdab(k, ih, 2)
1895  tdab(k, ih, 2) = 0.0_8
1896  CALL popreal8(rsa(k, ih, 1))
1897  tempb35 = rsa(k-1, ih, 1)*rsab(k, ih, 1)
1898  CALL popreal8(tta(k, ih, 1))
1899  tempb37 = denm*ttab(k, ih, 1)
1900  temp22 = rsa(k-1, ih, 1)
1901  temp21 = tda(k-1, ih, 1)*rr(k, ih)
1902  tdab(k-1, ih, 1) = tdab(k-1, ih, 1) + td(k, ih)*tdab(k, ih, 1)&
1903 & + (temp22*rr(k, ih)-1.0)*tempb37 + tt(k, ih)*ttab(k, ih, 1)
1904  ttb(k, ih) = ttb(k, ih) + tda(k-1, ih, 1)*ttab(k, ih, 1)
1905  rrb(k, ih) = rrb(k, ih) + temp22*tda(k-1, ih, 1)*tempb37
1906  ttab(k-1, ih, 1) = ttab(k-1, ih, 1) + tempb37
1907  denmb = (temp21*temp22+tta(k-1, ih, 1)-tda(k-1, ih, 1))*ttab(k&
1908 & , ih, 1) + ts(k, ih)*tempb35
1909  ttab(k, ih, 1) = 0.0_8
1910  CALL popreal8(tda(k, ih, 1))
1911  tdb(k, ih) = tdb(k, ih) + tda(k-1, ih, 1)*tdab(k, ih, 1)
1912  tdab(k, ih, 1) = 0.0_8
1913  temp20 = -(rsa(k-1, ih, 1)*rs(k, ih)) + 1.
1914  tsb(k, ih) = tsb(k, ih) + denmb/temp20 + denm*tempb35
1915  tempb36 = -(ts(k, ih)*denmb/temp20**2)
1916  rsb(k, ih) = rsb(k, ih) + rsab(k, ih, 1) - rsa(k-1, ih, 1)*&
1917 & tempb36
1918  rsab(k-1, ih, 1) = rsab(k-1, ih, 1) + ts(k, ih)*denm*rsab(k, &
1919 & ih, 1)
1920  rsab(k, ih, 1) = 0.0_8
1921  rsab(k-1, ih, 1) = rsab(k-1, ih, 1) + temp21*tempb37 - rs(k, &
1922 & ih)*tempb36
1923  CALL popreal8(denm)
1924  END DO
1925  CALL popreal8(rsa(0, ih, 2))
1926  rsb(0, ih) = rsb(0, ih) + rsab(0, ih, 2)
1927  rsab(0, ih, 2) = 0.0_8
1928  CALL popreal8(tta(0, ih, 2))
1929  ttb(0, ih) = ttb(0, ih) + ttab(0, ih, 2)
1930  ttab(0, ih, 2) = 0.0_8
1931  CALL popreal8(tda(0, ih, 2))
1932  tdb(0, ih) = tdb(0, ih) + tdab(0, ih, 2)
1933  tdab(0, ih, 2) = 0.0_8
1934  CALL popreal8(rsa(0, ih, 1))
1935  rsb(0, ih) = rsb(0, ih) + rsab(0, ih, 1)
1936  rsab(0, ih, 1) = 0.0_8
1937  CALL popreal8(tta(0, ih, 1))
1938  ttb(0, ih) = ttb(0, ih) + ttab(0, ih, 1)
1939  ttab(0, ih, 1) = 0.0_8
1940  CALL popreal8(tda(0, ih, 1))
1941  tdb(0, ih) = tdb(0, ih) + tdab(0, ih, 1)
1942  tdab(0, ih, 1) = 0.0_8
1943  END DO
1944  DO k=np+1,1,-1
1945  fallb(k) = 0.0_8
1946  END DO
1947  DO k=np,1,-1
1948  CALL popreal8(ts(k, 2))
1949  tstb = tsb(k, 2)
1950  tsb(k, 2) = 0.0_8
1951  CALL popreal8(rs(k, 2))
1952  rstb = rsb(k, 2)
1953  rsb(k, 2) = 0.0_8
1954  CALL popreal8(td(k, 2))
1955  tdtb = tdb(k, 2)
1956  tdb(k, 2) = 0.0_8
1957  CALL popreal8(tt(k, 2))
1958  tttb = ttb(k, 2)
1959  ttb(k, 2) = 0.0_8
1960  CALL popreal8(rr(k, 2))
1961  rrtb = rrb(k, 2)
1962  rrb(k, 2) = 0.0_8
1963  tauwv = xk_ir(ik)*wh(k)
1964  taurs = ry_ir(ib)*dp(k)
1965  tausto = taurs + tauwv + taua_dev(i, k, iv) + 1.0e-7
1966  tautof = tausto + tauclf(k)
1967  asysto = asya_dev(i, k, iv)
1968  asytof = (asysto+asycl(k)*ssacl(k)*tauclf(k))/(ssatof*tautof)
1969  tautofb = 0.0_8
1970  ssatofb = 0.0_8
1971  asytofb = 0.0_8
1972  dumb = 0.0_8
1973  CALL deledd_b(tautof, tautofb, ssatof, ssatofb, asytof, asytofb&
1974 & , dsm, rst, rstb, tst, tstb, dum, dumb)
1975  tautob = tausto + tauclb(k)
1976  asytob = (asysto+asycl(k)*ssacl(k)*tauclb(k))/(ssatob*tautob)
1977  tautobb = 0.0_8
1978  ssatobb = 0.0_8
1979  asytobb = 0.0_8
1980  CALL deledd_b(tautob, tautobb, ssatob, ssatobb, asytob, asytobb&
1981 & , cosz_dev(i), rrt, rrtb, ttt, tttb, tdt, tdtb)
1982  tempb33 = asytofb/(ssatof*tautof)
1983  temp19 = asycl(k)*ssacl(k)
1984  tempb34 = -((asysto+temp19*tauclf(k))*tempb33/(ssatof*tautof))
1985  asystob = tempb33
1986  asyclb(k) = asyclb(k) + tauclf(k)*ssacl(k)*tempb33
1987  ssaclb(k) = ssaclb(k) + tauclf(k)*asycl(k)*tempb33
1988  tauclfb(k) = tauclfb(k) + temp19*tempb33
1989  ssatofb = ssatofb + tautof*tempb34
1990  tautofb = tautofb + ssatof*tempb34
1991  CALL popcontrol1b(branch)
1992  IF (branch .EQ. 0) THEN
1993  ssatau = ssaa_dev(i, k, iv) + taurs + 1.0e-8
1994  ssatofb = 0.0_8
1995  ELSE
1996  ssatau = ssaa_dev(i, k, iv) + taurs + 1.0e-8
1997  END IF
1998  tempb31 = asytobb/(ssatob*tautob)
1999  CALL popreal8(ssatof)
2000  tempb30 = ssatofb/tautof
2001  ssataub = tempb30
2002  ssaclb(k) = ssaclb(k) + tauclb(k)*asycl(k)*tempb31 + tauclf(k)*&
2003 & tempb30
2004  tautofb = tautofb - (ssatau+ssacl(k)*tauclf(k))*tempb30/tautof
2005  tauclfb(k) = tauclfb(k) + tautofb + ssacl(k)*tempb30
2006  CALL popreal8(tautof)
2007  taustob = tautofb
2008  temp18 = asycl(k)*ssacl(k)
2009  tempb32 = -((asysto+temp18*tauclb(k))*tempb31/(ssatob*tautob))
2010  asystob = asystob + tempb31
2011  asyclb(k) = asyclb(k) + tauclb(k)*ssacl(k)*tempb31
2012  tauclbb(k) = tauclbb(k) + temp18*tempb31
2013  ssatobb = ssatobb + tautob*tempb32
2014  tautobb = tautobb + ssatob*tempb32
2015  CALL popcontrol1b(branch)
2016  IF (branch .EQ. 0) ssatobb = 0.0_8
2017  CALL popreal8(ssatob)
2018  tempb29 = ssatobb/tautob
2019  ssataub = ssataub + tempb29
2020  ssaclb(k) = ssaclb(k) + tauclb(k)*tempb29
2021  tautobb = tautobb - (ssatau+ssacl(k)*tauclb(k))*tempb29/tautob
2022  tauclbb(k) = tauclbb(k) + tautobb + ssacl(k)*tempb29
2023  taustob = taustob + tautobb
2024  CALL popreal8(ts(k, 1))
2025  tstb = tsb(k, 1)
2026  tsb(k, 1) = 0.0_8
2027  CALL popreal8(rs(k, 1))
2028  rstb = rsb(k, 1)
2029  rsb(k, 1) = 0.0_8
2030  CALL popreal8(td(k, 1))
2031  tdtb = tdb(k, 1)
2032  tdb(k, 1) = 0.0_8
2033  CALL popreal8(tt(k, 1))
2034  tttb = ttb(k, 1)
2035  ttb(k, 1) = 0.0_8
2036  CALL popreal8(rr(k, 1))
2037  rrtb = rrb(k, 1)
2038  rrb(k, 1) = 0.0_8
2039  tautob = tausto
2040  asytob = asysto/ssatau
2041  tautobb = 0.0_8
2042  ssatobb = 0.0_8
2043  asytobb = 0.0_8
2044  dumb = 0.0_8
2045  CALL deledd_b(tautob, tautobb, ssatob, ssatobb, asytob, asytobb&
2046 & , dsm, rst, rstb, tst, tstb, dum, dumb)
2047  CALL deledd_b(tautob, tautobb, ssatob, ssatobb, asytob, asytobb&
2048 & , cosz_dev(i), rrt, rrtb, ttt, tttb, tdt, tdtb)
2049  CALL popcontrol1b(branch)
2050  IF (branch .EQ. 0) ssatobb = 0.0_8
2051  CALL popreal8(ssatob)
2052  ssataub = ssataub + ssatobb/tautob - asysto*asytobb/ssatau**2
2053  tautobb = tautobb - ssatau*ssatobb/tautob**2
2054  asystob = asystob + asytobb/ssatau
2055  CALL popreal8(tautob)
2056  taustob = taustob + tautobb
2057  asya_devb(i, k, iv) = asya_devb(i, k, iv) + asystob
2058  ssaa_devb(i, k, iv) = ssaa_devb(i, k, iv) + ssataub
2059  tauwvb = taustob
2060  taua_devb(i, k, iv) = taua_devb(i, k, iv) + taustob
2061  whb(k) = whb(k) + xk_ir(ik)*tauwvb
2062  END DO
2063  CALL popreal8(td(0, 2))
2064  tdb(0, 1) = tdb(0, 1) + tdb(0, 2)
2065  tdb(0, 2) = 0.0_8
2066  CALL popreal8(td(0, 1))
2067  wvtoab = wvtoab - xk_ir(ik)*exp(-(xk_ir(ik)*(wvtoa/cosz_dev(i))))*&
2068 & tdb(0, 1)/cosz_dev(i)
2069  tdb(0, 1) = 0.0_8
2070  END DO
2071  taudiffb = 0.0_8
2072  taubeamb = 0.0_8
2073  DO k=np,1,-1
2074  CALL popreal8(tauclf(k))
2075  taudiffb(k, 1) = taudiffb(k, 1) + tauclfb(k)
2076  taudiffb(k, 2) = taudiffb(k, 2) + tauclfb(k)
2077  taudiffb(k, 3) = taudiffb(k, 3) + tauclfb(k)
2078  taudiffb(k, 4) = taudiffb(k, 4) + tauclfb(k)
2079  tauclfb(k) = 0.0_8
2080  CALL popreal8(tauclb(k))
2081  taubeamb(k, 1) = taubeamb(k, 1) + tauclbb(k)
2082  taubeamb(k, 2) = taubeamb(k, 2) + tauclbb(k)
2083  taubeamb(k, 3) = taubeamb(k, 3) + tauclbb(k)
2084  taubeamb(k, 4) = taubeamb(k, 4) + tauclbb(k)
2085  tauclbb(k) = 0.0_8
2086  END DO
2087  DO k=np,1,-1
2088  CALL popcontrol3b(branch)
2089  IF (branch .LT. 3) THEN
2090  IF (branch .EQ. 0) THEN
2091  CALL popreal8(cc3)
2092  fcld_devb(i, k) = fcld_devb(i, k) + cc3b
2093  cc3b = 0.0_8
2094  ELSE IF (branch .EQ. 1) THEN
2095  CALL popreal8(cc3)
2096  ELSE
2097  CALL popreal8(cc2)
2098  fcld_devb(i, k) = fcld_devb(i, k) + cc2b
2099  cc2b = 0.0_8
2100  END IF
2101  ELSE IF (branch .EQ. 3) THEN
2102  CALL popreal8(cc2)
2103  ELSE IF (branch .EQ. 4) THEN
2104  fcld_devb(i, k) = fcld_devb(i, k) + cc1b
2105  cc1b = 0.0_8
2106  END IF
2107  END DO
2108  CALL popreal8array(asycl, np)
2109  CALL popreal8array(ssacl, np)
2110  CALL getnirtau1_b(ib, np, cosz_dev(i), dp_pa, fcld_col, fcld_colb, &
2111 & reff_col, reff_colb, cwc_col, cwc_colb, ict, icb, &
2112 & taubeam, taubeamb, taudiff, taudiffb, asycl, asyclb, &
2113 & ssacl, ssaclb, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, &
2114 & arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, &
2115 & ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, &
2116 & cons_grav)
2117  CALL popinteger4(iv)
2118  END DO
2119  CALL popreal8(cc3)
2120  CALL popreal8(cc2)
2121  CALL popreal8(ts(0, 2))
2122  tsb(0, 2) = 0.0_8
2123  CALL popreal8(ts(0, 1))
2124  tsb(0, 1) = 0.0_8
2125  CALL popreal8(tt(0, 2))
2126  ttb(0, 2) = 0.0_8
2127  CALL popreal8(tt(0, 1))
2128  ttb(0, 1) = 0.0_8
2129  CALL popreal8(rs(0, 2))
2130  rsb(0, 2) = 0.0_8
2131  CALL popreal8(rs(0, 1))
2132  rsb(0, 1) = 0.0_8
2133  CALL popreal8(rr(0, 2))
2134  rrb(0, 2) = 0.0_8
2135  CALL popreal8(rr(0, 1))
2136  rrb(0, 1) = 0.0_8
2137  CALL popreal8(ts(np+1, 2))
2138  tsb(np+1, 2) = 0.0_8
2139  CALL popreal8(ts(np+1, 1))
2140  tsb(np+1, 1) = 0.0_8
2141  CALL popreal8(tt(np+1, 2))
2142  ttb(np+1, 2) = 0.0_8
2143  CALL popreal8(tt(np+1, 1))
2144  ttb(np+1, 1) = 0.0_8
2145  CALL popreal8(td(np+1, 2))
2146  tdb(np+1, 2) = 0.0_8
2147  CALL popreal8(td(np+1, 1))
2148  tdb(np+1, 1) = 0.0_8
2149  CALL popreal8(rs(np+1, 2))
2150  rsb(np+1, 2) = 0.0_8
2151  CALL popreal8(rs(np+1, 1))
2152  rsb(np+1, 1) = 0.0_8
2153  CALL popreal8(rr(np+1, 2))
2154  rrb(np+1, 2) = 0.0_8
2155  CALL popreal8(rr(np+1, 1))
2156  rrb(np+1, 1) = 0.0_8
2157  ohb = 0.0_8
2158  cc1b = 0.0_8
2159  cc2b = 0.0_8
2160  cc3b = 0.0_8
2161  o3toab = 0.0_8
2162  DO ib=nband_uv,1,-1
2163  DO k=np+1,1,-1
2164  fallb(k) = fallb(k) + hk_uv(ib)*flx_devb(i, k)
2165  END DO
2166  DO ih=2,1,-1
2167  chb = 0.0_8
2168  DO im=2,1,-1
2169  cmb = 0.0_8
2170  DO is=2,1,-1
2171  ctb = 0.0_8
2172  DO k=np+1,1,-1
2173  yy = tta(k-1, ih, im) - tda(k-1, ih, im)
2174  denm = 1./(1.-rsa(k-1, ih, im)*rxa(k, im, is))
2175  xx4 = tda(k-1, ih, im)*rra(k, im, is)
2176  fupdif = (xx4+yy*rxa(k, im, is))*denm
2177  fdndif = (xx4*rsa(k-1, ih, im)+yy)*denm
2178  fdndir = tda(k-1, ih, im)
2179  flxdn = fdndir + fdndif - fupdif
2180  flxdnb = ct*fallb(k)
2181  ctb = ctb + flxdn*fallb(k)
2182  fdndirb = flxdnb
2183  fdndifb = flxdnb
2184  fupdifb = -flxdnb
2185  tempb26 = denm*fupdifb
2186  denmb = (xx4*rsa(k-1, ih, im)+yy)*fdndifb + (xx4+yy*rxa(k, &
2187 & im, is))*fupdifb
2188  tempb27 = denm*fdndifb
2189  xx4b = rsa(k-1, ih, im)*tempb27 + tempb26
2190  yyb = tempb27 + rxa(k, im, is)*tempb26
2191  ttab(k-1, ih, im) = ttab(k-1, ih, im) + yyb
2192  tdab(k-1, ih, im) = tdab(k-1, ih, im) + rra(k, im, is)*xx4b &
2193 & + fdndirb - yyb
2194  rrab(k, im, is) = rrab(k, im, is) + tda(k-1, ih, im)*xx4b
2195  CALL popreal8(denm)
2196  temp17 = -(rsa(k-1, ih, im)*rxa(k, im, is)) + 1.
2197  tempb28 = -(denmb/temp17**2)
2198  rxab(k, im, is) = rxab(k, im, is) + yy*tempb26 - rsa(k-1, ih&
2199 & , im)*tempb28
2200  rsab(k-1, ih, im) = rsab(k-1, ih, im) + xx4*tempb27 - rxa(k&
2201 & , im, is)*tempb28
2202  END DO
2203  DO k=0,ict-1,1
2204  CALL popreal8(rxa(k, im, is))
2205  tempb23 = rxa(k+1, im, is)*rxab(k, im, is)
2206  CALL popreal8(rra(k, im, is))
2207  tempb25 = denm*rrab(k, im, is)
2208  temp16 = rxa(k+1, im, is)
2209  temp15 = tt(k, ih) - td(k, ih)
2210  rrb(k, ih) = rrb(k, ih) + rrab(k, im, is)
2211  tdb(k, ih) = tdb(k, ih) + (rra(k+1, im, is)-temp16)*tempb25
2212  rrab(k+1, im, is) = rrab(k+1, im, is) + td(k, ih)*tempb25
2213  denmb = (td(k, ih)*rra(k+1, im, is)+temp15*temp16)*rrab(k, &
2214 & im, is) + ts(k, ih)*tempb23
2215  ttb(k, ih) = ttb(k, ih) + temp16*tempb25
2216  rrab(k, im, is) = 0.0_8
2217  temp14 = -(rs(k, ih)*rxa(k+1, im, is)) + 1.
2218  tsb(k, ih) = tsb(k, ih) + denmb/temp14 + denm*tempb23
2219  tempb24 = -(ts(k, ih)*denmb/temp14**2)
2220  rsb(k, ih) = rsb(k, ih) + rxab(k, im, is) - rxa(k+1, im, is)&
2221 & *tempb24
2222  rxab(k+1, im, is) = rxab(k+1, im, is) + ts(k, ih)*denm*rxab(&
2223 & k, im, is)
2224  rxab(k, im, is) = 0.0_8
2225  rxab(k+1, im, is) = rxab(k+1, im, is) + temp15*tempb25 - rs(&
2226 & k, ih)*tempb24
2227  CALL popreal8(denm)
2228  END DO
2229  DO k=np,icb,-1
2230  CALL popreal8(rsa(k, ih, im))
2231  tempb20 = rsa(k-1, ih, im)*rsab(k, ih, im)
2232  CALL popreal8(tta(k, ih, im))
2233  tempb22 = denm*ttab(k, ih, im)
2234  temp13 = rsa(k-1, ih, im)
2235  temp12 = tda(k-1, ih, im)*rr(k, is)
2236  tdab(k-1, ih, im) = tdab(k-1, ih, im) + td(k, is)*tdab(k, ih&
2237 & , im) + (temp13*rr(k, is)-1.0)*tempb22 + tt(k, is)*ttab(k&
2238 & , ih, im)
2239  ttb(k, is) = ttb(k, is) + tda(k-1, ih, im)*ttab(k, ih, im)
2240  rrb(k, is) = rrb(k, is) + temp13*tda(k-1, ih, im)*tempb22
2241  ttab(k-1, ih, im) = ttab(k-1, ih, im) + tempb22
2242  denmb = (temp12*temp13+tta(k-1, ih, im)-tda(k-1, ih, im))*&
2243 & ttab(k, ih, im) + ts(k, is)*tempb20
2244  ttab(k, ih, im) = 0.0_8
2245  CALL popreal8(tda(k, ih, im))
2246  tdb(k, is) = tdb(k, is) + tda(k-1, ih, im)*tdab(k, ih, im)
2247  tdab(k, ih, im) = 0.0_8
2248  temp11 = -(rsa(k-1, ih, im)*rs(k, is)) + 1.
2249  tsb(k, is) = tsb(k, is) + denmb/temp11 + denm*tempb20
2250  tempb21 = -(ts(k, is)*denmb/temp11**2)
2251  rsb(k, is) = rsb(k, is) + rsab(k, ih, im) - rsa(k-1, ih, im)&
2252 & *tempb21
2253  rsab(k-1, ih, im) = rsab(k-1, ih, im) + ts(k, is)*denm*rsab(&
2254 & k, ih, im)
2255  rsab(k, ih, im) = 0.0_8
2256  rsab(k-1, ih, im) = rsab(k-1, ih, im) + temp12*tempb22 - rs(&
2257 & k, is)*tempb21
2258  CALL popreal8(denm)
2259  END DO
2260  CALL popcontrol1b(branch)
2261  IF (branch .EQ. 0) THEN
2262  CALL popreal8(ct)
2263  cmb = cmb + cc3*ctb
2264  cc3b = cc3b + cm*ctb
2265  ELSE
2266  CALL popreal8(ct)
2267  cmb = cmb + (1.0-cc3)*ctb
2268  cc3b = cc3b - cm*ctb
2269  END IF
2270  END DO
2271  CALL popcontrol1b(branch)
2272  IF (branch .EQ. 0) THEN
2273  CALL popreal8(cm)
2274  chb = chb + cc2*cmb
2275  cc2b = cc2b + ch*cmb
2276  ELSE
2277  CALL popreal8(cm)
2278  chb = chb + (1.0-cc2)*cmb
2279  cc2b = cc2b - ch*cmb
2280  END IF
2281  END DO
2282  CALL popcontrol1b(branch)
2283  IF (branch .EQ. 0) THEN
2284  CALL popreal8(ch)
2285  cc1b = cc1b + chb
2286  ELSE
2287  CALL popreal8(ch)
2288  cc1b = cc1b - chb
2289  END IF
2290  END DO
2291  DO is=2,1,-1
2292  DO k=ict,icb-1,1
2293  DO im=2,1,-1
2294  CALL popreal8(rxa(k, im, is))
2295  tempb17 = rxa(k+1, im, is)*rxab(k, im, is)
2296  CALL popreal8(rra(k, im, is))
2297  tempb19 = denm*rrab(k, im, is)
2298  temp10 = rxa(k+1, im, is)
2299  temp9 = tt(k, im) - td(k, im)
2300  rrb(k, im) = rrb(k, im) + rrab(k, im, is)
2301  tdb(k, im) = tdb(k, im) + (rra(k+1, im, is)-temp10)*tempb19
2302  rrab(k+1, im, is) = rrab(k+1, im, is) + td(k, im)*tempb19
2303  denmb = (td(k, im)*rra(k+1, im, is)+temp9*temp10)*rrab(k, im, &
2304 & is) + ts(k, im)*tempb17
2305  ttb(k, im) = ttb(k, im) + temp10*tempb19
2306  rrab(k, im, is) = 0.0_8
2307  temp8 = -(rs(k, im)*rxa(k+1, im, is)) + 1.
2308  tsb(k, im) = tsb(k, im) + denmb/temp8 + denm*tempb17
2309  tempb18 = -(ts(k, im)*denmb/temp8**2)
2310  rsb(k, im) = rsb(k, im) + rxab(k, im, is) - rxa(k+1, im, is)*&
2311 & tempb18
2312  rxab(k+1, im, is) = rxab(k+1, im, is) + ts(k, im)*denm*rxab(k&
2313 & , im, is)
2314  rxab(k, im, is) = 0.0_8
2315  rxab(k+1, im, is) = rxab(k+1, im, is) + temp9*tempb19 - rs(k, &
2316 & im)*tempb18
2317  CALL popreal8(denm)
2318  END DO
2319  END DO
2320  DO k=icb,np,1
2321  CALL popreal8(rxa(k, 2, is))
2322  rxab(k, 1, is) = rxab(k, 1, is) + rxab(k, 2, is)
2323  rxab(k, 2, is) = 0.0_8
2324  CALL popreal8(rra(k, 2, is))
2325  rrab(k, 1, is) = rrab(k, 1, is) + rrab(k, 2, is)
2326  rrab(k, 2, is) = 0.0_8
2327  CALL popreal8(rxa(k, 1, is))
2328  tempb14 = rxa(k+1, 1, is)*rxab(k, 1, is)
2329  CALL popreal8(rra(k, 1, is))
2330  tempb16 = denm*rrab(k, 1, is)
2331  temp7 = rxa(k+1, 1, is)
2332  temp6 = tt(k, is) - td(k, is)
2333  rrb(k, is) = rrb(k, is) + rrab(k, 1, is)
2334  tdb(k, is) = tdb(k, is) + (rra(k+1, 1, is)-temp7)*tempb16
2335  rrab(k+1, 1, is) = rrab(k+1, 1, is) + td(k, is)*tempb16
2336  denmb = (td(k, is)*rra(k+1, 1, is)+temp6*temp7)*rrab(k, 1, is) +&
2337 & ts(k, is)*tempb14
2338  ttb(k, is) = ttb(k, is) + temp7*tempb16
2339  rrab(k, 1, is) = 0.0_8
2340  temp5 = -(rs(k, is)*rxa(k+1, 1, is)) + 1.
2341  tsb(k, is) = tsb(k, is) + denmb/temp5 + denm*tempb14
2342  tempb15 = -(ts(k, is)*denmb/temp5**2)
2343  rsb(k, is) = rsb(k, is) + rxab(k, 1, is) - rxa(k+1, 1, is)*&
2344 & tempb15
2345  rxab(k+1, 1, is) = rxab(k+1, 1, is) + ts(k, is)*denm*rxab(k, 1, &
2346 & is)
2347  rxab(k, 1, is) = 0.0_8
2348  rxab(k+1, 1, is) = rxab(k+1, 1, is) + temp6*tempb16 - rs(k, is)*&
2349 & tempb15
2350  CALL popreal8(denm)
2351  END DO
2352  CALL popreal8(rxa(np+1, 2, is))
2353  rsb(np+1, is) = rsb(np+1, is) + rxab(np+1, 2, is)
2354  rxab(np+1, 2, is) = 0.0_8
2355  CALL popreal8(rra(np+1, 2, is))
2356  rrb(np+1, is) = rrb(np+1, is) + rrab(np+1, 2, is)
2357  rrab(np+1, 2, is) = 0.0_8
2358  CALL popreal8(rxa(np+1, 1, is))
2359  rsb(np+1, is) = rsb(np+1, is) + rxab(np+1, 1, is)
2360  rxab(np+1, 1, is) = 0.0_8
2361  CALL popreal8(rra(np+1, 1, is))
2362  rrb(np+1, is) = rrb(np+1, is) + rrab(np+1, 1, is)
2363  rrab(np+1, 1, is) = 0.0_8
2364  END DO
2365  DO ih=2,1,-1
2366  DO k=icb-1,ict,-1
2367  DO im=2,1,-1
2368  CALL popreal8(rsa(k, ih, im))
2369  tempb11 = rsa(k-1, ih, im)*rsab(k, ih, im)
2370  CALL popreal8(tta(k, ih, im))
2371  tempb13 = denm*ttab(k, ih, im)
2372  temp4 = rsa(k-1, ih, im)
2373  temp3 = tda(k-1, ih, im)*rr(k, im)
2374  tdab(k-1, ih, im) = tdab(k-1, ih, im) + td(k, im)*tdab(k, ih, &
2375 & im) + (temp4*rr(k, im)-1.0)*tempb13 + tt(k, im)*ttab(k, ih, &
2376 & im)
2377  ttb(k, im) = ttb(k, im) + tda(k-1, ih, im)*ttab(k, ih, im)
2378  rrb(k, im) = rrb(k, im) + temp4*tda(k-1, ih, im)*tempb13
2379  ttab(k-1, ih, im) = ttab(k-1, ih, im) + tempb13
2380  denmb = (temp3*temp4+tta(k-1, ih, im)-tda(k-1, ih, im))*ttab(k&
2381 & , ih, im) + ts(k, im)*tempb11
2382  ttab(k, ih, im) = 0.0_8
2383  CALL popreal8(tda(k, ih, im))
2384  tdb(k, im) = tdb(k, im) + tda(k-1, ih, im)*tdab(k, ih, im)
2385  tdab(k, ih, im) = 0.0_8
2386  temp2 = -(rsa(k-1, ih, im)*rs(k, im)) + 1.
2387  tsb(k, im) = tsb(k, im) + denmb/temp2 + denm*tempb11
2388  tempb12 = -(ts(k, im)*denmb/temp2**2)
2389  rsb(k, im) = rsb(k, im) + rsab(k, ih, im) - rsa(k-1, ih, im)*&
2390 & tempb12
2391  rsab(k-1, ih, im) = rsab(k-1, ih, im) + ts(k, im)*denm*rsab(k&
2392 & , ih, im)
2393  rsab(k, ih, im) = 0.0_8
2394  rsab(k-1, ih, im) = rsab(k-1, ih, im) + temp3*tempb13 - rs(k, &
2395 & im)*tempb12
2396  CALL popreal8(denm)
2397  END DO
2398  END DO
2399  DO k=ict-1,1,-1
2400  CALL popreal8(rsa(k, ih, 2))
2401  rsab(k, ih, 1) = rsab(k, ih, 1) + rsab(k, ih, 2)
2402  rsab(k, ih, 2) = 0.0_8
2403  CALL popreal8(tta(k, ih, 2))
2404  ttab(k, ih, 1) = ttab(k, ih, 1) + ttab(k, ih, 2)
2405  ttab(k, ih, 2) = 0.0_8
2406  CALL popreal8(tda(k, ih, 2))
2407  tdab(k, ih, 1) = tdab(k, ih, 1) + tdab(k, ih, 2)
2408  tdab(k, ih, 2) = 0.0_8
2409  CALL popreal8(rsa(k, ih, 1))
2410  tempb8 = rsa(k-1, ih, 1)*rsab(k, ih, 1)
2411  CALL popreal8(tta(k, ih, 1))
2412  tempb10 = denm*ttab(k, ih, 1)
2413  temp1 = rsa(k-1, ih, 1)
2414  temp0 = tda(k-1, ih, 1)*rr(k, ih)
2415  tdab(k-1, ih, 1) = tdab(k-1, ih, 1) + td(k, ih)*tdab(k, ih, 1) +&
2416 & (temp1*rr(k, ih)-1.0)*tempb10 + tt(k, ih)*ttab(k, ih, 1)
2417  ttb(k, ih) = ttb(k, ih) + tda(k-1, ih, 1)*ttab(k, ih, 1)
2418  rrb(k, ih) = rrb(k, ih) + temp1*tda(k-1, ih, 1)*tempb10
2419  ttab(k-1, ih, 1) = ttab(k-1, ih, 1) + tempb10
2420  denmb = (temp0*temp1+tta(k-1, ih, 1)-tda(k-1, ih, 1))*ttab(k, ih&
2421 & , 1) + ts(k, ih)*tempb8
2422  ttab(k, ih, 1) = 0.0_8
2423  CALL popreal8(tda(k, ih, 1))
2424  tdb(k, ih) = tdb(k, ih) + tda(k-1, ih, 1)*tdab(k, ih, 1)
2425  tdab(k, ih, 1) = 0.0_8
2426  temp = -(rsa(k-1, ih, 1)*rs(k, ih)) + 1.
2427  tsb(k, ih) = tsb(k, ih) + denmb/temp + denm*tempb8
2428  tempb9 = -(ts(k, ih)*denmb/temp**2)
2429  rsb(k, ih) = rsb(k, ih) + rsab(k, ih, 1) - rsa(k-1, ih, 1)*&
2430 & tempb9
2431  rsab(k-1, ih, 1) = rsab(k-1, ih, 1) + ts(k, ih)*denm*rsab(k, ih&
2432 & , 1)
2433  rsab(k, ih, 1) = 0.0_8
2434  rsab(k-1, ih, 1) = rsab(k-1, ih, 1) + temp0*tempb10 - rs(k, ih)*&
2435 & tempb9
2436  CALL popreal8(denm)
2437  END DO
2438  CALL popreal8(rsa(0, ih, 2))
2439  rsb(0, ih) = rsb(0, ih) + rsab(0, ih, 2)
2440  rsab(0, ih, 2) = 0.0_8
2441  CALL popreal8(tta(0, ih, 2))
2442  ttb(0, ih) = ttb(0, ih) + ttab(0, ih, 2)
2443  ttab(0, ih, 2) = 0.0_8
2444  CALL popreal8(tda(0, ih, 2))
2445  tdb(0, ih) = tdb(0, ih) + tdab(0, ih, 2)
2446  tdab(0, ih, 2) = 0.0_8
2447  CALL popreal8(rsa(0, ih, 1))
2448  rsb(0, ih) = rsb(0, ih) + rsab(0, ih, 1)
2449  rsab(0, ih, 1) = 0.0_8
2450  CALL popreal8(tta(0, ih, 1))
2451  ttb(0, ih) = ttb(0, ih) + ttab(0, ih, 1)
2452  ttab(0, ih, 1) = 0.0_8
2453  CALL popreal8(tda(0, ih, 1))
2454  tdb(0, ih) = tdb(0, ih) + tdab(0, ih, 1)
2455  tdab(0, ih, 1) = 0.0_8
2456  END DO
2457  DO k=np+1,1,-1
2458  fallb(k) = 0.0_8
2459  END DO
2460  DO k=np,1,-1
2461  CALL popreal8(ts(k, 2))
2462  tstb = tsb(k, 2)
2463  tsb(k, 2) = 0.0_8
2464  CALL popreal8(rs(k, 2))
2465  rstb = rsb(k, 2)
2466  rsb(k, 2) = 0.0_8
2467  CALL popreal8(td(k, 2))
2468  tdtb = tdb(k, 2)
2469  tdb(k, 2) = 0.0_8
2470  CALL popreal8(tt(k, 2))
2471  tttb = ttb(k, 2)
2472  ttb(k, 2) = 0.0_8
2473  CALL popreal8(rr(k, 2))
2474  rrtb = rrb(k, 2)
2475  rrb(k, 2) = 0.0_8
2476  tauwv = wk_uv(ib)*wh(k)
2477  taurs = ry_uv(ib)*dp(k)
2478  tauoz = zk_uv(ib)*oh(k)
2479  tausto = taurs + tauoz + tauwv + taua_dev(i, k, ib) + 1.0e-7
2480  tautof = tausto + tauclf(k)
2481  asysto = asya_dev(i, k, ib)
2482  asytof = (asysto+asycl(k)*tauclf(k))/(ssatof*tautof)
2483  tautofb = 0.0_8
2484  ssatofb = 0.0_8
2485  asytofb = 0.0_8
2486  dumb = 0.0_8
2487  CALL deledd_b(tautof, tautofb, ssatof, ssatofb, asytof, asytofb, &
2488 & dsm, rst, rstb, tst, tstb, dum, dumb)
2489  tautob = tausto + tauclb(k)
2490  asytob = (asysto+asycl(k)*tauclb(k))/(ssatob*tautob)
2491  tautobb = 0.0_8
2492  ssatobb = 0.0_8
2493  asytobb = 0.0_8
2494  CALL deledd_b(tautob, tautobb, ssatob, ssatobb, asytob, asytobb, &
2495 & cosz_dev(i), rrt, rrtb, ttt, tttb, tdt, tdtb)
2496  tempb6 = asytofb/(ssatof*tautof)
2497  tempb7 = -((asysto+asycl(k)*tauclf(k))*tempb6/(ssatof*tautof))
2498  asystob = tempb6
2499  asyclb(k) = asyclb(k) + tauclf(k)*tempb6
2500  tauclfb(k) = tauclfb(k) + asycl(k)*tempb6
2501  ssatofb = ssatofb + tautof*tempb7
2502  tautofb = tautofb + ssatof*tempb7
2503  CALL popcontrol1b(branch)
2504  IF (branch .EQ. 0) THEN
2505  ssatau = ssaa_dev(i, k, ib) + taurs
2506  ssatofb = 0.0_8
2507  ELSE
2508  ssatau = ssaa_dev(i, k, ib) + taurs
2509  END IF
2510  CALL popreal8(ssatof)
2511  tempb3 = ssatofb/tautof
2512  ssataub = tempb3
2513  tautofb = tautofb - (ssatau+tauclf(k))*tempb3/tautof
2514  tauclfb(k) = tauclfb(k) + tautofb + tempb3
2515  CALL popreal8(tautof)
2516  taustob = tautofb
2517  tempb4 = asytobb/(ssatob*tautob)
2518  tempb5 = -((asysto+asycl(k)*tauclb(k))*tempb4/(ssatob*tautob))
2519  asystob = asystob + tempb4
2520  asyclb(k) = asyclb(k) + tauclb(k)*tempb4
2521  tauclbb(k) = tauclbb(k) + asycl(k)*tempb4
2522  ssatobb = ssatobb + tautob*tempb5
2523  tautobb = tautobb + ssatob*tempb5
2524  CALL popcontrol1b(branch)
2525  IF (branch .EQ. 0) ssatobb = 0.0_8
2526  CALL popreal8(ssatob)
2527  tempb2 = ssatobb/tautob
2528  ssataub = ssataub + tempb2
2529  tautobb = tautobb - (ssatau+tauclb(k))*tempb2/tautob
2530  tauclbb(k) = tauclbb(k) + tautobb + tempb2
2531  taustob = taustob + tautobb
2532  CALL popreal8(ts(k, 1))
2533  tstb = tsb(k, 1)
2534  tsb(k, 1) = 0.0_8
2535  CALL popreal8(rs(k, 1))
2536  rstb = rsb(k, 1)
2537  rsb(k, 1) = 0.0_8
2538  CALL popreal8(td(k, 1))
2539  tdtb = tdb(k, 1)
2540  tdb(k, 1) = 0.0_8
2541  CALL popreal8(tt(k, 1))
2542  tttb = ttb(k, 1)
2543  ttb(k, 1) = 0.0_8
2544  CALL popreal8(rr(k, 1))
2545  rrtb = rrb(k, 1)
2546  rrb(k, 1) = 0.0_8
2547  tautob = tausto
2548  asytob = asysto/ssatau
2549  tautobb = 0.0_8
2550  ssatobb = 0.0_8
2551  asytobb = 0.0_8
2552  dumb = 0.0_8
2553  CALL deledd_b(tautob, tautobb, ssatob, ssatobb, asytob, asytobb, &
2554 & dsm, rst, rstb, tst, tstb, dum, dumb)
2555  CALL deledd_b(tautob, tautobb, ssatob, ssatobb, asytob, asytobb, &
2556 & cosz_dev(i), rrt, rrtb, ttt, tttb, tdt, tdtb)
2557  CALL popcontrol1b(branch)
2558  IF (branch .EQ. 0) ssatobb = 0.0_8
2559  CALL popreal8(ssatob)
2560  ssataub = ssataub + ssatobb/tautob - asysto*asytobb/ssatau**2
2561  tautobb = tautobb - ssatau*ssatobb/tautob**2
2562  asystob = asystob + asytobb/ssatau
2563  CALL popreal8(tautob)
2564  taustob = taustob + tautobb
2565  asya_devb(i, k, ib) = asya_devb(i, k, ib) + asystob
2566  ssaa_devb(i, k, ib) = ssaa_devb(i, k, ib) + ssataub
2567  tauozb = taustob
2568  tauwvb = taustob
2569  taua_devb(i, k, ib) = taua_devb(i, k, ib) + taustob
2570  whb(k) = whb(k) + wk_uv(ib)*tauwvb
2571  ohb(k) = ohb(k) + zk_uv(ib)*tauozb
2572  END DO
2573  CALL popreal8(td(0, 2))
2574  tdb(0, 1) = tdb(0, 1) + tdb(0, 2)
2575  tdb(0, 2) = 0.0_8
2576  CALL popreal8(td(0, 1))
2577  tempb1 = -(exp(-((wk_uv(ib)*wvtoa+zk_uv(ib)*o3toa)/cosz_dev(i)))*tdb&
2578 & (0, 1)/cosz_dev(i))
2579  wvtoab = wvtoab + wk_uv(ib)*tempb1
2580  o3toab = o3toab + zk_uv(ib)*tempb1
2581  tdb(0, 1) = 0.0_8
2582  END DO
2583  taudiffb = 0.0_8
2584  taubeamb = 0.0_8
2585  DO k=np,1,-1
2586  taudiffb(k, 1) = taudiffb(k, 1) + tauclfb(k)
2587  taudiffb(k, 2) = taudiffb(k, 2) + tauclfb(k)
2588  taudiffb(k, 3) = taudiffb(k, 3) + tauclfb(k)
2589  taudiffb(k, 4) = taudiffb(k, 4) + tauclfb(k)
2590  tauclfb(k) = 0.0_8
2591  taubeamb(k, 1) = taubeamb(k, 1) + tauclbb(k)
2592  taubeamb(k, 2) = taubeamb(k, 2) + tauclbb(k)
2593  taubeamb(k, 3) = taubeamb(k, 3) + tauclbb(k)
2594  taubeamb(k, 4) = taubeamb(k, 4) + tauclbb(k)
2595  tauclbb(k) = 0.0_8
2596  END DO
2597  DO k=np,1,-1
2598  CALL popcontrol3b(branch)
2599  IF (branch .LT. 3) THEN
2600  IF (branch .EQ. 0) THEN
2601  fcld_devb(i, k) = fcld_devb(i, k) + cc3b
2602  cc3b = 0.0_8
2603  ELSE IF (branch .NE. 1) THEN
2604  fcld_devb(i, k) = fcld_devb(i, k) + cc2b
2605  cc2b = 0.0_8
2606  END IF
2607  ELSE IF (branch .NE. 3) THEN
2608  IF (branch .EQ. 4) THEN
2609  fcld_devb(i, k) = fcld_devb(i, k) + cc1b
2610  cc1b = 0.0_8
2611  END IF
2612  END IF
2613  END DO
2614  CALL getvistau1_b(np, cosz_dev(i), dp_pa, fcld_col, fcld_colb, &
2615 & reff_col, reff_colb, cwc_col, cwc_colb, ict, icb, taubeam&
2616 & , taubeamb, taudiff, taudiffb, asycl, asyclb, aig_uv, &
2617 & awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, &
2618 & arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, &
2619 & arg_nir, caib, caif, cons_grav)
2620  DO k=np+1,1,-1
2621  flx_devb(i, k) = 0.0_8
2622  END DO
2623  DO k=np,1,-1
2624  DO l=4,1,-1
2625  cwc_devb(i, k, l) = cwc_devb(i, k, l) + cwc_colb(k, l)
2626  cwc_colb(k, l) = 0.0_8
2627  reff_devb(i, k, l) = reff_devb(i, k, l) + reff_colb(k, l)
2628  reff_colb(k, l) = 0.0_8
2629  END DO
2630  fcld_devb(i, k) = fcld_devb(i, k) + fcld_colb(k)
2631  fcld_colb(k) = 0.0_8
2632  oa_devb(i, k) = oa_devb(i, k) + dp(k)*466.7*1.02*ohb(k)
2633  ohb(k) = 0.0_8
2634  swhb(k) = swhb(k) + swhb(k+1)
2635  whb(k) = whb(k) + swhb(k+1)
2636  swhb(k+1) = 0.0_8
2637  tempb0 = scal(k)*1.02*whb(k)
2638  wa_devb(i, k) = wa_devb(i, k) + (0.00135*(ta_dev(i, k)-240.)+1.)*&
2639 & tempb0
2640  ta_devb(i, k) = ta_devb(i, k) + wa_dev(i, k)*0.00135*tempb0
2641  whb(k) = 0.0_8
2642  END DO
2643  wvtoab = wvtoab + swhb(1)
2644  tempb = scal0*1.02*wvtoab
2645  wa_devb(i, 1) = wa_devb(i, 1) + (0.00135*(ta_dev(i, 1)-240.)+1.0)*&
2646 & tempb
2647  ta_devb(i, 1) = ta_devb(i, 1) + wa_dev(i, 1)*0.00135*tempb
2648  oa_devb(i, 1) = oa_devb(i, 1) + xtoa*466.7*1.02*o3toab
2649 END SUBROUTINE sorad_b
2650 
2651 ! Differentiation of deledd in reverse (adjoint) mode:
2652 ! gradient of useful results: g01 tt1 td1 rr1 tau1 ssc1
2653 ! with respect to varying inputs: g01 tau1 ssc1
2654 !*********************************************************************
2655 SUBROUTINE deledd_b(tau1, tau1b, ssc1, ssc1b, g01, g01b, cza1, rr1, rr1b&
2656 & , tt1, tt1b, td1, td1b)
2657  IMPLICIT NONE
2658 ! 8 byte real
2659  INTEGER, PARAMETER :: real_de=8
2660 !integer,parameter :: REAL_SP = 4 ! 4 byte real
2661 !-----input parameters
2662  REAL*8, INTENT(IN) :: tau1, ssc1, g01, cza1
2663  REAL*8 :: tau1b, ssc1b, g01b
2664 !-----output parameters
2665  REAL*8 :: rr1, tt1, td1
2666  REAL*8 :: rr1b, tt1b, td1b
2667 !-----temporary parameters
2668  REAL*8, PARAMETER :: zero=0.0_real_de
2669  REAL*8, PARAMETER :: one=1.0_real_de
2670  REAL*8, PARAMETER :: two=2.0_real_de
2671  REAL*8, PARAMETER :: three=3.0_real_de
2672  REAL*8, PARAMETER :: four=4.0_real_de
2673  REAL*8, PARAMETER :: fourth=0.25_real_de
2674  REAL*8, PARAMETER :: seven=7.0_real_de
2675  REAL*8, PARAMETER :: thresh=1.e-8_real_de
2676  REAL*8 :: tau, ssc, g0, rr, tt, td
2677  REAL*8 :: taub, sscb, g0b, rrb, ttb, tdb
2678  REAL*8 :: zth, ff, xx, taup, sscp, gp, gm1, gm2, gm3, akk, alf1, alf2
2679  REAL*8 :: ffb, xxb, taupb, sscpb, gpb, gm1b, gm2b, gm3b, akkb, alf1b, &
2680 & alf2b
2681  REAL*8 :: all, bll, st7, st8, cll, dll, fll, ell, st1, st2, st3, st4
2682  REAL*8 :: allb, bllb, st7b, st8b, cllb, dllb, fllb, ellb, st1b, st2b, &
2683 & st3b, st4b
2684  INTRINSIC SQRT
2685  INTRINSIC abs
2686  INTRINSIC exp
2687  INTRINSIC max
2688  INTRINSIC real
2689  INTEGER :: branch
2690  REAL*8 :: tempb9
2691  REAL*8 :: tempb8
2692  REAL*8 :: tempb7
2693  REAL*8 :: tempb6
2694  REAL*8 :: tempb5
2695  REAL*8 :: tempb4
2696  REAL*8 :: tempb3
2697  REAL*8 :: tempb2
2698  REAL*8 :: tempb1
2699  REAL*8 :: tempb0
2700  REAL*8 :: tempb10
2701  REAL*8 :: tempb
2702  REAL*8 :: abs0
2703  REAL*8 :: temp
2704 !zth = real(cza1,kind=REAL_DE)
2705 !g0 = real(g01 ,kind=REAL_DE)
2706 !tau = real(tau1,kind=REAL_DE)
2707 !ssc = real(ssc1,kind=REAL_DE)
2708 !zth = dble(cza1)
2709 !g0 = dble(g01)
2710 !tau = dble(tau1)
2711 !ssc = dble(ssc1)
2712  zth = cza1
2713  g0 = g01
2714  tau = tau1
2715  ssc = ssc1
2716  ff = g0*g0
2717  xx = one - ff*ssc
2718  taup = tau*xx
2719  sscp = ssc*(one-ff)/xx
2720  gp = g0/(one+g0)
2721  xx = three*gp
2722  gm1 = (seven-sscp*(four+xx))*fourth
2723  gm2 = -((one-sscp*(four-xx))*fourth)
2724  akk = sqrt((gm1+gm2)*(gm1-gm2))
2725  xx = akk*zth
2726  st7 = one - xx
2727  st8 = one + xx
2728  st3 = st7*st8
2729  IF (st3 .GE. 0.) THEN
2730  abs0 = st3
2731  ELSE
2732  abs0 = -st3
2733  END IF
2734  IF (abs0 .LT. thresh) THEN
2735  CALL pushreal8(zth)
2736  zth = zth + 0.0010
2737  IF (zth .GT. 1.0) zth = zth - 0.0020
2738  xx = akk*zth
2739  CALL pushreal8(st7)
2740  st7 = one - xx
2741  CALL pushreal8(st8)
2742  st8 = one + xx
2743  st3 = st7*st8
2744  CALL pushcontrol1b(1)
2745  ELSE
2746  CALL pushcontrol1b(0)
2747  END IF
2748  td = exp(-(taup/zth))
2749  gm3 = (two-zth*three*gp)*fourth
2750  xx = gm1 - gm2
2751  alf1 = gm1 - gm3*xx
2752  alf2 = gm2 + gm3*xx
2753  xx = akk*two
2754  all = (gm3-alf2*zth)*xx*td
2755  bll = (one-gm3+alf1*zth)*xx
2756  xx = akk*gm3
2757  cll = (alf2+xx)*st7
2758  dll = (alf2-xx)*st8
2759  xx = akk*(one-gm3)
2760  fll = (alf1+xx)*st8
2761  ell = (alf1-xx)*st7
2762  st2 = exp(-(akk*taup))
2763  st4 = st2*st2
2764  st1 = sscp/((akk+gm1+(akk-gm1)*st4)*st3)
2765  rr = (cll-dll*st4-all*st2)*st1
2766  tt = -(((fll-ell*st4)*td-bll*st2)*st1)
2767  IF (rr .LT. zero) THEN
2768  rr = zero
2769  CALL pushcontrol1b(0)
2770  ELSE
2771  CALL pushcontrol1b(1)
2772  rr = rr
2773  END IF
2774  IF (tt .LT. zero) THEN
2775  tt = zero
2776  CALL pushcontrol1b(0)
2777  ELSE
2778  CALL pushcontrol1b(1)
2779  tt = tt
2780  END IF
2781  tt = tt + td
2782 !td1 = real(td,kind=REAL_SP)
2783 !rr1 = real(rr,kind=REAL_SP)
2784 !tt1 = real(tt,kind=REAL_SP)
2785  ttb = tt1b
2786  rrb = rr1b
2787  tdb = ttb + td1b
2788  CALL popcontrol1b(branch)
2789  IF (branch .EQ. 0) ttb = 0.0_8
2790  CALL popcontrol1b(branch)
2791  IF (branch .EQ. 0) rrb = 0.0_8
2792  st1b = (cll-dll*st4-all*st2)*rrb - ((fll-ell*st4)*td-bll*st2)*ttb
2793  tempb4 = st1*rrb
2794  temp = akk + gm1 + (akk-gm1)*st4
2795  tempb7 = st1b/(temp*st3)
2796  tempb8 = -(sscp*tempb7/(temp*st3))
2797  tempb5 = st3*tempb8
2798  tempb2 = -(st1*ttb)
2799  tempb3 = td*tempb2
2800  fllb = tempb3
2801  ellb = -(st4*tempb3)
2802  st4b = (akk-gm1)*tempb5 - dll*tempb4 - ell*tempb3
2803  bllb = -(st2*tempb2)
2804  st2b = 2*st2*st4b - all*tempb4 - bll*tempb2
2805  cllb = tempb4
2806  dllb = -(st4*tempb4)
2807  allb = -(st2*tempb4)
2808  sscpb = tempb7
2809  st3b = temp*tempb8
2810  tempb9 = exp(-(akk*taup))*st2b
2811  xxb = st8*fllb - st7*ellb
2812  akkb = (one-gm3)*xxb - taup*tempb9 + (st4+1.0_8)*tempb5
2813  st7b = (alf1-xx)*ellb
2814  st8b = (alf1+xx)*fllb
2815  gm3b = -(akk*xxb)
2816  xx = akk*gm3
2817  xxb = st7*cllb - st8*dllb
2818  st8b = st8b + (alf2-xx)*dllb
2819  st7b = st7b + (alf2+xx)*cllb
2820  akkb = akkb + gm3*xxb
2821  xx = akk*two
2822  alf1b = st8*fllb + xx*zth*bllb + st7*ellb
2823  gm3b = gm3b + akk*xxb - xx*bllb
2824  tempb10 = xx*td*allb
2825  alf2b = st7*cllb - zth*tempb10 + st8*dllb
2826  tempb6 = (gm3-zth*alf2)*allb
2827  tdb = tdb + xx*tempb6 + (fll-ell*st4)*tempb2
2828  taupb = -(exp(-(taup/zth))*tdb/zth) - akk*tempb9
2829  xxb = td*tempb6 + (one+zth*alf1-gm3)*bllb
2830  akkb = akkb + two*xxb
2831  xx = gm1 - gm2
2832  gm3b = gm3b + xx*alf2b - xx*alf1b + tempb10
2833  xxb = gm3*alf2b - gm3*alf1b
2834  gm1b = alf1b + xxb + (1.0_8-st4)*tempb5
2835  gm2b = alf2b - xxb
2836  gpb = -(fourth*zth*three*gm3b)
2837  CALL popcontrol1b(branch)
2838  IF (branch .NE. 0) THEN
2839  st7b = st7b + st8*st3b
2840  st8b = st8b + st7*st3b
2841  CALL popreal8(st8)
2842  xxb = st8b - st7b
2843  CALL popreal8(st7)
2844  akkb = akkb + zth*xxb
2845  CALL popreal8(zth)
2846  st3b = 0.0_8
2847  st7b = 0.0_8
2848  st8b = 0.0_8
2849  END IF
2850  st7b = st7b + st8*st3b
2851  st8b = st8b + st7*st3b
2852  xxb = st8b - st7b
2853  akkb = akkb + zth*xxb
2854  xx = three*gp
2855  IF ((gm1+gm2)*(gm1-gm2) .EQ. 0.0) THEN
2856  tempb = 0.0
2857  ELSE
2858  tempb = akkb/(2.0*sqrt((gm1+gm2)*(gm1-gm2)))
2859  END IF
2860  gm1b = gm1b + 2*gm1*tempb
2861  gm2b = gm2b - 2*gm2*tempb
2862  sscpb = sscpb + fourth*(four-xx)*gm2b - fourth*(four+xx)*gm1b
2863  xxb = -(fourth*sscp*gm1b) - sscp*fourth*gm2b
2864  gpb = gpb + three*xxb
2865  xx = one - ff*ssc
2866  tempb0 = gpb/(one+g0)
2867  tempb1 = (one-ff)*sscpb/xx
2868  xxb = tau*taupb - ssc*tempb1/xx
2869  ffb = -(ssc*xxb) - ssc*sscpb/xx
2870  g0b = 2*g0*ffb + (1.0_8-g0/(one+g0))*tempb0
2871  sscb = tempb1 - ff*xxb
2872  taub = xx*taupb
2873  ssc1b = ssc1b + sscb
2874  tau1b = tau1b + taub
2875  g01b = g01b + g0b
2876 END SUBROUTINE deledd_b
2877 
2878 ! Differentiation of getvistau1 in reverse (adjoint) mode:
2879 ! gradient of useful results: hydromets asycl taudiff fcld
2880 ! taubeam reff
2881 ! with respect to varying inputs: hydromets fcld reff
2882 SUBROUTINE getvistau1_b(nlevs, cosz, dp, fcld, fcldb, reff, reffb, &
2883 & hydromets, hydrometsb, ict, icb, taubeam, taubeamb, taudiff, taudiffb&
2884 & , asycl, asyclb, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, &
2885 & aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir&
2886 & , arg_nir, caib, caif, cons_grav)
2887  IMPLICIT NONE
2888 !EOC
2889 ! !INPUT PARAMETERS:
2890 ! Number of levels
2891  INTEGER, INTENT(IN) :: nlevs
2892 ! Cosine of solar zenith angle
2893  REAL*8, INTENT(IN) :: cosz
2894 ! Delta pressure (Pa)
2895  REAL*8, INTENT(IN) :: dp(nlevs)
2896 ! Cloud fraction (used sometimes)
2897  REAL*8, INTENT(IN) :: fcld(nlevs)
2898  REAL*8 :: fcldb(nlevs)
2899 ! Effective radius (microns)
2900  REAL*8, INTENT(IN) :: reff(nlevs, 4)
2901  REAL*8 :: reffb(nlevs, 4)
2902 ! Hydrometeors (kg/kg)
2903  REAL*8, INTENT(IN) :: hydromets(nlevs, 4)
2904  REAL*8 :: hydrometsb(nlevs, 4)
2905 ! Flags for various uses
2906  INTEGER, INTENT(IN) :: ict, icb
2907 ! ict = 0 Indicates that in-cloud values have been given
2908 ! and are expected
2909 ! ict != 0 Indicates that overlap computation is needed, and:
2910 ! ict is the level of the mid-high boundary
2911 ! icb is the level of the low-mid boundary
2912 !
2913 ! !OUTPUT PARAMETERS:
2914 ! Optical Depth for Beam Radiation
2915  REAL*8 :: taubeam(nlevs, 4)
2916  REAL*8 :: taubeamb(nlevs, 4)
2917 ! Optical Depth for Diffuse Radiation
2918  REAL*8 :: taudiff(nlevs, 4)
2919  REAL*8 :: taudiffb(nlevs, 4)
2920 ! Cloud Asymmetry Factor
2921  REAL*8 :: asycl(nlevs)
2922  REAL*8 :: asyclb(nlevs)
2923 ! !DESCRIPTION:
2924 ! Compute in-cloud or grid mean optical depths for visible wavelengths
2925 ! In general will compute in-cloud - will do grid mean when called
2926 ! for diagnostic use only. ict flag will indicate which to do.
2927 ! Slots for reff, hydrometeors, taubeam, taudiff, and asycl are as follows:
2928 ! 1 Cloud Ice
2929 ! 2 Cloud Liquid
2930 ! 3 Falling Liquid (Rain)
2931 ! 4 Falling Ice (Snow)
2932 !
2933 ! In the below calculations, the constants used in the tau calculation are in
2934 ! m$^2$ g$^{-1}$ and m$^2$ g$^{-1}$ $\mu$m. Thus, we must convert the kg contained in the
2935 ! pressure (Pa = kg m$^{-1}$ s$^{-2}$) to grams.
2936 !
2937 ! !REVISION HISTORY:
2938 ! 2011.10.27 Molod moved to Radiation_Shared and revised arg list, units
2939 ! 2011.11.16 MAT: Generalized to a call that is per-column
2940 !
2941 !EOP
2942 !------------------------------------------------------------------------------
2943 !BOC
2944  INTEGER :: k, in, im, it, ia, kk
2945  REAL*8 :: fm, ft, fa, xai, tauc, asyclt
2946  REAL*8 :: ftb, fab, xaib, taucb, asycltb
2947  REAL*8 :: cc(3)
2948  REAL*8 :: ccb(3)
2949  REAL*8 :: taucld1, taucld2, taucld3, taucld4
2950  REAL*8 :: taucld1b, taucld2b, taucld3b, taucld4b
2951  REAL*8 :: g1, g2, g3, g4
2952  REAL*8 :: g1b, g2b, g3b, g4b
2953  REAL*8 :: reff_snow
2954  REAL*8 :: reff_snowb
2955  INTEGER, PARAMETER :: nm=11, nt=9, na=11
2956  REAL*8, PARAMETER :: dm=0.1, dt=0.30103, da=0.1, t1=-0.9031
2957  REAL*8, INTENT(IN) :: aig_uv(3), awg_uv(3), arg_uv(3)
2958  REAL*8, INTENT(IN) :: aib_uv, awb_uv(2), arb_uv(2)
2959  REAL*8, INTENT(IN) :: aib_nir, awb_nir(3, 2), arb_nir(3, 2)
2960  REAL*8, INTENT(IN) :: aia_nir(3, 3), awa_nir(3, 3), ara_nir(3, 3)
2961  REAL*8, INTENT(IN) :: aig_nir(3, 3), awg_nir(3, 3), arg_nir(3, 3)
2962  REAL*8, INTENT(IN) :: caib(11, 9, 11), caif(9, 11)
2963  REAL*8, INTENT(IN) :: cons_grav
2964  INTRINSIC max
2965  INTRINSIC min
2966  INTRINSIC log10
2967  INTRINSIC int
2968  INTRINSIC real
2969  INTEGER :: branch
2970  REAL*8 :: temp3
2971  REAL*8 :: temp2
2972  REAL*8 :: temp1
2973  REAL*8 :: temp0
2974  REAL*8 :: tempb7
2975  REAL*8 :: tempb6
2976  REAL*8 :: tempb5
2977  REAL*8 :: tempb4
2978  REAL*8 :: tempb3
2979  REAL*8 :: tempb2
2980  REAL*8 :: tempb1
2981  REAL*8 :: tempb0
2982  REAL*8 :: tempb
2983  REAL*8 :: temp
2984  IF (ict .NE. 0) THEN
2985 !-----scale cloud optical thickness in each layer from taucld (with
2986 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
2987 ! taubeam is the scaled optical thickness for beam radiation and
2988 ! taudiff is for diffuse radiation (see section 7).
2989 !-----clouds within each of the high, middle, and low clouds are
2990 ! assumed to be maximally overlapped, and the cloud cover (cc)
2991 ! for a group (high, middle, or low) is the maximum cloud cover
2992 ! of all the layers within a group
2993  cc = 0.0
2994  DO k=1,ict-1
2995  IF (cc(1) .LT. fcld(k)) THEN
2996  cc(1) = fcld(k)
2997  CALL pushcontrol1b(0)
2998  ELSE
2999  CALL pushcontrol1b(1)
3000  cc(1) = cc(1)
3001  END IF
3002  END DO
3003  DO k=ict,icb-1
3004  IF (cc(2) .LT. fcld(k)) THEN
3005  cc(2) = fcld(k)
3006  CALL pushcontrol1b(0)
3007  ELSE
3008  CALL pushcontrol1b(1)
3009  cc(2) = cc(2)
3010  END IF
3011  END DO
3012  DO k=icb,nlevs
3013  IF (cc(3) .LT. fcld(k)) THEN
3014  cc(3) = fcld(k)
3015  CALL pushcontrol1b(0)
3016  ELSE
3017  CALL pushcontrol1b(1)
3018  cc(3) = cc(3)
3019  END IF
3020  END DO
3021  CALL pushcontrol1b(1)
3022  ELSE
3023  CALL pushcontrol1b(0)
3024  END IF
3025 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10)
3026 ! Note: the cloud optical properties are assumed to be independent
3027 ! of spectral bands in the UV and PAR regions.
3028 ! taucld1 is the optical thickness for ice particles
3029 ! taucld2 is the optical thickness for liquid particles
3030 ! taucld3 is the optical thickness for rain drops
3031 ! taucld4 is the optical thickness for snow
3032  DO k=1,nlevs
3033  IF (reff(k, 1) .LE. 0.) THEN
3034  CALL pushreal8(taucld1)
3035  taucld1 = 0.
3036  CALL pushcontrol1b(0)
3037  ELSE
3038  CALL pushreal8(taucld1)
3039  taucld1 = dp(k)*1.0e3/cons_grav*hydromets(k, 1)*aib_uv/reff(k, 1)
3040  CALL pushcontrol1b(1)
3041  END IF
3042  IF (reff(k, 2) .LE. 0.) THEN
3043  CALL pushreal8(taucld2)
3044  taucld2 = 0.
3045  CALL pushcontrol1b(1)
3046  ELSE
3047  CALL pushreal8(taucld2)
3048  taucld2 = dp(k)*1.0e3/cons_grav*hydromets(k, 2)*(awb_uv(1)+awb_uv(&
3049 & 2)/reff(k, 2))
3050  CALL pushcontrol1b(0)
3051  END IF
3052  CALL pushreal8(taucld3)
3053  taucld3 = dp(k)*1.0e3/cons_grav*hydromets(k, 3)*arb_uv(1)
3054  IF (reff(k, 4) .GT. 112.0) THEN
3055  CALL pushreal8(reff_snow)
3056  reff_snow = 112.0
3057  CALL pushcontrol1b(0)
3058  ELSE
3059  CALL pushreal8(reff_snow)
3060  reff_snow = reff(k, 4)
3061  CALL pushcontrol1b(1)
3062  END IF
3063  IF (reff_snow .LE. 0.) THEN
3064  CALL pushreal8(taucld4)
3065  taucld4 = 0.
3066  CALL pushcontrol1b(0)
3067  ELSE
3068  CALL pushreal8(taucld4)
3069  taucld4 = dp(k)*1.0e3/cons_grav*hydromets(k, 4)*aib_uv/reff_snow
3070  CALL pushcontrol1b(1)
3071  END IF
3072  IF (ict .NE. 0) THEN
3073 !-----scale cloud optical thickness in each layer from taucld (with
3074 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
3075 ! taubeam is the scaled optical thickness for beam radiation and
3076 ! taudiff is for diffuse radiation (see section 7).
3077 !-----clouds within each of the high, middle, and low clouds are
3078 ! assumed to be maximally overlapped, and the cloud cover (cc)
3079 ! for a group (high, middle, or low) is the maximum cloud cover
3080 ! of all the layers within a group
3081  IF (k .LT. ict) THEN
3082  CALL pushinteger4(kk)
3083  kk = 1
3084  CALL pushcontrol2b(0)
3085  ELSE IF (k .GE. ict .AND. k .LT. icb) THEN
3086  CALL pushinteger4(kk)
3087  kk = 2
3088  CALL pushcontrol2b(1)
3089  ELSE
3090  CALL pushinteger4(kk)
3091  kk = 3
3092  CALL pushcontrol2b(2)
3093  END IF
3094  tauc = taucld1 + taucld2 + taucld3 + taucld4
3095  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
3096 !-----normalize cloud cover following Eq. (7.8)
3097  CALL pushreal8(fa)
3098  fa = fcld(k)/cc(kk)
3099  IF (tauc .GT. 32.) THEN
3100  tauc = 32.
3101  CALL pushcontrol1b(0)
3102  ELSE
3103  CALL pushcontrol1b(1)
3104  tauc = tauc
3105  END IF
3106  fm = cosz/dm
3107  CALL pushreal8(ft)
3108  ft = (log10(tauc)-t1)/dt
3109  fa = fa/da
3110  CALL pushinteger4(im)
3111  im = int(fm + 1.5)
3112  CALL pushinteger4(it)
3113  it = int(ft + 1.5)
3114  CALL pushinteger4(ia)
3115  ia = int(fa + 1.5)
3116  IF (im .LT. 2) THEN
3117  im = 2
3118  ELSE
3119  im = im
3120  END IF
3121  IF (it .LT. 2) THEN
3122  it = 2
3123  ELSE
3124  it = it
3125  END IF
3126  IF (ia .LT. 2) THEN
3127  ia = 2
3128  ELSE
3129  ia = ia
3130  END IF
3131  IF (im .GT. nm - 1) THEN
3132  im = nm - 1
3133  ELSE
3134  im = im
3135  END IF
3136  IF (it .GT. nt - 1) THEN
3137  it = nt - 1
3138  ELSE
3139  it = it
3140  END IF
3141  IF (ia .GT. na - 1) THEN
3142  ia = na - 1
3143  ELSE
3144  ia = ia
3145  END IF
3146  fm = fm - REAL(im - 1)
3147  ft = ft - REAL(it - 1)
3148  fa = fa - REAL(ia - 1)
3149 !-----scale cloud optical thickness for beam radiation following
3150 ! Eq. (7.3).
3151 ! the scaling factor, xai, is a function of the solar zenith
3152 ! angle, optical thickness, and cloud cover.
3153  CALL pushreal8(xai)
3154  xai = (-(caib(im-1, it, ia)*(1.-fm))+caib(im+1, it, ia)*(1.+fm))&
3155 & *fm*.5 + caib(im, it, ia)*(1.-fm*fm)
3156  xai = xai + (-(caib(im, it-1, ia)*(1.-ft))+caib(im, it+1, ia)*(&
3157 & 1.+ft))*ft*.5 + caib(im, it, ia)*(1.-ft*ft)
3158  xai = xai + (-(caib(im, it, ia-1)*(1.-fa))+caib(im, it, ia+1)*(&
3159 & 1.+fa))*fa*.5 + caib(im, it, ia)*(1.-fa*fa)
3160  xai = xai - 2.*caib(im, it, ia)
3161  IF (xai .LT. 0.0) THEN
3162  xai = 0.0
3163  CALL pushcontrol1b(0)
3164  ELSE
3165  CALL pushcontrol1b(1)
3166  xai = xai
3167  END IF
3168  IF (xai .GT. 1.0) THEN
3169  xai = 1.0
3170  CALL pushcontrol1b(0)
3171  ELSE
3172  CALL pushcontrol1b(1)
3173  xai = xai
3174  END IF
3175 !-----scale cloud optical thickness for diffuse radiation following
3176 ! Eq. (7.4).
3177 ! the scaling factor, xai, is a function of the cloud optical
3178 ! thickness and cover but not the solar zenith angle.
3179  CALL pushreal8(xai)
3180  xai = (-(caif(it-1, ia)*(1.-ft))+caif(it+1, ia)*(1.+ft))*ft*.5 +&
3181 & caif(it, ia)*(1.-ft*ft)
3182  xai = xai + (-(caif(it, ia-1)*(1.-fa))+caif(it, ia+1)*(1.+fa))*&
3183 & fa*.5 + caif(it, ia)*(1.-fa*fa)
3184  xai = xai - caif(it, ia)
3185  IF (xai .LT. 0.0) THEN
3186  xai = 0.0
3187  CALL pushcontrol1b(0)
3188  ELSE
3189  CALL pushcontrol1b(1)
3190  xai = xai
3191  END IF
3192  IF (xai .GT. 1.0) THEN
3193  xai = 1.0
3194  CALL pushcontrol1b(0)
3195  ELSE
3196  CALL pushcontrol1b(1)
3197  xai = xai
3198  END IF
3199  CALL pushcontrol2b(0)
3200  ELSE
3201  CALL pushcontrol2b(1)
3202  END IF
3203  ELSE
3204  CALL pushcontrol2b(2)
3205  END IF
3206 !-----cloud asymmetry factor for a mixture of liquid and ice particles.
3207 ! unit of reff is micrometers. Eqs. (4.8) and (6.4)
3208  CALL pushreal8(tauc)
3209  tauc = taucld1 + taucld2 + taucld3 + taucld4
3210  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
3211  CALL pushcontrol1b(0)
3212  ELSE
3213  CALL pushcontrol1b(1)
3214  END IF
3215  END DO
3216  ccb = 0.0_8
3217  DO k=nlevs,1,-1
3218  asycltb = asyclb(k)
3219  asyclb(k) = 0.0_8
3220  CALL popcontrol1b(branch)
3221  IF (branch .EQ. 0) THEN
3222  g1 = (aig_uv(1)+(aig_uv(2)+aig_uv(3)*reff(k, 1))*reff(k, 1))*&
3223 & taucld1
3224  g2 = (awg_uv(1)+(awg_uv(2)+awg_uv(3)*reff(k, 2))*reff(k, 2))*&
3225 & taucld2
3226  taucld3 = dp(k)*1.0e3/cons_grav*hydromets(k, 3)*arb_uv(1)
3227  g3 = arg_uv(1)*taucld3
3228  g4 = (aig_uv(1)+(aig_uv(2)+aig_uv(3)*reff_snow)*reff_snow)*taucld4
3229  tauc = taucld1 + taucld2 + taucld3 + taucld4
3230  tempb7 = asycltb/tauc
3231  g1b = tempb7
3232  g2b = tempb7
3233  g3b = tempb7
3234  g4b = tempb7
3235  taucb = -((g1+g2+g3+g4)*tempb7/tauc)
3236  temp3 = aig_uv(2) + aig_uv(3)*reff_snow
3237  reff_snowb = (taucld4*temp3+reff_snow*taucld4*aig_uv(3))*g4b
3238  taucld4b = (aig_uv(1)+temp3*reff_snow)*g4b
3239  taucld3b = arg_uv(1)*g3b
3240  temp2 = awg_uv(2) + awg_uv(3)*reff(k, 2)
3241  reffb(k, 2) = reffb(k, 2) + (taucld2*temp2+reff(k, 2)*taucld2*&
3242 & awg_uv(3))*g2b
3243  taucld2b = (awg_uv(1)+temp2*reff(k, 2))*g2b
3244  temp1 = aig_uv(2) + aig_uv(3)*reff(k, 1)
3245  reffb(k, 1) = reffb(k, 1) + (taucld1*temp1+reff(k, 1)*taucld1*&
3246 & aig_uv(3))*g1b
3247  taucld1b = (aig_uv(1)+temp1*reff(k, 1))*g1b
3248  ELSE
3249  reff_snowb = 0.0_8
3250  taucld1b = 0.0_8
3251  taucld2b = 0.0_8
3252  taucld3b = 0.0_8
3253  taucld4b = 0.0_8
3254  taucb = 0.0_8
3255  END IF
3256  CALL popreal8(tauc)
3257  taucld1b = taucld1b + taucb
3258  taucld2b = taucld2b + taucb
3259  taucld3b = taucld3b + taucb
3260  taucld4b = taucld4b + taucb
3261  CALL popcontrol2b(branch)
3262  IF (branch .EQ. 0) THEN
3263  taucld4b = taucld4b + xai*taudiffb(k, 4)
3264  xaib = taucld4*taudiffb(k, 4)
3265  taudiffb(k, 4) = 0.0_8
3266  taucld3b = taucld3b + xai*taudiffb(k, 3)
3267  xaib = xaib + taucld3*taudiffb(k, 3)
3268  taudiffb(k, 3) = 0.0_8
3269  taucld2b = taucld2b + xai*taudiffb(k, 2)
3270  xaib = xaib + taucld2*taudiffb(k, 2)
3271  taudiffb(k, 2) = 0.0_8
3272  taucld1b = taucld1b + xai*taudiffb(k, 1)
3273  xaib = xaib + taucld1*taudiffb(k, 1)
3274  taudiffb(k, 1) = 0.0_8
3275  CALL popcontrol1b(branch)
3276  IF (branch .EQ. 0) xaib = 0.0_8
3277  CALL popcontrol1b(branch)
3278  IF (branch .EQ. 0) xaib = 0.0_8
3279  tempb5 = .5*fa*xaib
3280  fab = (.5*(caif(it, ia+1)*(fa+1.)-caif(it, ia-1)*(1.-fa))-caif(it&
3281 & , ia)*2*fa)*xaib + (caif(it, ia-1)+caif(it, ia+1))*tempb5
3282  CALL popreal8(xai)
3283  tempb6 = .5*ft*xaib
3284  ftb = (.5*(caif(it+1, ia)*(ft+1.)-caif(it-1, ia)*(1.-ft))-caif(it&
3285 & , ia)*2*ft)*xaib + (caif(it-1, ia)+caif(it+1, ia))*tempb6
3286  taucld4b = taucld4b + xai*taubeamb(k, 4)
3287  xaib = taucld4*taubeamb(k, 4)
3288  taubeamb(k, 4) = 0.0_8
3289  taucld3b = taucld3b + xai*taubeamb(k, 3)
3290  xaib = xaib + taucld3*taubeamb(k, 3)
3291  taubeamb(k, 3) = 0.0_8
3292  taucld2b = taucld2b + xai*taubeamb(k, 2)
3293  xaib = xaib + taucld2*taubeamb(k, 2)
3294  taubeamb(k, 2) = 0.0_8
3295  taucld1b = taucld1b + xai*taubeamb(k, 1)
3296  xaib = xaib + taucld1*taubeamb(k, 1)
3297  taubeamb(k, 1) = 0.0_8
3298  CALL popcontrol1b(branch)
3299  IF (branch .EQ. 0) xaib = 0.0_8
3300  CALL popcontrol1b(branch)
3301  IF (branch .EQ. 0) xaib = 0.0_8
3302  tempb3 = .5*fa*xaib
3303  fab = fab + (.5*(caib(im, it, ia+1)*(fa+1.)-caib(im, it, ia-1)*(1.&
3304 & -fa))-caib(im, it, ia)*2*fa)*xaib + (caib(im, it, ia-1)+caib(im&
3305 & , it, ia+1))*tempb3
3306  tempb4 = .5*ft*xaib
3307  ftb = ftb + (.5*(caib(im, it+1, ia)*(ft+1.)-caib(im, it-1, ia)*(1.&
3308 & -ft))-caib(im, it, ia)*2*ft)*xaib + (caib(im, it-1, ia)+caib(im&
3309 & , it+1, ia))*tempb4
3310  CALL popreal8(xai)
3311  CALL popinteger4(ia)
3312  CALL popinteger4(it)
3313  CALL popinteger4(im)
3314  fab = fab/da
3315  CALL popreal8(ft)
3316  taucb = ftb/(dt*tauc*log(10.0))
3317  CALL popcontrol1b(branch)
3318  IF (branch .EQ. 0) taucb = 0.0_8
3319  CALL popreal8(fa)
3320  tempb2 = fab/cc(kk)
3321  fcldb(k) = fcldb(k) + tempb2
3322  ccb(kk) = ccb(kk) - fcld(k)*tempb2/cc(kk)
3323  ELSE IF (branch .EQ. 1) THEN
3324  taucb = 0.0_8
3325  ELSE
3326  taucld4b = taucld4b + taubeamb(k, 4) + taudiffb(k, 4)
3327  taudiffb(k, 4) = 0.0_8
3328  taubeamb(k, 4) = 0.0_8
3329  taucld3b = taucld3b + taubeamb(k, 3) + taudiffb(k, 3)
3330  taudiffb(k, 3) = 0.0_8
3331  taubeamb(k, 3) = 0.0_8
3332  taucld2b = taucld2b + taubeamb(k, 2) + taudiffb(k, 2)
3333  taudiffb(k, 2) = 0.0_8
3334  taubeamb(k, 2) = 0.0_8
3335  taucld1b = taucld1b + taubeamb(k, 1) + taudiffb(k, 1)
3336  taudiffb(k, 1) = 0.0_8
3337  taubeamb(k, 1) = 0.0_8
3338  GOTO 100
3339  END IF
3340  taucld1b = taucld1b + taucb
3341  taucld2b = taucld2b + taucb
3342  taucld3b = taucld3b + taucb
3343  taucld4b = taucld4b + taucb
3344  CALL popcontrol2b(branch)
3345  IF (branch .EQ. 0) THEN
3346  CALL popinteger4(kk)
3347  ELSE IF (branch .EQ. 1) THEN
3348  CALL popinteger4(kk)
3349  ELSE
3350  CALL popinteger4(kk)
3351  END IF
3352  100 CALL popcontrol1b(branch)
3353  IF (branch .EQ. 0) THEN
3354  CALL popreal8(taucld4)
3355  ELSE
3356  CALL popreal8(taucld4)
3357  tempb1 = dp(k)*aib_uv*1.0e3*taucld4b/(cons_grav*reff_snow)
3358  hydrometsb(k, 4) = hydrometsb(k, 4) + tempb1
3359  reff_snowb = reff_snowb - hydromets(k, 4)*tempb1/reff_snow
3360  END IF
3361  CALL popcontrol1b(branch)
3362  IF (branch .EQ. 0) THEN
3363  CALL popreal8(reff_snow)
3364  ELSE
3365  CALL popreal8(reff_snow)
3366  reffb(k, 4) = reffb(k, 4) + reff_snowb
3367  END IF
3368  CALL popreal8(taucld3)
3369  hydrometsb(k, 3) = hydrometsb(k, 3) + dp(k)*arb_uv(1)*1.0e3*taucld3b&
3370 & /cons_grav
3371  CALL popcontrol1b(branch)
3372  IF (branch .EQ. 0) THEN
3373  CALL popreal8(taucld2)
3374  temp0 = awb_uv(2)/reff(k, 2)
3375  tempb0 = dp(k)*1.0e3*taucld2b
3376  hydrometsb(k, 2) = hydrometsb(k, 2) + (awb_uv(1)+temp0)*tempb0/&
3377 & cons_grav
3378  reffb(k, 2) = reffb(k, 2) - hydromets(k, 2)*temp0*tempb0/(reff(k, &
3379 & 2)*cons_grav)
3380  ELSE
3381  CALL popreal8(taucld2)
3382  END IF
3383  CALL popcontrol1b(branch)
3384  IF (branch .EQ. 0) THEN
3385  CALL popreal8(taucld1)
3386  ELSE
3387  CALL popreal8(taucld1)
3388  temp = cons_grav*reff(k, 1)
3389  tempb = dp(k)*aib_uv*1.0e3*taucld1b/temp
3390  hydrometsb(k, 1) = hydrometsb(k, 1) + tempb
3391  reffb(k, 1) = reffb(k, 1) - hydromets(k, 1)*cons_grav*tempb/temp
3392  END IF
3393  END DO
3394  CALL popcontrol1b(branch)
3395  IF (branch .NE. 0) THEN
3396  DO k=nlevs,icb,-1
3397  CALL popcontrol1b(branch)
3398  IF (branch .EQ. 0) THEN
3399  fcldb(k) = fcldb(k) + ccb(3)
3400  ccb(3) = 0.0_8
3401  END IF
3402  END DO
3403  DO k=icb-1,ict,-1
3404  CALL popcontrol1b(branch)
3405  IF (branch .EQ. 0) THEN
3406  fcldb(k) = fcldb(k) + ccb(2)
3407  ccb(2) = 0.0_8
3408  END IF
3409  END DO
3410  DO k=ict-1,1,-1
3411  CALL popcontrol1b(branch)
3412  IF (branch .EQ. 0) THEN
3413  fcldb(k) = fcldb(k) + ccb(1)
3414  ccb(1) = 0.0_8
3415  END IF
3416  END DO
3417  END IF
3418 END SUBROUTINE getvistau1_b
3419 
3420 ! Differentiation of getnirtau1 in reverse (adjoint) mode:
3421 ! gradient of useful results: hydromets asycl taudiff fcld
3422 ! ssacl taubeam reff
3423 ! with respect to varying inputs: hydromets asycl fcld ssacl
3424 ! reff
3425 SUBROUTINE getnirtau1_b(ib, nlevs, cosz, dp, fcld, fcldb, reff, reffb, &
3426 & hydromets, hydrometsb, ict, icb, taubeam, taubeamb, taudiff, taudiffb&
3427 & , asycl, asyclb, ssacl, ssaclb, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv&
3428 & , arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, &
3429 & aig_nir, awg_nir, arg_nir, caib, caif, cons_grav)
3430  IMPLICIT NONE
3431 ! !INPUT PARAMETERS:
3432 ! Band number
3433  INTEGER, INTENT(IN) :: ib
3434 ! Number of levels
3435  INTEGER, INTENT(IN) :: nlevs
3436 ! Cosine of solar zenith angle
3437  REAL*8, INTENT(IN) :: cosz
3438 ! Delta pressure in Pa
3439  REAL*8, INTENT(IN) :: dp(nlevs)
3440 ! Cloud fraction (used sometimes)
3441  REAL*8, INTENT(IN) :: fcld(nlevs)
3442  REAL*8 :: fcldb(nlevs)
3443 ! Effective radius (microns)
3444  REAL*8, INTENT(IN) :: reff(nlevs, 4)
3445  REAL*8 :: reffb(nlevs, 4)
3446 ! Hydrometeors (kg/kg)
3447  REAL*8, INTENT(IN) :: hydromets(nlevs, 4)
3448  REAL*8 :: hydrometsb(nlevs, 4)
3449 ! Flags for various uses
3450  INTEGER, INTENT(IN) :: ict, icb
3451  REAL*8, INTENT(IN) :: aig_uv(3), awg_uv(3), arg_uv(3)
3452  REAL*8, INTENT(IN) :: aib_uv, awb_uv(2), arb_uv(2)
3453  REAL*8, INTENT(IN) :: aib_nir, awb_nir(3, 2), arb_nir(3, 2)
3454  REAL*8, INTENT(IN) :: aia_nir(3, 3), awa_nir(3, 3), ara_nir(3, 3)
3455  REAL*8, INTENT(IN) :: aig_nir(3, 3), awg_nir(3, 3), arg_nir(3, 3)
3456  REAL*8, INTENT(IN) :: caib(11, 9, 11), caif(9, 11)
3457  REAL*8, INTENT(IN) :: cons_grav
3458 ! !OUTPUT PARAMETERS:
3459 ! Optical depth for beam radiation
3460  REAL*8 :: taubeam(nlevs, 4)
3461  REAL*8 :: taubeamb(nlevs, 4)
3462 ! Optical depth for diffuse radiation
3463  REAL*8 :: taudiff(nlevs, 4)
3464  REAL*8 :: taudiffb(nlevs, 4)
3465 ! Cloud single scattering albedo
3466  REAL*8 :: ssacl(nlevs)
3467  REAL*8 :: ssaclb(nlevs)
3468 ! Cloud asymmetry factor
3469  REAL*8 :: asycl(nlevs)
3470  REAL*8 :: asyclb(nlevs)
3471  INTEGER :: k, in, im, it, ia, kk
3472  REAL*8 :: fm, ft, fa, xai, tauc, asyclt, ssaclt
3473  REAL*8 :: ftb, fab, xaib, taucb, asycltb, ssacltb
3474  REAL*8 :: cc(3)
3475  REAL*8 :: ccb(3)
3476  REAL*8 :: taucld1, taucld2, taucld3, taucld4
3477  REAL*8 :: taucld1b, taucld2b, taucld3b, taucld4b
3478  REAL*8 :: g1, g2, g3, g4
3479  REAL*8 :: g1b, g2b, g3b, g4b
3480  REAL*8 :: w1, w2, w3, w4
3481  REAL*8 :: w1b, w2b, w3b, w4b
3482  REAL*8 :: reff_snow
3483  REAL*8 :: reff_snowb
3484  INTEGER, PARAMETER :: nm=11, nt=9, na=11
3485  REAL*8, PARAMETER :: dm=0.1, dt=0.30103, da=0.1, t1=-0.9031
3486  INTRINSIC max
3487  INTRINSIC min
3488  INTRINSIC log10
3489  INTRINSIC int
3490  INTRINSIC real
3491  INTEGER :: branch
3492  REAL*8 :: temp3
3493  REAL*8 :: temp2
3494  REAL*8 :: temp1
3495  REAL*8 :: temp0
3496  REAL*8 :: tempb9
3497  REAL*8 :: tempb8
3498  REAL*8 :: tempb7
3499  REAL*8 :: tempb6
3500  REAL*8 :: tempb5
3501  REAL*8 :: tempb4
3502  REAL*8 :: tempb3
3503  REAL*8 :: tempb2
3504  REAL*8 :: tempb1
3505  REAL*8 :: tempb0
3506  REAL*8 :: tempb
3507  REAL*8 :: temp
3508  REAL*8 :: temp6
3509  REAL*8 :: temp5
3510  REAL*8 :: temp4
3511  IF (ict .NE. 0) THEN
3512 !-----scale cloud optical thickness in each layer from taucld (with
3513 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
3514 ! taubeam is the scaled optical thickness for beam radiation and
3515 ! taudiff is for diffuse radiation (see section 7).
3516 !-----clouds within each of the high, middle, and low clouds are
3517 ! assumed to be maximally overlapped, and the cloud cover (cc)
3518 ! for a group (high, middle, or low) is the maximum cloud cover
3519 ! of all the layers within a group
3520  cc = 0.0
3521  DO k=1,ict-1
3522  IF (cc(1) .LT. fcld(k)) THEN
3523  cc(1) = fcld(k)
3524  CALL pushcontrol1b(0)
3525  ELSE
3526  CALL pushcontrol1b(1)
3527  cc(1) = cc(1)
3528  END IF
3529  END DO
3530  DO k=ict,icb-1
3531  IF (cc(2) .LT. fcld(k)) THEN
3532  cc(2) = fcld(k)
3533  CALL pushcontrol1b(0)
3534  ELSE
3535  CALL pushcontrol1b(1)
3536  cc(2) = cc(2)
3537  END IF
3538  END DO
3539  DO k=icb,nlevs
3540  IF (cc(3) .LT. fcld(k)) THEN
3541  cc(3) = fcld(k)
3542  CALL pushcontrol1b(0)
3543  ELSE
3544  CALL pushcontrol1b(1)
3545  cc(3) = cc(3)
3546  END IF
3547  END DO
3548  CALL pushcontrol1b(1)
3549  ELSE
3550  CALL pushcontrol1b(0)
3551  END IF
3552 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10)
3553 ! taucld1 is the optical thickness for ice particles
3554 ! taucld2 is the optical thickness for liquid particles
3555 ! taucld3 is the optical thickness for rain drops
3556 ! taucld4 is the optical thickness for snow
3557  DO k=1,nlevs
3558  IF (reff(k, 1) .LE. 0.) THEN
3559  CALL pushreal8(taucld1)
3560  taucld1 = 0.
3561  CALL pushcontrol1b(0)
3562  ELSE
3563  CALL pushreal8(taucld1)
3564  taucld1 = dp(k)*1.0e3/cons_grav*hydromets(k, 1)*aib_nir/reff(k, 1)
3565  CALL pushcontrol1b(1)
3566  END IF
3567  IF (reff(k, 2) .LE. 0.) THEN
3568  CALL pushreal8(taucld2)
3569  taucld2 = 0.
3570  CALL pushcontrol1b(1)
3571  ELSE
3572  CALL pushreal8(taucld2)
3573  taucld2 = dp(k)*1.0e3/cons_grav*hydromets(k, 2)*(awb_nir(ib, 1)+&
3574 & awb_nir(ib, 2)/reff(k, 2))
3575  CALL pushcontrol1b(0)
3576  END IF
3577  CALL pushreal8(taucld3)
3578  taucld3 = dp(k)*1.0e3/cons_grav*hydromets(k, 3)*arb_nir(ib, 1)
3579  IF (reff(k, 4) .GT. 112.0) THEN
3580  CALL pushreal8(reff_snow)
3581  reff_snow = 112.0
3582  CALL pushcontrol1b(0)
3583  ELSE
3584  CALL pushreal8(reff_snow)
3585  reff_snow = reff(k, 4)
3586  CALL pushcontrol1b(1)
3587  END IF
3588  IF (reff_snow .LE. 0.) THEN
3589  CALL pushreal8(taucld4)
3590  taucld4 = 0.
3591  CALL pushcontrol1b(0)
3592  ELSE
3593  CALL pushreal8(taucld4)
3594  taucld4 = dp(k)*1.0e3/cons_grav*hydromets(k, 4)*aib_nir/reff_snow
3595  CALL pushcontrol1b(1)
3596  END IF
3597  IF (ict .NE. 0) THEN
3598 !-----scale cloud optical thickness in each layer from taucld (with
3599 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
3600 ! taubeam is the scaled optical thickness for beam radiation and
3601 ! taudiff is for diffuse radiation (see section 7).
3602 !-----clouds within each of the high, middle, and low clouds are
3603 ! assumed to be maximally overlapped, and the cloud cover (cc)
3604 ! for a group (high, middle, or low) is the maximum cloud cover
3605 ! of all the layers within a group
3606  IF (k .LT. ict) THEN
3607  CALL pushinteger4(kk)
3608  kk = 1
3609  CALL pushcontrol2b(0)
3610  ELSE IF (k .GE. ict .AND. k .LT. icb) THEN
3611  CALL pushinteger4(kk)
3612  kk = 2
3613  CALL pushcontrol2b(1)
3614  ELSE
3615  CALL pushinteger4(kk)
3616  kk = 3
3617  CALL pushcontrol2b(2)
3618  END IF
3619  tauc = taucld1 + taucld2 + taucld3 + taucld4
3620  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
3621 !-----normalize cloud cover following Eq. (7.8)
3622  IF (cc(kk) .NE. 0.0) THEN
3623  CALL pushreal8(fa)
3624  fa = fcld(k)/cc(kk)
3625  CALL pushcontrol1b(1)
3626  ELSE
3627  CALL pushreal8(fa)
3628  fa = 0.0
3629  CALL pushcontrol1b(0)
3630  END IF
3631  IF (tauc .GT. 32.) THEN
3632  tauc = 32.
3633  CALL pushcontrol1b(0)
3634  ELSE
3635  CALL pushcontrol1b(1)
3636  tauc = tauc
3637  END IF
3638  fm = cosz/dm
3639  CALL pushreal8(ft)
3640  ft = (log10(tauc)-t1)/dt
3641  fa = fa/da
3642  CALL pushinteger4(im)
3643  im = int(fm + 1.5)
3644  CALL pushinteger4(it)
3645  it = int(ft + 1.5)
3646  CALL pushinteger4(ia)
3647  ia = int(fa + 1.5)
3648  IF (im .LT. 2) THEN
3649  im = 2
3650  ELSE
3651  im = im
3652  END IF
3653  IF (it .LT. 2) THEN
3654  it = 2
3655  ELSE
3656  it = it
3657  END IF
3658  IF (ia .LT. 2) THEN
3659  ia = 2
3660  ELSE
3661  ia = ia
3662  END IF
3663  IF (im .GT. nm - 1) THEN
3664  im = nm - 1
3665  ELSE
3666  im = im
3667  END IF
3668  IF (it .GT. nt - 1) THEN
3669  it = nt - 1
3670  ELSE
3671  it = it
3672  END IF
3673  IF (ia .GT. na - 1) THEN
3674  ia = na - 1
3675  ELSE
3676  ia = ia
3677  END IF
3678  fm = fm - REAL(im - 1)
3679  ft = ft - REAL(it - 1)
3680  fa = fa - REAL(ia - 1)
3681 !-----scale cloud optical thickness for beam radiation following
3682 ! Eq. (7.3).
3683 ! the scaling factor, xai, is a function of the solar zenith
3684 ! angle, optical thickness, and cloud cover.
3685  CALL pushreal8(xai)
3686  xai = (-(caib(im-1, it, ia)*(1.-fm))+caib(im+1, it, ia)*(1.+fm))&
3687 & *fm*.5 + caib(im, it, ia)*(1.-fm*fm)
3688  xai = xai + (-(caib(im, it-1, ia)*(1.-ft))+caib(im, it+1, ia)*(&
3689 & 1.+ft))*ft*.5 + caib(im, it, ia)*(1.-ft*ft)
3690  xai = xai + (-(caib(im, it, ia-1)*(1.-fa))+caib(im, it, ia+1)*(&
3691 & 1.+fa))*fa*.5 + caib(im, it, ia)*(1.-fa*fa)
3692  xai = xai - 2.*caib(im, it, ia)
3693  IF (xai .LT. 0.0) THEN
3694  xai = 0.0
3695  CALL pushcontrol1b(0)
3696  ELSE
3697  CALL pushcontrol1b(1)
3698  xai = xai
3699  END IF
3700  IF (xai .GT. 1.0) THEN
3701  xai = 1.0
3702  CALL pushcontrol1b(0)
3703  ELSE
3704  CALL pushcontrol1b(1)
3705  xai = xai
3706  END IF
3707 !-----scale cloud optical thickness for diffuse radiation following
3708 ! Eq. (7.4).
3709 ! the scaling factor, xai, is a function of the cloud optical
3710 ! thickness and cover but not the solar zenith angle.
3711  CALL pushreal8(xai)
3712  xai = (-(caif(it-1, ia)*(1.-ft))+caif(it+1, ia)*(1.+ft))*ft*.5 +&
3713 & caif(it, ia)*(1.-ft*ft)
3714  xai = xai + (-(caif(it, ia-1)*(1.-fa))+caif(it, ia+1)*(1.+fa))*&
3715 & fa*.5 + caif(it, ia)*(1.-fa*fa)
3716  xai = xai - caif(it, ia)
3717  IF (xai .LT. 0.0) THEN
3718  xai = 0.0
3719  CALL pushcontrol1b(0)
3720  ELSE
3721  CALL pushcontrol1b(1)
3722  xai = xai
3723  END IF
3724  IF (xai .GT. 1.0) THEN
3725  xai = 1.0
3726  CALL pushcontrol1b(0)
3727  ELSE
3728  CALL pushcontrol1b(1)
3729  xai = xai
3730  END IF
3731  CALL pushcontrol2b(0)
3732  ELSE
3733  CALL pushcontrol2b(1)
3734  END IF
3735  ELSE
3736  CALL pushcontrol2b(2)
3737  END IF
3738 !-----compute cloud single scattering albedo and asymmetry factor
3739 ! for a mixture of ice and liquid particles.
3740 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
3741  CALL pushreal8(tauc)
3742  tauc = taucld1 + taucld2 + taucld3 + taucld4
3743  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
3744  CALL pushreal8(w1)
3745  w1 = (1.-(aia_nir(ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff(k, 1)&
3746 & )*reff(k, 1)))*taucld1
3747  CALL pushreal8(w2)
3748  w2 = (1.-(awa_nir(ib, 1)+(awa_nir(ib, 2)+awa_nir(ib, 3)*reff(k, 2)&
3749 & )*reff(k, 2)))*taucld2
3750  CALL pushreal8(w3)
3751  w3 = (1.-ara_nir(ib, 1))*taucld3
3752  CALL pushreal8(w4)
3753  w4 = (1.-(aia_nir(ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff_snow)&
3754 & *reff_snow))*taucld4
3755  g1 = (aig_nir(ib, 1)+(aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 1))*&
3756 & reff(k, 1))*w1
3757  g2 = (awg_nir(ib, 1)+(awg_nir(ib, 2)+awg_nir(ib, 3)*reff(k, 2))*&
3758 & reff(k, 2))*w2
3759  g3 = arg_nir(ib, 1)*w3
3760  g4 = (aig_nir(ib, 1)+(aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 4))*&
3761 & reff(k, 4))*w4
3762  IF (w1 + w2 + w3 + w4 .NE. 0.0) THEN
3763  CALL pushcontrol2b(0)
3764  ELSE
3765  CALL pushcontrol2b(1)
3766  END IF
3767  ELSE
3768  CALL pushcontrol2b(2)
3769  END IF
3770  END DO
3771  ccb = 0.0_8
3772  DO k=nlevs,1,-1
3773  asycltb = asyclb(k)
3774  asyclb(k) = 0.0_8
3775  ssacltb = ssaclb(k)
3776  ssaclb(k) = 0.0_8
3777  CALL popcontrol2b(branch)
3778  IF (branch .EQ. 0) THEN
3779  w1 = (1.-(aia_nir(ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff(k, 1)&
3780 & )*reff(k, 1)))*taucld1
3781  w2 = (1.-(awa_nir(ib, 1)+(awa_nir(ib, 2)+awa_nir(ib, 3)*reff(k, 2)&
3782 & )*reff(k, 2)))*taucld2
3783  taucld3 = dp(k)*1.0e3/cons_grav*hydromets(k, 3)*arb_nir(ib, 1)
3784  w3 = (1.-ara_nir(ib, 1))*taucld3
3785  w4 = (1.-(aia_nir(ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff_snow)&
3786 & *reff_snow))*taucld4
3787  g1 = (aig_nir(ib, 1)+(aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 1))*&
3788 & reff(k, 1))*w1
3789  g2 = (awg_nir(ib, 1)+(awg_nir(ib, 2)+awg_nir(ib, 3)*reff(k, 2))*&
3790 & reff(k, 2))*w2
3791  g3 = arg_nir(ib, 1)*w3
3792  g4 = (aig_nir(ib, 1)+(aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 4))*&
3793 & reff(k, 4))*w4
3794  tempb8 = asycltb/(w1+w2+w3+w4)
3795  tempb9 = -((g1+g2+g3+g4)*tempb8/(w1+w2+w3+w4))
3796  g1b = tempb8
3797  g2b = tempb8
3798  g3b = tempb8
3799  g4b = tempb8
3800  w1b = tempb9
3801  w2b = tempb9
3802  w3b = tempb9
3803  w4b = tempb9
3804  ELSE IF (branch .EQ. 1) THEN
3805  w1b = 0.0_8
3806  w2b = 0.0_8
3807  w3b = 0.0_8
3808  w4b = 0.0_8
3809  g1b = 0.0_8
3810  g2b = 0.0_8
3811  g3b = 0.0_8
3812  g4b = 0.0_8
3813  ELSE
3814  reff_snowb = 0.0_8
3815  taucld1b = 0.0_8
3816  taucld2b = 0.0_8
3817  taucld3b = 0.0_8
3818  taucld4b = 0.0_8
3819  taucb = 0.0_8
3820  GOTO 100
3821  END IF
3822  tauc = taucld1 + taucld2 + taucld3 + taucld4
3823  tempb7 = ssacltb/tauc
3824  temp6 = aig_nir(ib, 2) + aig_nir(ib, 3)*reff(k, 4)
3825  reffb(k, 4) = reffb(k, 4) + (w4*temp6+reff(k, 4)*w4*aig_nir(ib, 3))*&
3826 & g4b
3827  w4b = w4b + tempb7 + (aig_nir(ib, 1)+temp6*reff(k, 4))*g4b
3828  w3b = w3b + tempb7 + arg_nir(ib, 1)*g3b
3829  temp5 = awg_nir(ib, 2) + awg_nir(ib, 3)*reff(k, 2)
3830  reffb(k, 2) = reffb(k, 2) + (w2*temp5+reff(k, 2)*w2*awg_nir(ib, 3))*&
3831 & g2b
3832  w2b = w2b + tempb7 + (awg_nir(ib, 1)+temp5*reff(k, 2))*g2b
3833  temp4 = aig_nir(ib, 2) + aig_nir(ib, 3)*reff(k, 1)
3834  reffb(k, 1) = reffb(k, 1) + (w1*temp4+reff(k, 1)*w1*aig_nir(ib, 3))*&
3835 & g1b
3836  w1b = w1b + tempb7 + (aig_nir(ib, 1)+temp4*reff(k, 1))*g1b
3837  taucb = -((w1+w2+w3+w4)*tempb7/tauc)
3838  CALL popreal8(w4)
3839  temp3 = aia_nir(ib, 2) + aia_nir(ib, 3)*reff_snow
3840  reff_snowb = (-(taucld4*temp3)-reff_snow*taucld4*aia_nir(ib, 3))*w4b
3841  taucld4b = (1.-temp3*reff_snow-aia_nir(ib, 1))*w4b
3842  CALL popreal8(w3)
3843  taucld3b = (1.-ara_nir(ib, 1))*w3b
3844  CALL popreal8(w2)
3845  temp2 = awa_nir(ib, 2) + awa_nir(ib, 3)*reff(k, 2)
3846  reffb(k, 2) = reffb(k, 2) + (-(taucld2*temp2)-reff(k, 2)*taucld2*&
3847 & awa_nir(ib, 3))*w2b
3848  taucld2b = (1.-temp2*reff(k, 2)-awa_nir(ib, 1))*w2b
3849  CALL popreal8(w1)
3850  temp1 = aia_nir(ib, 2) + aia_nir(ib, 3)*reff(k, 1)
3851  reffb(k, 1) = reffb(k, 1) + (-(taucld1*temp1)-reff(k, 1)*taucld1*&
3852 & aia_nir(ib, 3))*w1b
3853  taucld1b = (1.-temp1*reff(k, 1)-aia_nir(ib, 1))*w1b
3854  100 CALL popreal8(tauc)
3855  taucld1b = taucld1b + taucb
3856  taucld2b = taucld2b + taucb
3857  taucld3b = taucld3b + taucb
3858  taucld4b = taucld4b + taucb
3859  CALL popcontrol2b(branch)
3860  IF (branch .EQ. 0) THEN
3861  taucld4b = taucld4b + xai*taudiffb(k, 4)
3862  xaib = taucld4*taudiffb(k, 4)
3863  taudiffb(k, 4) = 0.0_8
3864  taucld3b = taucld3b + xai*taudiffb(k, 3)
3865  xaib = xaib + taucld3*taudiffb(k, 3)
3866  taudiffb(k, 3) = 0.0_8
3867  taucld2b = taucld2b + xai*taudiffb(k, 2)
3868  xaib = xaib + taucld2*taudiffb(k, 2)
3869  taudiffb(k, 2) = 0.0_8
3870  taucld1b = taucld1b + xai*taudiffb(k, 1)
3871  xaib = xaib + taucld1*taudiffb(k, 1)
3872  taudiffb(k, 1) = 0.0_8
3873  CALL popcontrol1b(branch)
3874  IF (branch .EQ. 0) xaib = 0.0_8
3875  CALL popcontrol1b(branch)
3876  IF (branch .EQ. 0) xaib = 0.0_8
3877  tempb5 = .5*fa*xaib
3878  fab = (.5*(caif(it, ia+1)*(fa+1.)-caif(it, ia-1)*(1.-fa))-caif(it&
3879 & , ia)*2*fa)*xaib + (caif(it, ia-1)+caif(it, ia+1))*tempb5
3880  CALL popreal8(xai)
3881  tempb6 = .5*ft*xaib
3882  ftb = (.5*(caif(it+1, ia)*(ft+1.)-caif(it-1, ia)*(1.-ft))-caif(it&
3883 & , ia)*2*ft)*xaib + (caif(it-1, ia)+caif(it+1, ia))*tempb6
3884  taucld4b = taucld4b + xai*taubeamb(k, 4)
3885  xaib = taucld4*taubeamb(k, 4)
3886  taubeamb(k, 4) = 0.0_8
3887  taucld3b = taucld3b + xai*taubeamb(k, 3)
3888  xaib = xaib + taucld3*taubeamb(k, 3)
3889  taubeamb(k, 3) = 0.0_8
3890  taucld2b = taucld2b + xai*taubeamb(k, 2)
3891  xaib = xaib + taucld2*taubeamb(k, 2)
3892  taubeamb(k, 2) = 0.0_8
3893  taucld1b = taucld1b + xai*taubeamb(k, 1)
3894  xaib = xaib + taucld1*taubeamb(k, 1)
3895  taubeamb(k, 1) = 0.0_8
3896  CALL popcontrol1b(branch)
3897  IF (branch .EQ. 0) xaib = 0.0_8
3898  CALL popcontrol1b(branch)
3899  IF (branch .EQ. 0) xaib = 0.0_8
3900  tempb3 = .5*fa*xaib
3901  fab = fab + (.5*(caib(im, it, ia+1)*(fa+1.)-caib(im, it, ia-1)*(1.&
3902 & -fa))-caib(im, it, ia)*2*fa)*xaib + (caib(im, it, ia-1)+caib(im&
3903 & , it, ia+1))*tempb3
3904  tempb4 = .5*ft*xaib
3905  ftb = ftb + (.5*(caib(im, it+1, ia)*(ft+1.)-caib(im, it-1, ia)*(1.&
3906 & -ft))-caib(im, it, ia)*2*ft)*xaib + (caib(im, it-1, ia)+caib(im&
3907 & , it+1, ia))*tempb4
3908  CALL popreal8(xai)
3909  CALL popinteger4(ia)
3910  CALL popinteger4(it)
3911  CALL popinteger4(im)
3912  fab = fab/da
3913  CALL popreal8(ft)
3914  taucb = ftb/(dt*tauc*log(10.0))
3915  CALL popcontrol1b(branch)
3916  IF (branch .EQ. 0) taucb = 0.0_8
3917  CALL popcontrol1b(branch)
3918  IF (branch .EQ. 0) THEN
3919  CALL popreal8(fa)
3920  ELSE
3921  CALL popreal8(fa)
3922  tempb2 = fab/cc(kk)
3923  fcldb(k) = fcldb(k) + tempb2
3924  ccb(kk) = ccb(kk) - fcld(k)*tempb2/cc(kk)
3925  END IF
3926  ELSE IF (branch .EQ. 1) THEN
3927  taucb = 0.0_8
3928  ELSE
3929  taucld4b = taucld4b + taubeamb(k, 4) + taudiffb(k, 4)
3930  taudiffb(k, 4) = 0.0_8
3931  taubeamb(k, 4) = 0.0_8
3932  taucld3b = taucld3b + taubeamb(k, 3) + taudiffb(k, 3)
3933  taudiffb(k, 3) = 0.0_8
3934  taubeamb(k, 3) = 0.0_8
3935  taucld2b = taucld2b + taubeamb(k, 2) + taudiffb(k, 2)
3936  taudiffb(k, 2) = 0.0_8
3937  taubeamb(k, 2) = 0.0_8
3938  taucld1b = taucld1b + taubeamb(k, 1) + taudiffb(k, 1)
3939  taudiffb(k, 1) = 0.0_8
3940  taubeamb(k, 1) = 0.0_8
3941  GOTO 110
3942  END IF
3943  taucld1b = taucld1b + taucb
3944  taucld2b = taucld2b + taucb
3945  taucld3b = taucld3b + taucb
3946  taucld4b = taucld4b + taucb
3947  CALL popcontrol2b(branch)
3948  IF (branch .EQ. 0) THEN
3949  CALL popinteger4(kk)
3950  ELSE IF (branch .EQ. 1) THEN
3951  CALL popinteger4(kk)
3952  ELSE
3953  CALL popinteger4(kk)
3954  END IF
3955  110 CALL popcontrol1b(branch)
3956  IF (branch .EQ. 0) THEN
3957  CALL popreal8(taucld4)
3958  ELSE
3959  CALL popreal8(taucld4)
3960  tempb1 = dp(k)*aib_nir*1.0e3*taucld4b/(cons_grav*reff_snow)
3961  hydrometsb(k, 4) = hydrometsb(k, 4) + tempb1
3962  reff_snowb = reff_snowb - hydromets(k, 4)*tempb1/reff_snow
3963  END IF
3964  CALL popcontrol1b(branch)
3965  IF (branch .EQ. 0) THEN
3966  CALL popreal8(reff_snow)
3967  ELSE
3968  CALL popreal8(reff_snow)
3969  reffb(k, 4) = reffb(k, 4) + reff_snowb
3970  END IF
3971  CALL popreal8(taucld3)
3972  hydrometsb(k, 3) = hydrometsb(k, 3) + dp(k)*1.0e3*arb_nir(ib, 1)*&
3973 & taucld3b/cons_grav
3974  CALL popcontrol1b(branch)
3975  IF (branch .EQ. 0) THEN
3976  CALL popreal8(taucld2)
3977  temp0 = awb_nir(ib, 2)/reff(k, 2)
3978  tempb0 = dp(k)*1.0e3*taucld2b
3979  hydrometsb(k, 2) = hydrometsb(k, 2) + (awb_nir(ib, 1)+temp0)*&
3980 & tempb0/cons_grav
3981  reffb(k, 2) = reffb(k, 2) - hydromets(k, 2)*temp0*tempb0/(reff(k, &
3982 & 2)*cons_grav)
3983  ELSE
3984  CALL popreal8(taucld2)
3985  END IF
3986  CALL popcontrol1b(branch)
3987  IF (branch .EQ. 0) THEN
3988  CALL popreal8(taucld1)
3989  ELSE
3990  CALL popreal8(taucld1)
3991  temp = cons_grav*reff(k, 1)
3992  tempb = dp(k)*aib_nir*1.0e3*taucld1b/temp
3993  hydrometsb(k, 1) = hydrometsb(k, 1) + tempb
3994  reffb(k, 1) = reffb(k, 1) - hydromets(k, 1)*cons_grav*tempb/temp
3995  END IF
3996  END DO
3997  CALL popcontrol1b(branch)
3998  IF (branch .NE. 0) THEN
3999  DO k=nlevs,icb,-1
4000  CALL popcontrol1b(branch)
4001  IF (branch .EQ. 0) THEN
4002  fcldb(k) = fcldb(k) + ccb(3)
4003  ccb(3) = 0.0_8
4004  END IF
4005  END DO
4006  DO k=icb-1,ict,-1
4007  CALL popcontrol1b(branch)
4008  IF (branch .EQ. 0) THEN
4009  fcldb(k) = fcldb(k) + ccb(2)
4010  ccb(2) = 0.0_8
4011  END IF
4012  END DO
4013  DO k=ict-1,1,-1
4014  CALL popcontrol1b(branch)
4015  IF (branch .EQ. 0) THEN
4016  fcldb(k) = fcldb(k) + ccb(1)
4017  ccb(1) = 0.0_8
4018  END IF
4019  END DO
4020  END IF
4021 END SUBROUTINE getnirtau1_b
4022 
4023 end module sorad_ad
soradmod::deledd
subroutine, public deledd(tau1, ssc1, g01, cza1, rr1, tt1, td1)
Definition: sorad.F90:1495
popinteger4
subroutine popinteger4(x)
Definition: adBuffer.f:541
sorad_ad::getnirtau1_b
subroutine getnirtau1_b(ib, nlevs, cosz, dp, fcld, fcldb, reff, reffb, hydromets, hydrometsb, ict, icb, taubeam, taubeamb, taudiff, taudiffb, asycl, asyclb, ssacl, ssaclb, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, cons_grav)
Definition: sorad_ad.F90:3430
popcontrol2b
subroutine popcontrol2b(cc)
Definition: adBuffer.f:146
popreal8array
void popreal8array(double *x, int n)
Definition: adStack.c:375
pushcontrol1b
subroutine pushcontrol1b(cc)
Definition: adBuffer.f:115
sorad_ad::deledd_b
subroutine deledd_b(tau1, tau1b, ssc1, ssc1b, g01, g01b, cza1, rr1, rr1b, tt1, tt1b, td1, td1b)
Definition: sorad_ad.F90:2657
soradmod::getnirtau1
subroutine, public getnirtau1(ib, nlevs, cosz, dp, fcld, reff, hydromets, ict, icb, taubeam, taudiff, asycl, ssacl, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, CONS_GRAV)
Definition: sorad.F90:1908
pushcontrol2b
subroutine pushcontrol2b(cc)
Definition: adBuffer.f:140
soradmod::getvistau1
subroutine, public getvistau1(nlevs, cosz, dp, fcld, reff, hydromets, ict, icb, taubeam, taudiff, asycl, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, CONS_GRAV)
Definition: sorad.F90:1631
pushreal8array
void pushreal8array(double *x, int n)
Definition: adStack.c:372
popreal8
subroutine popreal8(x)
Definition: adBuffer.f:820
qg_constants::u1
real(kind=kind_real), parameter u1
Definition: qg_constants.f90:35
sorad_ad::sorad_b
subroutine, public sorad_b(m, np, nb, cosz_dev, pl_dev, ta_dev, ta_devb, wa_dev, wa_devb, oa_dev, oa_devb, co2, cwc_dev, cwc_devb, fcld_dev, fcld_devb, ict, icb, reff_dev, reff_devb, hk_uv, hk_ir, taua_dev, taua_devb, ssaa_dev, ssaa_devb, asya_dev, asya_devb, rsuvbm_dev, rsuvdf_dev, rsirbm_dev, rsirdf_dev, flx_dev, flx_devb, cons_grav, wk_uv, zk_uv, ry_uv, xk_ir, ry_ir, cah, coa, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, arg_nir, caib, caif)
Definition: sorad_ad.F90:20
popcontrol3b
subroutine popcontrol3b(cc)
Definition: adBuffer.f:175
pushreal8
subroutine pushreal8(x)
Definition: adBuffer.f:763
sorad_ad::getvistau1_b
subroutine getvistau1_b(nlevs, cosz, dp, fcld, fcldb, reff, reffb, hydromets, hydrometsb, ict, icb, taubeam, taubeamb, taudiff, taudiffb, asycl, asyclb, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, cons_grav)
Definition: sorad_ad.F90:2887
popcontrol1b
subroutine popcontrol1b(cc)
Definition: adBuffer.f:120
max
#define max(a, b)
Definition: mosaic_util.h:33
pushcontrol3b
subroutine pushcontrol3b(cc)
Definition: adBuffer.f:168
min
#define min(a, b)
Definition: mosaic_util.h:32
pushinteger4
subroutine pushinteger4(x)
Definition: adBuffer.f:484
sorad_ad
Definition: sorad_ad.F90:1