FV3 Bundle
sorad_tl.F90
Go to the documentation of this file.
1 module sorad_tl
2 
3 use soradmod
4 
5 IMPLICIT NONE
6 
7 PRIVATE
8 PUBLIC :: sorad_d
9 
10 contains
11 
12 SUBROUTINE sorad_d(m, np, nb, cosz_dev, pl_dev, ta_dev, ta_devd, wa_dev&
13 & , wa_devd, oa_dev, oa_devd, co2, cwc_dev, cwc_devd, fcld_dev, &
14 & fcld_devd, ict, icb, reff_dev, reff_devd, hk_uv, hk_ir, taua_dev, &
15 & taua_devd, ssaa_dev, ssaa_devd, asya_dev, asya_devd, rsuvbm_dev, &
16 & rsuvdf_dev, rsirbm_dev, rsirdf_dev, flx_dev, flx_devd, 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_devd(m, np), wa_devd(m, np), oa_devd(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_devd(m, np, 4), fcld_devd(m, np), reff_devd(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_devd(m, np, nb)
45  REAL*8 :: ssaa_dev(m, np, nb)
46  REAL*8 :: ssaa_devd(m, np, nb)
47  REAL*8 :: asya_dev(m, np, nb)
48  REAL*8 :: asya_devd(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_devd(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 :: dpd(np), whd(np), ohd(np)
73  REAL*8 :: scal(np)
74  REAL*8 :: scald(np)
75  REAL*8 :: swh(np+1), so2(np+1), df(0:np+1)
76  REAL*8 :: swhd(np+1), so2d(np+1), dfd(0:np+1)
77  REAL*8 :: scal0, wvtoa, o3toa, pa
78  REAL*8 :: wvtoad, o3toad
79  REAL*8 :: snt, cnt, x, xx4, xtoa
80  REAL*8 :: xx4d
81  REAL*8 :: dp_pa(np)
82  REAL*8 :: dp_pad(np)
83 !-----parameters for co2 transmission tables
84  REAL*8 :: w1, dw, u1, du
85  INTEGER :: ib, rc
86  REAL*8 :: tauclb(np), tauclf(np), asycl(np)
87  REAL*8 :: tauclbd(np), tauclfd(np), asycld(np)
88  REAL*8 :: taubeam(np, 4), taudiff(np, 4)
89  REAL*8 :: taubeamd(np, 4), taudiffd(np, 4)
90  REAL*8 :: fcld_col(np)
91  REAL*8 :: fcld_cold(np)
92  REAL*8 :: cwc_col(np, 4)
93  REAL*8 :: cwc_cold(np, 4)
94  REAL*8 :: reff_col(np, 4)
95  REAL*8 :: reff_cold(np, 4)
96  REAL*8 :: taurs, tauoz, tauwv
97  REAL*8 :: tauozd, tauwvd
98  REAL*8 :: tausto, ssatau, asysto
99  REAL*8 :: taustod, ssataud, asystod
100  REAL*8 :: tautob, ssatob, asytob
101  REAL*8 :: tautobd, ssatobd, asytobd
102  REAL*8 :: tautof, ssatof, asytof
103  REAL*8 :: tautofd, ssatofd, asytofd
104  REAL*8 :: rr(0:np+1, 2), tt(0:np+1, 2), td(0:np+1, 2)
105  REAL*8 :: rrd(0:np+1, 2), ttd(0:np+1, 2), tdd(0:np+1, 2)
106  REAL*8 :: rs(0:np+1, 2), ts(0:np+1, 2)
107  REAL*8 :: rsd(0:np+1, 2), tsd(0:np+1, 2)
108  REAL*8 :: fall(np+1), fclr(np+1), fsdir, fsdif
109  REAL*8 :: falld(np+1)
110  REAL*8 :: fupa(np+1), fupc(np+1)
111  REAL*8 :: cc1, cc2, cc3
112  REAL*8 :: cc1d, cc2d, cc3d
113  REAL*8 :: rrt, ttt, tdt, rst, tst
114  REAL*8 :: rrtd, tttd, tdtd, rstd, tstd
115  INTEGER :: iv, ik
116  REAL*8 :: ssacl(np)
117  REAL*8 :: ssacld(np)
118  INTEGER :: im
119  INTEGER :: ic, iw
120  REAL*8 :: ulog, wlog, dc, dd, x0, x1, x2, y0, y1, y2, du2, dw2
121  REAL*8 :: wlogd, ddd, x2d, y2d
122  INTEGER :: ih
123 !if (overcast == true) then
124 !real(8) :: rra(0:np+1),rxa(0:np+1)
125 !real(8) :: ttaold,tdaold,rsaold
126 !real(8) :: ttanew,tdanew,rsanew
127 !else
128  REAL*8 :: rra(0:np+1, 2, 2), tta(0:np, 2, 2)
129  REAL*8 :: rrad(0:np+1, 2, 2), ttad(0:np, 2, 2)
130  REAL*8 :: tda(0:np, 2, 2)
131  REAL*8 :: tdad(0:np, 2, 2)
132  REAL*8 :: rsa(0:np, 2, 2), rxa(0:np+1, 2, 2)
133  REAL*8 :: rsad(0:np, 2, 2), rxad(0:np+1, 2, 2)
134 !endif
135  REAL*8 :: flxdn
136  REAL*8 :: flxdnd
137  REAL*8 :: fdndir, fdndif, fupdif
138  REAL*8 :: fdndird, fdndifd, fupdifd
139  REAL*8 :: denm, yy
140  REAL*8 :: denmd, yyd
141  INTEGER :: is
142  REAL*8 :: ch, cm, ct
143  REAL*8 :: chd, cmd, ctd
144  INTEGER :: foundtop
145  REAL*8 :: dftop
146  REAL*8 :: dftopd
147 !-----Variables for aerosols
148  INTEGER :: ii, jj, irhp1, an
149  REAL*8 :: dum
150  REAL*8 :: dumd
151  INTRINSIC max
152  INTRINSIC exp
153  INTRINSIC min
154  INTRINSIC sqrt
155  INTRINSIC real
156  INTRINSIC log10
157  INTRINSIC int
158  INTRINSIC abs
159  INTRINSIC epsilon
160  REAL*8 :: arg1
161  REAL*8 :: arg1d
162  REAL*8 :: result1
163  REAL :: result10
164  REAL*8 :: x6
165  REAL*8 :: x5
166  REAL*8 :: x4
167  REAL*8 :: x3
168  REAL*8 :: x4d
169  REAL*8 :: abs0
170  i = 1
171 !RUN_LOOP: do i=1,m
172  ntop = 0
173  fdndir = 0.0
174  fdndif = 0.0
175 !-----Beginning of sorad code
176 !-----wvtoa and o3toa are the water vapor and o3 amounts of the region
177 ! above the pl(1) level.
178 ! snt is the secant of the solar zenith angle
179  snt = 1.0/cosz_dev(i)
180  IF (pl_dev(i, 1) .LT. 1.e-3) THEN
181  xtoa = 1.e-3
182  ELSE
183  xtoa = pl_dev(i, 1)
184  END IF
185  scal0 = xtoa*(0.5*xtoa/300.)**.8
186  o3toad = 1.02*xtoa*466.7*oa_devd(i, 1)
187  o3toa = 1.02*oa_dev(i, 1)*xtoa*466.7 + 1.0e-8
188  wvtoad = 1.02*scal0*(wa_devd(i, 1)*(1.0+0.00135*(ta_dev(i, 1)-240.))+&
189 & wa_dev(i, 1)*0.00135*ta_devd(i, 1))
190  wvtoa = 1.02*wa_dev(i, 1)*scal0*(1.0+0.00135*(ta_dev(i, 1)-240.)) + &
191 & 1.0e-9
192  swhd = 0.0_8
193  swhd(1) = wvtoad
194  swh(1) = wvtoa
195  reff_cold = 0.0_8
196  ohd = 0.0_8
197  fcld_cold = 0.0_8
198  cwc_cold = 0.0_8
199  whd = 0.0_8
200  DO k=1,np
201 !-----compute layer thickness. indices for the surface level and
202 ! surface layer are np+1 and np, respectively.
203  dpd(k) = 0.0_8
204  dp(k) = pl_dev(i, k+1) - pl_dev(i, k)
205 ! dp in pascals
206  dp_pad(k) = 0.0_8
207  dp_pa(k) = dp(k)*100.
208 !-----compute scaled water vapor amount following Eqs. (3.3) and (3.5)
209 ! unit is g/cm**2
210 !
211  pa = 0.5*(pl_dev(i, k)+pl_dev(i, k+1))
212  scald(k) = 0.0_8
213  scal(k) = dp(k)*(pa/300.)**.8
214  whd(k) = 1.02*scal(k)*(wa_devd(i, k)*(1.+0.00135*(ta_dev(i, k)-240.)&
215 & )+wa_dev(i, k)*0.00135*ta_devd(i, k))
216  wh(k) = 1.02*wa_dev(i, k)*scal(k)*(1.+0.00135*(ta_dev(i, k)-240.)) +&
217 & 1.e-9
218  swhd(k+1) = swhd(k) + whd(k)
219  swh(k+1) = swh(k) + wh(k)
220 !-----compute ozone amount, unit is (cm-atm)stp
221 ! the number 466.7 is the unit conversion factor
222 ! from g/cm**2 to (cm-atm)stp
223  ohd(k) = 1.02*dp(k)*466.7*oa_devd(i, k)
224  oh(k) = 1.02*oa_dev(i, k)*dp(k)*466.7 + 1.e-8
225 !-----Fill the reff, cwc, and fcld for the column
226  fcld_cold(k) = fcld_devd(i, k)
227  fcld_col(k) = fcld_dev(i, k)
228  DO l=1,4
229  reff_cold(k, l) = reff_devd(i, k, l)
230  reff_col(k, l) = reff_dev(i, k, l)
231  cwc_cold(k, l) = cwc_devd(i, k, l)
232  cwc_col(k, l) = cwc_dev(i, k, l)
233  END DO
234  END DO
235 !-----Initialize temporary arrays to zero to avoid UMR
236  rr = 0.0
237  tt = 0.0
238  td = 0.0
239  rs = 0.0
240  ts = 0.0
241  rra = 0.0
242  rxa = 0.0
243 !if( OVERCAST == .false. ) then
244  tta = 0.0
245  tda = 0.0
246  rsa = 0.0
247 !endif
248 !-----initialize fluxes for all-sky (flx), clear-sky (flc), and
249 ! flux reduction (df)
250 !
251  DO k=1,np+1
252  flx_devd(i, k) = 0.0_8
253  flx_dev(i, k) = 0.
254  flc_dev(i, k) = 0.
255  flxu_dev(i, k) = 0.
256  flcu_dev(i, k) = 0.
257  END DO
258 !-----Initialize new per-band surface fluxes
259  DO ib=1,nband
260  flx_sfc_band_dev(i, ib) = 0.
261  END DO
262 !-----Begin inline of SOLUV
263 !-----compute solar uv and par fluxes
264 !-----initialize fdiruv, fdifuv, surface reflectances and transmittances.
265 ! the reflectance and transmittance of the clear and cloudy portions
266 ! of a layer are denoted by 1 and 2, respectively.
267 ! cc is the maximum cloud cover in each of the high, middle, and low
268 ! cloud groups.
269 ! 1/dsm=1/cos(53) = 1.66
270  fdiruv_dev(i) = 0.0
271  fdifuv_dev(i) = 0.0
272  rrd(np+1, 1) = 0.0_8
273  rr(np+1, 1) = rsuvbm_dev
274  rrd(np+1, 2) = 0.0_8
275  rr(np+1, 2) = rsuvbm_dev
276  rsd(np+1, 1) = 0.0_8
277  rs(np+1, 1) = rsuvdf_dev
278  rsd(np+1, 2) = 0.0_8
279  rs(np+1, 2) = rsuvdf_dev
280  tdd(np+1, 1) = 0.0_8
281  td(np+1, 1) = 0.0
282  tdd(np+1, 2) = 0.0_8
283  td(np+1, 2) = 0.0
284  ttd(np+1, 1) = 0.0_8
285  tt(np+1, 1) = 0.0
286  ttd(np+1, 2) = 0.0_8
287  tt(np+1, 2) = 0.0
288  tsd(np+1, 1) = 0.0_8
289  ts(np+1, 1) = 0.0
290  tsd(np+1, 2) = 0.0_8
291  ts(np+1, 2) = 0.0
292  rrd(0, 1) = 0.0_8
293  rr(0, 1) = 0.0
294  rrd(0, 2) = 0.0_8
295  rr(0, 2) = 0.0
296  rsd(0, 1) = 0.0_8
297  rs(0, 1) = 0.0
298  rsd(0, 2) = 0.0_8
299  rs(0, 2) = 0.0
300 ! td(0,1)=1.0
301 ! td(0,2)=1.0
302  ttd(0, 1) = 0.0_8
303  tt(0, 1) = 1.0
304  ttd(0, 2) = 0.0_8
305  tt(0, 2) = 1.0
306  tsd(0, 1) = 0.0_8
307  ts(0, 1) = 1.0
308  tsd(0, 2) = 0.0_8
309  ts(0, 2) = 1.0
310  cc1 = 0.0
311  cc2 = 0.0
312  cc3 = 0.0
313 !-----options for scaling cloud optical thickness
314 !if ( OVERCAST == .true. ) then
315 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10) and
316 ! compute cloud asymmetry factor
317 ! Note: the cloud optical properties are assumed to be independent
318 ! of spectral bands in the UV and PAR regions.
319 ! for a mixture of ice and liquid particles.
320 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
321 !
322 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
323 ! respectively.
324 ! call getvistau1(np,cosz_dev(i),dp_pa,fcld_col,reff_col,cwc_col,0,0,&
325 ! taubeam,taudiff,asycl, &
326 ! aig_uv, awg_uv, arg_uv, &
327 ! aib_uv, awb_uv, arb_uv, &
328 ! aib_nir, awb_nir, arb_nir, &
329 ! aia_nir, awa_nir, ara_nir, &
330 ! aig_nir, awg_nir, arg_nir, &
331 ! caib, caif, &
332 ! CONS_GRAV )
333 !else
334 !-----Compute scaled cloud optical thickness. Eqs. (4.6) and (4.10) and
335 ! compute cloud asymmetry factor
336 ! Note: the cloud optical properties are assumed to be independent
337 ! of spectral bands in the UV and PAR regions.
338 ! for a mixture of ice and liquid particles.
339 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
340 !
341 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
342 ! respectively.
343  CALL getvistau1_d(np, cosz_dev(i), dp_pa, fcld_col, fcld_cold, &
344 & reff_col, reff_cold, cwc_col, cwc_cold, ict, icb, taubeam&
345 & , taubeamd, taudiff, taudiffd, asycl, asycld, aig_uv, &
346 & awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, aib_nir, awb_nir, &
347 & arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir, &
348 & arg_nir, caib, caif, cons_grav)
349  cc1d = 0.0_8
350  cc2d = 0.0_8
351  cc3d = 0.0_8
352 !-----clouds within each of the high, middle, and low clouds are
353 ! assumed to be maximally overlapped, and the cloud cover (cc)
354 ! for a group (high, middle, or low) is the maximum cloud cover
355 ! of all the layers within a group
356 ! The cc1,2,3 are still needed in the flux calculations below
357 !MAT---DO NOT FUSE THIS LOOP
358 !MAT---Loop must run to completion so that cc[1,2,3] are correct.
359  DO k=1,np
360  IF (k .LT. ict) THEN
361  IF (cc1 .LT. fcld_dev(i, k)) THEN
362  cc1d = fcld_devd(i, k)
363  cc1 = fcld_dev(i, k)
364  ELSE
365  cc1 = cc1
366  END IF
367  ELSE IF (k .LT. icb) THEN
368  IF (cc2 .LT. fcld_dev(i, k)) THEN
369  cc2d = fcld_devd(i, k)
370  cc2 = fcld_dev(i, k)
371  ELSE
372  cc2 = cc2
373  END IF
374  ELSE IF (cc3 .LT. fcld_dev(i, k)) THEN
375  cc3d = fcld_devd(i, k)
376  cc3 = fcld_dev(i, k)
377  ELSE
378  cc3 = cc3
379  END IF
380  END DO
381  tauclbd = 0.0_8
382  tauclfd = 0.0_8
383 !MAT---DO NOT FUSE THIS LOOP
384 !endif !overcast
385  DO k=1,np
386  tauclbd(k) = taubeamd(k, 1) + taubeamd(k, 2) + taubeamd(k, 3) + &
387 & taubeamd(k, 4)
388  tauclb(k) = taubeam(k, 1) + taubeam(k, 2) + taubeam(k, 3) + taubeam(&
389 & k, 4)
390  tauclfd(k) = taudiffd(k, 1) + taudiffd(k, 2) + taudiffd(k, 3) + &
391 & taudiffd(k, 4)
392  tauclf(k) = taudiff(k, 1) + taudiff(k, 2) + taudiff(k, 3) + taudiff(&
393 & k, 4)
394  END DO
395  tsd = 0.0_8
396  ttd = 0.0_8
397  rsad = 0.0_8
398  rrd = 0.0_8
399  rsd = 0.0_8
400  falld = 0.0_8
401  ttad = 0.0_8
402  rxad = 0.0_8
403  tdad = 0.0_8
404  rrad = 0.0_8
405  tdd = 0.0_8
406 !-----integration over spectral bands
407 !-----Compute optical thickness, single-scattering albedo and asymmetry
408 ! factor for a mixture of "na" aerosol types. [Eqs. (4.16)-(4.18)]
409  DO ib=1,nband_uv
410 !-----compute direct beam transmittances of the layer above pl(1)
411  arg1d = -((wk_uv(ib)*wvtoad+zk_uv(ib)*o3toad)/cosz_dev(i))
412  arg1 = -((wvtoa*wk_uv(ib)+o3toa*zk_uv(ib))/cosz_dev(i))
413  tdd(0, 1) = arg1d*exp(arg1)
414  td(0, 1) = exp(arg1)
415  tdd(0, 2) = tdd(0, 1)
416  td(0, 2) = td(0, 1)
417  DO k=1,np
418 !-----compute clear-sky optical thickness, single scattering albedo,
419 ! and asymmetry factor (Eqs. 6.2-6.4)
420  taurs = ry_uv(ib)*dp(k)
421  tauozd = zk_uv(ib)*ohd(k)
422  tauoz = zk_uv(ib)*oh(k)
423  tauwvd = wk_uv(ib)*whd(k)
424  tauwv = wk_uv(ib)*wh(k)
425  taustod = tauozd + tauwvd + taua_devd(i, k, ib)
426  tausto = taurs + tauoz + tauwv + taua_dev(i, k, ib) + 1.0e-7
427  ssataud = ssaa_devd(i, k, ib)
428  ssatau = ssaa_dev(i, k, ib) + taurs
429  asystod = asya_devd(i, k, ib)
430  asysto = asya_dev(i, k, ib)
431  tautobd = taustod
432  tautob = tausto
433  asytobd = (asystod*ssatau-asysto*ssataud)/ssatau**2
434  asytob = asysto/ssatau
435  ssatobd = (ssataud*tautob-ssatau*tautobd)/tautob**2
436  ssatob = ssatau/tautob + 1.0e-8
437  IF (ssatob .GT. 0.999999) THEN
438  ssatob = 0.999999
439  ssatobd = 0.0_8
440  ELSE
441  ssatob = ssatob
442  END IF
443 !-----for direct incident radiation
444  CALL deledd_d(tautob, tautobd, ssatob, ssatobd, asytob, asytobd, &
445 & cosz_dev(i), rrt, rrtd, ttt, tttd, tdt, tdtd)
446 !-----diffuse incident radiation is approximated by beam radiation with
447 ! an incident angle of 53 degrees, Eqs. (6.5) and (6.6)
448  CALL deledd_d(tautob, tautobd, ssatob, ssatobd, asytob, asytobd, &
449 & dsm, rst, rstd, tst, tstd, dum, dumd)
450  rrd(k, 1) = rrtd
451  rr(k, 1) = rrt
452  ttd(k, 1) = tttd
453  tt(k, 1) = ttt
454  tdd(k, 1) = tdtd
455  td(k, 1) = tdt
456  rsd(k, 1) = rstd
457  rs(k, 1) = rst
458  tsd(k, 1) = tstd
459  ts(k, 1) = tst
460 !-----compute reflectance and transmittance of the cloudy portion
461 ! of a layer
462 !-----for direct incident radiation
463 ! The effective layer optical properties. Eqs. (6.2)-(6.4)
464  tautobd = taustod + tauclbd(k)
465  tautob = tausto + tauclb(k)
466  ssatobd = ((ssataud+tauclbd(k))*tautob-(ssatau+tauclb(k))*tautobd)&
467 & /tautob**2
468  ssatob = (ssatau+tauclb(k))/tautob + 1.0e-8
469  IF (ssatob .GT. 0.999999) THEN
470  ssatob = 0.999999
471  ssatobd = 0.0_8
472  ELSE
473  ssatob = ssatob
474  END IF
475  asytobd = ((asystod+asycld(k)*tauclb(k)+asycl(k)*tauclbd(k))*&
476 & ssatob*tautob-(asysto+asycl(k)*tauclb(k))*(ssatobd*tautob+ssatob&
477 & *tautobd))/(ssatob*tautob)**2
478  asytob = (asysto+asycl(k)*tauclb(k))/(ssatob*tautob)
479 !-----for diffuse incident radiation
480  tautofd = taustod + tauclfd(k)
481  tautof = tausto + tauclf(k)
482  ssatofd = ((ssataud+tauclfd(k))*tautof-(ssatau+tauclf(k))*tautofd)&
483 & /tautof**2
484  ssatof = (ssatau+tauclf(k))/tautof + 1.0e-8
485  IF (ssatof .GT. 0.999999) THEN
486  ssatof = 0.999999
487  ssatofd = 0.0_8
488  ELSE
489  ssatof = ssatof
490  END IF
491  asytofd = ((asystod+asycld(k)*tauclf(k)+asycl(k)*tauclfd(k))*&
492 & ssatof*tautof-(asysto+asycl(k)*tauclf(k))*(ssatofd*tautof+ssatof&
493 & *tautofd))/(ssatof*tautof)**2
494  asytof = (asysto+asycl(k)*tauclf(k))/(ssatof*tautof)
495 !-----for direct incident radiation
496 ! note that the cloud optical thickness is scaled differently
497 ! for direct and diffuse insolation, Eqs. (7.3) and (7.4).
498  CALL deledd_d(tautob, tautobd, ssatob, ssatobd, asytob, asytobd, &
499 & cosz_dev(i), rrt, rrtd, ttt, tttd, tdt, tdtd)
500 !-----diffuse incident radiation is approximated by beam radiation
501 ! with an incident angle of 53 degrees, Eqs. (6.5) and (6.6)
502  CALL deledd_d(tautof, tautofd, ssatof, ssatofd, asytof, asytofd, &
503 & dsm, rst, rstd, tst, tstd, dum, dumd)
504  rrd(k, 2) = rrtd
505  rr(k, 2) = rrt
506  ttd(k, 2) = tttd
507  tt(k, 2) = ttt
508  tdd(k, 2) = tdtd
509  td(k, 2) = tdt
510  rsd(k, 2) = rstd
511  rs(k, 2) = rst
512  tsd(k, 2) = tstd
513  ts(k, 2) = tst
514  END DO
515 !-----flux calculations
516 ! initialize clear-sky flux (fclr), all-sky flux (fall),
517 ! and surface downward fluxes (fsdir and fsdif)
518  DO k=1,np+1
519  fclr(k) = 0.0
520  falld(k) = 0.0_8
521  fall(k) = 0.0
522  fupa(k) = 0.0
523  fupc(k) = 0.0
524  END DO
525  fsdir = 0.0
526  fsdif = 0.0
527 !if ( OVERCAST == .true. ) then
528 !-----Inline CLDFLXY
529 !-----for clear- and all-sky flux calculations when fractional
530 ! cloud cover is either 0 or 1.
531 !-----layers are added one at a time, going up from the surface.
532 ! rra is the composite reflectance illuminated by beam radiation
533 ! rxa is the composite reflectance illuminated from above
534 ! by diffuse radiation
535 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
536 !-----compute transmittances and reflectances for a composite of
537 ! layers. layers are added one at a time, going down from the top.
538 ! tda is the composite direct transmittance illuminated by
539 ! beam radiation
540 ! tta is the composite total transmittance illuminated by
541 ! beam radiation
542 ! rsa is the composite reflectance illuminated from below
543 ! by diffuse radiation
544 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
545 !-----compute fluxes following Eq. (6.15) for fupdif and
546 ! Eq. (6.16) for (fdndir+fdndif)
547 ! fdndir is the direct downward flux
548 ! fdndif is the diffuse downward flux
549 ! fupdif is the diffuse upward flux
550 !-----ih=1 for clear sky; ih=2 for cloudy sky.
551 !-----First set is ih = 1
552 ! rra(np+1)=rr(np+1,1)
553 ! rxa(np+1)=rs(np+1,1)
554 !
555 ! do k=np,0,-1
556 ! denm=ts(k,1)/(1.-rs(k,1)*rxa(k+1))
557 ! rra(k)=rr(k,1)+(td(k,1)*rra(k+1)+(tt(k,1)-td(k,1))*rxa(k+1))*denm
558 ! rxa(k)=rs(k,1)+ts(k,1)*rxa(k+1)*denm
559 ! end do
560 !
561 ! do k=1,np+1
562 ! if (k <= np) then
563 ! if (k == 1) then
564 ! tdaold = td(0,1)
565 ! ttaold = tt(0,1)
566 ! rsaold = rs(0,1)
567 !
568 ! tdanew = 0.0
569 ! ttanew = 0.0
570 ! rsanew = 0.0
571 ! end if
572 ! denm=ts(k,1)/(1.-rsaold*rs(k,1))
573 ! tdanew=tdaold*td(k,1)
574 ! ttanew=tdaold*tt(k,1)+(tdaold*rsaold*rr(k,1)+ttaold-tdaold)*denm
575 ! rsanew=rs(k,1)+ts(k,1)*rsaold*denm
576 ! end if
577 !
578 ! denm=1./(1.-rsaold*rxa(k))
579 ! fdndir=tdaold
580 ! xx4=tdaold*rra(k)
581 ! yy=ttaold-tdaold
582 ! fdndif=(xx4*rsaold+yy)*denm
583 ! fupdif=(xx4+yy*rxa(k))*denm
584 ! flxdn=fdndir+fdndif-fupdif
585 ! fupc(k)=fupdif
586 ! fclr(k)=flxdn
587 !
588 ! tdaold = tdanew
589 ! ttaold = ttanew
590 ! rsaold = rsanew
591 !
592 ! tdanew = 0.0
593 ! ttanew = 0.0
594 ! rsanew = 0.0
595 ! end do
596 !
597 !!-----Second set is ih = 2
598 !
599 ! rra(np+1)=rr(np+1,2)
600 ! rxa(np+1)=rs(np+1,2)
601 !
602 ! do k=np,0,-1
603 ! denm=ts(k,2)/(1.-rs(k,2)*rxa(k+1))
604 ! rra(k)=rr(k,2)+(td(k,2)*rra(k+1)+(tt(k,2)-td(k,2))*rxa(k+1))*denm
605 ! rxa(k)=rs(k,2)+ts(k,2)*rxa(k+1)*denm
606 ! end do
607 !
608 ! do k=1,np+1
609 ! if (k <= np) then
610 ! if (k == 1) then
611 ! tdaold = td(0,2)
612 ! ttaold = tt(0,2)
613 ! rsaold = rs(0,2)
614 ! tdanew = 0.0
615 ! ttanew = 0.0
616 ! rsanew = 0.0
617 ! end if
618 ! denm=ts(k,2)/(1.-rsaold*rs(k,2))
619 ! tdanew=tdaold*td(k,2)
620 ! ttanew=tdaold*tt(k,2)+(tdaold*rsaold*rr(k,2)+ttaold-tdaold)*denm
621 ! rsanew=rs(k,2)+ts(k,2)*rsaold*denm
622 ! end if
623 !
624 ! denm=1./(1.-rsaold*rxa(k))
625 ! fdndir=tdaold
626 ! xx4=tdaold*rra(k)
627 ! yy=ttaold-tdaold
628 ! fdndif=(xx4*rsaold+yy)*denm
629 ! fupdif=(xx4+yy*rxa(k))*denm
630 ! flxdn=fdndir+fdndif-fupdif
631 !
632 ! fupa(k)=fupdif
633 ! fall(k)=flxdn
634 !
635 ! tdaold = tdanew
636 ! ttaold = ttanew
637 ! rsaold = rsanew
638 !
639 ! tdanew = 0.0
640 ! ttanew = 0.0
641 ! rsanew = 0.0
642 ! end do
643 !
644 ! fsdir=fdndir
645 ! fsdif=fdndif
646 !
647 !!-----End CLDFLXY inline
648 !
649 !else
650 !-----for clear- and all-sky flux calculations when fractional
651 ! cloud cover is allowed to be between 0 and 1.
652 ! the all-sky flux, fall is the summation inside the brackets
653 ! of Eq. (7.11)
654 !-----Inline CLDFLX
655 !-----compute transmittances and reflectances for a composite of
656 ! layers. layers are added one at a time, going down from the top.
657 ! tda is the composite direct transmittance illuminated
658 ! by beam radiation
659 ! tta is the composite total transmittance illuminated by
660 ! beam radiation
661 ! rsa is the composite reflectance illuminated from below
662 ! by diffuse radiation
663 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
664 !-----To save memory space, tda, tta, and rsa are pre-computed
665 ! for k<icb. The dimension of these parameters is (m,np,2,2).
666 ! It would have been (m,np,2,2,2) if these parameters were
667 ! computed for all k's.
668 !-----for high clouds
669 ! ih=1 for clear-sky condition, ih=2 for cloudy-sky condition
670  DO ih=1,2
671  tdad(0, ih, 1) = tdd(0, ih)
672  tda(0, ih, 1) = td(0, ih)
673  ttad(0, ih, 1) = ttd(0, ih)
674  tta(0, ih, 1) = tt(0, ih)
675  rsad(0, ih, 1) = rsd(0, ih)
676  rsa(0, ih, 1) = rs(0, ih)
677  tdad(0, ih, 2) = tdd(0, ih)
678  tda(0, ih, 2) = td(0, ih)
679  ttad(0, ih, 2) = ttd(0, ih)
680  tta(0, ih, 2) = tt(0, ih)
681  rsad(0, ih, 2) = rsd(0, ih)
682  rsa(0, ih, 2) = rs(0, ih)
683  DO k=1,ict-1
684  denmd = (tsd(k, ih)*(1.-rsa(k-1, ih, 1)*rs(k, ih))-ts(k, ih)*(-(&
685 & rsad(k-1, ih, 1)*rs(k, ih))-rsa(k-1, ih, 1)*rsd(k, ih)))/(1.-&
686 & rsa(k-1, ih, 1)*rs(k, ih))**2
687  denm = ts(k, ih)/(1.-rsa(k-1, ih, 1)*rs(k, ih))
688  tdad(k, ih, 1) = tdad(k-1, ih, 1)*td(k, ih) + tda(k-1, ih, 1)*&
689 & tdd(k, ih)
690  tda(k, ih, 1) = tda(k-1, ih, 1)*td(k, ih)
691  ttad(k, ih, 1) = tdad(k-1, ih, 1)*tt(k, ih) + tda(k-1, ih, 1)*&
692 & ttd(k, ih) + ((tdad(k-1, ih, 1)*rr(k, ih)+tda(k-1, ih, 1)*rrd(&
693 & k, ih))*rsa(k-1, ih, 1)+tda(k-1, ih, 1)*rr(k, ih)*rsad(k-1, ih&
694 & , 1)+ttad(k-1, ih, 1)-tdad(k-1, ih, 1))*denm + (tda(k-1, ih, 1&
695 & )*rsa(k-1, ih, 1)*rr(k, ih)+tta(k-1, ih, 1)-tda(k-1, ih, 1))*&
696 & denmd
697  tta(k, ih, 1) = tda(k-1, ih, 1)*tt(k, ih) + (tda(k-1, ih, 1)*rsa&
698 & (k-1, ih, 1)*rr(k, ih)+tta(k-1, ih, 1)-tda(k-1, ih, 1))*denm
699  rsad(k, ih, 1) = rsd(k, ih) + (tsd(k, ih)*denm+ts(k, ih)*denmd)*&
700 & rsa(k-1, ih, 1) + ts(k, ih)*denm*rsad(k-1, ih, 1)
701  rsa(k, ih, 1) = rs(k, ih) + ts(k, ih)*rsa(k-1, ih, 1)*denm
702  tdad(k, ih, 2) = tdad(k, ih, 1)
703  tda(k, ih, 2) = tda(k, ih, 1)
704  ttad(k, ih, 2) = ttad(k, ih, 1)
705  tta(k, ih, 2) = tta(k, ih, 1)
706  rsad(k, ih, 2) = rsad(k, ih, 1)
707  rsa(k, ih, 2) = rsa(k, ih, 1)
708  END DO
709 ! k loop
710 !-----for middle clouds
711 ! im=1 for clear-sky condition, im=2 for cloudy-sky condition
712  DO k=ict,icb-1
713  DO im=1,2
714  denmd = (tsd(k, im)*(1.-rsa(k-1, ih, im)*rs(k, im))-ts(k, im)*&
715 & (-(rsad(k-1, ih, im)*rs(k, im))-rsa(k-1, ih, im)*rsd(k, im))&
716 & )/(1.-rsa(k-1, ih, im)*rs(k, im))**2
717  denm = ts(k, im)/(1.-rsa(k-1, ih, im)*rs(k, im))
718  tdad(k, ih, im) = tdad(k-1, ih, im)*td(k, im) + tda(k-1, ih, &
719 & im)*tdd(k, im)
720  tda(k, ih, im) = tda(k-1, ih, im)*td(k, im)
721  ttad(k, ih, im) = tdad(k-1, ih, im)*tt(k, im) + tda(k-1, ih, &
722 & im)*ttd(k, im) + ((tdad(k-1, ih, im)*rr(k, im)+tda(k-1, ih, &
723 & im)*rrd(k, im))*rsa(k-1, ih, im)+tda(k-1, ih, im)*rr(k, im)*&
724 & rsad(k-1, ih, im)+ttad(k-1, ih, im)-tdad(k-1, ih, im))*denm &
725 & + (tda(k-1, ih, im)*rsa(k-1, ih, im)*rr(k, im)+tta(k-1, ih, &
726 & im)-tda(k-1, ih, im))*denmd
727  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, im) + (tda(k-1, ih, im&
728 & )*rsa(k-1, ih, im)*rr(k, im)+tta(k-1, ih, im)-tda(k-1, ih, &
729 & im))*denm
730  rsad(k, ih, im) = rsd(k, im) + (tsd(k, im)*denm+ts(k, im)*&
731 & denmd)*rsa(k-1, ih, im) + ts(k, im)*denm*rsad(k-1, ih, im)
732  rsa(k, ih, im) = rs(k, im) + ts(k, im)*rsa(k-1, ih, im)*denm
733  END DO
734  END DO
735  END DO
736 ! im loop
737 ! k loop
738 ! ih loop
739 !-----layers are added one at a time, going up from the surface.
740 ! rra is the composite reflectance illuminated by beam radiation
741 ! rxa is the composite reflectance illuminated from above
742 ! by diffuse radiation
743 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
744 !-----To save memory space, rra and rxa are pre-computed for k>=icb.
745 ! the dimension of these parameters is (m,np,2,2). It would have
746 ! been (m,np,2,2,2) if these parameters were computed for all k's.
747 !-----for the low clouds
748 ! is=1 for clear-sky condition, is=2 for cloudy-sky condition
749  DO is=1,2
750  rrad(np+1, 1, is) = rrd(np+1, is)
751  rra(np+1, 1, is) = rr(np+1, is)
752  rxad(np+1, 1, is) = rsd(np+1, is)
753  rxa(np+1, 1, is) = rs(np+1, is)
754  rrad(np+1, 2, is) = rrd(np+1, is)
755  rra(np+1, 2, is) = rr(np+1, is)
756  rxad(np+1, 2, is) = rsd(np+1, is)
757  rxa(np+1, 2, is) = rs(np+1, is)
758  DO k=np,icb,-1
759  denmd = (tsd(k, is)*(1.-rs(k, is)*rxa(k+1, 1, is))-ts(k, is)*(-(&
760 & rsd(k, is)*rxa(k+1, 1, is))-rs(k, is)*rxad(k+1, 1, is)))/(1.-&
761 & rs(k, is)*rxa(k+1, 1, is))**2
762  denm = ts(k, is)/(1.-rs(k, is)*rxa(k+1, 1, is))
763  rrad(k, 1, is) = rrd(k, is) + (tdd(k, is)*rra(k+1, 1, is)+td(k, &
764 & is)*rrad(k+1, 1, is)+(ttd(k, is)-tdd(k, is))*rxa(k+1, 1, is)+(&
765 & tt(k, is)-td(k, is))*rxad(k+1, 1, is))*denm + (td(k, is)*rra(k&
766 & +1, 1, is)+(tt(k, is)-td(k, is))*rxa(k+1, 1, is))*denmd
767  rra(k, 1, is) = rr(k, is) + (td(k, is)*rra(k+1, 1, is)+(tt(k, is&
768 & )-td(k, is))*rxa(k+1, 1, is))*denm
769  rxad(k, 1, is) = rsd(k, is) + (tsd(k, is)*denm+ts(k, is)*denmd)*&
770 & rxa(k+1, 1, is) + ts(k, is)*denm*rxad(k+1, 1, is)
771  rxa(k, 1, is) = rs(k, is) + ts(k, is)*rxa(k+1, 1, is)*denm
772  rrad(k, 2, is) = rrad(k, 1, is)
773  rra(k, 2, is) = rra(k, 1, is)
774  rxad(k, 2, is) = rxad(k, 1, is)
775  rxa(k, 2, is) = rxa(k, 1, is)
776  END DO
777 ! k loop
778 !-----for middle clouds
779  DO k=icb-1,ict,-1
780  DO im=1,2
781  denmd = (tsd(k, im)*(1.-rs(k, im)*rxa(k+1, im, is))-ts(k, im)*&
782 & (-(rsd(k, im)*rxa(k+1, im, is))-rs(k, im)*rxad(k+1, im, is))&
783 & )/(1.-rs(k, im)*rxa(k+1, im, is))**2
784  denm = ts(k, im)/(1.-rs(k, im)*rxa(k+1, im, is))
785  rrad(k, im, is) = rrd(k, im) + (tdd(k, im)*rra(k+1, im, is)+td&
786 & (k, im)*rrad(k+1, im, is)+(ttd(k, im)-tdd(k, im))*rxa(k+1, &
787 & im, is)+(tt(k, im)-td(k, im))*rxad(k+1, im, is))*denm + (td(&
788 & k, im)*rra(k+1, im, is)+(tt(k, im)-td(k, im))*rxa(k+1, im, &
789 & is))*denmd
790  rra(k, im, is) = rr(k, im) + (td(k, im)*rra(k+1, im, is)+(tt(k&
791 & , im)-td(k, im))*rxa(k+1, im, is))*denm
792  rxad(k, im, is) = rsd(k, im) + (tsd(k, im)*denm+ts(k, im)*&
793 & denmd)*rxa(k+1, im, is) + ts(k, im)*denm*rxad(k+1, im, is)
794  rxa(k, im, is) = rs(k, im) + ts(k, im)*rxa(k+1, im, is)*denm
795  END DO
796  END DO
797  END DO
798 ! im loop
799 ! k loop
800 ! is loop
801 !-----integration over eight sky situations.
802 ! ih, im, is denote high, middle and low cloud groups.
803  DO ih=1,2
804 !-----clear portion
805  IF (ih .EQ. 1) THEN
806  chd = -cc1d
807  ch = 1.0 - cc1
808 !-----cloudy portion
809  ELSE
810  chd = cc1d
811  ch = cc1
812  END IF
813  DO im=1,2
814 !-----clear portion
815  IF (im .EQ. 1) THEN
816  cmd = chd*(1.0-cc2) - ch*cc2d
817  cm = ch*(1.0-cc2)
818 !-----cloudy portion
819  ELSE
820  cmd = chd*cc2 + ch*cc2d
821  cm = ch*cc2
822  END IF
823  DO is=1,2
824 !-----clear portion
825  IF (is .EQ. 1) THEN
826  ctd = cmd*(1.0-cc3) - cm*cc3d
827  ct = cm*(1.0-cc3)
828 !-----cloudy portion
829  ELSE
830  ctd = cmd*cc3 + cm*cc3d
831  ct = cm*cc3
832  END IF
833 !-----add one layer at a time, going down.
834  DO k=icb,np
835  denmd = (tsd(k, is)*(1.-rsa(k-1, ih, im)*rs(k, is))-ts(k, is&
836 & )*(-(rsad(k-1, ih, im)*rs(k, is))-rsa(k-1, ih, im)*rsd(k, &
837 & is)))/(1.-rsa(k-1, ih, im)*rs(k, is))**2
838  denm = ts(k, is)/(1.-rsa(k-1, ih, im)*rs(k, is))
839  tdad(k, ih, im) = tdad(k-1, ih, im)*td(k, is) + tda(k-1, ih&
840 & , im)*tdd(k, is)
841  tda(k, ih, im) = tda(k-1, ih, im)*td(k, is)
842  ttad(k, ih, im) = tdad(k-1, ih, im)*tt(k, is) + tda(k-1, ih&
843 & , im)*ttd(k, is) + ((tdad(k-1, ih, im)*rr(k, is)+tda(k-1, &
844 & ih, im)*rrd(k, is))*rsa(k-1, ih, im)+tda(k-1, ih, im)*rr(k&
845 & , is)*rsad(k-1, ih, im)+ttad(k-1, ih, im)-tdad(k-1, ih, im&
846 & ))*denm + (tda(k-1, ih, im)*rr(k, is)*rsa(k-1, ih, im)+tta&
847 & (k-1, ih, im)-tda(k-1, ih, im))*denmd
848  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, is) + (tda(k-1, ih, &
849 & im)*rr(k, is)*rsa(k-1, ih, im)+tta(k-1, ih, im)-tda(k-1, &
850 & ih, im))*denm
851  rsad(k, ih, im) = rsd(k, is) + (tsd(k, is)*denm+ts(k, is)*&
852 & denmd)*rsa(k-1, ih, im) + ts(k, is)*denm*rsad(k-1, ih, im)
853  rsa(k, ih, im) = rs(k, is) + ts(k, is)*rsa(k-1, ih, im)*denm
854  END DO
855 ! k loop
856 !-----add one layer at a time, going up.
857  DO k=ict-1,0,-1
858  denmd = (tsd(k, ih)*(1.-rs(k, ih)*rxa(k+1, im, is))-ts(k, ih&
859 & )*(-(rsd(k, ih)*rxa(k+1, im, is))-rs(k, ih)*rxad(k+1, im, &
860 & is)))/(1.-rs(k, ih)*rxa(k+1, im, is))**2
861  denm = ts(k, ih)/(1.-rs(k, ih)*rxa(k+1, im, is))
862  rrad(k, im, is) = rrd(k, ih) + (tdd(k, ih)*rra(k+1, im, is)+&
863 & td(k, ih)*rrad(k+1, im, is)+(ttd(k, ih)-tdd(k, ih))*rxa(k+&
864 & 1, im, is)+(tt(k, ih)-td(k, ih))*rxad(k+1, im, is))*denm +&
865 & (td(k, ih)*rra(k+1, im, is)+(tt(k, ih)-td(k, ih))*rxa(k+1&
866 & , im, is))*denmd
867  rra(k, im, is) = rr(k, ih) + (td(k, ih)*rra(k+1, im, is)+(tt&
868 & (k, ih)-td(k, ih))*rxa(k+1, im, is))*denm
869  rxad(k, im, is) = rsd(k, ih) + (tsd(k, ih)*denm+ts(k, ih)*&
870 & denmd)*rxa(k+1, im, is) + ts(k, ih)*denm*rxad(k+1, im, is)
871  rxa(k, im, is) = rs(k, ih) + ts(k, ih)*rxa(k+1, im, is)*denm
872  END DO
873 ! k loop
874 !-----compute fluxes following Eq. (6.15) for fupdif and
875 ! Eq. (6.16) for (fdndir+fdndif)
876 ! fdndir is the direct downward flux
877 ! fdndif is the diffuse downward flux
878 ! fupdif is the diffuse upward flux
879  DO k=1,np+1
880  denmd = -((-(rsad(k-1, ih, im)*rxa(k, im, is))-rsa(k-1, ih, &
881 & im)*rxad(k, im, is))/(1.-rsa(k-1, ih, im)*rxa(k, im, is))&
882 & **2)
883  denm = 1./(1.-rsa(k-1, ih, im)*rxa(k, im, is))
884  fdndird = tdad(k-1, ih, im)
885  fdndir = tda(k-1, ih, im)
886  xx4d = tdad(k-1, ih, im)*rra(k, im, is) + tda(k-1, ih, im)*&
887 & rrad(k, im, is)
888  xx4 = tda(k-1, ih, im)*rra(k, im, is)
889  yyd = ttad(k-1, ih, im) - tdad(k-1, ih, im)
890  yy = tta(k-1, ih, im) - tda(k-1, ih, im)
891  fdndifd = (xx4d*rsa(k-1, ih, im)+xx4*rsad(k-1, ih, im)+yyd)*&
892 & denm + (xx4*rsa(k-1, ih, im)+yy)*denmd
893  fdndif = (xx4*rsa(k-1, ih, im)+yy)*denm
894  fupdifd = (xx4d+yyd*rxa(k, im, is)+yy*rxad(k, im, is))*denm &
895 & + (xx4+yy*rxa(k, im, is))*denmd
896  fupdif = (xx4+yy*rxa(k, im, is))*denm
897  flxdnd = fdndird + fdndifd - fupdifd
898  flxdn = fdndir + fdndif - fupdif
899 !-----summation of fluxes over all sky situations;
900 ! the term in the brackets of Eq. (7.11)
901  IF (ih .EQ. 1 .AND. im .EQ. 1 .AND. is .EQ. 1) THEN
902  fupc(k) = fupdif
903  fclr(k) = flxdn
904  END IF
905  fupa(k) = fupa(k) + fupdif*ct
906  falld(k) = falld(k) + flxdnd*ct + flxdn*ctd
907  fall(k) = fall(k) + flxdn*ct
908  END DO
909 ! k loop
910  fsdir = fsdir + fdndir*ct
911  fsdif = fsdif + fdndif*ct
912  END DO
913  END DO
914  END DO
915 ! is loop
916 ! im loop
917 ! ih loop
918 !-----End CLDFLX inline
919 !endif !overcast
920 !-----flux integration, Eq. (6.1)
921  DO k=1,np+1
922  flx_devd(i, k) = flx_devd(i, k) + hk_uv(ib)*falld(k)
923  flx_dev(i, k) = flx_dev(i, k) + fall(k)*hk_uv(ib)
924  flc_dev(i, k) = flc_dev(i, k) + fclr(k)*hk_uv(ib)
925  flxu_dev(i, k) = flxu_dev(i, k) + fupa(k)*hk_uv(ib)
926  flcu_dev(i, k) = flcu_dev(i, k) + fupc(k)*hk_uv(ib)
927  END DO
928 !-----get surface flux for each band
929  flx_sfc_band_dev(i, ib) = flx_sfc_band_dev(i, ib) + fall(np+1)*hk_uv&
930 & (ib)
931 !-----compute direct and diffuse downward surface fluxes in the UV
932 ! and par regions
933  IF (ib .LT. 5) THEN
934  fdiruv_dev(i) = fdiruv_dev(i) + fsdir*hk_uv(ib)
935  fdifuv_dev(i) = fdifuv_dev(i) + fsdif*hk_uv(ib)
936  ELSE
937  fdirpar_dev(i) = fsdir*hk_uv(ib)
938  fdifpar_dev(i) = fsdif*hk_uv(ib)
939  END IF
940  END DO
941 !-----Inline SOLIR
942 !-----compute and update solar ir fluxes
943  fdirir_dev(i) = 0.0
944  fdifir_dev(i) = 0.0
945  rrd(np+1, 1) = 0.0_8
946  rr(np+1, 1) = rsirbm_dev
947  rrd(np+1, 2) = 0.0_8
948  rr(np+1, 2) = rsirbm_dev
949  rsd(np+1, 1) = 0.0_8
950  rs(np+1, 1) = rsirdf_dev
951  rsd(np+1, 2) = 0.0_8
952  rs(np+1, 2) = rsirdf_dev
953  tdd(np+1, 1) = 0.0_8
954  td(np+1, 1) = 0.0
955  tdd(np+1, 2) = 0.0_8
956  td(np+1, 2) = 0.0
957  ttd(np+1, 1) = 0.0_8
958  tt(np+1, 1) = 0.0
959  ttd(np+1, 2) = 0.0_8
960  tt(np+1, 2) = 0.0
961  tsd(np+1, 1) = 0.0_8
962  ts(np+1, 1) = 0.0
963  tsd(np+1, 2) = 0.0_8
964  ts(np+1, 2) = 0.0
965  rrd(0, 1) = 0.0_8
966  rr(0, 1) = 0.0
967  rrd(0, 2) = 0.0_8
968  rr(0, 2) = 0.0
969  rsd(0, 1) = 0.0_8
970  rs(0, 1) = 0.0
971  rsd(0, 2) = 0.0_8
972  rs(0, 2) = 0.0
973 ! td(0,1)=1.0
974 ! td(0,2)=1.0
975  ttd(0, 1) = 0.0_8
976  tt(0, 1) = 1.0
977  ttd(0, 2) = 0.0_8
978  tt(0, 2) = 1.0
979  tsd(0, 1) = 0.0_8
980  ts(0, 1) = 1.0
981  tsd(0, 2) = 0.0_8
982  ts(0, 2) = 1.0
983  cc1 = 0.0
984  cc2 = 0.0
985  cc3 = 0.0
986  cc1d = 0.0_8
987  cc2d = 0.0_8
988  cc3d = 0.0_8
989  ssacld = 0.0_8
990 !-----integration over spectral bands
991 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10)
992 ! The indices 1, 2, 3 are for ice, water, rain particles,
993 ! respectively.
994  DO ib=1,nband_ir
995  iv = ib + 5
996 !-----options for scaling cloud optical thickness
997 !if ( OVERCAST == .true. ) then
998 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10) and
999 ! compute cloud single scattering albedo and asymmetry factor
1000 ! for a mixture of ice and liquid particles.
1001 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
1002 !
1003 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
1004 ! respectively.
1005 ! call getnirtau1(ib,np,cosz_dev(i),dp_pa,fcld_col,reff_col,cwc_col,0,0,&
1006 ! taubeam,taudiff,asycl,ssacl, &
1007 ! aig_uv, awg_uv, arg_uv, &
1008 ! aib_uv, awb_uv, arb_uv, &
1009 ! aib_nir, awb_nir, arb_nir, &
1010 ! aia_nir, awa_nir, ara_nir, &
1011 ! aig_nir, awg_nir, arg_nir, &
1012 ! caib, caif, &
1013 ! CONS_GRAV )
1014 !else
1015 !-----Compute scaled cloud optical thickness. Eqs. (4.6) and (4.10) and
1016 ! compute cloud single scattering albedo and asymmetry factor
1017 ! for a mixture of ice and liquid particles.
1018 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
1019 !
1020 ! The indices 1, 2, 3, 4 are for ice, water, rain, snow particles,
1021 ! respectively.
1022  CALL getnirtau1_d(ib, np, cosz_dev(i), dp_pa, fcld_col, fcld_cold, &
1023 & reff_col, reff_cold, cwc_col, cwc_cold, ict, icb, &
1024 & taubeam, taubeamd, taudiff, taudiffd, asycl, asycld, &
1025 & ssacl, ssacld, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, &
1026 & arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, &
1027 & ara_nir, aig_nir, awg_nir, arg_nir, caib, caif, &
1028 & cons_grav)
1029 !-----clouds within each of the high, middle, and low clouds are
1030 ! assumed to be maximally overlapped, and the cloud cover (cc)
1031 ! for a group (high, middle, or low) is the maximum cloud cover
1032 ! of all the layers within a group
1033 !MAT--DO NOT FUSE THIS LOOP
1034 !MAT Loop must run to completion so that cc[1,2,3] are correct.
1035  DO k=1,np
1036  IF (k .LT. ict) THEN
1037  IF (cc1 .LT. fcld_dev(i, k)) THEN
1038  cc1d = fcld_devd(i, k)
1039  cc1 = fcld_dev(i, k)
1040  ELSE
1041  cc1 = cc1
1042  END IF
1043  ELSE IF (k .LT. icb) THEN
1044  IF (cc2 .LT. fcld_dev(i, k)) THEN
1045  cc2d = fcld_devd(i, k)
1046  cc2 = fcld_dev(i, k)
1047  ELSE
1048  cc2 = cc2
1049  END IF
1050  ELSE IF (cc3 .LT. fcld_dev(i, k)) THEN
1051  cc3d = fcld_devd(i, k)
1052  cc3 = fcld_dev(i, k)
1053  ELSE
1054  cc3 = cc3
1055  END IF
1056  END DO
1057 !MAT--DO NOT FUSE THIS LOOP
1058 !endif !overcast
1059  DO k=1,np
1060  tauclbd(k) = taubeamd(k, 1) + taubeamd(k, 2) + taubeamd(k, 3) + &
1061 & taubeamd(k, 4)
1062  tauclb(k) = taubeam(k, 1) + taubeam(k, 2) + taubeam(k, 3) + &
1063 & taubeam(k, 4)
1064  tauclfd(k) = taudiffd(k, 1) + taudiffd(k, 2) + taudiffd(k, 3) + &
1065 & taudiffd(k, 4)
1066  tauclf(k) = taudiff(k, 1) + taudiff(k, 2) + taudiff(k, 3) + &
1067 & taudiff(k, 4)
1068  END DO
1069 !-----integration over the k-distribution function
1070  DO ik=1,nk_ir
1071 !-----compute direct beam transmittances of the layer above pl(1)
1072  tdd(0, 1) = -(xk_ir(ik)*wvtoad*exp(-(wvtoa*xk_ir(ik)/cosz_dev(i)))&
1073 & /cosz_dev(i))
1074  td(0, 1) = exp(-(wvtoa*xk_ir(ik)/cosz_dev(i)))
1075  tdd(0, 2) = tdd(0, 1)
1076  td(0, 2) = td(0, 1)
1077  DO k=1,np
1078  taurs = ry_ir(ib)*dp(k)
1079  tauwvd = xk_ir(ik)*whd(k)
1080  tauwv = xk_ir(ik)*wh(k)
1081 !-----compute clear-sky optical thickness, single scattering albedo,
1082 ! and asymmetry factor. Eqs.(6.2)-(6.4)
1083  taustod = tauwvd + taua_devd(i, k, iv)
1084  tausto = taurs + tauwv + taua_dev(i, k, iv) + 1.0e-7
1085  ssataud = ssaa_devd(i, k, iv)
1086  ssatau = ssaa_dev(i, k, iv) + taurs + 1.0e-8
1087  asystod = asya_devd(i, k, iv)
1088  asysto = asya_dev(i, k, iv)
1089  tautobd = taustod
1090  tautob = tausto
1091  asytobd = (asystod*ssatau-asysto*ssataud)/ssatau**2
1092  asytob = asysto/ssatau
1093  ssatobd = (ssataud*tautob-ssatau*tautobd)/tautob**2
1094  ssatob = ssatau/tautob + 1.0e-8
1095  IF (ssatob .GT. 0.999999) THEN
1096  ssatob = 0.999999
1097  ssatobd = 0.0_8
1098  ELSE
1099  ssatob = ssatob
1100  END IF
1101 !-----Compute reflectance and transmittance of the clear portion
1102 ! of a layer
1103 !-----for direct incident radiation
1104  CALL deledd_d(tautob, tautobd, ssatob, ssatobd, asytob, asytobd&
1105 & , cosz_dev(i), rrt, rrtd, ttt, tttd, tdt, tdtd)
1106 !-----diffuse incident radiation is approximated by beam radiation with
1107 ! an incident angle of 53 degrees, Eqs. (6.5) and (6.6)
1108  CALL deledd_d(tautob, tautobd, ssatob, ssatobd, asytob, asytobd&
1109 & , dsm, rst, rstd, tst, tstd, dum, dumd)
1110  rrd(k, 1) = rrtd
1111  rr(k, 1) = rrt
1112  ttd(k, 1) = tttd
1113  tt(k, 1) = ttt
1114  tdd(k, 1) = tdtd
1115  td(k, 1) = tdt
1116  rsd(k, 1) = rstd
1117  rs(k, 1) = rst
1118  tsd(k, 1) = tstd
1119  ts(k, 1) = tst
1120 !-----compute reflectance and transmittance of the cloudy portion
1121 ! of a layer
1122 !-----for direct incident radiation. Eqs.(6.2)-(6.4)
1123  tautobd = taustod + tauclbd(k)
1124  tautob = tausto + tauclb(k)
1125  ssatobd = ((ssataud+ssacld(k)*tauclb(k)+ssacl(k)*tauclbd(k))*&
1126 & tautob-(ssatau+ssacl(k)*tauclb(k))*tautobd)/tautob**2
1127  ssatob = (ssatau+ssacl(k)*tauclb(k))/tautob + 1.0e-8
1128  IF (ssatob .GT. 0.999999) THEN
1129  ssatob = 0.999999
1130  ssatobd = 0.0_8
1131  ELSE
1132  ssatob = ssatob
1133  END IF
1134  asytobd = ((asystod+(asycld(k)*ssacl(k)+asycl(k)*ssacld(k))*&
1135 & tauclb(k)+asycl(k)*ssacl(k)*tauclbd(k))*ssatob*tautob-(asysto+&
1136 & asycl(k)*ssacl(k)*tauclb(k))*(ssatobd*tautob+ssatob*tautobd))/&
1137 & (ssatob*tautob)**2
1138  asytob = (asysto+asycl(k)*ssacl(k)*tauclb(k))/(ssatob*tautob)
1139 !-----for diffuse incident radiation
1140  tautofd = taustod + tauclfd(k)
1141  tautof = tausto + tauclf(k)
1142  ssatofd = ((ssataud+ssacld(k)*tauclf(k)+ssacl(k)*tauclfd(k))*&
1143 & tautof-(ssatau+ssacl(k)*tauclf(k))*tautofd)/tautof**2
1144  ssatof = (ssatau+ssacl(k)*tauclf(k))/tautof + 1.0e-8
1145  IF (ssatof .GT. 0.999999) THEN
1146  ssatof = 0.999999
1147  ssatofd = 0.0_8
1148  ELSE
1149  ssatof = ssatof
1150  END IF
1151  asytofd = ((asystod+(asycld(k)*ssacl(k)+asycl(k)*ssacld(k))*&
1152 & tauclf(k)+asycl(k)*ssacl(k)*tauclfd(k))*ssatof*tautof-(asysto+&
1153 & asycl(k)*ssacl(k)*tauclf(k))*(ssatofd*tautof+ssatof*tautofd))/&
1154 & (ssatof*tautof)**2
1155  asytof = (asysto+asycl(k)*ssacl(k)*tauclf(k))/(ssatof*tautof)
1156 !-----for direct incident radiation
1157  CALL deledd_d(tautob, tautobd, ssatob, ssatobd, asytob, asytobd&
1158 & , cosz_dev(i), rrt, rrtd, ttt, tttd, tdt, tdtd)
1159 !-----diffuse incident radiation is approximated by beam radiation with
1160 ! an incident angle of 53 degrees, Eqs.(6.5) and (6.6)
1161  CALL deledd_d(tautof, tautofd, ssatof, ssatofd, asytof, asytofd&
1162 & , dsm, rst, rstd, tst, tstd, dum, dumd)
1163  rrd(k, 2) = rrtd
1164  rr(k, 2) = rrt
1165  ttd(k, 2) = tttd
1166  tt(k, 2) = ttt
1167  tdd(k, 2) = tdtd
1168  td(k, 2) = tdt
1169  rsd(k, 2) = rstd
1170  rs(k, 2) = rst
1171  tsd(k, 2) = tstd
1172  ts(k, 2) = tst
1173  END DO
1174 !-----FLUX CALCULATIONS
1175 ! initialize clear-sky flux (fclr), all-sky flux (fall),
1176 ! and surface downward fluxes (fsdir and fsdif)
1177  DO k=1,np+1
1178  fclr(k) = 0.0
1179  falld(k) = 0.0_8
1180  fall(k) = 0.0
1181  fupc(k) = 0.0
1182  fupa(k) = 0.0
1183  END DO
1184  fsdir = 0.0
1185  fsdif = 0.0
1186 !-----for clear- and all-sky flux calculations when fractional
1187 ! cloud cover is either 0 or 1.
1188 !if ( OVERCAST == .true. ) then
1189 !-----Inline CLDFLXY
1190 !-----layers are added one at a time, going up from the surface.
1191 ! rra is the composite reflectance illuminated by beam radiation
1192 ! rxa is the composite reflectance illuminated from above
1193 ! by diffuse radiation
1194 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
1195 !-----compute transmittances and reflectances for a composite of
1196 ! layers. layers are added one at a time, going down from the top.
1197 ! tda is the composite direct transmittance illuminated by
1198 ! beam radiation
1199 ! tta is the composite total transmittance illuminated by
1200 ! beam radiation
1201 ! rsa is the composite reflectance illuminated from below
1202 ! by diffuse radiation
1203 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
1204 !-----compute fluxes following Eq. (6.15) for fupdif and
1205 ! Eq. (6.16) for (fdndir+fdndif)
1206 ! fdndir is the direct downward flux
1207 ! fdndif is the diffuse downward flux
1208 ! fupdif is the diffuse upward flux
1209 !-----ih=1 for clear sky; ih=2 for cloudy sky.
1210 !-----First set is ih = 1
1211 ! rra(np+1)=rr(np+1,1)
1212 ! rxa(np+1)=rs(np+1,1)
1213 !
1214 ! do k=np,0,-1
1215 ! denm=ts(k,1)/(1.-rs(k,1)*rxa(k+1))
1216 ! rra(k)=rr(k,1)+(td(k,1)*rra(k+1)+(tt(k,1)-td(k,1))*rxa(k+1))*denm
1217 ! rxa(k)=rs(k,1)+ts(k,1)*rxa(k+1)*denm
1218 ! end do
1219 !
1220 ! do k=1,np+1
1221 ! if (k <= np) then
1222 ! if (k == 1) then
1223 ! tdaold = td(0,1)
1224 ! ttaold = tt(0,1)
1225 ! rsaold = rs(0,1)
1226 !
1227 ! tdanew = 0.0
1228 ! ttanew = 0.0
1229 ! rsanew = 0.0
1230 ! end if
1231 ! denm=ts(k,1)/(1.-rsaold*rs(k,1))
1232 ! tdanew=tdaold*td(k,1)
1233 ! ttanew=tdaold*tt(k,1)+(tdaold*rsaold*rr(k,1)+ttaold-tdaold)*denm
1234 ! rsanew=rs(k,1)+ts(k,1)*rsaold*denm
1235 ! end if
1236 !
1237 ! denm=1./(1.-rsaold*rxa(k))
1238 ! fdndir=tdaold
1239 ! xx4=tdaold*rra(k)
1240 ! yy=ttaold-tdaold
1241 ! fdndif=(xx4*rsaold+yy)*denm
1242 ! fupdif=(xx4+yy*rxa(k))*denm
1243 ! flxdn=fdndir+fdndif-fupdif
1244 !
1245 ! fupc(k)=fupdif
1246 ! fclr(k)=flxdn
1247 !
1248 ! tdaold = tdanew
1249 ! ttaold = ttanew
1250 ! rsaold = rsanew
1251 !
1252 ! tdanew = 0.0
1253 ! ttanew = 0.0
1254 ! rsanew = 0.0
1255 ! end do
1256 !
1257 !!-----Second set is ih = 2
1258 !
1259 ! rra(np+1)=rr(np+1,2)
1260 ! rxa(np+1)=rs(np+1,2)
1261 !
1262 ! do k=np,0,-1
1263 ! denm=ts(k,2)/(1.-rs(k,2)*rxa(k+1))
1264 ! rra(k)=rr(k,2)+(td(k,2)*rra(k+1)+(tt(k,2)-td(k,2))*rxa(k+1))*denm
1265 ! rxa(k)=rs(k,2)+ts(k,2)*rxa(k+1)*denm
1266 ! end do
1267 !
1268 ! do k=1,np+1
1269 ! if (k <= np) then
1270 ! if (k == 1) then
1271 ! tdaold = td(0,2)
1272 ! ttaold = tt(0,2)
1273 ! rsaold = rs(0,2)
1274 !
1275 ! tdanew = 0.0
1276 ! ttanew = 0.0
1277 ! rsanew = 0.0
1278 ! end if
1279 ! denm=ts(k,2)/(1.-rsaold*rs(k,2))
1280 ! tdanew=tdaold*td(k,2)
1281 ! ttanew=tdaold*tt(k,2)+(tdaold*rsaold*rr(k,2)+ttaold-tdaold)*denm
1282 ! rsanew=rs(k,2)+ts(k,2)*rsaold*denm
1283 ! end if
1284 !
1285 ! denm=1./(1.-rsaold*rxa(k))
1286 ! fdndir=tdaold
1287 ! xx4=tdaold*rra(k)
1288 ! yy=ttaold-tdaold
1289 ! fdndif=(xx4*rsaold+yy)*denm
1290 ! fupdif=(xx4+yy*rxa(k))*denm
1291 ! flxdn=fdndir+fdndif-fupdif
1292 !
1293 ! fupa(k)=fupdif
1294 ! fall(k)=flxdn
1295 !
1296 ! tdaold = tdanew
1297 ! ttaold = ttanew
1298 ! rsaold = rsanew
1299 !
1300 ! tdanew = 0.0
1301 ! ttanew = 0.0
1302 ! rsanew = 0.0
1303 ! end do
1304 !
1305 ! fsdir=fdndir
1306 ! fsdif=fdndif
1307 !
1308 !!-----End CLDFLXY inline
1309 !
1310 !else
1311 !-----for clear- and all-sky flux calculations when fractional
1312 ! cloud cover is allowed to be between 0 and 1.
1313 ! the all-sky flux, fall is the summation inside the brackets
1314 ! of Eq. (7.11)
1315 !-----Inline CLDFLX
1316 !-----compute transmittances and reflectances for a composite of
1317 ! layers. layers are added one at a time, going down from the top.
1318 ! tda is the composite direct transmittance illuminated
1319 ! by beam radiation
1320 ! tta is the composite total transmittance illuminated by
1321 ! beam radiation
1322 ! rsa is the composite reflectance illuminated from below
1323 ! by diffuse radiation
1324 ! tta and rsa are computed from Eqs. (6.10) and (6.12)
1325 !-----To save memory space, tda, tta, and rsa are pre-computed
1326 ! for k<icb. The dimension of these parameters is (m,np,2,2).
1327 ! It would have been (m,np,2,2,2) if these parameters were
1328 ! computed for all k's.
1329 !-----for high clouds
1330 ! ih=1 for clear-sky condition, ih=2 for cloudy-sky condition
1331  DO ih=1,2
1332  tdad(0, ih, 1) = tdd(0, ih)
1333  tda(0, ih, 1) = td(0, ih)
1334  ttad(0, ih, 1) = ttd(0, ih)
1335  tta(0, ih, 1) = tt(0, ih)
1336  rsad(0, ih, 1) = rsd(0, ih)
1337  rsa(0, ih, 1) = rs(0, ih)
1338  tdad(0, ih, 2) = tdd(0, ih)
1339  tda(0, ih, 2) = td(0, ih)
1340  ttad(0, ih, 2) = ttd(0, ih)
1341  tta(0, ih, 2) = tt(0, ih)
1342  rsad(0, ih, 2) = rsd(0, ih)
1343  rsa(0, ih, 2) = rs(0, ih)
1344  DO k=1,ict-1
1345  denmd = (tsd(k, ih)*(1.-rsa(k-1, ih, 1)*rs(k, ih))-ts(k, ih)*(&
1346 & -(rsad(k-1, ih, 1)*rs(k, ih))-rsa(k-1, ih, 1)*rsd(k, ih)))/(&
1347 & 1.-rsa(k-1, ih, 1)*rs(k, ih))**2
1348  denm = ts(k, ih)/(1.-rsa(k-1, ih, 1)*rs(k, ih))
1349  tdad(k, ih, 1) = tdad(k-1, ih, 1)*td(k, ih) + tda(k-1, ih, 1)*&
1350 & tdd(k, ih)
1351  tda(k, ih, 1) = tda(k-1, ih, 1)*td(k, ih)
1352  ttad(k, ih, 1) = tdad(k-1, ih, 1)*tt(k, ih) + tda(k-1, ih, 1)*&
1353 & ttd(k, ih) + ((tdad(k-1, ih, 1)*rr(k, ih)+tda(k-1, ih, 1)*&
1354 & rrd(k, ih))*rsa(k-1, ih, 1)+tda(k-1, ih, 1)*rr(k, ih)*rsad(k&
1355 & -1, ih, 1)+ttad(k-1, ih, 1)-tdad(k-1, ih, 1))*denm + (tda(k-&
1356 & 1, ih, 1)*rsa(k-1, ih, 1)*rr(k, ih)+tta(k-1, ih, 1)-tda(k-1&
1357 & , ih, 1))*denmd
1358  tta(k, ih, 1) = tda(k-1, ih, 1)*tt(k, ih) + (tda(k-1, ih, 1)*&
1359 & rsa(k-1, ih, 1)*rr(k, ih)+tta(k-1, ih, 1)-tda(k-1, ih, 1))*&
1360 & denm
1361  rsad(k, ih, 1) = rsd(k, ih) + (tsd(k, ih)*denm+ts(k, ih)*denmd&
1362 & )*rsa(k-1, ih, 1) + ts(k, ih)*denm*rsad(k-1, ih, 1)
1363  rsa(k, ih, 1) = rs(k, ih) + ts(k, ih)*rsa(k-1, ih, 1)*denm
1364  tdad(k, ih, 2) = tdad(k, ih, 1)
1365  tda(k, ih, 2) = tda(k, ih, 1)
1366  ttad(k, ih, 2) = ttad(k, ih, 1)
1367  tta(k, ih, 2) = tta(k, ih, 1)
1368  rsad(k, ih, 2) = rsad(k, ih, 1)
1369  rsa(k, ih, 2) = rsa(k, ih, 1)
1370  END DO
1371 ! k loop
1372 !-----for middle clouds
1373 ! im=1 for clear-sky condition, im=2 for cloudy-sky condition
1374  DO k=ict,icb-1
1375  DO im=1,2
1376  denmd = (tsd(k, im)*(1.-rsa(k-1, ih, im)*rs(k, im))-ts(k, im&
1377 & )*(-(rsad(k-1, ih, im)*rs(k, im))-rsa(k-1, ih, im)*rsd(k, &
1378 & im)))/(1.-rsa(k-1, ih, im)*rs(k, im))**2
1379  denm = ts(k, im)/(1.-rsa(k-1, ih, im)*rs(k, im))
1380  tdad(k, ih, im) = tdad(k-1, ih, im)*td(k, im) + tda(k-1, ih&
1381 & , im)*tdd(k, im)
1382  tda(k, ih, im) = tda(k-1, ih, im)*td(k, im)
1383  ttad(k, ih, im) = tdad(k-1, ih, im)*tt(k, im) + tda(k-1, ih&
1384 & , im)*ttd(k, im) + ((tdad(k-1, ih, im)*rr(k, im)+tda(k-1, &
1385 & ih, im)*rrd(k, im))*rsa(k-1, ih, im)+tda(k-1, ih, im)*rr(k&
1386 & , im)*rsad(k-1, ih, im)+ttad(k-1, ih, im)-tdad(k-1, ih, im&
1387 & ))*denm + (tda(k-1, ih, im)*rsa(k-1, ih, im)*rr(k, im)+tta&
1388 & (k-1, ih, im)-tda(k-1, ih, im))*denmd
1389  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, im) + (tda(k-1, ih, &
1390 & im)*rsa(k-1, ih, im)*rr(k, im)+tta(k-1, ih, im)-tda(k-1, &
1391 & ih, im))*denm
1392  rsad(k, ih, im) = rsd(k, im) + (tsd(k, im)*denm+ts(k, im)*&
1393 & denmd)*rsa(k-1, ih, im) + ts(k, im)*denm*rsad(k-1, ih, im)
1394  rsa(k, ih, im) = rs(k, im) + ts(k, im)*rsa(k-1, ih, im)*denm
1395  END DO
1396  END DO
1397  END DO
1398 ! im loop
1399 ! k loop
1400 ! ih loop
1401 !-----layers are added one at a time, going up from the surface.
1402 ! rra is the composite reflectance illuminated by beam radiation
1403 ! rxa is the composite reflectance illuminated from above
1404 ! by diffuse radiation
1405 ! rra and rxa are computed from Eqs. (6.9) and (6.11)
1406 !-----To save memory space, rra and rxa are pre-computed for k>=icb.
1407 ! the dimension of these parameters is (m,np,2,2). It would have
1408 ! been (m,np,2,2,2) if these parameters were computed for all k's.
1409 !-----for the low clouds
1410 ! is=1 for clear-sky condition, is=2 for cloudy-sky condition
1411  DO is=1,2
1412  rrad(np+1, 1, is) = rrd(np+1, is)
1413  rra(np+1, 1, is) = rr(np+1, is)
1414  rxad(np+1, 1, is) = rsd(np+1, is)
1415  rxa(np+1, 1, is) = rs(np+1, is)
1416  rrad(np+1, 2, is) = rrd(np+1, is)
1417  rra(np+1, 2, is) = rr(np+1, is)
1418  rxad(np+1, 2, is) = rsd(np+1, is)
1419  rxa(np+1, 2, is) = rs(np+1, is)
1420  DO k=np,icb,-1
1421  denmd = (tsd(k, is)*(1.-rs(k, is)*rxa(k+1, 1, is))-ts(k, is)*(&
1422 & -(rsd(k, is)*rxa(k+1, 1, is))-rs(k, is)*rxad(k+1, 1, is)))/(&
1423 & 1.-rs(k, is)*rxa(k+1, 1, is))**2
1424  denm = ts(k, is)/(1.-rs(k, is)*rxa(k+1, 1, is))
1425  rrad(k, 1, is) = rrd(k, is) + (tdd(k, is)*rra(k+1, 1, is)+td(k&
1426 & , is)*rrad(k+1, 1, is)+(ttd(k, is)-tdd(k, is))*rxa(k+1, 1, &
1427 & is)+(tt(k, is)-td(k, is))*rxad(k+1, 1, is))*denm + (td(k, is&
1428 & )*rra(k+1, 1, is)+(tt(k, is)-td(k, is))*rxa(k+1, 1, is))*&
1429 & denmd
1430  rra(k, 1, is) = rr(k, is) + (td(k, is)*rra(k+1, 1, is)+(tt(k, &
1431 & is)-td(k, is))*rxa(k+1, 1, is))*denm
1432  rxad(k, 1, is) = rsd(k, is) + (tsd(k, is)*denm+ts(k, is)*denmd&
1433 & )*rxa(k+1, 1, is) + ts(k, is)*denm*rxad(k+1, 1, is)
1434  rxa(k, 1, is) = rs(k, is) + ts(k, is)*rxa(k+1, 1, is)*denm
1435  rrad(k, 2, is) = rrad(k, 1, is)
1436  rra(k, 2, is) = rra(k, 1, is)
1437  rxad(k, 2, is) = rxad(k, 1, is)
1438  rxa(k, 2, is) = rxa(k, 1, is)
1439  END DO
1440 ! k loop
1441 !-----for middle clouds
1442  DO k=icb-1,ict,-1
1443  DO im=1,2
1444  denmd = (tsd(k, im)*(1.-rs(k, im)*rxa(k+1, im, is))-ts(k, im&
1445 & )*(-(rsd(k, im)*rxa(k+1, im, is))-rs(k, im)*rxad(k+1, im, &
1446 & is)))/(1.-rs(k, im)*rxa(k+1, im, is))**2
1447  denm = ts(k, im)/(1.-rs(k, im)*rxa(k+1, im, is))
1448  rrad(k, im, is) = rrd(k, im) + (tdd(k, im)*rra(k+1, im, is)+&
1449 & td(k, im)*rrad(k+1, im, is)+(ttd(k, im)-tdd(k, im))*rxa(k+&
1450 & 1, im, is)+(tt(k, im)-td(k, im))*rxad(k+1, im, is))*denm +&
1451 & (td(k, im)*rra(k+1, im, is)+(tt(k, im)-td(k, im))*rxa(k+1&
1452 & , im, is))*denmd
1453  rra(k, im, is) = rr(k, im) + (td(k, im)*rra(k+1, im, is)+(tt&
1454 & (k, im)-td(k, im))*rxa(k+1, im, is))*denm
1455  rxad(k, im, is) = rsd(k, im) + (tsd(k, im)*denm+ts(k, im)*&
1456 & denmd)*rxa(k+1, im, is) + ts(k, im)*denm*rxad(k+1, im, is)
1457  rxa(k, im, is) = rs(k, im) + ts(k, im)*rxa(k+1, im, is)*denm
1458  END DO
1459  END DO
1460  END DO
1461 ! im loop
1462 ! k loop
1463 ! is loop
1464 !-----integration over eight sky situations.
1465 ! ih, im, is denote high, middle and low cloud groups.
1466  DO ih=1,2
1467 !-----clear portion
1468  IF (ih .EQ. 1) THEN
1469  chd = -cc1d
1470  ch = 1.0 - cc1
1471 !-----cloudy portion
1472  ELSE
1473  chd = cc1d
1474  ch = cc1
1475  END IF
1476  DO im=1,2
1477 !-----clear portion
1478  IF (im .EQ. 1) THEN
1479  cmd = chd*(1.0-cc2) - ch*cc2d
1480  cm = ch*(1.0-cc2)
1481 !-----cloudy portion
1482  ELSE
1483  cmd = chd*cc2 + ch*cc2d
1484  cm = ch*cc2
1485  END IF
1486  DO is=1,2
1487 !-----clear portion
1488  IF (is .EQ. 1) THEN
1489  ctd = cmd*(1.0-cc3) - cm*cc3d
1490  ct = cm*(1.0-cc3)
1491 !-----cloudy portion
1492  ELSE
1493  ctd = cmd*cc3 + cm*cc3d
1494  ct = cm*cc3
1495  END IF
1496 !-----add one layer at a time, going down.
1497  DO k=icb,np
1498  denmd = (tsd(k, is)*(1.-rsa(k-1, ih, im)*rs(k, is))-ts(k, &
1499 & is)*(-(rsad(k-1, ih, im)*rs(k, is))-rsa(k-1, ih, im)*rsd&
1500 & (k, is)))/(1.-rsa(k-1, ih, im)*rs(k, is))**2
1501  denm = ts(k, is)/(1.-rsa(k-1, ih, im)*rs(k, is))
1502  tdad(k, ih, im) = tdad(k-1, ih, im)*td(k, is) + tda(k-1, &
1503 & ih, im)*tdd(k, is)
1504  tda(k, ih, im) = tda(k-1, ih, im)*td(k, is)
1505  ttad(k, ih, im) = tdad(k-1, ih, im)*tt(k, is) + tda(k-1, &
1506 & ih, im)*ttd(k, is) + ((tdad(k-1, ih, im)*rr(k, is)+tda(k&
1507 & -1, ih, im)*rrd(k, is))*rsa(k-1, ih, im)+tda(k-1, ih, im&
1508 & )*rr(k, is)*rsad(k-1, ih, im)+ttad(k-1, ih, im)-tdad(k-1&
1509 & , ih, im))*denm + (tda(k-1, ih, im)*rr(k, is)*rsa(k-1, &
1510 & ih, im)+tta(k-1, ih, im)-tda(k-1, ih, im))*denmd
1511  tta(k, ih, im) = tda(k-1, ih, im)*tt(k, is) + (tda(k-1, ih&
1512 & , im)*rr(k, is)*rsa(k-1, ih, im)+tta(k-1, ih, im)-tda(k-&
1513 & 1, ih, im))*denm
1514  rsad(k, ih, im) = rsd(k, is) + (tsd(k, is)*denm+ts(k, is)*&
1515 & denmd)*rsa(k-1, ih, im) + ts(k, is)*denm*rsad(k-1, ih, &
1516 & im)
1517  rsa(k, ih, im) = rs(k, is) + ts(k, is)*rsa(k-1, ih, im)*&
1518 & denm
1519  END DO
1520 ! k loop
1521 !-----add one layer at a time, going up.
1522  DO k=ict-1,0,-1
1523  denmd = (tsd(k, ih)*(1.-rs(k, ih)*rxa(k+1, im, is))-ts(k, &
1524 & ih)*(-(rsd(k, ih)*rxa(k+1, im, is))-rs(k, ih)*rxad(k+1, &
1525 & im, is)))/(1.-rs(k, ih)*rxa(k+1, im, is))**2
1526  denm = ts(k, ih)/(1.-rs(k, ih)*rxa(k+1, im, is))
1527  rrad(k, im, is) = rrd(k, ih) + (tdd(k, ih)*rra(k+1, im, is&
1528 & )+td(k, ih)*rrad(k+1, im, is)+(ttd(k, ih)-tdd(k, ih))*&
1529 & rxa(k+1, im, is)+(tt(k, ih)-td(k, ih))*rxad(k+1, im, is)&
1530 & )*denm + (td(k, ih)*rra(k+1, im, is)+(tt(k, ih)-td(k, ih&
1531 & ))*rxa(k+1, im, is))*denmd
1532  rra(k, im, is) = rr(k, ih) + (td(k, ih)*rra(k+1, im, is)+(&
1533 & tt(k, ih)-td(k, ih))*rxa(k+1, im, is))*denm
1534  rxad(k, im, is) = rsd(k, ih) + (tsd(k, ih)*denm+ts(k, ih)*&
1535 & denmd)*rxa(k+1, im, is) + ts(k, ih)*denm*rxad(k+1, im, &
1536 & is)
1537  rxa(k, im, is) = rs(k, ih) + ts(k, ih)*rxa(k+1, im, is)*&
1538 & denm
1539  END DO
1540 ! k loop
1541 !-----compute fluxes following Eq. (6.15) for fupdif and
1542 ! Eq. (6.16) for (fdndir+fdndif)
1543 ! fdndir is the direct downward flux
1544 ! fdndif is the diffuse downward flux
1545 ! fupdif is the diffuse upward flux
1546  DO k=1,np+1
1547  denmd = -((-(rsad(k-1, ih, im)*rxa(k, im, is))-rsa(k-1, ih&
1548 & , im)*rxad(k, im, is))/(1.-rsa(k-1, ih, im)*rxa(k, im, &
1549 & is))**2)
1550  denm = 1./(1.-rsa(k-1, ih, im)*rxa(k, im, is))
1551  fdndird = tdad(k-1, ih, im)
1552  fdndir = tda(k-1, ih, im)
1553  xx4d = tdad(k-1, ih, im)*rra(k, im, is) + tda(k-1, ih, im)&
1554 & *rrad(k, im, is)
1555  xx4 = tda(k-1, ih, im)*rra(k, im, is)
1556  yyd = ttad(k-1, ih, im) - tdad(k-1, ih, im)
1557  yy = tta(k-1, ih, im) - tda(k-1, ih, im)
1558  fdndifd = (xx4d*rsa(k-1, ih, im)+xx4*rsad(k-1, ih, im)+yyd&
1559 & )*denm + (xx4*rsa(k-1, ih, im)+yy)*denmd
1560  fdndif = (xx4*rsa(k-1, ih, im)+yy)*denm
1561  fupdifd = (xx4d+yyd*rxa(k, im, is)+yy*rxad(k, im, is))*&
1562 & denm + (xx4+yy*rxa(k, im, is))*denmd
1563  fupdif = (xx4+yy*rxa(k, im, is))*denm
1564  flxdnd = fdndird + fdndifd - fupdifd
1565  flxdn = fdndir + fdndif - fupdif
1566 !-----summation of fluxes over all sky situations;
1567 ! the term in the brackets of Eq. (7.11)
1568  IF (ih .EQ. 1 .AND. im .EQ. 1 .AND. is .EQ. 1) THEN
1569  fupc(k) = fupdif
1570  fclr(k) = flxdn
1571  END IF
1572  fupa(k) = fupa(k) + fupdif*ct
1573  falld(k) = falld(k) + flxdnd*ct + flxdn*ctd
1574  fall(k) = fall(k) + flxdn*ct
1575  END DO
1576 ! k loop
1577  fsdir = fsdir + fdndir*ct
1578  fsdif = fsdif + fdndif*ct
1579  END DO
1580  END DO
1581  END DO
1582 ! is loop
1583 ! im loop
1584 ! ih loop
1585 !-----End CLDFLX inline
1586 !endif !overcast
1587 !-----flux integration following Eq. (6.1)
1588  DO k=1,np+1
1589  flx_devd(i, k) = flx_devd(i, k) + hk_ir(ib, ik)*falld(k)
1590  flx_dev(i, k) = flx_dev(i, k) + fall(k)*hk_ir(ib, ik)
1591  flc_dev(i, k) = flc_dev(i, k) + fclr(k)*hk_ir(ib, ik)
1592  flxu_dev(i, k) = flxu_dev(i, k) + fupa(k)*hk_ir(ib, ik)
1593  flcu_dev(i, k) = flcu_dev(i, k) + fupc(k)*hk_ir(ib, ik)
1594  END DO
1595 !-----compute downward surface fluxes in the ir region
1596  fdirir_dev(i) = fdirir_dev(i) + fsdir*hk_ir(ib, ik)
1597  fdifir_dev(i) = fdifir_dev(i) + fsdif*hk_ir(ib, ik)
1598 !-----tabulate surface flux at ir bands
1599  flx_sfc_band_dev(i, iv) = flx_sfc_band_dev(i, iv) + fall(np+1)*&
1600 & hk_ir(ib, ik)
1601  END DO
1602  END DO
1603 ! ik loop
1604 !-----compute pressure-scaled o2 amount following Eq. (3.5) with f=1.
1605 ! unit is (cm-atm)stp. 165.22 = (1000/980)*23.14%*(22400/32)
1606 ! compute flux reduction due to oxygen following Eq. (3.18). 0.0633 is the
1607 ! fraction of insolation contained in the oxygen bands
1608  dfd(0) = 0.0_8
1609  df(0) = 0.0
1610  cnt = 165.22*snt
1611  so2d(1) = 0.0_8
1612  so2(1) = scal0*cnt
1613 ! LLT increased parameter 145 to 155 to enhance effect
1614  result1 = sqrt(so2(1))
1615  dfd(1) = 0.0_8
1616  df(1) = 0.0633*(1.-exp(-(0.000155*result1)))
1617  DO k=1,np
1618  so2d(k+1) = 0.0_8
1619  so2(k+1) = so2(k) + scal(k)*cnt
1620 ! LLT increased parameter 145 to 155 to enhance effect
1621  result1 = sqrt(so2(k+1))
1622  dfd(k+1) = 0.0_8
1623  df(k+1) = 0.0633*(1.0-exp(-(0.000155*result1)))
1624  END DO
1625 !-----for solar heating due to co2 scaling follows Eq(3.5) with f=1.
1626 ! unit is (cm-atm)stp. 789 = (1000/980)*(44/28.97)*(22400/44)
1627  so2d(1) = 0.0_8
1628  so2(1) = 789.*co2*scal0
1629  DO k=1,np
1630  so2d(k+1) = 0.0_8
1631  so2(k+1) = so2(k) + 789.*co2*scal(k)
1632  END DO
1633 !-----The updated flux reduction for co2 absorption in Band 7 where absorption due to
1634 ! water vapor and co2 are both moderate. df is given by the second term on the
1635 ! right-hand-side of Eq. (3.24) divided by So. so2 and swh are the co2 and
1636 ! water vapor amounts integrated from the top of the atmosphere
1637  u1 = -3.0
1638  du = 0.15
1639  w1 = -4.0
1640  dw = 0.15
1641 !-----Inline RFLX
1642  du2 = du*du
1643  dw2 = dw*dw
1644  x0 = u1 + REAL(nu)*du
1645  y0 = w1 + REAL(nw)*dw
1646  x1 = u1 - 0.5*du
1647  y1 = w1 - 0.5*dw
1648  dfd = 0.0_8
1649  DO k=1,np+1
1650  x3 = log10(so2(k)*snt)
1651  IF (x3 .GT. x0) THEN
1652  ulog = x0
1653  ELSE
1654  ulog = x3
1655  END IF
1656  x4d = swhd(k)/(swh(k)*log(10.0))
1657  x4 = log10(swh(k)*snt)
1658  IF (x4 .GT. y0) THEN
1659  wlog = y0
1660  wlogd = 0.0_8
1661  ELSE
1662  wlogd = x4d
1663  wlog = x4
1664  END IF
1665  ic = int((ulog-x1)/du + 1.)
1666  iw = int((wlog-y1)/dw + 1.)
1667  IF (ic .LT. 2) ic = 2
1668  IF (iw .LT. 2) iw = 2
1669  IF (ic .GT. nu) ic = nu
1670  IF (iw .GT. nw) iw = nw
1671  dc = ulog - REAL(ic-2)*du - u1
1672  ddd = wlogd
1673  dd = wlog - REAL(iw-2)*dw - w1
1674  x2d = (cah(ic-1, iw)-cah(ic-1, iw-1))*ddd/dw
1675  x2 = cah(ic-1, iw-1) + (cah(ic-1, iw)-cah(ic-1, iw-1))/dw*dd
1676  y2d = x2d
1677  y2 = x2 + (cah(ic, iw-1)-cah(ic-1, iw-1))/du*dc
1678  IF (y2 .LT. 0.0) THEN
1679  y2 = 0.0
1680  y2d = 0.0_8
1681  ELSE
1682  y2 = y2
1683  END IF
1684 ! LLT increase CO2 effect to help reduce cold tropopause bias
1685  dfd(k) = dfd(k) + 1.5*y2d
1686  df(k) = df(k) + 1.5*y2
1687  END DO
1688 !-----df is the updated flux reduction for co2 absorption
1689 ! in Band 8 where the co2 absorption has a large impact
1690 ! on the heating of middle atmosphere. From the table
1691 ! given by Eq. (3.19)
1692  u1 = 0.000250
1693  du = 0.000050
1694  w1 = -2.0
1695  dw = 0.05
1696 !-----Inline RFLX
1697  du2 = du*du
1698  dw2 = dw*dw
1699  x0 = u1 + REAL(nx)*du
1700  y0 = w1 + REAL(ny)*dw
1701  x1 = u1 - 0.5*du
1702  y1 = w1 - 0.5*dw
1703  DO k=1,np+1
1704  IF (co2*snt .GT. x0) THEN
1705  ulog = x0
1706  ELSE
1707  ulog = co2*snt
1708  END IF
1709  x5 = log10(pl_dev(i, k))
1710  IF (x5 .GT. y0) THEN
1711  wlog = y0
1712  ELSE
1713  wlog = x5
1714  END IF
1715  ic = int((ulog-x1)/du + 1.)
1716  iw = int((wlog-y1)/dw + 1.)
1717  IF (ic .LT. 2) ic = 2
1718  IF (iw .LT. 2) iw = 2
1719  IF (ic .GT. nx) ic = nx
1720  IF (iw .GT. ny) iw = ny
1721  dc = ulog - REAL(ic-2)*du - u1
1722  dd = wlog - REAL(iw-2)*dw - w1
1723  x2 = coa(ic-1, iw-1) + (coa(ic-1, iw)-coa(ic-1, iw-1))/dw*dd
1724  y2 = x2 + (coa(ic, iw-1)-coa(ic-1, iw-1))/du*dc
1725  IF (y2 .LT. 0.0) THEN
1726  y2 = 0.0
1727  ELSE
1728  y2 = y2
1729  END IF
1730 ! LLT increase CO2 effect to help reduce cold tropopause bias
1731  df(k) = df(k) + 1.5*y2
1732  END DO
1733 !-----adjust the o2-co2 reduction below cloud top following Eq. (6.18)
1734  foundtop = 0
1735  DO k=1,np
1736  IF (fcld_dev(i, k) .GT. 0.02 .AND. foundtop .EQ. 0) THEN
1737  foundtop = 1
1738  ntop = k
1739  END IF
1740  END DO
1741  IF (foundtop .EQ. 0) ntop = np + 1
1742  dftopd = dfd(ntop)
1743  dftop = df(ntop)
1744  DO k=1,np+1
1745  IF (k .GT. ntop) THEN
1746  xx4d = (flx_devd(i, k)*flx_dev(i, ntop)-flx_dev(i, k)*flx_devd(i, &
1747 & ntop))/flx_dev(i, ntop)**2
1748  xx4 = flx_dev(i, k)/flx_dev(i, ntop)
1749  dfd(k) = dftopd + xx4d*(df(k)-dftop) + xx4*(dfd(k)-dftopd)
1750  df(k) = dftop + xx4*(df(k)-dftop)
1751  END IF
1752  END DO
1753 !-----update the net fluxes
1754  DO k=1,np+1
1755  IF (df(k) .GT. flx_dev(i, k) - 1.0e-8) THEN
1756  dfd(k) = flx_devd(i, k)
1757  df(k) = flx_dev(i, k) - 1.0e-8
1758  ELSE
1759  df(k) = df(k)
1760  END IF
1761 ! df(k) = 0.0
1762  flx_devd(i, k) = flx_devd(i, k) - dfd(k)
1763  flx_dev(i, k) = flx_dev(i, k) - df(k)
1764  flc_dev(i, k) = flc_dev(i, k) - df(k)
1765  END DO
1766 !-----update the downward surface fluxes
1767 ! xx4 = fdirir (i) + fdifir (i) +&
1768 ! fdiruv (i) + fdifuv (i) +&
1769 ! fdirpar(i) + fdifpar(i)
1770  xx4 = flx_dev(i, np+1) + df(np+1)
1771  IF (xx4 .GE. 0.) THEN
1772  abs0 = xx4
1773  ELSE
1774  abs0 = -xx4
1775  END IF
1776  result10 = epsilon(1.0)
1777  IF (abs0 .GT. result10) THEN
1778  IF (1.0 - df(np+1)/xx4 .GT. 1.) THEN
1779  x6 = 1.
1780  ELSE
1781  x6 = 1.0 - df(np+1)/xx4
1782  END IF
1783  IF (x6 .LT. 0.) THEN
1784  xx4 = 0.
1785  ELSE
1786  xx4 = x6
1787  END IF
1788  ELSE
1789  xx4 = 0.0
1790  END IF
1791  fdirir_dev(i) = xx4*fdirir_dev(i)
1792  fdifir_dev(i) = xx4*fdifir_dev(i)
1793  fdiruv_dev(i) = xx4*fdiruv_dev(i)
1794  fdifuv_dev(i) = xx4*fdifuv_dev(i)
1795  fdirpar_dev(i) = xx4*fdirpar_dev(i)
1796  fdifpar_dev(i) = xx4*fdifpar_dev(i)
1797  DO ib=1,nband
1798  flx_sfc_band_dev(i, ib) = xx4*flx_sfc_band_dev(i, ib)
1799  END DO
1800 END SUBROUTINE sorad_d
1801 
1802 ! Differentiation of deledd in forward (tangent) mode:
1803 ! variations of useful results: tt1 td1 rr1
1804 ! with respect to varying inputs: g01 tau1 ssc1
1805 !*********************************************************************
1806 SUBROUTINE deledd_d(tau1, tau1d, ssc1, ssc1d, g01, g01d, cza1, rr1, rr1d&
1807 & , tt1, tt1d, td1, td1d)
1808  IMPLICIT NONE
1809 ! 8 byte real
1810  INTEGER, PARAMETER :: real_de=8
1811 !integer,parameter :: REAL_SP = 4 ! 4 byte real
1812 !-----input parameters
1813  REAL*8, INTENT(IN) :: tau1, ssc1, g01, cza1
1814  REAL*8, INTENT(IN) :: tau1d, ssc1d, g01d
1815 !-----output parameters
1816  REAL*8, INTENT(OUT) :: rr1, tt1, td1
1817  REAL*8, INTENT(OUT) :: rr1d, tt1d, td1d
1818 !-----temporary parameters
1819  REAL*8, PARAMETER :: zero=0.0_real_de
1820  REAL*8, PARAMETER :: one=1.0_real_de
1821  REAL*8, PARAMETER :: two=2.0_real_de
1822  REAL*8, PARAMETER :: three=3.0_real_de
1823  REAL*8, PARAMETER :: four=4.0_real_de
1824  REAL*8, PARAMETER :: fourth=0.25_real_de
1825  REAL*8, PARAMETER :: seven=7.0_real_de
1826  REAL*8, PARAMETER :: thresh=1.e-8_real_de
1827  REAL*8 :: tau, ssc, g0, rr, tt, td
1828  REAL*8 :: taud, sscd, g0d, rrd, ttd, tdd
1829  REAL*8 :: zth, ff, xx, taup, sscp, gp, gm1, gm2, gm3, akk, alf1, alf2
1830  REAL*8 :: ffd, xxd, taupd, sscpd, gpd, gm1d, gm2d, gm3d, akkd, alf1d, &
1831 & alf2d
1832  REAL*8 :: all, bll, st7, st8, cll, dll, fll, ell, st1, st2, st3, st4
1833  REAL*8 :: alld, blld, st7d, st8d, clld, dlld, flld, elld, st1d, st2d, &
1834 & st3d, st4d
1835  INTRINSIC SQRT
1836  INTRINSIC abs
1837  INTRINSIC exp
1838  INTRINSIC max
1839  INTRINSIC real
1840  REAL*8 :: arg1
1841  REAL*8 :: arg1d
1842  REAL*8 :: abs0
1843 !zth = real(cza1,kind=REAL_DE)
1844 !g0 = real(g01 ,kind=REAL_DE)
1845 !tau = real(tau1,kind=REAL_DE)
1846 !ssc = real(ssc1,kind=REAL_DE)
1847 !zth = dble(cza1)
1848 !g0 = dble(g01)
1849 !tau = dble(tau1)
1850 !ssc = dble(ssc1)
1851  zth = cza1
1852  g0d = g01d
1853  g0 = g01
1854  taud = tau1d
1855  tau = tau1
1856  sscd = ssc1d
1857  ssc = ssc1
1858  ffd = g0d*g0 + g0*g0d
1859  ff = g0*g0
1860  xxd = -(ffd*ssc) - ff*sscd
1861  xx = one - ff*ssc
1862  taupd = taud*xx + tau*xxd
1863  taup = tau*xx
1864  sscpd = ((sscd*(one-ff)-ssc*ffd)*xx-ssc*(one-ff)*xxd)/xx**2
1865  sscp = ssc*(one-ff)/xx
1866  gpd = (g0d*(one+g0)-g0*g0d)/(one+g0)**2
1867  gp = g0/(one+g0)
1868  xxd = three*gpd
1869  xx = three*gp
1870  gm1d = fourth*(-(sscpd*(four+xx))-sscp*xxd)
1871  gm1 = (seven-sscp*(four+xx))*fourth
1872  gm2d = -(fourth*(sscp*xxd-sscpd*(four-xx)))
1873  gm2 = -((one-sscp*(four-xx))*fourth)
1874  arg1d = (gm1d+gm2d)*(gm1-gm2) + (gm1+gm2)*(gm1d-gm2d)
1875  arg1 = (gm1+gm2)*(gm1-gm2)
1876  IF (arg1 .EQ. 0.0) THEN
1877  akkd = 0.0_8
1878  ELSE
1879  akkd = arg1d/(2.0*sqrt(arg1))
1880  END IF
1881  akk = sqrt(arg1)
1882  xxd = zth*akkd
1883  xx = akk*zth
1884  st7d = -xxd
1885  st7 = one - xx
1886  st8d = xxd
1887  st8 = one + xx
1888  st3d = st7d*st8 + st7*st8d
1889  st3 = st7*st8
1890  IF (st3 .GE. 0.) THEN
1891  abs0 = st3
1892  ELSE
1893  abs0 = -st3
1894  END IF
1895  IF (abs0 .LT. thresh) THEN
1896  zth = zth + 0.0010
1897  IF (zth .GT. 1.0) zth = zth - 0.0020
1898  xxd = zth*akkd
1899  xx = akk*zth
1900  st7d = -xxd
1901  st7 = one - xx
1902  st8d = xxd
1903  st8 = one + xx
1904  st3d = st7d*st8 + st7*st8d
1905  st3 = st7*st8
1906  END IF
1907  tdd = -(taupd*exp(-(taup/zth))/zth)
1908  td = exp(-(taup/zth))
1909  gm3d = -(fourth*zth*three*gpd)
1910  gm3 = (two-zth*three*gp)*fourth
1911  xxd = gm1d - gm2d
1912  xx = gm1 - gm2
1913  alf1d = gm1d - gm3d*xx - gm3*xxd
1914  alf1 = gm1 - gm3*xx
1915  alf2d = gm2d + gm3d*xx + gm3*xxd
1916  alf2 = gm2 + gm3*xx
1917  xxd = two*akkd
1918  xx = akk*two
1919  alld = (gm3d-zth*alf2d)*xx*td + (gm3-alf2*zth)*(xxd*td+xx*tdd)
1920  all = (gm3-alf2*zth)*xx*td
1921  blld = (zth*alf1d-gm3d)*xx + (one-gm3+alf1*zth)*xxd
1922  bll = (one-gm3+alf1*zth)*xx
1923  xxd = akkd*gm3 + akk*gm3d
1924  xx = akk*gm3
1925  clld = (alf2d+xxd)*st7 + (alf2+xx)*st7d
1926  cll = (alf2+xx)*st7
1927  dlld = (alf2d-xxd)*st8 + (alf2-xx)*st8d
1928  dll = (alf2-xx)*st8
1929  xxd = akkd*(one-gm3) - akk*gm3d
1930  xx = akk*(one-gm3)
1931  flld = (alf1d+xxd)*st8 + (alf1+xx)*st8d
1932  fll = (alf1+xx)*st8
1933  elld = (alf1d-xxd)*st7 + (alf1-xx)*st7d
1934  ell = (alf1-xx)*st7
1935  st2d = -((akkd*taup+akk*taupd)*exp(-(akk*taup)))
1936  st2 = exp(-(akk*taup))
1937  st4d = st2d*st2 + st2*st2d
1938  st4 = st2*st2
1939  st1d = (sscpd*(akk+gm1+(akk-gm1)*st4)*st3-sscp*((akkd+gm1d+(akkd-gm1d)&
1940 & *st4+(akk-gm1)*st4d)*st3+(akk+gm1+(akk-gm1)*st4)*st3d))/((akk+gm1+(&
1941 & akk-gm1)*st4)*st3)**2
1942  st1 = sscp/((akk+gm1+(akk-gm1)*st4)*st3)
1943  rrd = (clld-dlld*st4-dll*st4d-alld*st2-all*st2d)*st1 + (cll-dll*st4-&
1944 & all*st2)*st1d
1945  rr = (cll-dll*st4-all*st2)*st1
1946  ttd = -(((flld-elld*st4-ell*st4d)*td+(fll-ell*st4)*tdd-blld*st2-bll*&
1947 & st2d)*st1+((fll-ell*st4)*td-bll*st2)*st1d)
1948  tt = -(((fll-ell*st4)*td-bll*st2)*st1)
1949  IF (rr .LT. zero) THEN
1950  rr = zero
1951  rrd = 0.0_8
1952  ELSE
1953  rr = rr
1954  END IF
1955  IF (tt .LT. zero) THEN
1956  tt = zero
1957  ttd = 0.0_8
1958  ELSE
1959  tt = tt
1960  END IF
1961  ttd = ttd + tdd
1962  tt = tt + td
1963 !td1 = real(td,kind=REAL_SP)
1964 !rr1 = real(rr,kind=REAL_SP)
1965 !tt1 = real(tt,kind=REAL_SP)
1966  td1d = tdd
1967  td1 = REAL(td)
1968  rr1d = rrd
1969  rr1 = REAL(rr)
1970  tt1d = ttd
1971  tt1 = REAL(tt)
1972 END SUBROUTINE deledd_d
1973 
1974 ! Differentiation of getvistau1 in forward (tangent) mode:
1975 ! variations of useful results: asycl taudiff taubeam
1976 ! with respect to varying inputs: hydromets fcld reff
1977 SUBROUTINE getvistau1_d(nlevs, cosz, dp, fcld, fcldd, reff, reffd, &
1978 & hydromets, hydrometsd, ict, icb, taubeam, taubeamd, taudiff, taudiffd&
1979 & , asycl, asycld, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv, arb_uv, &
1980 & aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, aig_nir, awg_nir&
1981 & , arg_nir, caib, caif, cons_grav)
1982  IMPLICIT NONE
1983 !EOC
1984 ! !INPUT PARAMETERS:
1985 ! Number of levels
1986  INTEGER, INTENT(IN) :: nlevs
1987 ! Cosine of solar zenith angle
1988  REAL*8, INTENT(IN) :: cosz
1989 ! Delta pressure (Pa)
1990  REAL*8, INTENT(IN) :: dp(nlevs)
1991 ! Cloud fraction (used sometimes)
1992  REAL*8, INTENT(IN) :: fcld(nlevs)
1993  REAL*8, INTENT(IN) :: fcldd(nlevs)
1994 ! Effective radius (microns)
1995  REAL*8, INTENT(IN) :: reff(nlevs, 4)
1996  REAL*8, INTENT(IN) :: reffd(nlevs, 4)
1997 ! Hydrometeors (kg/kg)
1998  REAL*8, INTENT(IN) :: hydromets(nlevs, 4)
1999  REAL*8, INTENT(IN) :: hydrometsd(nlevs, 4)
2000 ! Flags for various uses
2001  INTEGER, INTENT(IN) :: ict, icb
2002 ! ict = 0 Indicates that in-cloud values have been given
2003 ! and are expected
2004 ! ict != 0 Indicates that overlap computation is needed, and:
2005 ! ict is the level of the mid-high boundary
2006 ! icb is the level of the low-mid boundary
2007 !
2008 ! !OUTPUT PARAMETERS:
2009 ! Optical Depth for Beam Radiation
2010  REAL*8, INTENT(OUT) :: taubeam(nlevs, 4)
2011  REAL*8, INTENT(OUT) :: taubeamd(nlevs, 4)
2012 ! Optical Depth for Diffuse Radiation
2013  REAL*8, INTENT(OUT) :: taudiff(nlevs, 4)
2014  REAL*8, INTENT(OUT) :: taudiffd(nlevs, 4)
2015 ! Cloud Asymmetry Factor
2016  REAL*8, INTENT(OUT) :: asycl(nlevs)
2017  REAL*8, INTENT(OUT) :: asycld(nlevs)
2018 ! !DESCRIPTION:
2019 ! Compute in-cloud or grid mean optical depths for visible wavelengths
2020 ! In general will compute in-cloud - will do grid mean when called
2021 ! for diagnostic use only. ict flag will indicate which to do.
2022 ! Slots for reff, hydrometeors, taubeam, taudiff, and asycl are as follows:
2023 ! 1 Cloud Ice
2024 ! 2 Cloud Liquid
2025 ! 3 Falling Liquid (Rain)
2026 ! 4 Falling Ice (Snow)
2027 !
2028 ! In the below calculations, the constants used in the tau calculation are in
2029 ! m$^2$ g$^{-1}$ and m$^2$ g$^{-1}$ $\mu$m. Thus, we must convert the kg contained in the
2030 ! pressure (Pa = kg m$^{-1}$ s$^{-2}$) to grams.
2031 !
2032 ! !REVISION HISTORY:
2033 ! 2011.10.27 Molod moved to Radiation_Shared and revised arg list, units
2034 ! 2011.11.16 MAT: Generalized to a call that is per-column
2035 !
2036 !EOP
2037 !------------------------------------------------------------------------------
2038 !BOC
2039  INTEGER :: k, in, im, it, ia, kk
2040  REAL*8 :: fm, ft, fa, xai, tauc, asyclt
2041  REAL*8 :: ftd, fad, xaid, taucd, asycltd
2042  REAL*8 :: cc(3)
2043  REAL*8 :: ccd(3)
2044  REAL*8 :: taucld1, taucld2, taucld3, taucld4
2045  REAL*8 :: taucld1d, taucld2d, taucld3d, taucld4d
2046  REAL*8 :: g1, g2, g3, g4
2047  REAL*8 :: g1d, g2d, g3d, g4d
2048  REAL*8 :: reff_snow
2049  REAL*8 :: reff_snowd
2050  INTEGER, PARAMETER :: nm=11, nt=9, na=11
2051  REAL*8, PARAMETER :: dm=0.1, dt=0.30103, da=0.1, t1=-0.9031
2052  REAL*8, INTENT(IN) :: aig_uv(3), awg_uv(3), arg_uv(3)
2053  REAL*8, INTENT(IN) :: aib_uv, awb_uv(2), arb_uv(2)
2054  REAL*8, INTENT(IN) :: aib_nir, awb_nir(3, 2), arb_nir(3, 2)
2055  REAL*8, INTENT(IN) :: aia_nir(3, 3), awa_nir(3, 3), ara_nir(3, 3)
2056  REAL*8, INTENT(IN) :: aig_nir(3, 3), awg_nir(3, 3), arg_nir(3, 3)
2057  REAL*8, INTENT(IN) :: caib(11, 9, 11), caif(9, 11)
2058  REAL*8, INTENT(IN) :: cons_grav
2059  INTRINSIC max
2060  INTRINSIC min
2061  INTRINSIC log10
2062  INTRINSIC int
2063  INTRINSIC real
2064  taubeam = 0.0
2065  taudiff = 0.0
2066  IF (ict .NE. 0) THEN
2067 !-----scale cloud optical thickness in each layer from taucld (with
2068 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
2069 ! taubeam is the scaled optical thickness for beam radiation and
2070 ! taudiff is for diffuse radiation (see section 7).
2071 !-----clouds within each of the high, middle, and low clouds are
2072 ! assumed to be maximally overlapped, and the cloud cover (cc)
2073 ! for a group (high, middle, or low) is the maximum cloud cover
2074 ! of all the layers within a group
2075  cc = 0.0
2076  ccd = 0.0_8
2077  DO k=1,ict-1
2078  IF (cc(1) .LT. fcld(k)) THEN
2079  ccd(1) = fcldd(k)
2080  cc(1) = fcld(k)
2081  ELSE
2082  cc(1) = cc(1)
2083  END IF
2084  END DO
2085  DO k=ict,icb-1
2086  IF (cc(2) .LT. fcld(k)) THEN
2087  ccd(2) = fcldd(k)
2088  cc(2) = fcld(k)
2089  ELSE
2090  cc(2) = cc(2)
2091  END IF
2092  END DO
2093  DO k=icb,nlevs
2094  IF (cc(3) .LT. fcld(k)) THEN
2095  ccd(3) = fcldd(k)
2096  cc(3) = fcld(k)
2097  ELSE
2098  cc(3) = cc(3)
2099  END IF
2100  END DO
2101  asycld = 0.0_8
2102  taudiffd = 0.0_8
2103  taubeamd = 0.0_8
2104  ELSE
2105  asycld = 0.0_8
2106  taudiffd = 0.0_8
2107  taubeamd = 0.0_8
2108  ccd = 0.0_8
2109  END IF
2110 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10)
2111 ! Note: the cloud optical properties are assumed to be independent
2112 ! of spectral bands in the UV and PAR regions.
2113 ! taucld1 is the optical thickness for ice particles
2114 ! taucld2 is the optical thickness for liquid particles
2115 ! taucld3 is the optical thickness for rain drops
2116 ! taucld4 is the optical thickness for snow
2117  DO k=1,nlevs
2118  IF (reff(k, 1) .LE. 0.) THEN
2119  taucld1 = 0.
2120  taucld1d = 0.0_8
2121  ELSE
2122  taucld1d = (dp(k)*1.0e3*aib_uv*hydrometsd(k, 1)*reff(k, 1)/&
2123 & cons_grav-dp(k)*1.0e3*hydromets(k, 1)*aib_uv*reffd(k, 1)/&
2124 & cons_grav)/reff(k, 1)**2
2125  taucld1 = dp(k)*1.0e3/cons_grav*hydromets(k, 1)*aib_uv/reff(k, 1)
2126  END IF
2127  IF (reff(k, 2) .LE. 0.) THEN
2128  taucld2 = 0.
2129  taucld2d = 0.0_8
2130  ELSE
2131  taucld2d = dp(k)*1.0e3*(hydrometsd(k, 2)*(awb_uv(1)+awb_uv(2)/reff&
2132 & (k, 2))-hydromets(k, 2)*awb_uv(2)*reffd(k, 2)/reff(k, 2)**2)/&
2133 & cons_grav
2134  taucld2 = dp(k)*1.0e3/cons_grav*hydromets(k, 2)*(awb_uv(1)+awb_uv(&
2135 & 2)/reff(k, 2))
2136  END IF
2137  taucld3d = dp(k)*1.0e3*arb_uv(1)*hydrometsd(k, 3)/cons_grav
2138  taucld3 = dp(k)*1.0e3/cons_grav*hydromets(k, 3)*arb_uv(1)
2139  IF (reff(k, 4) .GT. 112.0) THEN
2140  reff_snow = 112.0
2141  reff_snowd = 0.0_8
2142  ELSE
2143  reff_snowd = reffd(k, 4)
2144  reff_snow = reff(k, 4)
2145  END IF
2146  IF (reff_snow .LE. 0.) THEN
2147  taucld4 = 0.
2148  taucld4d = 0.0_8
2149  ELSE
2150  taucld4d = (dp(k)*1.0e3*aib_uv*hydrometsd(k, 4)*reff_snow/&
2151 & cons_grav-dp(k)*1.0e3*hydromets(k, 4)*aib_uv*reff_snowd/&
2152 & cons_grav)/reff_snow**2
2153  taucld4 = dp(k)*1.0e3/cons_grav*hydromets(k, 4)*aib_uv/reff_snow
2154  END IF
2155  IF (ict .NE. 0) THEN
2156 !-----scale cloud optical thickness in each layer from taucld (with
2157 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
2158 ! taubeam is the scaled optical thickness for beam radiation and
2159 ! taudiff is for diffuse radiation (see section 7).
2160 !-----clouds within each of the high, middle, and low clouds are
2161 ! assumed to be maximally overlapped, and the cloud cover (cc)
2162 ! for a group (high, middle, or low) is the maximum cloud cover
2163 ! of all the layers within a group
2164  IF (k .LT. ict) THEN
2165  kk = 1
2166  ELSE IF (k .GE. ict .AND. k .LT. icb) THEN
2167  kk = 2
2168  ELSE
2169  kk = 3
2170  END IF
2171  taucd = taucld1d + taucld2d + taucld3d + taucld4d
2172  tauc = taucld1 + taucld2 + taucld3 + taucld4
2173  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
2174 !-----normalize cloud cover following Eq. (7.8)
2175  fad = (fcldd(k)*cc(kk)-fcld(k)*ccd(kk))/cc(kk)**2
2176  fa = fcld(k)/cc(kk)
2177  IF (tauc .GT. 32.) THEN
2178  tauc = 32.
2179  taucd = 0.0_8
2180  ELSE
2181  tauc = tauc
2182  END IF
2183  fm = cosz/dm
2184  ftd = taucd/(tauc*log(10.0))/dt
2185  ft = (log10(tauc)-t1)/dt
2186  fad = fad/da
2187  fa = fa/da
2188  im = int(fm + 1.5)
2189  it = int(ft + 1.5)
2190  ia = int(fa + 1.5)
2191  IF (im .LT. 2) THEN
2192  im = 2
2193  ELSE
2194  im = im
2195  END IF
2196  IF (it .LT. 2) THEN
2197  it = 2
2198  ELSE
2199  it = it
2200  END IF
2201  IF (ia .LT. 2) THEN
2202  ia = 2
2203  ELSE
2204  ia = ia
2205  END IF
2206  IF (im .GT. nm - 1) THEN
2207  im = nm - 1
2208  ELSE
2209  im = im
2210  END IF
2211  IF (it .GT. nt - 1) THEN
2212  it = nt - 1
2213  ELSE
2214  it = it
2215  END IF
2216  IF (ia .GT. na - 1) THEN
2217  ia = na - 1
2218  ELSE
2219  ia = ia
2220  END IF
2221  fm = fm - REAL(im - 1)
2222  ft = ft - REAL(it - 1)
2223  fa = fa - REAL(ia - 1)
2224 !-----scale cloud optical thickness for beam radiation following
2225 ! Eq. (7.3).
2226 ! the scaling factor, xai, is a function of the solar zenith
2227 ! angle, optical thickness, and cloud cover.
2228  xai = (-(caib(im-1, it, ia)*(1.-fm))+caib(im+1, it, ia)*(1.+fm))&
2229 & *fm*.5 + caib(im, it, ia)*(1.-fm*fm)
2230  xaid = .5*((caib(im, it-1, ia)*ftd+caib(im, it+1, ia)*ftd)*ft+(-&
2231 & (caib(im, it-1, ia)*(1.-ft))+caib(im, it+1, ia)*(1.+ft))*ftd) &
2232 & + caib(im, it, ia)*(-(ftd*ft)-ft*ftd)
2233  xai = xai + (-(caib(im, it-1, ia)*(1.-ft))+caib(im, it+1, ia)*(&
2234 & 1.+ft))*ft*.5 + caib(im, it, ia)*(1.-ft*ft)
2235  xaid = xaid + .5*((caib(im, it, ia-1)*fad+caib(im, it, ia+1)*fad&
2236 & )*fa+(-(caib(im, it, ia-1)*(1.-fa))+caib(im, it, ia+1)*(1.+fa)&
2237 & )*fad) + caib(im, it, ia)*(-(fad*fa)-fa*fad)
2238  xai = xai + (-(caib(im, it, ia-1)*(1.-fa))+caib(im, it, ia+1)*(&
2239 & 1.+fa))*fa*.5 + caib(im, it, ia)*(1.-fa*fa)
2240  xai = xai - 2.*caib(im, it, ia)
2241  IF (xai .LT. 0.0) THEN
2242  xai = 0.0
2243  xaid = 0.0_8
2244  ELSE
2245  xai = xai
2246  END IF
2247  IF (xai .GT. 1.0) THEN
2248  xai = 1.0
2249  xaid = 0.0_8
2250  ELSE
2251  xai = xai
2252  END IF
2253  taubeamd(k, 1) = taucld1d*xai + taucld1*xaid
2254  taubeam(k, 1) = taucld1*xai
2255  taubeamd(k, 2) = taucld2d*xai + taucld2*xaid
2256  taubeam(k, 2) = taucld2*xai
2257  taubeamd(k, 3) = taucld3d*xai + taucld3*xaid
2258  taubeam(k, 3) = taucld3*xai
2259  taubeamd(k, 4) = taucld4d*xai + taucld4*xaid
2260  taubeam(k, 4) = taucld4*xai
2261 !-----scale cloud optical thickness for diffuse radiation following
2262 ! Eq. (7.4).
2263 ! the scaling factor, xai, is a function of the cloud optical
2264 ! thickness and cover but not the solar zenith angle.
2265  xaid = .5*((caif(it-1, ia)*ftd+caif(it+1, ia)*ftd)*ft+(-(caif(it&
2266 & -1, ia)*(1.-ft))+caif(it+1, ia)*(1.+ft))*ftd) + caif(it, ia)*(&
2267 & -(ftd*ft)-ft*ftd)
2268  xai = (-(caif(it-1, ia)*(1.-ft))+caif(it+1, ia)*(1.+ft))*ft*.5 +&
2269 & caif(it, ia)*(1.-ft*ft)
2270  xaid = xaid + .5*((caif(it, ia-1)*fad+caif(it, ia+1)*fad)*fa+(-(&
2271 & caif(it, ia-1)*(1.-fa))+caif(it, ia+1)*(1.+fa))*fad) + caif(it&
2272 & , ia)*(-(fad*fa)-fa*fad)
2273  xai = xai + (-(caif(it, ia-1)*(1.-fa))+caif(it, ia+1)*(1.+fa))*&
2274 & fa*.5 + caif(it, ia)*(1.-fa*fa)
2275  xai = xai - caif(it, ia)
2276  IF (xai .LT. 0.0) THEN
2277  xai = 0.0
2278  xaid = 0.0_8
2279  ELSE
2280  xai = xai
2281  END IF
2282  IF (xai .GT. 1.0) THEN
2283  xai = 1.0
2284  xaid = 0.0_8
2285  ELSE
2286  xai = xai
2287  END IF
2288  taudiffd(k, 1) = taucld1d*xai + taucld1*xaid
2289  taudiff(k, 1) = taucld1*xai
2290  taudiffd(k, 2) = taucld2d*xai + taucld2*xaid
2291  taudiff(k, 2) = taucld2*xai
2292  taudiffd(k, 3) = taucld3d*xai + taucld3*xaid
2293  taudiff(k, 3) = taucld3*xai
2294  taudiffd(k, 4) = taucld4d*xai + taucld4*xaid
2295  taudiff(k, 4) = taucld4*xai
2296  END IF
2297  ELSE
2298 ! Overlap calculation scaling not needed
2299  taubeamd(k, 1) = taucld1d
2300  taubeam(k, 1) = taucld1
2301  taubeamd(k, 2) = taucld2d
2302  taubeam(k, 2) = taucld2
2303  taubeamd(k, 3) = taucld3d
2304  taubeam(k, 3) = taucld3
2305  taubeamd(k, 4) = taucld4d
2306  taubeam(k, 4) = taucld4
2307  taudiffd(k, 1) = taucld1d
2308  taudiff(k, 1) = taucld1
2309  taudiffd(k, 2) = taucld2d
2310  taudiff(k, 2) = taucld2
2311  taudiffd(k, 3) = taucld3d
2312  taudiff(k, 3) = taucld3
2313  taudiffd(k, 4) = taucld4d
2314  taudiff(k, 4) = taucld4
2315  END IF
2316 !-----cloud asymmetry factor for a mixture of liquid and ice particles.
2317 ! unit of reff is micrometers. Eqs. (4.8) and (6.4)
2318  asyclt = 1.0
2319  taucd = taucld1d + taucld2d + taucld3d + taucld4d
2320  tauc = taucld1 + taucld2 + taucld3 + taucld4
2321  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
2322  g1d = (aig_uv(3)*reffd(k, 1)*reff(k, 1)+(aig_uv(2)+aig_uv(3)*reff(&
2323 & k, 1))*reffd(k, 1))*taucld1 + (aig_uv(1)+(aig_uv(2)+aig_uv(3)*&
2324 & reff(k, 1))*reff(k, 1))*taucld1d
2325  g1 = (aig_uv(1)+(aig_uv(2)+aig_uv(3)*reff(k, 1))*reff(k, 1))*&
2326 & taucld1
2327  g2d = (awg_uv(3)*reffd(k, 2)*reff(k, 2)+(awg_uv(2)+awg_uv(3)*reff(&
2328 & k, 2))*reffd(k, 2))*taucld2 + (awg_uv(1)+(awg_uv(2)+awg_uv(3)*&
2329 & reff(k, 2))*reff(k, 2))*taucld2d
2330  g2 = (awg_uv(1)+(awg_uv(2)+awg_uv(3)*reff(k, 2))*reff(k, 2))*&
2331 & taucld2
2332  g3d = arg_uv(1)*taucld3d
2333  g3 = arg_uv(1)*taucld3
2334  g4d = (aig_uv(3)*reff_snowd*reff_snow+(aig_uv(2)+aig_uv(3)*&
2335 & reff_snow)*reff_snowd)*taucld4 + (aig_uv(1)+(aig_uv(2)+aig_uv(3)&
2336 & *reff_snow)*reff_snow)*taucld4d
2337  g4 = (aig_uv(1)+(aig_uv(2)+aig_uv(3)*reff_snow)*reff_snow)*taucld4
2338  asycltd = ((g1d+g2d+g3d+g4d)*tauc-(g1+g2+g3+g4)*taucd)/tauc**2
2339  asyclt = (g1+g2+g3+g4)/tauc
2340  ELSE
2341  asycltd = 0.0_8
2342  END IF
2343  asycld(k) = asycltd
2344  asycl(k) = asyclt
2345  END DO
2346  RETURN
2347 END SUBROUTINE getvistau1_d
2348 
2349 ! Differentiation of getnirtau1 in forward (tangent) mode:
2350 ! variations of useful results: asycl taudiff ssacl taubeam
2351 ! with respect to varying inputs: hydromets asycl fcld ssacl
2352 ! reff
2353 SUBROUTINE getnirtau1_d(ib, nlevs, cosz, dp, fcld, fcldd, reff, reffd, &
2354 & hydromets, hydrometsd, ict, icb, taubeam, taubeamd, taudiff, taudiffd&
2355 & , asycl, asycld, ssacl, ssacld, aig_uv, awg_uv, arg_uv, aib_uv, awb_uv&
2356 & , arb_uv, aib_nir, awb_nir, arb_nir, aia_nir, awa_nir, ara_nir, &
2357 & aig_nir, awg_nir, arg_nir, caib, caif, cons_grav)
2358  IMPLICIT NONE
2359 ! !INPUT PARAMETERS:
2360 ! Band number
2361  INTEGER, INTENT(IN) :: ib
2362 ! Number of levels
2363  INTEGER, INTENT(IN) :: nlevs
2364 ! Cosine of solar zenith angle
2365  REAL*8, INTENT(IN) :: cosz
2366 ! Delta pressure in Pa
2367  REAL*8, INTENT(IN) :: dp(nlevs)
2368 ! Cloud fraction (used sometimes)
2369  REAL*8, INTENT(IN) :: fcld(nlevs)
2370  REAL*8, INTENT(IN) :: fcldd(nlevs)
2371 ! Effective radius (microns)
2372  REAL*8, INTENT(IN) :: reff(nlevs, 4)
2373  REAL*8, INTENT(IN) :: reffd(nlevs, 4)
2374 ! Hydrometeors (kg/kg)
2375  REAL*8, INTENT(IN) :: hydromets(nlevs, 4)
2376  REAL*8, INTENT(IN) :: hydrometsd(nlevs, 4)
2377 ! Flags for various uses
2378  INTEGER, INTENT(IN) :: ict, icb
2379  REAL*8, INTENT(IN) :: aig_uv(3), awg_uv(3), arg_uv(3)
2380  REAL*8, INTENT(IN) :: aib_uv, awb_uv(2), arb_uv(2)
2381  REAL*8, INTENT(IN) :: aib_nir, awb_nir(3, 2), arb_nir(3, 2)
2382  REAL*8, INTENT(IN) :: aia_nir(3, 3), awa_nir(3, 3), ara_nir(3, 3)
2383  REAL*8, INTENT(IN) :: aig_nir(3, 3), awg_nir(3, 3), arg_nir(3, 3)
2384  REAL*8, INTENT(IN) :: caib(11, 9, 11), caif(9, 11)
2385  REAL*8, INTENT(IN) :: cons_grav
2386 ! !OUTPUT PARAMETERS:
2387 ! Optical depth for beam radiation
2388  REAL*8, INTENT(OUT) :: taubeam(nlevs, 4)
2389  REAL*8, INTENT(OUT) :: taubeamd(nlevs, 4)
2390 ! Optical depth for diffuse radiation
2391  REAL*8, INTENT(OUT) :: taudiff(nlevs, 4)
2392  REAL*8, INTENT(OUT) :: taudiffd(nlevs, 4)
2393 ! Cloud single scattering albedo
2394  REAL*8, INTENT(OUT) :: ssacl(nlevs)
2395  REAL*8, INTENT(OUT) :: ssacld(nlevs)
2396 ! Cloud asymmetry factor
2397  REAL*8, INTENT(OUT) :: asycl(nlevs)
2398  REAL*8, INTENT(OUT) :: asycld(nlevs)
2399  INTEGER :: k, in, im, it, ia, kk
2400  REAL*8 :: fm, ft, fa, xai, tauc, asyclt, ssaclt
2401  REAL*8 :: ftd, fad, xaid, taucd, asycltd, ssacltd
2402  REAL*8 :: cc(3)
2403  REAL*8 :: ccd(3)
2404  REAL*8 :: taucld1, taucld2, taucld3, taucld4
2405  REAL*8 :: taucld1d, taucld2d, taucld3d, taucld4d
2406  REAL*8 :: g1, g2, g3, g4
2407  REAL*8 :: g1d, g2d, g3d, g4d
2408  REAL*8 :: w1, w2, w3, w4
2409  REAL*8 :: w1d, w2d, w3d, w4d
2410  REAL*8 :: reff_snow
2411  REAL*8 :: reff_snowd
2412  INTEGER, PARAMETER :: nm=11, nt=9, na=11
2413  REAL*8, PARAMETER :: dm=0.1, dt=0.30103, da=0.1, t1=-0.9031
2414  INTRINSIC max
2415  INTRINSIC min
2416  INTRINSIC log10
2417  INTRINSIC int
2418  INTRINSIC real
2419  taubeam = 0.0
2420  taudiff = 0.0
2421  IF (ict .NE. 0) THEN
2422 !-----scale cloud optical thickness in each layer from taucld (with
2423 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
2424 ! taubeam is the scaled optical thickness for beam radiation and
2425 ! taudiff is for diffuse radiation (see section 7).
2426 !-----clouds within each of the high, middle, and low clouds are
2427 ! assumed to be maximally overlapped, and the cloud cover (cc)
2428 ! for a group (high, middle, or low) is the maximum cloud cover
2429 ! of all the layers within a group
2430  cc = 0.0
2431  ccd = 0.0_8
2432  DO k=1,ict-1
2433  IF (cc(1) .LT. fcld(k)) THEN
2434  ccd(1) = fcldd(k)
2435  cc(1) = fcld(k)
2436  ELSE
2437  cc(1) = cc(1)
2438  END IF
2439  END DO
2440  DO k=ict,icb-1
2441  IF (cc(2) .LT. fcld(k)) THEN
2442  ccd(2) = fcldd(k)
2443  cc(2) = fcld(k)
2444  ELSE
2445  cc(2) = cc(2)
2446  END IF
2447  END DO
2448  DO k=icb,nlevs
2449  IF (cc(3) .LT. fcld(k)) THEN
2450  ccd(3) = fcldd(k)
2451  cc(3) = fcld(k)
2452  ELSE
2453  cc(3) = cc(3)
2454  END IF
2455  END DO
2456  taudiffd = 0.0_8
2457  taubeamd = 0.0_8
2458  ELSE
2459  taudiffd = 0.0_8
2460  taubeamd = 0.0_8
2461  ccd = 0.0_8
2462  END IF
2463 !-----Compute cloud optical thickness. Eqs. (4.6) and (4.10)
2464 ! taucld1 is the optical thickness for ice particles
2465 ! taucld2 is the optical thickness for liquid particles
2466 ! taucld3 is the optical thickness for rain drops
2467 ! taucld4 is the optical thickness for snow
2468  DO k=1,nlevs
2469  IF (reff(k, 1) .LE. 0.) THEN
2470  taucld1 = 0.
2471  taucld1d = 0.0_8
2472  ELSE
2473  taucld1d = (dp(k)*1.0e3*aib_nir*hydrometsd(k, 1)*reff(k, 1)/&
2474 & cons_grav-dp(k)*1.0e3*hydromets(k, 1)*aib_nir*reffd(k, 1)/&
2475 & cons_grav)/reff(k, 1)**2
2476  taucld1 = dp(k)*1.0e3/cons_grav*hydromets(k, 1)*aib_nir/reff(k, 1)
2477  END IF
2478  IF (reff(k, 2) .LE. 0.) THEN
2479  taucld2 = 0.
2480  taucld2d = 0.0_8
2481  ELSE
2482  taucld2d = dp(k)*1.0e3*(hydrometsd(k, 2)*(awb_nir(ib, 1)+awb_nir(&
2483 & ib, 2)/reff(k, 2))-hydromets(k, 2)*awb_nir(ib, 2)*reffd(k, 2)/&
2484 & reff(k, 2)**2)/cons_grav
2485  taucld2 = dp(k)*1.0e3/cons_grav*hydromets(k, 2)*(awb_nir(ib, 1)+&
2486 & awb_nir(ib, 2)/reff(k, 2))
2487  END IF
2488  taucld3d = dp(k)*1.0e3*arb_nir(ib, 1)*hydrometsd(k, 3)/cons_grav
2489  taucld3 = dp(k)*1.0e3/cons_grav*hydromets(k, 3)*arb_nir(ib, 1)
2490  IF (reff(k, 4) .GT. 112.0) THEN
2491  reff_snow = 112.0
2492  reff_snowd = 0.0_8
2493  ELSE
2494  reff_snowd = reffd(k, 4)
2495  reff_snow = reff(k, 4)
2496  END IF
2497  IF (reff_snow .LE. 0.) THEN
2498  taucld4 = 0.
2499  taucld4d = 0.0_8
2500  ELSE
2501  taucld4d = (dp(k)*1.0e3*aib_nir*hydrometsd(k, 4)*reff_snow/&
2502 & cons_grav-dp(k)*1.0e3*hydromets(k, 4)*aib_nir*reff_snowd/&
2503 & cons_grav)/reff_snow**2
2504  taucld4 = dp(k)*1.0e3/cons_grav*hydromets(k, 4)*aib_nir/reff_snow
2505  END IF
2506  IF (ict .NE. 0) THEN
2507 !-----scale cloud optical thickness in each layer from taucld (with
2508 ! cloud amount fcld) to taubeam and taudiff (with cloud amount cc).
2509 ! taubeam is the scaled optical thickness for beam radiation and
2510 ! taudiff is for diffuse radiation (see section 7).
2511 !-----clouds within each of the high, middle, and low clouds are
2512 ! assumed to be maximally overlapped, and the cloud cover (cc)
2513 ! for a group (high, middle, or low) is the maximum cloud cover
2514 ! of all the layers within a group
2515  IF (k .LT. ict) THEN
2516  kk = 1
2517  ELSE IF (k .GE. ict .AND. k .LT. icb) THEN
2518  kk = 2
2519  ELSE
2520  kk = 3
2521  END IF
2522  taucd = taucld1d + taucld2d + taucld3d + taucld4d
2523  tauc = taucld1 + taucld2 + taucld3 + taucld4
2524  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
2525 !-----normalize cloud cover following Eq. (7.8)
2526  IF (cc(kk) .NE. 0.0) THEN
2527  fad = (fcldd(k)*cc(kk)-fcld(k)*ccd(kk))/cc(kk)**2
2528  fa = fcld(k)/cc(kk)
2529  ELSE
2530  fa = 0.0
2531  fad = 0.0_8
2532  END IF
2533  IF (tauc .GT. 32.) THEN
2534  tauc = 32.
2535  taucd = 0.0_8
2536  ELSE
2537  tauc = tauc
2538  END IF
2539  fm = cosz/dm
2540  ftd = taucd/(tauc*log(10.0))/dt
2541  ft = (log10(tauc)-t1)/dt
2542  fad = fad/da
2543  fa = fa/da
2544  im = int(fm + 1.5)
2545  it = int(ft + 1.5)
2546  ia = int(fa + 1.5)
2547  IF (im .LT. 2) THEN
2548  im = 2
2549  ELSE
2550  im = im
2551  END IF
2552  IF (it .LT. 2) THEN
2553  it = 2
2554  ELSE
2555  it = it
2556  END IF
2557  IF (ia .LT. 2) THEN
2558  ia = 2
2559  ELSE
2560  ia = ia
2561  END IF
2562  IF (im .GT. nm - 1) THEN
2563  im = nm - 1
2564  ELSE
2565  im = im
2566  END IF
2567  IF (it .GT. nt - 1) THEN
2568  it = nt - 1
2569  ELSE
2570  it = it
2571  END IF
2572  IF (ia .GT. na - 1) THEN
2573  ia = na - 1
2574  ELSE
2575  ia = ia
2576  END IF
2577  fm = fm - REAL(im - 1)
2578  ft = ft - REAL(it - 1)
2579  fa = fa - REAL(ia - 1)
2580 !-----scale cloud optical thickness for beam radiation following
2581 ! Eq. (7.3).
2582 ! the scaling factor, xai, is a function of the solar zenith
2583 ! angle, optical thickness, and cloud cover.
2584  xai = (-(caib(im-1, it, ia)*(1.-fm))+caib(im+1, it, ia)*(1.+fm))&
2585 & *fm*.5 + caib(im, it, ia)*(1.-fm*fm)
2586  xaid = .5*((caib(im, it-1, ia)*ftd+caib(im, it+1, ia)*ftd)*ft+(-&
2587 & (caib(im, it-1, ia)*(1.-ft))+caib(im, it+1, ia)*(1.+ft))*ftd) &
2588 & + caib(im, it, ia)*(-(ftd*ft)-ft*ftd)
2589  xai = xai + (-(caib(im, it-1, ia)*(1.-ft))+caib(im, it+1, ia)*(&
2590 & 1.+ft))*ft*.5 + caib(im, it, ia)*(1.-ft*ft)
2591  xaid = xaid + .5*((caib(im, it, ia-1)*fad+caib(im, it, ia+1)*fad&
2592 & )*fa+(-(caib(im, it, ia-1)*(1.-fa))+caib(im, it, ia+1)*(1.+fa)&
2593 & )*fad) + caib(im, it, ia)*(-(fad*fa)-fa*fad)
2594  xai = xai + (-(caib(im, it, ia-1)*(1.-fa))+caib(im, it, ia+1)*(&
2595 & 1.+fa))*fa*.5 + caib(im, it, ia)*(1.-fa*fa)
2596  xai = xai - 2.*caib(im, it, ia)
2597  IF (xai .LT. 0.0) THEN
2598  xai = 0.0
2599  xaid = 0.0_8
2600  ELSE
2601  xai = xai
2602  END IF
2603  IF (xai .GT. 1.0) THEN
2604  xai = 1.0
2605  xaid = 0.0_8
2606  ELSE
2607  xai = xai
2608  END IF
2609  taubeamd(k, 1) = taucld1d*xai + taucld1*xaid
2610  taubeam(k, 1) = taucld1*xai
2611  taubeamd(k, 2) = taucld2d*xai + taucld2*xaid
2612  taubeam(k, 2) = taucld2*xai
2613  taubeamd(k, 3) = taucld3d*xai + taucld3*xaid
2614  taubeam(k, 3) = taucld3*xai
2615  taubeamd(k, 4) = taucld4d*xai + taucld4*xaid
2616  taubeam(k, 4) = taucld4*xai
2617 !-----scale cloud optical thickness for diffuse radiation following
2618 ! Eq. (7.4).
2619 ! the scaling factor, xai, is a function of the cloud optical
2620 ! thickness and cover but not the solar zenith angle.
2621  xaid = .5*((caif(it-1, ia)*ftd+caif(it+1, ia)*ftd)*ft+(-(caif(it&
2622 & -1, ia)*(1.-ft))+caif(it+1, ia)*(1.+ft))*ftd) + caif(it, ia)*(&
2623 & -(ftd*ft)-ft*ftd)
2624  xai = (-(caif(it-1, ia)*(1.-ft))+caif(it+1, ia)*(1.+ft))*ft*.5 +&
2625 & caif(it, ia)*(1.-ft*ft)
2626  xaid = xaid + .5*((caif(it, ia-1)*fad+caif(it, ia+1)*fad)*fa+(-(&
2627 & caif(it, ia-1)*(1.-fa))+caif(it, ia+1)*(1.+fa))*fad) + caif(it&
2628 & , ia)*(-(fad*fa)-fa*fad)
2629  xai = xai + (-(caif(it, ia-1)*(1.-fa))+caif(it, ia+1)*(1.+fa))*&
2630 & fa*.5 + caif(it, ia)*(1.-fa*fa)
2631  xai = xai - caif(it, ia)
2632  IF (xai .LT. 0.0) THEN
2633  xai = 0.0
2634  xaid = 0.0_8
2635  ELSE
2636  xai = xai
2637  END IF
2638  IF (xai .GT. 1.0) THEN
2639  xai = 1.0
2640  xaid = 0.0_8
2641  ELSE
2642  xai = xai
2643  END IF
2644  taudiffd(k, 1) = taucld1d*xai + taucld1*xaid
2645  taudiff(k, 1) = taucld1*xai
2646  taudiffd(k, 2) = taucld2d*xai + taucld2*xaid
2647  taudiff(k, 2) = taucld2*xai
2648  taudiffd(k, 3) = taucld3d*xai + taucld3*xaid
2649  taudiff(k, 3) = taucld3*xai
2650  taudiffd(k, 4) = taucld4d*xai + taucld4*xaid
2651  taudiff(k, 4) = taucld4*xai
2652  END IF
2653  ELSE
2654 ! Overlap calculation scaling not needed
2655  taubeamd(k, 1) = taucld1d
2656  taubeam(k, 1) = taucld1
2657  taubeamd(k, 2) = taucld2d
2658  taubeam(k, 2) = taucld2
2659  taubeamd(k, 3) = taucld3d
2660  taubeam(k, 3) = taucld3
2661  taubeamd(k, 4) = taucld4d
2662  taubeam(k, 4) = taucld4
2663  taudiffd(k, 1) = taucld1d
2664  taudiff(k, 1) = taucld1
2665  taudiffd(k, 2) = taucld2d
2666  taudiff(k, 2) = taucld2
2667  taudiffd(k, 3) = taucld3d
2668  taudiff(k, 3) = taucld3
2669  taudiffd(k, 4) = taucld4d
2670  taudiff(k, 4) = taucld4
2671  END IF
2672 !-----compute cloud single scattering albedo and asymmetry factor
2673 ! for a mixture of ice and liquid particles.
2674 ! Eqs.(4.6)-(4.8), (6.2)-(6.4)
2675  ssaclt = 0.99999
2676  asyclt = 1.0
2677  taucd = taucld1d + taucld2d + taucld3d + taucld4d
2678  tauc = taucld1 + taucld2 + taucld3 + taucld4
2679  IF (tauc .GT. 0.02 .AND. fcld(k) .GT. 0.01) THEN
2680  w1d = (-(aia_nir(ib, 3)*reffd(k, 1)*reff(k, 1))-(aia_nir(ib, 2)+&
2681 & aia_nir(ib, 3)*reff(k, 1))*reffd(k, 1))*taucld1 + (1.-(aia_nir(&
2682 & ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff(k, 1))*reff(k, 1)))*&
2683 & taucld1d
2684  w1 = (1.-(aia_nir(ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff(k, 1)&
2685 & )*reff(k, 1)))*taucld1
2686  w2d = (-(awa_nir(ib, 3)*reffd(k, 2)*reff(k, 2))-(awa_nir(ib, 2)+&
2687 & awa_nir(ib, 3)*reff(k, 2))*reffd(k, 2))*taucld2 + (1.-(awa_nir(&
2688 & ib, 1)+(awa_nir(ib, 2)+awa_nir(ib, 3)*reff(k, 2))*reff(k, 2)))*&
2689 & taucld2d
2690  w2 = (1.-(awa_nir(ib, 1)+(awa_nir(ib, 2)+awa_nir(ib, 3)*reff(k, 2)&
2691 & )*reff(k, 2)))*taucld2
2692  w3d = (1.-ara_nir(ib, 1))*taucld3d
2693  w3 = (1.-ara_nir(ib, 1))*taucld3
2694  w4d = (-(aia_nir(ib, 3)*reff_snowd*reff_snow)-(aia_nir(ib, 2)+&
2695 & aia_nir(ib, 3)*reff_snow)*reff_snowd)*taucld4 + (1.-(aia_nir(ib&
2696 & , 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff_snow)*reff_snow))*&
2697 & taucld4d
2698  w4 = (1.-(aia_nir(ib, 1)+(aia_nir(ib, 2)+aia_nir(ib, 3)*reff_snow)&
2699 & *reff_snow))*taucld4
2700  ssacltd = ((w1d+w2d+w3d+w4d)*tauc-(w1+w2+w3+w4)*taucd)/tauc**2
2701  ssaclt = (w1+w2+w3+w4)/tauc
2702  g1d = (aig_nir(ib, 3)*reffd(k, 1)*reff(k, 1)+(aig_nir(ib, 2)+&
2703 & aig_nir(ib, 3)*reff(k, 1))*reffd(k, 1))*w1 + (aig_nir(ib, 1)+(&
2704 & aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 1))*reff(k, 1))*w1d
2705  g1 = (aig_nir(ib, 1)+(aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 1))*&
2706 & reff(k, 1))*w1
2707  g2d = (awg_nir(ib, 3)*reffd(k, 2)*reff(k, 2)+(awg_nir(ib, 2)+&
2708 & awg_nir(ib, 3)*reff(k, 2))*reffd(k, 2))*w2 + (awg_nir(ib, 1)+(&
2709 & awg_nir(ib, 2)+awg_nir(ib, 3)*reff(k, 2))*reff(k, 2))*w2d
2710  g2 = (awg_nir(ib, 1)+(awg_nir(ib, 2)+awg_nir(ib, 3)*reff(k, 2))*&
2711 & reff(k, 2))*w2
2712  g3d = arg_nir(ib, 1)*w3d
2713  g3 = arg_nir(ib, 1)*w3
2714  g4d = (aig_nir(ib, 3)*reffd(k, 4)*reff(k, 4)+(aig_nir(ib, 2)+&
2715 & aig_nir(ib, 3)*reff(k, 4))*reffd(k, 4))*w4 + (aig_nir(ib, 1)+(&
2716 & aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 4))*reff(k, 4))*w4d
2717  g4 = (aig_nir(ib, 1)+(aig_nir(ib, 2)+aig_nir(ib, 3)*reff(k, 4))*&
2718 & reff(k, 4))*w4
2719  IF (w1 + w2 + w3 + w4 .NE. 0.0) THEN
2720  asycltd = ((g1d+g2d+g3d+g4d)*(w1+w2+w3+w4)-(g1+g2+g3+g4)*(w1d+&
2721 & w2d+w3d+w4d))/(w1+w2+w3+w4)**2
2722  asyclt = (g1+g2+g3+g4)/(w1+w2+w3+w4)
2723  ELSE
2724  asycltd = 0.0_8
2725  END IF
2726  ELSE
2727  ssacltd = 0.0_8
2728  asycltd = 0.0_8
2729  END IF
2730  ssacld(k) = ssacltd
2731  ssacl(k) = ssaclt
2732  asycld(k) = asycltd
2733  asycl(k) = asyclt
2734  END DO
2735  RETURN
2736 END SUBROUTINE getnirtau1_d
2737 
2738 end module sorad_tl
sorad_tl
Definition: sorad_tl.F90:1
sorad_tl::sorad_d
subroutine, public sorad_d(m, np, nb, cosz_dev, pl_dev, ta_dev, ta_devd, wa_dev, wa_devd, oa_dev, oa_devd, co2, cwc_dev, cwc_devd, fcld_dev, fcld_devd, ict, icb, reff_dev, reff_devd, hk_uv, hk_ir, taua_dev, taua_devd, ssaa_dev, ssaa_devd, asya_dev, asya_devd, rsuvbm_dev, rsuvdf_dev, rsirbm_dev, rsirdf_dev, flx_dev, flx_devd, 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_tl.F90:20
qg_constants::u1
real(kind=kind_real), parameter u1
Definition: qg_constants.f90:35
sorad_tl::getnirtau1_d
subroutine getnirtau1_d(ib, nlevs, cosz, dp, fcld, fcldd, reff, reffd, hydromets, hydrometsd, ict, icb, taubeam, taubeamd, taudiff, taudiffd, asycl, asycld, ssacl, ssacld, 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_tl.F90:2358
max
#define max(a, b)
Definition: mosaic_util.h:33
min
#define min(a, b)
Definition: mosaic_util.h:32
sorad_tl::deledd_d
subroutine deledd_d(tau1, tau1d, ssc1, ssc1d, g01, g01d, cza1, rr1, rr1d, tt1, tt1d, td1, td1d)
Definition: sorad_tl.F90:1808
sorad_tl::getvistau1_d
subroutine getvistau1_d(nlevs, cosz, dp, fcld, fcldd, reff, reffd, hydromets, hydrometsd, ict, icb, taubeam, taubeamd, taudiff, taudiffd, asycl, asycld, 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_tl.F90:1982