Radiation Scheme in CCPP
radsw_main.f
Go to the documentation of this file.
1 ! ============================================================== !!!!!
2 ! sw-rrtm3 radiation package description !!!!!
3 ! ============================================================== !!!!!
4 ! !
5 ! this package includes ncep's modifications of the rrtm-sw radiation !
6 ! code from aer inc. !
7 ! !
8 ! the sw-rrtm3 package includes these parts: !
9 ! !
10 ! 'radsw_rrtm3_param.f' !
11 ! 'radsw_rrtm3_datatb.f' !
12 ! 'radsw_rrtm3_main.f' !
13 ! !
14 ! the 'radsw_rrtm3_param.f' contains: !
15 ! !
16 ! 'module_radsw_parameters' -- band parameters set up !
17 ! !
18 ! the 'radsw_rrtm3_datatb.f' contains: !
19 ! !
20 ! 'module_radsw_ref' -- reference temperature and pressure !
21 ! 'module_radsw_cldprtb' -- cloud property coefficients table !
22 ! 'module_radsw_sflux' -- spectral distribution of solar flux !
23 ! 'module_radsw_kgbnn' -- absorption coeffients for 14 !
24 ! bands, where nn = 16-29 !
25 ! !
26 ! the 'radsw_rrtm3_main.f' contains: !
27 ! !
28 ! 'module_radsw_main' -- main sw radiation transfer !
29 ! !
30 ! in the main module 'module_radsw_main' there are only two !
31 ! externally callable subroutines: !
32 ! !
33 ! 'swrad' -- main sw radiation routine !
34 ! inputs: !
35 ! (plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, !
36 ! clouds,icseed,aerosols,sfcalb, !
37 ! cosz,solcon,NDAY,idxday, !
38 ! npts, nlay, nlp1, lprnt, !
39 ! outputs: !
40 ! hswc,topflx,sfcflx, !
41 !! optional outputs: !
42 ! HSW0,HSWB,FLXPRF,FDNCMP) !
43 ! ) !
44 ! !
45 ! 'rswinit' -- initialization routine !
46 ! inputs: !
47 ! ( me ) !
48 ! outputs: !
49 ! (none) !
50 ! !
51 ! all the sw radiation subprograms become contained subprograms !
52 ! in module 'module_radsw_main' and many of them are not directly !
53 ! accessable from places outside the module. !
54 ! !
55 ! derived data type constructs used: !
56 ! !
57 ! 1. radiation flux at toa: (from module 'module_radsw_parameters') !
58 ! topfsw_type - derived data type for toa rad fluxes !
59 ! upfxc total sky upward flux at toa !
60 ! dnfxc total sky downward flux at toa !
61 ! upfx0 clear sky upward flux at toa !
62 ! !
63 ! 2. radiation flux at sfc: (from module 'module_radsw_parameters') !
64 ! sfcfsw_type - derived data type for sfc rad fluxes !
65 ! upfxc total sky upward flux at sfc !
66 ! dnfxc total sky downward flux at sfc !
67 ! upfx0 clear sky upward flux at sfc !
68 ! dnfx0 clear sky downward flux at sfc !
69 ! !
70 ! 3. radiation flux profiles(from module 'module_radsw_parameters') !
71 ! profsw_type - derived data type for rad vertical prof !
72 ! upfxc level upward flux for total sky !
73 ! dnfxc level downward flux for total sky !
74 ! upfx0 level upward flux for clear sky !
75 ! dnfx0 level downward flux for clear sky !
76 ! !
77 ! 4. surface component fluxes(from module 'module_radsw_parameters' !
78 ! cmpfsw_type - derived data type for component sfc flux !
79 ! uvbfc total sky downward uv-b flux at sfc !
80 ! uvbf0 clear sky downward uv-b flux at sfc !
81 ! nirbm surface downward nir direct beam flux !
82 ! nirdf surface downward nir diffused flux !
83 ! visbm surface downward uv+vis direct beam flx !
84 ! visdf surface downward uv+vis diffused flux !
85 ! !
86 ! external modules referenced: !
87 ! !
88 ! 'module physparam' !
89 ! 'module physcons' !
90 ! 'mersenne_twister' !
91 ! !
92 ! compilation sequence is: !
93 ! !
94 ! 'radsw_rrtm3_param.f' !
95 ! 'radsw_rrtm3_datatb.f' !
96 ! 'radsw_rrtm3_main.f' !
97 ! !
98 ! and all should be put in front of routines that use sw modules !
99 ! !
100 !==========================================================================!
101 ! !
102 ! the original program declarations: !
103 ! !
104 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
105 ! !
106 ! Copyright 2002-2007, Atmospheric & Environmental Research, Inc. (AER). !
107 ! This software may be used, copied, or redistributed as long as it is !
108 ! not sold and this copyright notice is reproduced on each copy made. !
109 ! This model is provided as is without any express or implied warranties. !
110 ! (http://www.rtweb.aer.com/) !
111 ! !
112 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
113 ! !
114 ! ************************************************************************ !
115 ! !
116 ! rrtmg_sw !
117 ! !
118 ! !
119 ! a rapid radiative transfer model !
120 ! for the solar spectral region !
121 ! atmospheric and environmental research, inc. !
122 ! 131 hartwell avenue !
123 ! lexington, ma 02421 !
124 ! !
125 ! eli j. mlawer !
126 ! jennifer s. delamere !
127 ! michael j. iacono !
128 ! shepard a. clough !
129 ! !
130 ! !
131 ! email: miacono@aer.com !
132 ! email: emlawer@aer.com !
133 ! email: jdelamer@aer.com !
134 ! !
135 ! the authors wish to acknowledge the contributions of the !
136 ! following people: steven j. taubman, patrick d. brown, !
137 ! ronald e. farren, luke chen, robert bergstrom. !
138 ! !
139 ! ************************************************************************ !
140 ! !
141 ! references: !
142 ! (rrtm_sw/rrtmg_sw): !
143 ! clough, s.a., m.w. shephard, e.j. mlawer, j.s. delamere, !
144 ! m.j. iacono, k. cady-pereira, s. boukabara, and p.d. brown: !
145 ! atmospheric radiative transfer modeling: a summary of the aer !
146 ! codes, j. quant. spectrosc. radiat. transfer, 91, 233-244, 2005. !
147 ! !
148 ! (mcica): !
149 ! pincus, r., h. w. barker, and j.-j. morcrette: a fast, flexible, !
150 ! approximation technique for computing radiative transfer in !
151 ! inhomogeneous cloud fields, j. geophys. res., 108(d13), 4376, !
152 ! doi:10.1029/2002jd003322, 2003. !
153 ! !
154 ! ************************************************************************ !
155 ! !
156 ! aer's revision history: !
157 ! this version of rrtmg_sw has been modified from rrtm_sw to use a !
158 ! reduced set of g-point intervals and a two-stream model for !
159 ! application to gcms. !
160 ! !
161 ! -- original version (derived from rrtm_sw) !
162 ! 2002: aer. inc. !
163 ! -- conversion to f90 formatting; addition of 2-stream radiative transfer!
164 ! feb 2003: j.-j. morcrette, ecmwf !
165 ! -- additional modifications for gcm application !
166 ! aug 2003: m. j. iacono, aer inc. !
167 ! -- total number of g-points reduced from 224 to 112. original !
168 ! set of 224 can be restored by exchanging code in module parrrsw.f90 !
169 ! and in file rrtmg_sw_init.f90. !
170 ! apr 2004: m. j. iacono, aer, inc. !
171 ! -- modifications to include output for direct and diffuse !
172 ! downward fluxes. there are output as "true" fluxes without !
173 ! any delta scaling applied. code can be commented to exclude !
174 ! this calculation in source file rrtmg_sw_spcvrt.f90. !
175 ! jan 2005: e. j. mlawer, m. j. iacono, aer, inc. !
176 ! -- revised to add mcica capability. !
177 ! nov 2005: m. j. iacono, aer, inc. !
178 ! -- reformatted for consistency with rrtmg_lw. !
179 ! feb 2007: m. j. iacono, aer, inc. !
180 ! -- modifications to formatting to use assumed-shape arrays. !
181 ! aug 2007: m. j. iacono, aer, inc. !
182 ! !
183 ! ************************************************************************ !
184 ! !
185 ! ncep modifications history log: !
186 ! !
187 ! sep 2003, yu-tai hou -- received aer's rrtm-sw gcm version !
188 ! code (v224) !
189 ! nov 2003, yu-tai hou -- corrected errors in direct/diffuse !
190 ! surface alabedo components. !
191 ! jan 2004, yu-tai hou -- modified code into standard modular!
192 ! f9x code for ncep models. the original three cloud !
193 ! control flags are simplified into two: iflagliq and !
194 ! iflagice. combined the org subr sw_224 and setcoef !
195 ! into radsw (the main program); put all kgb##together !
196 ! and reformat into a separated data module; combine !
197 ! reftra and vrtqdr as swflux; optimized taumol and all !
198 ! taubgs to form a contained subroutines. !
199 ! jun 2004, yu-tai hou -- modified code based on aer's faster!
200 ! version rrtmg_sw (v2.0) with 112 g-points. !
201 ! mar 2005, yu-tai hou -- modified to aer v2.3, correct cloud!
202 ! scaling error, total sky properties are delta scaled !
203 ! after combining clear and cloudy parts. the testing !
204 ! criterion of ssa is saved before scaling. added cloud !
205 ! layer rain and snow contributions. all cloud water !
206 ! partical contents are treated the same way as other !
207 ! atmos particles. !
208 ! apr 2005, yu-tai hou -- modified on module structures (this!
209 ! version of code was given back to aer in jun 2006) !
210 ! nov 2006, yu-tai hou -- modified code to include the !
211 ! generallized aerosol optical property scheme for gcms.!
212 ! apr 2007, yu-tai hou -- added spectral band heating as an !
213 ! optional output to support the 500km model's upper !
214 ! stratospheric radiation calculations. restructure !
215 ! optional outputs for easy access by different models. !
216 ! oct 2008, yu-tai hou -- modified to include new features !
217 ! from aer's newer release v3.5-v3.61, including mcica !
218 ! sub-grid cloud option and true direct/diffuse fluxes !
219 ! without delta scaling. added rain/snow opt properties !
220 ! support to cloudy sky calculations. simplified and !
221 ! unified sw and lw sub-column cloud subroutines into !
222 ! one module by using optional parameters. !
223 ! mar 2009, yu-tai hou -- replaced the original random number!
224 ! generator coming with the original code with ncep w3 !
225 ! library to simplify the program and moved sub-column !
226 ! cloud subroutines inside the main module. added !
227 ! option of user provided permutation seeds that could !
228 ! be randomly generated from forecast time stamp. !
229 ! mar 2009, yu-tai hou -- replaced random number generator !
230 ! programs coming from the original code with the ncep !
231 ! w3 library to simplify the program and moved sub-col !
232 ! cloud subroutines inside the main module. added !
233 ! option of user provided permutation seeds that could !
234 ! be randomly generated from forecast time stamp. !
235 ! nov 2009, yu-tai hou -- updated to aer v3.7-v3.8 version. !
236 ! notice the input cloud ice/liquid are assumed as !
237 ! in-cloud quantities, not grid average quantities. !
238 ! aug 2010, yu-tai hou -- uptimized code to improve efficiency
239 ! splited subroutine spcvrt into two subs, spcvrc and !
240 ! spcvrm, to handling non-mcica and mcica type of calls.!
241 ! apr 2012, b. ferrier and y. hou -- added conversion factor to fu's!
242 ! cloud-snow optical property scheme. !
243 ! jul 2012, s. moorthi and Y. hou -- eliminated the pointer array !
244 ! in subr 'spcvrt' for multi-threading issue running !
245 ! under intel's fortran compiler. !
246 ! nov 2012, yu-tai hou -- modified control parameters thru !
247 ! module 'physparam'. !
248 ! jun 2013, yu-tai hou -- moving band 9 surface treatment !
249 ! back as in the rrtm2 version, spliting surface flux !
250 ! into two spectral regions (vis & nir), instead of !
251 ! designated it in nir region only. !
252 ! !
253 !!!!! ============================================================== !!!!!
254 !!!!! end descriptions !!!!!
255 !!!!! ============================================================== !!!!!
256 
259 !========================================!
261 !........................................!
262 !
263  use physparam, only : iswrate, iswrgas, iswcliq, iswcice, &
264  & isubcsw, icldflg, iovrsw, ivflip, &
265  & iswmode, kind_phys
266  use physcons, only : con_g, con_cp, con_avgd, con_amd, &
267  & con_amw, con_amo3
268 
270  use mersenne_twister, only : random_setseed, random_number, &
271  & random_stat
272  use module_radsw_ref, only : preflog, tref
274 !
275  implicit none
276 !
277  private
278 !
279 ! --- version tag and last revision date
280  character(40), parameter :: &
281  & VTAGSW='NCEP SW v5.1 Nov 2012 -RRTMG-SW v3.8 '
282 ! & VTAGSW='NCEP SW v5.0 Aug 2012 -RRTMG-SW v3.8 '
283 ! & VTAGSW='RRTMG-SW v3.8 Nov 2009'
284 ! & VTAGSW='RRTMG-SW v3.7 Nov 2009'
285 ! & VTAGSW='RRTMG-SW v3.61 Oct 2008'
286 ! & VTAGSW='RRTMG-SW v3.5 Oct 2008'
287 ! & VTAGSW='RRTM-SW 112v2.3 Apr 2007'
288 ! & VTAGSW='RRTM-SW 112v2.3 Mar 2005'
289 ! & VTAGSW='RRTM-SW 112v2.0 Jul 2004'
290 
291 ! --- constant values
292  real (kind=kind_phys), parameter :: eps = 1.0e-6
293  real (kind=kind_phys), parameter :: oneminus= 1.0 - eps
294  real (kind=kind_phys), parameter :: bpade = 1.0/0.278 ! pade approx constant
295  real (kind=kind_phys), parameter :: stpfac = 296.0/1013.0
296  real (kind=kind_phys), parameter :: ftiny = 1.0e-12
297  real (kind=kind_phys), parameter :: s0 = 1368.22 ! internal solar const
298  ! adj through input
299  real (kind=kind_phys), parameter :: f_zero = 0.0
300  real (kind=kind_phys), parameter :: f_one = 1.0
301 
302 ! --- atomic weights for conversion from mass to volume mixing ratios
303  real (kind=kind_phys), parameter :: amdw = con_amd/con_amw
304  real (kind=kind_phys), parameter :: amdo3 = con_amd/con_amo3
305 
306 ! --- band indices
307  integer, dimension(nblow:nbhgh) :: nspa, nspb, idxebc, idxsfc
308 
309  data nspa(:) / 9, 9, 9, 9, 1, 9, 9, 1, 9, 1, 0, 1, 9, 1 /
310  data nspb(:) / 1, 5, 1, 1, 1, 5, 1, 0, 1, 0, 0, 1, 5, 1 /
311 
312 ! data idxsfc(:) / 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1 / ! band index for sfc flux
313  data idxsfc(:) / 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1 / ! band index for sfc flux
314  data idxebc(:) / 5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 1, 1, 1, 5 / ! band index for cld prop
315 
316 ! --- band wavenumber intervals
317 ! real (kind=kind_phys), dimension(nblow:nbhgh):: wavenum1,wavenum2
318 ! data wavenum1(:) / &
319 ! & 2600.0, 3250.0, 4000.0, 4650.0, 5150.0, 6150.0, 7700.0, &
320 ! & 8050.0,12850.0,16000.0,22650.0,29000.0,38000.0, 820.0 /
321 ! data wavenum2(:) / &
322 ! 3250.0, 4000.0, 4650.0, 5150.0, 6150.0, 7700.0, 8050.0, &
323 ! & 12850.0,16000.0,22650.0,29000.0,38000.0,50000.0, 2600.0 /
324 ! real (kind=kind_phys), dimension(nblow:nbhgh) :: delwave
325 ! data delwave(:) / &
326 ! & 650.0, 750.0, 650.0, 500.0, 1000.0, 1550.0, 350.0, &
327 ! & 4800.0, 3150.0, 6650.0, 6350.0, 9000.0,12000.0, 1780.0 /
328 
329  integer, parameter :: nuvb = 27 !uv-b band index
330 
331 !! --- logical flags for optional output fields
332 
333  logical :: lhswb = .false.
334  logical :: lhsw0 = .false.
335  logical :: lflxprf= .false.
336  logical :: lfdncmp= .false.
337 
338 ! --- those data will be set up only once by "rswinit"
339 
340  real (kind=kind_phys) :: exp_tbl(0:ntbmx)
341 
342 ! ... heatfac is the factor for heating rates
343 ! (in k/day, or k/sec set by subroutine 'rswinit')
344 
345  real (kind=kind_phys) :: heatfac
346 
347 ! --- the following variables are used for sub-column cloud scheme
348 
349  integer, parameter :: ipsdsw0 = 1 ! initial permutation seed
350 
351 ! --- public accessable subprograms
352 
353  public swrad, rswinit
354 
355 
356 ! =================
357  contains
358 ! =================
359 
437 !-----------------------------------
438  subroutine swrad &
439 !...................................
441 ! --- inputs:
442  & ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, &
443  & clouds,icseed,aerosols,sfcalb, &
444  & cosz,solcon,nday,idxday, &
445  & npts, nlay, nlp1, lprnt, &
446 ! --- outputs:
447  & hswc,topflx,sfcflx &
448 !! --- optional:
449  &, hsw0,hswb,flxprf,fdncmp &
450  & )
451 
452 ! ==================== defination of variables ==================== !
453 ! !
454 ! input variables: !
455 ! plyr (npts,nlay) : model layer mean pressure in mb !
456 ! plvl (npts,nlp1) : model level pressure in mb !
457 ! tlyr (npts,nlay) : model layer mean temperature in k !
458 ! tlvl (npts,nlp1) : model level temperature in k (not in use) !
459 ! qlyr (npts,nlay) : layer specific humidity in gm/gm *see inside !
460 ! olyr (npts,nlay) : layer ozone concentration in gm/gm !
461 ! gasvmr(npts,nlay,:): atmospheric constent gases: !
462 ! (check module_radiation_gases for definition) !
463 ! gasvmr(:,:,1) - co2 volume mixing ratio !
464 ! gasvmr(:,:,2) - n2o volume mixing ratio !
465 ! gasvmr(:,:,3) - ch4 volume mixing ratio !
466 ! gasvmr(:,:,4) - o2 volume mixing ratio !
467 ! gasvmr(:,:,5) - co volume mixing ratio (not used) !
468 ! gasvmr(:,:,6) - cfc11 volume mixing ratio (not used) !
469 ! gasvmr(:,:,7) - cfc12 volume mixing ratio (not used) !
470 ! gasvmr(:,:,8) - cfc22 volume mixing ratio (not used) !
471 ! gasvmr(:,:,9) - ccl4 volume mixing ratio (not used) !
472 ! clouds(npts,nlay,:): cloud profile !
473 ! (check module_radiation_clouds for definition) !
474 ! --- for iswcliq > 0 --- !
475 ! clouds(:,:,1) - layer total cloud fraction !
476 ! clouds(:,:,2) - layer in-cloud liq water path (g/m**2) !
477 ! clouds(:,:,3) - mean eff radius for liq cloud (micron) !
478 ! clouds(:,:,4) - layer in-cloud ice water path (g/m**2) !
479 ! clouds(:,:,5) - mean eff radius for ice cloud (micron) !
480 ! clouds(:,:,6) - layer rain drop water path (g/m**2) !
481 ! clouds(:,:,7) - mean eff radius for rain drop (micron) !
482 ! clouds(:,:,8) - layer snow flake water path (g/m**2) !
483 ! clouds(:,:,9) - mean eff radius for snow flake (micron) !
484 ! --- for iswcliq = 0 --- !
485 ! clouds(:,:,1) - layer total cloud fraction !
486 ! clouds(:,:,2) - layer cloud optical depth !
487 ! clouds(:,:,3) - layer cloud single scattering albedo !
488 ! clouds(:,:,4) - layer cloud asymmetry factor !
489 ! icseed(npts) : auxiliary special cloud related array !
490 ! when module variable isubcsw=2, it provides !
491 ! permutation seed for each column profile that !
492 ! are used for generating random numbers. !
493 ! when isubcsw /=2, it will not be used. !
494 ! aerosols(npts,nlay,nbdsw,:) : aerosol optical properties !
495 ! (check module_radiation_aerosols for definition) !
496 ! (:,:,:,1) - optical depth !
497 ! (:,:,:,2) - single scattering albedo !
498 ! (:,:,:,3) - asymmetry parameter !
499 ! sfcalb(npts, : ) : surface albedo in fraction !
500 ! (check module_radiation_surface for definition) !
501 ! ( :, 1 ) - near ir direct beam albedo !
502 ! ( :, 2 ) - near ir diffused albedo !
503 ! ( :, 3 ) - uv+vis direct beam albedo !
504 ! ( :, 4 ) - uv+vis diffused albedo !
505 ! cosz (npts) : cosine of solar zenith angle !
506 ! solcon : solar constant (w/m**2) !
507 ! NDAY : num of daytime points !
508 ! idxday(npts) : index array for daytime points !
509 ! npts : number of horizontal points !
510 ! nlay,nlp1 : vertical layer/lavel numbers !
511 ! lprnt : logical check print flag !
512 ! !
513 ! output variables: !
514 ! hswc (npts,nlay): total sky heating rates (k/sec or k/day) !
515 ! topflx(npts) : radiation fluxes at toa (w/m**2), components: !
516 ! (check module_radsw_parameters for definition) !
517 ! upfxc - total sky upward flux at toa !
518 ! dnflx - total sky downward flux at toa !
519 ! upfx0 - clear sky upward flux at toa !
520 ! sfcflx(npts) : radiation fluxes at sfc (w/m**2), components: !
521 ! (check module_radsw_parameters for definition) !
522 ! upfxc - total sky upward flux at sfc !
523 ! dnfxc - total sky downward flux at sfc !
524 ! upfx0 - clear sky upward flux at sfc !
525 ! dnfx0 - clear sky downward flux at sfc !
526 ! !
527 !!optional outputs variables: !
528 ! hswb(npts,nlay,nbdsw): spectral band total sky heating rates !
529 ! hsw0 (npts,nlay): clear sky heating rates (k/sec or k/day) !
530 ! flxprf(npts,nlp1): level radiation fluxes (w/m**2), components: !
531 ! (check module_radsw_parameters for definition) !
532 ! dnfxc - total sky downward flux at interface !
533 ! upfxc - total sky upward flux at interface !
534 ! dnfx0 - clear sky downward flux at interface !
535 ! upfx0 - clear sky upward flux at interface !
536 ! fdncmp(npts) : component surface downward fluxes (w/m**2): !
537 ! (check module_radsw_parameters for definition) !
538 ! uvbfc - total sky downward uv-b flux at sfc !
539 ! uvbf0 - clear sky downward uv-b flux at sfc !
540 ! nirbm - downward surface nir direct beam flux !
541 ! nirdf - downward surface nir diffused flux !
542 ! visbm - downward surface uv+vis direct beam flux !
543 ! visdf - downward surface uv+vis diffused flux !
544 ! !
545 ! external module variables: (in physparam) !
546 ! iswrgas - control flag for rare gases (ch4,n2o,o2, etc.) !
547 ! =0: do not include rare gases !
548 ! >0: include all rare gases !
549 ! iswcliq - control flag for liq-cloud optical properties !
550 ! =0: input cloud optical depth, fixed ssa, asy !
551 ! =1: use hu and stamnes(1993) method for liq cld !
552 ! =2: not used !
553 ! iswcice - control flag for ice-cloud optical properties !
554 ! *** if iswcliq==0, iswcice is ignored !
555 ! =1: use ebert and curry (1992) scheme for ice clouds !
556 ! =2: use streamer v3.0 (2001) method for ice clouds !
557 ! =3: use fu's method (1996) for ice clouds !
558 ! iswmode - control flag for 2-stream transfer scheme !
559 ! =1; delta-eddington (joseph et al., 1976) !
560 ! =2: pifm (zdunkowski et al., 1980) !
561 ! =3: discrete ordinates (liou, 1973) !
562 ! isubcsw - sub-column cloud approximation control flag !
563 ! =0: no sub-col cld treatment, use grid-mean cld quantities !
564 ! =1: mcica sub-col, prescribed seeds to get random numbers !
565 ! =2: mcica sub-col, providing array icseed for random numbers!
566 ! iovrsw - cloud overlapping control flag !
567 ! =0: random overlapping clouds !
568 ! =1: maximum/random overlapping clouds !
569 ! =2: maximum overlap cloud !
570 ! ivflip - control flg for direction of vertical index !
571 ! =0: index from toa to surface !
572 ! =1: index from surface to toa !
573 ! !
574 ! module parameters, control variables: !
575 ! nblow,nbhgh - lower and upper limits of spectral bands !
576 ! maxgas - maximum number of absorbing gaseous !
577 ! ngptsw - total number of g-point subintervals !
578 ! ng## - number of g-points in band (##=16-29) !
579 ! ngb(ngptsw) - band indices for each g-point !
580 ! bpade - pade approximation constant (1/0.278) !
581 ! nspa,nspb(nblow:nbhgh) !
582 ! - number of lower/upper ref atm's per band !
583 ! ipsdsw0 - permutation seed for mcica sub-col clds !
584 ! !
585 ! major local variables: !
586 ! pavel (nlay) - layer pressures (mb) !
587 ! delp (nlay) - layer pressure thickness (mb) !
588 ! tavel (nlay) - layer temperatures (k) !
589 ! coldry (nlay) - dry air column amount !
590 ! (1.e-20*molecules/cm**2) !
591 ! cldfrc (nlay) - layer cloud fraction (norm by tot cld) !
592 ! cldfmc (nlay,ngptsw) - layer cloud fraction for g-point !
593 ! taucw (nlay,nbdsw) - cloud optical depth !
594 ! ssacw (nlay,nbdsw) - cloud single scattering albedo (weighted) !
595 ! asycw (nlay,nbdsw) - cloud asymmetry factor (weighted) !
596 ! tauaer (nlay,nbdsw) - aerosol optical depths !
597 ! ssaaer (nlay,nbdsw) - aerosol single scattering albedo !
598 ! asyaer (nlay,nbdsw) - aerosol asymmetry factor !
599 ! colamt (nlay,maxgas) - column amounts of absorbing gases !
600 ! 1 to maxgas are for h2o, co2, o3, n2o, !
601 ! ch4, o2, co, respectively (mol/cm**2) !
602 ! facij (nlay) - indicator of interpolation factors !
603 ! =0/1: indicate lower/higher temp & height !
604 ! selffac(nlay) - scale factor for self-continuum, equals !
605 ! (w.v. density)/(atm density at 296K,1013 mb) !
606 ! selffrac(nlay) - factor for temp interpolation of ref !
607 ! self-continuum data !
608 ! indself(nlay) - index of the lower two appropriate ref !
609 ! temp for the self-continuum interpolation !
610 ! forfac (nlay) - scale factor for w.v. foreign-continuum !
611 ! forfrac(nlay) - factor for temp interpolation of ref !
612 ! w.v. foreign-continuum data !
613 ! indfor (nlay) - index of the lower two appropriate ref !
614 ! temp for the foreign-continuum interp !
615 ! laytrop - layer at which switch is made from one !
616 ! combination of key species to another !
617 ! jp(nlay),jt(nlay),jt1(nlay) !
618 ! - lookup table indexes !
619 ! flxucb(nlp1,nbdsw) - spectral bnd total-sky upward flx (w/m2) !
620 ! flxdcb(nlp1,nbdsw) - spectral bnd total-sky downward flx (w/m2)!
621 ! flxu0b(nlp1,nbdsw) - spectral bnd clear-sky upward flx (w/m2) !
622 ! flxd0b(nlp1,nbdsw) - spectral b d clear-sky downward flx (w/m2)!
623 ! !
624 ! !
625 ! ===================== end of definitions ==================== !
626 
627 ! --- inputs:
628  integer, intent(in) :: npts, nlay, nlp1, NDAY
629 
630  integer, dimension(:), intent(in) :: idxday, icseed
631 
632  logical, intent(in) :: lprnt
633 
634  real (kind=kind_phys), dimension(npts,nlp1), intent(in) :: &
635  & plvl, tlvl
636  real (kind=kind_phys), dimension(npts,nlay), intent(in) :: &
637  & plyr, tlyr, qlyr, olyr
638  real (kind=kind_phys), dimension(npts,4), intent(in) :: sfcalb
639 
640  real (kind=kind_phys), dimension(npts,nlay,9),intent(in):: gasvmr
641  real (kind=kind_phys), dimension(npts,nlay,9),intent(in):: clouds
642  real (kind=kind_phys), dimension(npts,nlay,nbdsw,3),intent(in):: &
643  & aerosols
644 
645  real (kind=kind_phys), intent(in) :: cosz(npts), solcon
646 
647 ! --- outputs:
648  real (kind=kind_phys), dimension(npts,nlay), intent(out) :: hswc
649 
650  type(topfsw_type), dimension(npts), intent(out) :: topflx
651  type(sfcfsw_type), dimension(npts), intent(out) :: sfcflx
652 
653 !! --- optional outputs:
654  real (kind=kind_phys), dimension(npts,nlay,nbdsw), optional, &
655  & intent(out) :: hswb
656 
657  real (kind=kind_phys), dimension(npts,nlay), optional, &
658  & intent(out) :: hsw0
659  type (profsw_type), dimension(npts,nlp1), optional, &
660  & intent(out) :: flxprf
661  type (cmpfsw_type), dimension(npts), optional, &
662  & intent(out) :: fdncmp
663 
664 ! --- locals:
665  real (kind=kind_phys), dimension(nlay,ngptsw) :: cldfmc, &
666  & taug, taur
667  real (kind=kind_phys), dimension(nlp1,nbdsw):: fxupc, fxdnc, &
668  & fxup0, fxdn0
669 
670  real (kind=kind_phys), dimension(nlay,nbdsw) :: &
671  & tauae, ssaae, asyae, taucw, ssacw, asycw
672 
673  real (kind=kind_phys), dimension(ngptsw) :: sfluxzen
674 
675  real (kind=kind_phys), dimension(nlay) :: cldfrc, delp, &
676  & pavel, tavel, coldry, colmol, h2ovmr, o3vmr, temcol, &
677  & cliqp, reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, &
678  & cfrac, fac00, fac01, fac10, fac11, forfac, forfrac, &
679  & selffac, selffrac, rfdelp
680 
681  real (kind=kind_phys), dimension(nlp1) :: fnet, flxdc, flxuc, &
682  & flxd0, flxu0
683 
684  real (kind=kind_phys), dimension(2) :: albbm, albdf, sfbmc, &
685  & sfbm0, sfdfc, sfdf0
686 
687  real (kind=kind_phys) :: cosz1, sntz1, tem0, tem1, tem2, s0fac, &
688  & ssolar, zcf0, zcf1, ftoau0, ftoauc, ftoadc, &
689  & fsfcu0, fsfcuc, fsfcd0, fsfcdc, suvbfc, suvbf0
690 
691 ! --- column amount of absorbing gases:
692 ! (:,m) m = 1-h2o, 2-co2, 3-o3, 4-n2o, 5-ch4, 6-o2, 7-co
693  real (kind=kind_phys) :: colamt(nlay,maxgas)
694 
695  integer, dimension(npts) :: ipseed
696  integer, dimension(nlay) :: indfor, indself, jp, jt, jt1
697 
698  integer :: i, ib, ipt, j1, k, kk, laytrop, mb
699 !
700 !===> ... begin here
701 !
702 
703  lhswb = present ( hswb )
704  lhsw0 = present ( hsw0 )
705  lflxprf= present ( flxprf )
706  lfdncmp= present ( fdncmp )
707 
712 ! --- ... compute solar constant adjustment factor according to solcon.
713 ! *** s0, the solar constant at toa in w/m**2, is hard-coded with
714 ! each spectra band, the total flux is about 1368.22 w/m**2.
715 
716  s0fac = solcon / s0
717 
719 ! --- ... initial output arrays
720 
721  hswc(:,:) = f_zero
722  topflx = topfsw_type( f_zero, f_zero, f_zero )
723  sfcflx = sfcfsw_type( f_zero, f_zero, f_zero, f_zero )
724 
725 !! --- ... initial optional outputs
726  if ( lflxprf ) then
727  flxprf = profsw_type( f_zero, f_zero, f_zero, f_zero )
728  endif
729 
730  if ( lfdncmp ) then
732  endif
733 
734  if ( lhsw0 ) then
735  hsw0(:,:) = f_zero
736  endif
737 
738  if ( lhswb ) then
739  hswb(:,:,:) = f_zero
740  endif
741 
743 ! --- ... change random number seed value for each radiation invocation
744 
745  if ( isubcsw == 1 ) then ! advance prescribed permutation seed
746  do i = 1, npts
747  ipseed(i) = ipsdsw0 + i
748  enddo
749  elseif ( isubcsw == 2 ) then ! use input array of permutaion seeds
750  do i = 1, npts
751  ipseed(i) = icseed(i)
752  enddo
753  endif
754 
755  if ( lprnt ) then
756  write(0,*)' In radsw, isubcsw, ipsdsw0,ipseed =', &
757  & isubcsw, ipsdsw0, ipseed
758  endif
759 
760 ! --- ... loop over each daytime grid point
761 
762  lab_do_ipt : do ipt = 1, nday
763 
764  j1 = idxday(ipt)
765 
766  cosz1 = cosz(j1)
767  sntz1 = f_one / cosz(j1)
768  ssolar = s0fac * cosz(j1)
769 
771 ! --- ... surface albedo: bm,df - dir,dif; 1,2 - nir,uvv
772  albbm(1) = sfcalb(j1,1)
773  albdf(1) = sfcalb(j1,2)
774  albbm(2) = sfcalb(j1,3)
775  albdf(2) = sfcalb(j1,4)
776 
783 
784 ! --- ... prepare atmospheric profile for use in rrtm
785 ! the vertical index of internal array is from surface to top
786 
787  if (ivflip == 0) then ! input from toa to sfc
788 
789  tem1 = 100.0 * con_g
790  tem2 = 1.0e-20 * 1.0e3 * con_avgd
791 
792  do k = 1, nlay
793  kk = nlp1 - k
794  pavel(k) = plyr(j1,kk)
795  tavel(k) = tlyr(j1,kk)
796  delp(k) = plvl(j1,kk+1) - plvl(j1,kk)
797 
798 ! --- ... set absorber amount
799 !test use
800 ! h2ovmr(k)= max(f_zero,qlyr(j1,kk)*amdw) ! input mass mixing ratio
801 ! h2ovmr(k)= max(f_zero,qlyr(j1,kk)) ! input vol mixing ratio
802 ! o3vmr (k)= max(f_zero,olyr(j1,kk)) ! input vol mixing ratio
803 !ncep model use
804  h2ovmr(k)= max(f_zero,qlyr(j1,kk)*amdw/(f_one-qlyr(j1,kk))) ! input specific humidity
805  o3vmr(k)= max(f_zero,olyr(j1,kk)*amdo3) ! input mass mixing ratio
806 
807  tem0 = (f_one - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
808  coldry(k) = tem2 * delp(k) / (tem1*tem0*(f_one + h2ovmr(k)))
809  temcol(k) = 1.0e-12 * coldry(k)
810 
811  colamt(k,1) = max(f_zero, coldry(k)*h2ovmr(k)) ! h2o
812  colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(j1,kk,1)) ! co2
813  colamt(k,3) = max(f_zero, coldry(k)*o3vmr(k)) ! o3
814  colmol(k) = coldry(k) + colamt(k,1)
815  enddo
816 
817 ! --- ... set up gas column amount, convert from volume mixing ratio
818 ! to molec/cm2 based on coldry (scaled to 1.0e-20)
819 
820  if (iswrgas > 0) then
821  do k = 1, nlay
822  kk = nlp1 - k
823  colamt(k,4) = max(temcol(k), coldry(k)*gasvmr(j1,kk,2)) ! n2o
824  colamt(k,5) = max(temcol(k), coldry(k)*gasvmr(j1,kk,3)) ! ch4
825  colamt(k,6) = max(temcol(k), coldry(k)*gasvmr(j1,kk,4)) ! o2
826 ! colamt(k,7) = max(temcol(k), coldry(k)*gasvmr(j1,kk,5)) ! co - notused
827  enddo
828  else
829  do k = 1, nlay
830  colamt(k,4) = temcol(k) ! n2o
831  colamt(k,5) = temcol(k) ! ch4
832  colamt(k,6) = temcol(k) ! o2
833 ! colamt(k,7) = temcol(k) ! co - notused
834  enddo
835  endif
836 
838 ! --- ... set aerosol optical properties
839 
840  do k = 1, nlay
841  kk = nlp1 - k
842  do ib = 1, nbdsw
843  tauae(k,ib) = aerosols(j1,kk,ib,1)
844  ssaae(k,ib) = aerosols(j1,kk,ib,2)
845  asyae(k,ib) = aerosols(j1,kk,ib,3)
846  enddo
847  enddo
848 
850  if (iswcliq > 0) then ! use prognostic cloud method
851  do k = 1, nlay
852  kk = nlp1 - k
853  cfrac(k) = clouds(j1,kk,1) ! cloud fraction
854  cliqp(k) = clouds(j1,kk,2) ! cloud liq path
855  reliq(k) = clouds(j1,kk,3) ! liq partical effctive radius
856  cicep(k) = clouds(j1,kk,4) ! cloud ice path
857  reice(k) = clouds(j1,kk,5) ! ice partical effctive radius
858  cdat1(k) = clouds(j1,kk,6) ! cloud rain drop path
859  cdat2(k) = clouds(j1,kk,7) ! rain partical effctive radius
860  cdat3(k) = clouds(j1,kk,8) ! cloud snow path
861  cdat4(k) = clouds(j1,kk,9) ! snow partical effctive radius
862  enddo
863  else ! use diagnostic cloud method
864  do k = 1, nlay
865  kk = nlp1 - k
866  cfrac(k) = clouds(j1,kk,1) ! cloud fraction
867  cdat1(k) = clouds(j1,kk,2) ! cloud optical depth
868  cdat2(k) = clouds(j1,kk,3) ! cloud single scattering albedo
869  cdat3(k) = clouds(j1,kk,4) ! cloud asymmetry factor
870  enddo
871  endif ! end if_iswcliq
872 
873  else ! input from sfc to toa
874 
875  tem1 = 100.0 * con_g
876  tem2 = 1.0e-20 * 1.0e3 * con_avgd
877 
878  do k = 1, nlay
879  pavel(k) = plyr(j1,k)
880  tavel(k) = tlyr(j1,k)
881  delp(k) = plvl(j1,k) - plvl(j1,k+1)
882 
883 ! --- ... set absorber amount
884 !test use
885 ! h2ovmr(k)= max(f_zero,qlyr(j1,k)*amdw) ! input mass mixing ratio
886 ! h2ovmr(k)= max(f_zero,qlyr(j1,k)) ! input vol mixing ratio
887 ! o3vmr (k)= max(f_zero,olyr(j1,k)) ! input vol mixing ratio
888 !ncep model use
889  h2ovmr(k)= max(f_zero,qlyr(j1,k)*amdw/(f_one-qlyr(j1,k))) ! input specific humidity
890  o3vmr(k)= max(f_zero,olyr(j1,k)*amdo3) ! input mass mixing ratio
891 
892  tem0 = (f_one - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
893  coldry(k) = tem2 * delp(k) / (tem1*tem0*(f_one + h2ovmr(k)))
894  temcol(k) = 1.0e-12 * coldry(k)
895 
896  colamt(k,1) = max(f_zero, coldry(k)*h2ovmr(k)) ! h2o
897  colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(j1,k,1)) ! co2
898  colamt(k,3) = max(f_zero, coldry(k)*o3vmr(k)) ! o3
899  colmol(k) = coldry(k) + colamt(k,1)
900  enddo
901 
902 
903  if (lprnt) then
904  if (ipt == 1) then
905  write(0,*)' pavel=',pavel
906  write(0,*)' tavel=',tavel
907  write(0,*)' delp=',delp
908  write(0,*)' h2ovmr=',h2ovmr*1000
909  write(0,*)' o3vmr=',o3vmr*1000000
910  endif
911  endif
912 
913 ! --- ... set up gas column amount, convert from volume mixing ratio
914 ! to molec/cm2 based on coldry (scaled to 1.0e-20)
915 
916  if (iswrgas > 0) then
917  do k = 1, nlay
918  colamt(k,4) = max(temcol(k), coldry(k)*gasvmr(j1,k,2)) ! n2o
919  colamt(k,5) = max(temcol(k), coldry(k)*gasvmr(j1,k,3)) ! ch4
920  colamt(k,6) = max(temcol(k), coldry(k)*gasvmr(j1,k,4)) ! o2
921 ! colamt(k,7) = max(temcol(k), coldry(k)*gasvmr(j1,k,5)) ! co - notused
922  enddo
923  else
924  do k = 1, nlay
925  colamt(k,4) = temcol(k) ! n2o
926  colamt(k,5) = temcol(k) ! ch4
927  colamt(k,6) = temcol(k) ! o2
928 ! colamt(k,7) = temcol(k) ! co - notused
929  enddo
930  endif
931 
932 ! --- ... set aerosol optical properties
933 
934  do ib = 1, nbdsw
935  do k = 1, nlay
936  tauae(k,ib) = aerosols(j1,k,ib,1)
937  ssaae(k,ib) = aerosols(j1,k,ib,2)
938  asyae(k,ib) = aerosols(j1,k,ib,3)
939  enddo
940  enddo
941 
942  if (iswcliq > 0) then ! use prognostic cloud method
943  do k = 1, nlay
944  cfrac(k) = clouds(j1,k,1) ! cloud fraction
945  cliqp(k) = clouds(j1,k,2) ! cloud liq path
946  reliq(k) = clouds(j1,k,3) ! liq partical effctive radius
947  cicep(k) = clouds(j1,k,4) ! cloud ice path
948  reice(k) = clouds(j1,k,5) ! ice partical effctive radius
949  cdat1(k) = clouds(j1,k,6) ! cloud rain drop path
950  cdat2(k) = clouds(j1,k,7) ! rain partical effctive radius
951  cdat3(k) = clouds(j1,k,8) ! cloud snow path
952  cdat4(k) = clouds(j1,k,9) ! snow partical effctive radius
953  enddo
954  else ! use diagnostic cloud method
955  do k = 1, nlay
956  cfrac(k) = clouds(j1,k,1) ! cloud fraction
957  cdat1(k) = clouds(j1,k,2) ! cloud optical depth
958  cdat2(k) = clouds(j1,k,3) ! cloud single scattering albedo
959  cdat3(k) = clouds(j1,k,4) ! cloud asymmetry factor
960  enddo
961  endif ! end if_iswcliq
962 
963  endif ! if_ivflip
964 
969 ! --- ... compute fractions of clear sky view
970 
971  zcf0 = f_one
972  zcf1 = f_one
973  if (iovrsw == 0) then ! random overlapping
974  do k = 1, nlay
975  zcf0 = zcf0 * (f_one - cfrac(k))
976  enddo
977  else if (iovrsw == 1) then ! max/ran overlapping
978  do k = 1, nlay
979  if (cfrac(k) > ftiny) then ! cloudy layer
980  zcf1 = min( zcf1, f_one-cfrac(k) )
981  elseif (zcf1 < f_one) then ! clear layer
982  zcf0 = zcf0 * zcf1
983  zcf1 = f_one
984  endif
985  enddo
986  zcf0 = zcf0 * zcf1
987  else if (iovrsw == 2) then ! maximum overlapping
988  do k = 1, nlay
989  zcf0 = min( zcf0, f_one-cfrac(k) )
990  enddo
991  endif
992 
993  if (zcf0 <= ftiny) zcf0 = f_zero
994  if (zcf0 > oneminus) zcf0 = f_one
995  zcf1 = f_one - zcf0
996 
998 
999 ! --- ... compute cloud optical properties
1000 
1001  if (zcf1 > f_zero) then ! cloudy sky column
1002 
1003  call cldprop &
1004 ! --- inputs:
1005  & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, &
1006  & zcf1, nlay, ipseed(j1), &
1007 ! --- outputs:
1008  & taucw, ssacw, asycw, cldfrc, cldfmc &
1009  & )
1010 
1011  else ! clear sky column
1012  cldfrc(:) = f_zero
1013  cldfmc(:,:)= f_zero
1014  do i = 1, nbdsw
1015  do k = 1, nlay
1016  taucw(k,i) = f_zero
1017  ssacw(k,i) = f_zero
1018  asycw(k,i) = f_zero
1019  enddo
1020  enddo
1021  endif ! end if_zcf1_block
1022 
1024  call setcoef &
1025 ! --- inputs:
1026  & ( pavel,tavel,h2ovmr, nlay,nlp1, &
1027 ! --- outputs:
1028  & laytrop,jp,jt,jt1,fac00,fac01,fac10,fac11, &
1029  & selffac,selffrac,indself,forfac,forfrac,indfor &
1030  & )
1031 
1036 ! --- ... calculate optical depths for gaseous absorption and Rayleigh
1037 ! scattering
1038 
1039  call taumol &
1040 ! --- inputs:
1041  & ( colamt,colmol,fac00,fac01,fac10,fac11,jp,jt,jt1,laytrop, &
1042  & forfac,forfrac,indfor,selffac,selffrac,indself, nlay, &
1043 ! --- outputs:
1044  & sfluxzen, taug, taur &
1045  & )
1046 
1050 ! --- ... call the 2-stream radiation transfer model
1051 
1052  if ( isubcsw <= 0 ) then ! use standard cloud scheme
1053 
1054  call spcvrtc &
1055 ! --- inputs:
1056  & ( ssolar,cosz1,sntz1,albbm,albdf,sfluxzen,cldfrc, &
1057  & zcf1,zcf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
1058  & nlay, nlp1, &
1059 ! --- outputs:
1060  & fxupc,fxdnc,fxup0,fxdn0, &
1061  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
1062  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
1063  & )
1064 
1065  else ! use mcica cloud scheme
1066 
1067  call spcvrtm &
1068 ! --- inputs:
1069  & ( ssolar,cosz1,sntz1,albbm,albdf,sfluxzen,cldfmc, &
1070  & zcf1,zcf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
1071  & nlay, nlp1, &
1072 ! --- outputs:
1073  & fxupc,fxdnc,fxup0,fxdn0, &
1074  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
1075  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
1076  & )
1077 
1078  endif
1080 ! --- ... sum up total spectral fluxes for total-sky
1081 
1082  do k = 1, nlp1
1083  flxuc(k) = f_zero
1084  flxdc(k) = f_zero
1085 
1086  do ib = 1, nbdsw
1087  flxuc(k) = flxuc(k) + fxupc(k,ib)
1088  flxdc(k) = flxdc(k) + fxdnc(k,ib)
1089  enddo
1090  enddo
1091 
1092 !! --- ... optional clear sky fluxes
1093 
1094  if ( lhsw0 .or. lflxprf ) then
1095  do k = 1, nlp1
1096  flxu0(k) = f_zero
1097  flxd0(k) = f_zero
1098 
1099  do ib = 1, nbdsw
1100  flxu0(k) = flxu0(k) + fxup0(k,ib)
1101  flxd0(k) = flxd0(k) + fxdn0(k,ib)
1102  enddo
1103  enddo
1104  endif
1105 
1106 ! --- ... prepare for final outputs
1107 
1108  do k = 1, nlay
1109  rfdelp(k) = heatfac / delp(k)
1110  enddo
1111 
1112  if ( lfdncmp ) then
1113 !! --- ... optional uv-b surface downward flux
1114  fdncmp(j1)%uvbf0 = suvbf0
1115  fdncmp(j1)%uvbfc = suvbfc
1116 
1117 !! --- ... optional beam and diffuse sfc fluxes
1118  fdncmp(j1)%nirbm = sfbmc(1)
1119  fdncmp(j1)%nirdf = sfdfc(1)
1120  fdncmp(j1)%visbm = sfbmc(2)
1121  fdncmp(j1)%visdf = sfdfc(2)
1122  endif ! end if_lfdncmp
1123 
1124 ! --- ... toa and sfc fluxes
1125 
1126  topflx(j1)%upfxc = ftoauc
1127  topflx(j1)%dnfxc = ftoadc
1128  topflx(j1)%upfx0 = ftoau0
1129 
1130  sfcflx(j1)%upfxc = fsfcuc
1131  sfcflx(j1)%dnfxc = fsfcdc
1132  sfcflx(j1)%upfx0 = fsfcu0
1133  sfcflx(j1)%dnfx0 = fsfcd0
1134 
1135  if (ivflip == 0) then ! output from toa to sfc
1136 
1137 ! --- ... compute heating rates
1138 
1139  fnet(1) = flxdc(1) - flxuc(1)
1140 
1141  do k = 2, nlp1
1142  kk = nlp1 - k + 1
1143  fnet(k) = flxdc(k) - flxuc(k)
1144  hswc(j1,kk) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1145  enddo
1146 
1147 !! --- ... optional flux profiles
1148 
1149  if ( lflxprf ) then
1150  do k = 1, nlp1
1151  kk = nlp1 - k + 1
1152  flxprf(j1,kk)%upfxc = flxuc(k)
1153  flxprf(j1,kk)%dnfxc = flxdc(k)
1154  flxprf(j1,kk)%upfx0 = flxu0(k)
1155  flxprf(j1,kk)%dnfx0 = flxd0(k)
1156  enddo
1157  endif
1158 
1159 !! --- ... optional clear sky heating rates
1160 
1161  if ( lhsw0 ) then
1162  fnet(1) = flxd0(1) - flxu0(1)
1163 
1164  do k = 2, nlp1
1165  kk = nlp1 - k + 1
1166  fnet(k) = flxd0(k) - flxu0(k)
1167  hsw0(j1,kk) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1168  enddo
1169  endif
1170 
1171 !! --- ... optional spectral band heating rates
1172 
1173  if ( lhswb ) then
1174  do mb = 1, nbdsw
1175  fnet(1) = fxdnc(1,mb) - fxupc(1,mb)
1176 
1177  do k = 2, nlp1
1178  kk = nlp1 - k + 1
1179  fnet(k) = fxdnc(k,mb) - fxupc(k,mb)
1180  hswb(j1,kk,mb) = (fnet(k) - fnet(k-1)) * rfdelp(k-1)
1181  enddo
1182  enddo
1183  endif
1184 
1185  else ! output from sfc to toa
1186 
1187 ! --- ... compute heating rates
1188 
1189  fnet(1) = flxdc(1) - flxuc(1)
1190 
1191  do k = 2, nlp1
1192  fnet(k) = flxdc(k) - flxuc(k)
1193  hswc(j1,k-1) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1194  enddo
1195 
1196 !! --- ... optional flux profiles
1197 
1198  if ( lflxprf ) then
1199  do k = 1, nlp1
1200  flxprf(j1,k)%upfxc = flxuc(k)
1201  flxprf(j1,k)%dnfxc = flxdc(k)
1202  flxprf(j1,k)%upfx0 = flxu0(k)
1203  flxprf(j1,k)%dnfx0 = flxd0(k)
1204  enddo
1205  endif
1206 
1207 !! --- ... optional clear sky heating rates
1208 
1209  if ( lhsw0 ) then
1210  fnet(1) = flxd0(1) - flxu0(1)
1211 
1212  do k = 2, nlp1
1213  fnet(k) = flxd0(k) - flxu0(k)
1214  hsw0(j1,k-1) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1215  enddo
1216  endif
1217 
1218 !! --- ... optional spectral band heating rates
1219 
1220  if ( lhswb ) then
1221  do mb = 1, nbdsw
1222  fnet(1) = fxdnc(1,mb) - fxupc(1,mb)
1223 
1224  do k = 1, nlay
1225  fnet(k+1) = fxdnc(k+1,mb) - fxupc(k+1,mb)
1226  hswb(j1,k,mb) = (fnet(k+1) - fnet(k)) * rfdelp(k)
1227  enddo
1228  enddo
1229  endif
1230 
1231  endif ! if_ivflip
1232 
1233  enddo lab_do_ipt
1234 
1235  return
1236 !...................................
1237  end subroutine swrad
1238 !-----------------------------------
1240 
1269 !-----------------------------------
1270  subroutine rswinit &
1271 !...................................
1273 ! --- inputs:
1274  & ( me )
1275 ! --- outputs: (none)
1276 
1277 ! =================== program usage description =================== !
1278 ! !
1279 ! purpose: initialize non-varying module variables, conversion factors,!
1280 ! and look-up tables. !
1281 ! !
1282 ! subprograms called: none !
1283 ! !
1284 ! ==================== defination of variables ==================== !
1285 ! !
1286 ! inputs: !
1287 ! me - print control for parallel process !
1288 ! !
1289 ! outputs: (none) !
1290 ! !
1291 ! external module variables: (in physparam) !
1292 ! iswrate - heating rate unit selections !
1293 ! =1: output in k/day !
1294 ! =2: output in k/second !
1295 ! iswrgas - control flag for rare gases (ch4,n2o,o2, etc.) !
1296 ! =0: do not include rare gases !
1297 ! >0: include all rare gases !
1298 ! iswcliq - liquid cloud optical properties contrl flag !
1299 ! =0: input cloud opt depth from diagnostic scheme !
1300 ! >0: input cwp,rew, and other cloud content parameters !
1301 ! isubcsw - sub-column cloud approximation control flag !
1302 ! =0: no sub-col cld treatment, use grid-mean cld quantities !
1303 ! =1: mcica sub-col, prescribed seeds to get random numbers !
1304 ! =2: mcica sub-col, providing array icseed for random numbers!
1305 ! icldflg - cloud scheme control flag !
1306 ! =0: diagnostic scheme gives cloud tau, omiga, and g. !
1307 ! =1: prognostic scheme gives cloud liq/ice path, etc. !
1308 ! iovrsw - clouds vertical overlapping control flag !
1309 ! =0: random overlapping clouds !
1310 ! =1: maximum/random overlapping clouds !
1311 ! =2: maximum overlap cloud !
1312 ! iswmode - control flag for 2-stream transfer scheme !
1313 ! =1; delta-eddington (joseph et al., 1976) !
1314 ! =2: pifm (zdunkowski et al., 1980) !
1315 ! =3: discrete ordinates (liou, 1973) !
1316 ! !
1317 ! ******************************************************************* !
1318 ! !
1319 ! definitions: !
1320 ! arrays for 10000-point look-up tables: !
1321 ! tau_tbl clear-sky optical depth !
1322 ! exp_tbl exponential lookup table for transmittance !
1323 ! !
1324 ! ******************************************************************* !
1325 ! !
1326 ! ====================== end of description block ================= !
1327 
1328 ! --- inputs:
1329  integer, intent(in) :: me
1330 
1331 ! --- outputs: none
1332 
1333 ! --- locals:
1334  real (kind=kind_phys), parameter :: expeps = 1.e-20
1335 
1336  integer :: i
1337 
1338  real (kind=kind_phys) :: tfn, tau
1339 
1340 !
1341 !===> ... begin here
1342 !
1343  if ( iovrsw<0 .or. iovrsw>2 ) then
1344  print *,' *** Error in specification of cloud overlap flag', &
1345  & ' IOVRSW=',iovrsw,' in RSWINIT !!'
1346  stop
1347  endif
1348 
1349  if (me == 0) then
1350  print *,' - Using AER Shortwave Radiation, Version: ',vtagsw
1351 
1352  if (iswmode == 1) then
1353  print *,' --- Delta-eddington 2-stream transfer scheme'
1354  else if (iswmode == 2) then
1355  print *,' --- PIFM 2-stream transfer scheme'
1356  else if (iswmode == 3) then
1357  print *,' --- Discrete ordinates 2-stream transfer scheme'
1358  endif
1359 
1360  if (iswrgas <= 0) then
1361  print *,' --- Rare gases absorption is NOT included in SW'
1362  else
1363  print *,' --- Include rare gases N2O, CH4, O2, absorptions',&
1364  & ' in SW'
1365  endif
1366 
1367  if ( isubcsw == 0 ) then
1368  print *,' --- Using standard grid average clouds, no ', &
1369  & 'sub-column clouds approximation applied'
1370  elseif ( isubcsw == 1 ) then
1371  print *,' --- Using MCICA sub-colum clouds approximation ', &
1372  & 'with a prescribed sequence of permutation seeds'
1373  elseif ( isubcsw == 2 ) then
1374  print *,' --- Using MCICA sub-colum clouds approximation ', &
1375  & 'with provided input array of permutation seeds'
1376  else
1377  print *,' *** Error in specification of sub-column cloud ', &
1378  & ' control flag isubcsw =',isubcsw,' !!'
1379  stop
1380  endif
1381  endif
1382 
1383 ! --- ... check cloud flags for consistency
1384 
1385  if ((icldflg == 0 .and. iswcliq /= 0) .or. &
1386  & (icldflg == 1 .and. iswcliq == 0)) then
1387  print *,' *** Model cloud scheme inconsistent with SW', &
1388  & ' radiation cloud radiative property setup !!'
1389  stop
1390  endif
1391 
1392 ! --- ... setup constant factors for heating rate
1393 ! the 1.0e-2 is to convert pressure from mb to N/m**2
1394 
1395  if (iswrate == 1) then
1396 ! heatfac = 8.4391
1397 ! heatfac = con_g * 86400. * 1.0e-2 / con_cp ! (in k/day)
1398  heatfac = con_g * 864.0 / con_cp ! (in k/day)
1399  else
1400  heatfac = con_g * 1.0e-2 / con_cp ! (in k/second)
1401  endif
1402 
1403 ! --- ... define exponential lookup tables for transmittance. tau is
1404 ! computed as a function of the tau transition function, and
1405 ! transmittance is calculated as a function of tau. all tables
1406 ! are computed at intervals of 0.0001. the inverse of the
1407 ! constant used in the Pade approximation to the tau transition
1408 ! function is set to bpade.
1409 
1410  exp_tbl(0) = 1.0
1411  exp_tbl(ntbmx) = expeps
1412 
1413  do i = 1, ntbmx-1
1414  tfn = float(i) / float(ntbmx-i)
1415  tau = bpade * tfn
1416  exp_tbl(i) = exp( -tau )
1417  enddo
1418 
1419  return
1420 !...................................
1421  end subroutine rswinit
1422 !-----------------------------------
1423 
1426 !-----------------------------------
1427  subroutine cldprop &
1428 !...................................
1429 ! --- inputs:
1430  & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, &
1431  & cf1, nlay, ipseed, &
1432 ! --- output:
1433  & taucw, ssacw, asycw, cldfrc, cldfmc &
1434  & )
1435 
1436 ! =================== program usage description =================== !
1437 ! !
1438 ! Purpose: Compute the cloud optical properties for each cloudy layer !
1439 ! and g-point interval. !
1440 ! !
1441 ! subprograms called: none !
1442 ! !
1443 ! ==================== defination of variables ==================== !
1444 ! !
1445 ! inputs: size !
1446 ! cfrac - real, layer cloud fraction nlay !
1447 ! ..... for iswcliq > 0 (prognostic cloud sckeme) - - - !
1448 ! cliqp - real, layer in-cloud liq water path (g/m**2) nlay !
1449 ! reliq - real, mean eff radius for liq cloud (micron) nlay !
1450 ! cicep - real, layer in-cloud ice water path (g/m**2) nlay !
1451 ! reice - real, mean eff radius for ice cloud (micron) nlay !
1452 ! cdat1 - real, layer rain drop water path (g/m**2) nlay !
1453 ! cdat2 - real, effective radius for rain drop (micron) nlay !
1454 ! cdat3 - real, layer snow flake water path(g/m**2) nlay !
1455 ! cdat4 - real, mean eff radius for snow flake(micron) nlay !
1456 ! ..... for iswcliq = 0 (diagnostic cloud sckeme) - - - !
1457 ! cdat1 - real, layer cloud optical depth nlay !
1458 ! cdat2 - real, layer cloud single scattering albedo nlay !
1459 ! cdat3 - real, layer cloud asymmetry factor nlay !
1460 ! cdat4 - real, optional use nlay !
1461 ! cliqp - real, not used nlay !
1462 ! cicep - real, not used nlay !
1463 ! reliq - real, not used nlay !
1464 ! reice - real, not used nlay !
1465 ! !
1466 ! cf1 - real, effective total cloud cover at surface 1 !
1467 ! nlay - integer, vertical layer number 1 !
1468 ! ipseed- permutation seed for generating random numbers (isubcsw>0) !
1469 ! !
1470 ! outputs: !
1471 ! taucw - real, cloud optical depth, w/o delta scaled nlay*nbdsw !
1472 ! ssacw - real, weighted cloud single scattering albedo nlay*nbdsw !
1473 ! (ssa = ssacw / taucw) !
1474 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
1475 ! (asy = asycw / ssacw) !
1476 ! cldfrc - real, cloud fraction of grid mean value nlay !
1477 ! cldfmc - real, cloud fraction for each sub-column nlay*ngptsw!
1478 ! !
1479 ! !
1480 ! explanation of the method for each value of iswcliq, and iswcice. !
1481 ! set up in module "physparam" !
1482 ! !
1483 ! iswcliq=0 : input cloud optical property (tau, ssa, asy). !
1484 ! (used for diagnostic cloud method) !
1485 ! iswcliq>0 : input cloud liq/ice path and effective radius, also !
1486 ! require the user of 'iswcice' to specify the method !
1487 ! used to compute aborption due to water/ice parts. !
1488 ! ................................................................... !
1489 ! !
1490 ! iswcliq=1 : liquid water cloud optical properties are computed !
1491 ! as in hu and stamnes (1993), j. clim., 6, 728-742. !
1492 ! !
1493 ! iswcice used only when iswcliq > 0 !
1494 ! the cloud ice path (g/m2) and ice effective radius !
1495 ! (microns) are inputs. !
1496 ! iswcice=1 : ice cloud optical properties are computed as in !
1497 ! ebert and curry (1992), jgr, 97, 3831-3836. !
1498 ! iswcice=2 : ice cloud optical properties are computed as in !
1499 ! streamer v3.0 (2001), key, streamer user's guide, !
1500 ! cooperative institude for meteorological studies,95pp!
1501 ! iswcice=3 : ice cloud optical properties are computed as in !
1502 ! fu (1996), j. clim., 9. !
1503 ! !
1504 ! other cloud control module variables: !
1505 ! isubcsw =0: standard cloud scheme, no sub-col cloud approximation !
1506 ! >0: mcica sub-col cloud scheme using ipseed as permutation!
1507 ! seed for generating rundom numbers !
1508 ! !
1509 ! ====================== end of description block ================= !
1510 !
1512 
1513 ! --- inputs:
1514  integer, intent(in) :: nlay, ipseed
1515  real (kind=kind_phys), intent(in) :: cf1
1516 
1517  real (kind=kind_phys), dimension(nlay), intent(in) :: cliqp, &
1518  & reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, cfrac
1519 
1520 ! --- outputs:
1521  real (kind=kind_phys), dimension(nlay,ngptsw), intent(out) :: &
1522  & cldfmc
1523  real (kind=kind_phys), dimension(nlay,nbdsw), intent(out) :: &
1524  & taucw, ssacw, asycw
1525  real (kind=kind_phys), dimension(nlay), intent(out) :: cldfrc
1526 
1527 ! --- locals:
1528  real (kind=kind_phys), dimension(nblow:nbhgh) :: tauliq, tauice, &
1529  & ssaliq, ssaice, ssaran, ssasnw, asyliq, asyice, &
1530  & asyran, asysnw
1531  real (kind=kind_phys), dimension(nlay) :: cldf
1532 
1533  real (kind=kind_phys) :: dgeice, factor, fint, tauran, tausnw, &
1534  & cldliq, refliq, cldice, refice, cldran, cldsnw, refsnw, &
1535  & extcoliq, ssacoliq, asycoliq, extcoice, ssacoice, asycoice,&
1536  & dgesnw
1537 
1538  logical :: lcloudy(nlay,ngptsw)
1539  integer :: ia, ib, ig, jb, k, index
1540 
1541 !
1542 !===> ... begin here
1543 !
1544  do ib = 1, nbdsw
1545  do k = 1, nlay
1546  taucw(k,ib) = f_zero
1547  ssacw(k,ib) = f_one
1548  asycw(k,ib) = f_zero
1549  enddo
1550  enddo
1551 
1552 ! --- ... compute cloud radiative properties for a cloudy column
1553 
1554  lab_if_iswcliq : if (iswcliq > 0) then
1555 
1556  lab_do_k : do k = 1, nlay
1557  lab_if_cld : if (cfrac(k) > ftiny) then
1558 
1559 ! --- ... optical properties for rain and snow
1560  cldran = cdat1(k)
1561 ! refran = cdat2(k)
1562  cldsnw = cdat3(k)
1563  refsnw = cdat4(k)
1564  dgesnw = 1.0315 * refsnw ! for fu's snow formula
1565 
1566  tauran = cldran * a0r
1567 
1568 
1569 ! --- if use fu's formula it needs to be normalized by snow/ice density
1570 ! !not use snow density = 0.1 g/cm**3 = 0.1 g/(mu * m**2)
1571 ! use ice density = 0.9167 g/cm**3 = 0.9167 g/(mu * m**2)
1572 ! 1/0.9167 = 1.09087
1573 ! factor 1.5396=8/(3*sqrt(3)) converts reff to generalized ice particle size
1574 ! use newer factor value 1.0315
1575  if (cldsnw>f_zero .and. refsnw>10.0_kind_phys) then
1576 ! tausnw = cldsnw * (a0s + a1s/refsnw)
1577  tausnw = cldsnw*1.09087*(a0s + a1s/dgesnw) ! fu's formula
1578  else
1579  tausnw = f_zero
1580  endif
1581 
1582  do ib = nblow, nbhgh
1583  ssaran(ib) = tauran * (f_one - b0r(ib))
1584  ssasnw(ib) = tausnw * (f_one - (b0s(ib)+b1s(ib)*dgesnw))
1585  asyran(ib) = ssaran(ib) * c0r(ib)
1586  asysnw(ib) = ssasnw(ib) * c0s(ib)
1587  enddo
1588 
1589  cldliq = cliqp(k)
1590  cldice = cicep(k)
1591  refliq = reliq(k)
1592  refice = reice(k)
1593 
1594 ! --- ... calculation of absorption coefficients due to water clouds.
1595 
1596  if ( cldliq <= f_zero ) then
1597  do ib = nblow, nbhgh
1598  tauliq(ib) = f_zero
1599  ssaliq(ib) = f_zero
1600  asyliq(ib) = f_zero
1601  enddo
1602  else
1603  if ( iswcliq == 1 ) then
1604  factor = refliq - 1.5
1605  index = max( 1, min( 57, int( factor ) ))
1606  fint = factor - float(index)
1607 
1608  do ib = nblow, nbhgh
1609  extcoliq = max(f_zero, extliq1(index,ib) &
1610  & + fint*(extliq1(index+1,ib)-extliq1(index,ib)) )
1611  ssacoliq = max(f_zero, min(f_one, ssaliq1(index,ib) &
1612  & + fint*(ssaliq1(index+1,ib)-ssaliq1(index,ib)) ))
1613 
1614  asycoliq = max(f_zero, min(f_one, asyliq1(index,ib) &
1615  & + fint*(asyliq1(index+1,ib)-asyliq1(index,ib)) ))
1616 ! forcoliq = asycoliq * asycoliq
1617 
1618  tauliq(ib) = cldliq * extcoliq
1619  ssaliq(ib) = tauliq(ib) * ssacoliq
1620  asyliq(ib) = ssaliq(ib) * asycoliq
1621  enddo
1622  endif ! end if_iswcliq_block
1623  endif ! end if_cldliq_block
1624 
1625 ! --- ... calculation of absorption coefficients due to ice clouds.
1626 
1627  if ( cldice <= f_zero ) then
1628  do ib = nblow, nbhgh
1629  tauice(ib) = f_zero
1630  ssaice(ib) = f_zero
1631  asyice(ib) = f_zero
1632  enddo
1633  else
1634 
1635 ! --- ... ebert and curry approach for all particle sizes though somewhat
1636 ! unjustified for large ice particles
1637 
1638  if ( iswcice == 1 ) then
1639  refice = min(130.0_kind_phys,max(13.0_kind_phys,refice))
1640 
1641  do ib = nblow, nbhgh
1642  ia = idxebc(ib) ! eb_&_c band index for ice cloud coeff
1643 
1644  extcoice = max(f_zero, abari(ia)+bbari(ia)/refice )
1645  ssacoice = max(f_zero, min(f_one, &
1646  & f_one-cbari(ia)-dbari(ia)*refice ))
1647  asycoice = max(f_zero, min(f_one, &
1648  & ebari(ia)+fbari(ia)*refice ))
1649 ! forcoice = asycoice * asycoice
1650 
1651  tauice(ib) = cldice * extcoice
1652  ssaice(ib) = tauice(ib) * ssacoice
1653  asyice(ib) = ssaice(ib) * asycoice
1654  enddo
1655 
1656 ! --- ... streamer approach for ice effective radius between 5.0 and 131.0 microns
1657 
1658  elseif ( iswcice == 2 ) then
1659  refice = min(131.0_kind_phys,max(5.0_kind_phys,refice))
1660 
1661  factor = (refice - 2.0) / 3.0
1662  index = max( 1, min( 42, int( factor ) ))
1663  fint = factor - float(index)
1664 
1665  do ib = nblow, nbhgh
1666  extcoice = max(f_zero, extice2(index,ib) &
1667  & + fint*(extice2(index+1,ib)-extice2(index,ib)) )
1668  ssacoice = max(f_zero, min(f_one, ssaice2(index,ib) &
1669  & + fint*(ssaice2(index+1,ib)-ssaice2(index,ib)) ))
1670  asycoice = max(f_zero, min(f_one, asyice2(index,ib) &
1671  & + fint*(asyice2(index+1,ib)-asyice2(index,ib)) ))
1672 ! forcoice = asycoice * asycoice
1673 
1674  tauice(ib) = cldice * extcoice
1675  ssaice(ib) = tauice(ib) * ssacoice
1676  asyice(ib) = ssaice(ib) * asycoice
1677  enddo
1678 
1679 ! --- ... fu's approach for ice effective radius between 4.8 and 135 microns
1680 ! (generalized effective size from 5 to 140 microns)
1681 
1682  elseif ( iswcice == 3 ) then
1683  dgeice = max( 5.0, min( 140.0, 1.0315*refice ))
1684 
1685  factor = (dgeice - 2.0) / 3.0
1686  index = max( 1, min( 45, int( factor ) ))
1687  fint = factor - float(index)
1688 
1689  do ib = nblow, nbhgh
1690  extcoice = max(f_zero, extice3(index,ib) &
1691  & + fint*(extice3(index+1,ib)-extice3(index,ib)) )
1692  ssacoice = max(f_zero, min(f_one, ssaice3(index,ib) &
1693  & + fint*(ssaice3(index+1,ib)-ssaice3(index,ib)) ))
1694  asycoice = max(f_zero, min(f_one, asyice3(index,ib) &
1695  & + fint*(asyice3(index+1,ib)-asyice3(index,ib)) ))
1696 ! fdelta = max(f_zero, min(f_one, fdlice3(index,ib) &
1697 ! & + fint*(fdlice3(index+1,ib)-fdlice3(index,ib)) ))
1698 ! forcoice = min( asycoice, fdelta+0.5/ssacoice ) ! see fu 1996 p. 2067
1699 
1700  tauice(ib) = cldice * extcoice
1701  ssaice(ib) = tauice(ib) * ssacoice
1702  asyice(ib) = ssaice(ib) * asycoice
1703  enddo
1704 
1705  endif ! end if_iswcice_block
1706  endif ! end if_cldice_block
1707 
1708  do ib = 1, nbdsw
1709  jb = nblow + ib - 1
1710  taucw(k,ib) = tauliq(jb)+tauice(jb)+tauran+tausnw
1711  ssacw(k,ib) = ssaliq(jb)+ssaice(jb)+ssaran(jb)+ssasnw(jb)
1712  asycw(k,ib) = asyliq(jb)+asyice(jb)+asyran(jb)+asysnw(jb)
1713  enddo
1714 
1715  endif lab_if_cld
1716  enddo lab_do_k
1717 
1718  else lab_if_iswcliq
1719 
1720  do k = 1, nlay
1721  if (cfrac(k) > ftiny) then
1722  do ib = 1, nbdsw
1723  taucw(k,ib) = cdat1(k)
1724  ssacw(k,ib) = cdat1(k) * cdat2(k)
1725  asycw(k,ib) = ssacw(k,ib) * cdat3(k)
1726  enddo
1727  endif
1728  enddo
1729 
1730  endif lab_if_iswcliq
1731 
1732 ! --- distribute cloud properties to each g-point
1733 
1734  if ( isubcsw > 0 ) then ! mcica sub-col clouds approx
1735 
1736  cldf(:) = cfrac(:)
1737  where (cldf(:) < ftiny)
1738  cldf(:) = f_zero
1739  end where
1740 
1741 ! --- ... call sub-column cloud generator
1742 
1743  call mcica_subcol &
1744 ! --- inputs:
1745  & ( cldf, nlay, ipseed, &
1746 ! --- outputs:
1747  & lcloudy &
1748  & )
1749 
1750  do ig = 1, ngptsw
1751  do k = 1, nlay
1752  if ( lcloudy(k,ig) ) then
1753  cldfmc(k,ig) = f_one
1754  else
1755  cldfmc(k,ig) = f_zero
1756  endif
1757  enddo
1758  enddo
1759 
1760  else ! non-mcica, normalize cloud
1761 
1762  do k = 1, nlay
1763  cldfrc(k) = cfrac(k) / cf1
1764  enddo
1765  endif ! end if_isubcsw_block
1766 
1767  return
1768 !...................................
1769  end subroutine cldprop
1770 !-----------------------------------
1771 
1780 ! ----------------------------------
1781  subroutine mcica_subcol &
1782 ! ..................................
1783 ! --- inputs:
1784  & ( cldf, nlay, ipseed, &
1785 ! --- outputs:
1786  & lcloudy &
1787  & )
1788 
1789 ! ==================== defination of variables ==================== !
1790 ! !
1791 ! input variables: size !
1792 ! cldf - real, layer cloud fraction nlay !
1793 ! nlay - integer, number of model vertical layers 1 !
1794 ! ipseed - integer, permute seed for random num generator 1 !
1795 ! ** note : if the cloud generator is called multiple times, need !
1796 ! to permute the seed between each call; if between calls !
1797 ! for lw and sw, use values differ by the number of g-pts. !
1798 ! !
1799 ! output variables: !
1800 ! lcloudy - logical, sub-colum cloud profile flag array nlay*ngptsw!
1801 ! !
1802 ! other control flags from module variables: !
1803 ! iovrsw : control flag for cloud overlapping method !
1804 ! =0:random; =1:maximum/random; =2:maximum !
1805 ! !
1806 ! !
1807 ! ===================== end of definitions ==================== !
1808 
1809  implicit none
1810 
1811 ! --- inputs:
1812  integer, intent(in) :: nlay, ipseed
1813 
1814  real (kind=kind_phys), dimension(nlay), intent(in) :: cldf
1815 
1816 ! --- outputs:
1817  logical, dimension(nlay,ngptsw), intent(out):: lcloudy
1818 
1819 ! --- locals:
1820  real (kind=kind_phys) :: cdfunc(nlay,ngptsw), tem1, &
1821  & rand2d(nlay*ngptsw), rand1d(ngptsw)
1822 
1823  type(random_stat) :: stat ! for thread safe random generator
1824 
1825  integer :: k, n, k1
1826 !
1827 !===> ... begin here
1828 !
1829 ! --- ... advance randum number generator by ipseed values
1830 
1831  call random_setseed &
1832 ! --- inputs:
1833  & ( ipseed, &
1834 ! --- outputs:
1835  & stat &
1836  & )
1837 
1838 ! --- ... sub-column set up according to overlapping assumption
1839 
1840  select case ( iovrsw )
1841 
1842  case( 0 ) ! random overlap, pick a random value at every level
1843 
1844  call random_number &
1845 ! --- inputs: ( none )
1846 ! --- outputs:
1847  & ( rand2d, stat )
1848 
1849  k1 = 0
1850  do n = 1, ngptsw
1851  do k = 1, nlay
1852  k1 = k1 + 1
1853  cdfunc(k,n) = rand2d(k1)
1854  enddo
1855  enddo
1856 
1857  case( 1 ) ! max-ran overlap
1858 
1859  call random_number &
1860 ! --- inputs: ( none )
1861 ! --- outputs:
1862  & ( rand2d, stat )
1863 
1864  k1 = 0
1865  do n = 1, ngptsw
1866  do k = 1, nlay
1867  k1 = k1 + 1
1868  cdfunc(k,n) = rand2d(k1)
1869  enddo
1870  enddo
1871 
1872 ! --- first pick a random number for bottom/top layer.
1873 ! then walk up the column: (aer's code)
1874 ! if layer below is cloudy, use the same rand num in the layer below
1875 ! if layer below is clear, use a new random number
1876 
1877 ! --- from bottom up
1878  do k = 2, nlay
1879  k1 = k - 1
1880  tem1 = f_one - cldf(k1)
1881 
1882  do n = 1, ngptsw
1883  if ( cdfunc(k1,n) > tem1 ) then
1884  cdfunc(k,n) = cdfunc(k1,n)
1885  else
1886  cdfunc(k,n) = cdfunc(k,n) * tem1
1887  endif
1888  enddo
1889  enddo
1890 
1891 ! --- then walk down the column: (if use original author's method)
1892 ! if layer above is cloudy, use the same rand num in the layer above
1893 ! if layer above is clear, use a new random number
1894 
1895 ! --- from top down
1896 ! do k = nlay-1, 1, -1
1897 ! k1 = k + 1
1898 ! tem1 = f_one - cldf(k1)
1899 
1900 ! do n = 1, ngptsw
1901 ! if ( cdfunc(k1,n) > tem1 ) then
1902 ! cdfunc(k,n) = cdfunc(k1,n)
1903 ! else
1904 ! cdfunc(k,n) = cdfunc(k,n) * tem1
1905 ! endif
1906 ! enddo
1907 ! enddo
1908 
1909  case( 2 ) ! maximum overlap, pick same random numebr at every level
1910 
1911  call random_number &
1912 ! --- inputs: ( none )
1913 ! --- outputs:
1914  & ( rand1d, stat )
1915 
1916  do n = 1, ngptsw
1917  tem1 = rand1d(n)
1918 
1919  do k = 1, nlay
1920  cdfunc(k,n) = tem1
1921  enddo
1922  enddo
1923 
1924  end select
1925 
1926 ! --- ... generate subcolumns for homogeneous clouds
1927 
1928  do k = 1, nlay
1929  tem1 = f_one - cldf(k)
1930 
1931  do n = 1, ngptsw
1932  lcloudy(k,n) = cdfunc(k,n) >= tem1
1933  enddo
1934  enddo
1935 
1936  return
1937 ! ..................................
1938  end subroutine mcica_subcol
1939 ! ----------------------------------
1940 
1942 ! ----------------------------------
1943  subroutine setcoef &
1944 ! ..................................
1945 ! --- inputs:
1946  & ( pavel,tavel,h2ovmr, nlay,nlp1, &
1947 ! --- outputs:
1948  & laytrop,jp,jt,jt1,fac00,fac01,fac10,fac11, &
1949  & selffac,selffrac,indself,forfac,forfrac,indfor &
1950  & )
1951 
1952 ! =================== program usage description =================== !
1953 ! !
1954 ! purpose: compute various coefficients needed in radiative transfer !
1955 ! calculations. !
1956 ! !
1957 ! subprograms called: none !
1958 ! !
1959 ! ==================== defination of variables ==================== !
1960 ! !
1961 ! inputs: -size- !
1962 ! pavel - real, layer pressures (mb) nlay !
1963 ! tavel - real, layer temperatures (k) nlay !
1964 ! h2ovmr - real, layer w.v. volum mixing ratio (kg/kg) nlay !
1965 ! nlay/nlp1 - integer, total number of vertical layers, levels 1 !
1966 ! !
1967 ! outputs: !
1968 ! laytrop - integer, tropopause layer index (unitless) 1 !
1969 ! jp - real, indices of lower reference pressure nlay !
1970 ! jt, jt1 - real, indices of lower reference temperatures nlay !
1971 ! at levels of jp and jp+1 !
1972 ! facij - real, factors multiply the reference ks, nlay !
1973 ! i,j=0/1 for lower/higher of the 2 appropriate !
1974 ! temperatures and altitudes. !
1975 ! selffac - real, scale factor for w. v. self-continuum nlay !
1976 ! equals (w. v. density)/(atmospheric density !
1977 ! at 296k and 1013 mb) !
1978 ! selffrac - real, factor for temperature interpolation of nlay !
1979 ! reference w. v. self-continuum data !
1980 ! indself - integer, index of lower ref temp for selffac nlay !
1981 ! forfac - real, scale factor for w. v. foreign-continuum nlay !
1982 ! forfrac - real, factor for temperature interpolation of nlay !
1983 ! reference w.v. foreign-continuum data !
1984 ! indfor - integer, index of lower ref temp for forfac nlay !
1985 ! !
1986 ! ====================== end of definitions =================== !
1987 
1988 ! --- inputs:
1989  integer, intent(in) :: nlay, nlp1
1990 
1991  real (kind=kind_phys), dimension(:), intent(in) :: pavel, tavel, &
1992  & h2ovmr
1993 
1994 ! --- outputs:
1995  integer, dimension(nlay), intent(out) :: indself, indfor, &
1996  & jp, jt, jt1
1997  integer, intent(out) :: laytrop
1998 
1999  real (kind=kind_phys), dimension(nlay), intent(out) :: fac00, &
2000  & fac01, fac10, fac11, selffac, selffrac, forfac, forfrac
2001 
2002 ! --- locals:
2003  real (kind=kind_phys) :: plog, fp, fp1, ft, ft1, tem1, tem2
2004 
2005  integer :: i, k, jp1
2006 !
2007 !===> ... begin here
2008 !
2009  laytrop= nlay
2010 
2011  do k = 1, nlay
2012 
2013  forfac(k) = pavel(k)*stpfac / (tavel(k)*(f_one + h2ovmr(k)))
2014 
2015 ! --- ... find the two reference pressures on either side of the
2016 ! layer pressure. store them in jp and jp1. store in fp the
2017 ! fraction of the difference (in ln(pressure)) between these
2018 ! two values that the layer pressure lies.
2019 
2020  plog = log(pavel(k))
2021  jp(k) = max(1, min(58, int(36.0 - 5.0*(plog+0.04)) ))
2022  jp1 = jp(k) + 1
2023  fp = 5.0 * (preflog(jp(k)) - plog)
2024 
2025 ! --- ... determine, for each reference pressure (jp and jp1), which
2026 ! reference temperature (these are different for each reference
2027 ! pressure) is nearest the layer temperature but does not exceed it.
2028 ! store these indices in jt and jt1, resp. store in ft (resp. ft1)
2029 ! the fraction of the way between jt (jt1) and the next highest
2030 ! reference temperature that the layer temperature falls.
2031 
2032  tem1 = (tavel(k) - tref(jp(k))) / 15.0
2033  tem2 = (tavel(k) - tref(jp1 )) / 15.0
2034  jt(k) = max(1, min(4, int(3.0 + tem1) ))
2035  jt1(k) = max(1, min(4, int(3.0 + tem2) ))
2036  ft = tem1 - float(jt(k) - 3)
2037  ft1 = tem2 - float(jt1(k) - 3)
2038 
2039 ! --- ... we have now isolated the layer ln pressure and temperature,
2040 ! between two reference pressures and two reference temperatures
2041 ! (for each reference pressure). we multiply the pressure
2042 ! fraction fp with the appropriate temperature fractions to get
2043 ! the factors that will be needed for the interpolation that yields
2044 ! the optical depths (performed in routines taugbn for band n).
2045 
2046  fp1 = f_one - fp
2047  fac10(k) = fp1 * ft
2048  fac00(k) = fp1 * (f_one - ft)
2049  fac11(k) = fp * ft1
2050  fac01(k) = fp * (f_one - ft1)
2051 
2052 ! --- ... if the pressure is less than ~100mb, perform a different
2053 ! set of species interpolations.
2054 
2055  if ( plog > 4.56 ) then
2056 
2057  laytrop = k
2058 
2059 ! --- ... set up factors needed to separately include the water vapor
2060 ! foreign-continuum in the calculation of absorption coefficient.
2061 
2062  tem1 = (332.0 - tavel(k)) / 36.0
2063  indfor(k) = min(2, max(1, int(tem1)))
2064  forfrac(k) = tem1 - float(indfor(k))
2065 
2066 ! --- ... set up factors needed to separately include the water vapor
2067 ! self-continuum in the calculation of absorption coefficient.
2068 
2069  tem2 = (tavel(k) - 188.0) / 7.2
2070  indself(k) = min(9, max(1, int(tem2)-7))
2071  selffrac(k) = tem2 - float(indself(k) + 7)
2072  selffac(k) = h2ovmr(k) * forfac(k)
2073 
2074  else
2075 
2076 ! --- ... set up factors needed to separately include the water vapor
2077 ! foreign-continuum in the calculation of absorption coefficient.
2078 
2079  tem1 = (tavel(k) - 188.0) / 36.0
2080  indfor(k) = 3
2081  forfrac(k) = tem1 - f_one
2082 
2083  indself(k) = 0
2084  selffrac(k) = f_zero
2085  selffac(k) = f_zero
2086 
2087  endif
2088 
2089  enddo ! end_do_k_loop
2090 
2091  return
2092 ! ..................................
2093  end subroutine setcoef
2094 ! ----------------------------------
2095 
2097 !-----------------------------------
2098  subroutine spcvrtc &
2099 !...................................
2100 ! --- inputs:
2101  & ( ssolar,cosz,sntz,albbm,albdf,sfluxzen,cldfrc, &
2102  & cf1,cf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
2103  & nlay, nlp1, &
2104 ! --- outputs:
2105  & fxupc,fxdnc,fxup0,fxdn0, &
2106  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
2107  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
2108  & )
2109 
2110 ! =================== program usage description =================== !
2111 ! !
2112 ! purpose: computes the shortwave radiative fluxes using two-stream !
2113 ! method !
2114 ! !
2115 ! subprograms called: swflux !
2116 ! !
2117 ! ==================== defination of variables ==================== !
2118 ! !
2119 ! inputs: size !
2120 ! ssolar - real, incoming solar flux at top 1 !
2121 ! cosz - real, cosine solar zenith angle 1 !
2122 ! sntz - real, secant solar zenith angle 1 !
2123 ! albbm - real, surface albedo for direct beam radiation 2 !
2124 ! albdf - real, surface albedo for diffused radiation 2 !
2125 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
2126 ! cldfrc - real, layer cloud fraction nlay !
2127 ! cf1 - real, >0: cloudy sky, otherwise: clear sky 1 !
2128 ! cf0 - real, =1-cf1 1 !
2129 ! taug - real, spectral optical depth for gases nlay*ngptsw!
2130 ! taur - real, optical depth for rayleigh scattering nlay*ngptsw!
2131 ! tauae - real, aerosols optical depth nlay*nbdsw !
2132 ! ssaae - real, aerosols single scattering albedo nlay*nbdsw !
2133 ! asyae - real, aerosols asymmetry factor nlay*nbdsw !
2134 ! taucw - real, weighted cloud optical depth nlay*nbdsw !
2135 ! ssacw - real, weighted cloud single scat albedo nlay*nbdsw !
2136 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
2137 ! nlay,nlp1 - integer, number of layers/levels 1 !
2138 ! !
2139 ! output variables: !
2140 ! fxupc - real, tot sky upward flux nlp1*nbdsw !
2141 ! fxdnc - real, tot sky downward flux nlp1*nbdsw !
2142 ! fxup0 - real, clr sky upward flux nlp1*nbdsw !
2143 ! fxdn0 - real, clr sky downward flux nlp1*nbdsw !
2144 ! ftoauc - real, tot sky toa upwd flux 1 !
2145 ! ftoau0 - real, clr sky toa upwd flux 1 !
2146 ! ftoadc - real, toa downward (incoming) solar flux 1 !
2147 ! fsfcuc - real, tot sky sfc upwd flux 1 !
2148 ! fsfcu0 - real, clr sky sfc upwd flux 1 !
2149 ! fsfcdc - real, tot sky sfc dnwd flux 1 !
2150 ! fsfcd0 - real, clr sky sfc dnwd flux 1 !
2151 ! sfbmc - real, tot sky sfc dnwd beam flux (nir/uv+vis) 2 !
2152 ! sfdfc - real, tot sky sfc dnwd diff flux (nir/uv+vis) 2 !
2153 ! sfbm0 - real, clr sky sfc dnwd beam flux (nir/uv+vis) 2 !
2154 ! sfdf0 - real, clr sky sfc dnwd diff flux (nir/uv+vis) 2 !
2155 ! suvbfc - real, tot sky sfc dnwd uv-b flux 1 !
2156 ! suvbf0 - real, clr sky sfc dnwd uv-b flux 1 !
2157 ! !
2158 ! internal variables: !
2159 ! zrefb - real, direct beam reflectivity for clear/cloudy nlp1 !
2160 ! zrefd - real, diffuse reflectivity for clear/cloudy nlp1 !
2161 ! ztrab - real, direct beam transmissivity for clear/cloudy nlp1 !
2162 ! ztrad - real, diffuse transmissivity for clear/cloudy nlp1 !
2163 ! zldbt - real, layer beam transmittance for clear/cloudy nlp1 !
2164 ! ztdbt - real, lev total beam transmittance for clr/cld nlp1 !
2165 ! !
2166 ! control parameters in module "physparam" !
2167 ! iswmode - control flag for 2-stream transfer schemes !
2168 ! = 1 delta-eddington (joseph et al., 1976) !
2169 ! = 2 pifm (zdunkowski et al., 1980) !
2170 ! = 3 discrete ordinates (liou, 1973) !
2171 ! !
2172 ! ******************************************************************* !
2173 ! original code description !
2174 ! !
2175 ! method: !
2176 ! ------- !
2177 ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations. !
2178 ! kmodts = 1 eddington (joseph et al., 1976) !
2179 ! = 2 pifm (zdunkowski et al., 1980) !
2180 ! = 3 discrete ordinates (liou, 1973) !
2181 ! !
2182 ! modifications: !
2183 ! -------------- !
2184 ! original: h. barker !
2185 ! revision: merge with rrtmg_sw: j.-j.morcrette, ecmwf, feb 2003 !
2186 ! revision: add adjustment for earth/sun distance:mjiacono,aer,oct2003!
2187 ! revision: bug fix for use of palbp and palbd: mjiacono, aer, nov2003!
2188 ! revision: bug fix to apply delta scaling to clear sky: aer, dec2004 !
2189 ! revision: code modified so that delta scaling is not done in cloudy !
2190 ! profiles if routine cldprop is used; delta scaling can be !
2191 ! applied by swithcing code below if cldprop is not used to !
2192 ! get cloud properties. aer, jan 2005 !
2193 ! revision: uniform formatting for rrtmg: mjiacono, aer, jul 2006 !
2194 ! revision: use exponential lookup table for transmittance: mjiacono, !
2195 ! aer, aug 2007 !
2196 ! !
2197 ! ******************************************************************* !
2198 ! ====================== end of description block ================= !
2199 
2200 ! --- constant parameters:
2201  real (kind=kind_phys), parameter :: zcrit = 0.9999995 ! thresold for conservative scattering
2202  real (kind=kind_phys), parameter :: zsr3 = sqrt(3.0)
2203  real (kind=kind_phys), parameter :: od_lo = 0.06
2204  real (kind=kind_phys), parameter :: eps1 = 1.0e-8
2205 
2206 ! --- inputs:
2207  integer, intent(in) :: nlay, nlp1
2208 
2209  real (kind=kind_phys), dimension(nlay,ngptsw), intent(in) :: &
2210  & taug, taur
2211  real (kind=kind_phys), dimension(nlay,nbdsw), intent(in) :: &
2212  & taucw, ssacw, asycw, tauae, ssaae, asyae
2213 
2214  real (kind=kind_phys), dimension(ngptsw), intent(in) :: sfluxzen
2215  real (kind=kind_phys), dimension(nlay), intent(in) :: cldfrc
2216 
2217  real (kind=kind_phys), dimension(2), intent(in) :: albbm, albdf
2218 
2219  real (kind=kind_phys), intent(in) :: cosz, sntz, cf1, cf0, ssolar
2220 
2221 ! --- outputs:
2222  real (kind=kind_phys), dimension(nlp1,nbdsw), intent(out) :: &
2223  & fxupc, fxdnc, fxup0, fxdn0
2224 
2225  real (kind=kind_phys), dimension(2), intent(out) :: sfbmc, sfdfc, &
2226  & sfbm0, sfdf0
2227 
2228  real (kind=kind_phys), intent(out) :: suvbfc, suvbf0, ftoadc, &
2229  & ftoauc, ftoau0, fsfcuc, fsfcu0, fsfcdc, fsfcd0
2230 
2231 ! --- locals:
2232  real (kind=kind_phys), dimension(nlay) :: ztaus, zssas, zasys, &
2233  & zldbt0
2234 
2235  real (kind=kind_phys), dimension(nlp1) :: zrefb, zrefd, ztrab, &
2236  & ztrad, ztdbt, zldbt, zfu, zfd
2237 
2238  real (kind=kind_phys) :: ztau1, zssa1, zasy1, ztau0, zssa0, &
2239  & zasy0, zasy3, zssaw, zasyw, zgam1, zgam2, zgam3, zgam4, &
2240  & zc0, zc1, za1, za2, zb1, zb2, zrk, zrk2, zrp, zrp1, zrm1, &
2241  & zrpp, zrkg1, zrkg3, zrkg4, zexp1, zexm1, zexp2, zexm2, &
2242  & zexp3, zexp4, zden1, ze1r45, ftind, zsolar, zrefb1, &
2243  & zrefd1, ztrab1, ztrad1, ztdbt0, zr1, zr2, zr3, zr4, zr5, &
2244  & zt1, zt2, zt3, zf1, zf2
2245 
2246  integer :: ib, ibd, jb, jg, k, kp, itind
2247 !
2248 !===> ... begin here
2249 !
2250 ! --- ... initialization of output fluxes
2251 
2252  do ib = 1, nbdsw
2253  do k = 1, nlp1
2254  fxdnc(k,ib) = f_zero
2255  fxupc(k,ib) = f_zero
2256  fxdn0(k,ib) = f_zero
2257  fxup0(k,ib) = f_zero
2258  enddo
2259  enddo
2260 
2261  ftoadc = f_zero
2262  ftoauc = f_zero
2263  ftoau0 = f_zero
2264  fsfcuc = f_zero
2265  fsfcu0 = f_zero
2266  fsfcdc = f_zero
2267  fsfcd0 = f_zero
2268 
2269 !! --- ... uv-b surface downward fluxes
2270  suvbfc = f_zero
2271  suvbf0 = f_zero
2272 
2273 !! --- ... output surface flux components
2274  sfbmc(1) = f_zero
2275  sfbmc(2) = f_zero
2276  sfdfc(1) = f_zero
2277  sfdfc(2) = f_zero
2278  sfbm0(1) = f_zero
2279  sfbm0(2) = f_zero
2280  sfdf0(1) = f_zero
2281  sfdf0(2) = f_zero
2282 
2283 ! --- ... loop over all g-points in each band
2284 
2285  lab_do_jg : do jg = 1, ngptsw
2286 
2287  jb = ngb(jg)
2288  ib = jb + 1 - nblow
2289  ibd = idxsfc(jb)
2290 
2291  zsolar = ssolar * sfluxzen(jg)
2292 
2293 ! --- ... set up toa direct beam and surface values (beam and diff)
2294 
2295  ztdbt(nlp1) = f_one
2296  ztdbt0 = f_one
2297 
2298  zldbt(1) = f_zero
2299  if (ibd /= 0) then
2300  zrefb(1) = albbm(ibd)
2301  zrefd(1) = albdf(ibd)
2302  else
2303  zrefb(1) = 0.5 * (albbm(1) + albbm(2))
2304  zrefd(1) = 0.5 * (albdf(1) + albdf(2))
2305  endif
2306  ztrab(1) = f_zero
2307  ztrad(1) = f_zero
2308 
2309 ! --- ... compute clear-sky optical parameters, layer reflectance and transmittance
2310 
2311  do k = nlay, 1, -1
2312  kp = k + 1
2313 
2314  ztau0 = max( ftiny, taur(k,jg)+taug(k,jg)+tauae(k,ib) )
2315  zssa0 = taur(k,jg) + tauae(k,ib)*ssaae(k,ib)
2316  zasy0 = asyae(k,ib)*ssaae(k,ib)*tauae(k,ib)
2317  zssaw = min( oneminus, zssa0 / ztau0 )
2318  zasyw = zasy0 / max( ftiny, zssa0 )
2319 
2320 ! --- ... saving clear-sky quantities for later total-sky usage
2321  ztaus(k) = ztau0
2322  zssas(k) = zssa0
2323  zasys(k) = zasy0
2324 
2325 ! --- ... delta scaling for clear-sky condition
2326  za1 = zasyw * zasyw
2327  za2 = zssaw * za1
2328 
2329  ztau1 = (f_one - za2) * ztau0
2330  zssa1 = (zssaw - za2) / (f_one - za2)
2331 !org zasy1 = (zasyw - za1) / (f_one - za1) ! this line is replaced by the next
2332  zasy1 = zasyw / (f_one + zasyw) ! to reduce truncation error
2333  zasy3 = 0.75 * zasy1
2334 
2335 ! --- ... general two-stream expressions
2336  if ( iswmode == 1 ) then
2337  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
2338  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
2339  zgam3 = 0.5 - zasy3 * cosz
2340  elseif ( iswmode == 2 ) then ! pifm
2341  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
2342  zgam2 = 0.75* zssa1 * (f_one- zasy1)
2343  zgam3 = 0.5 - zasy3 * cosz
2344  elseif ( iswmode == 3 ) then ! discrete ordinates
2345  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
2346  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
2347  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
2348  endif
2349  zgam4 = f_one - zgam3
2350 
2351 ! --- ... compute homogeneous reflectance and transmittance
2352 
2353  if ( zssaw >= zcrit ) then ! for conservative scattering
2354  za1 = zgam1 * cosz - zgam3
2355  za2 = zgam1 * ztau1
2356 
2357 ! --- ... use exponential lookup table for transmittance, or expansion
2358 ! of exponential for low optical depth
2359 
2360  zb1 = min( ztau1*sntz , 500.0 )
2361  if ( zb1 <= od_lo ) then
2362  zb2 = f_one - zb1 + 0.5*zb1*zb1
2363  else
2364  ftind = zb1 / (bpade + zb1)
2365  itind = ftind*ntbmx + 0.5
2366  zb2 = exp_tbl(itind)
2367  endif
2368 
2369 ! ... collimated beam
2370  zrefb(kp) = max(f_zero, min(f_one, &
2371  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
2372  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp) ))
2373 
2374 ! ... isotropic incidence
2375  zrefd(kp) = max(f_zero, min(f_one, za2/(f_one + za2) ))
2376  ztrad(kp) = max(f_zero, min(f_one, f_one-zrefd(kp) ))
2377 
2378  else ! for non-conservative scattering
2379  za1 = zgam1*zgam4 + zgam2*zgam3
2380  za2 = zgam1*zgam3 + zgam2*zgam4
2381  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
2382  zrk2= 2.0 * zrk
2383 
2384  zrp = zrk * cosz
2385  zrp1 = f_one + zrp
2386  zrm1 = f_one - zrp
2387  zrpp = f_one - zrp*zrp
2388  zrkg1= zrk + zgam1
2389  zrkg3= zrk * zgam3
2390  zrkg4= zrk * zgam4
2391 
2392  zr1 = zrm1 * (za2 + zrkg3)
2393  zr2 = zrp1 * (za2 - zrkg3)
2394  zr3 = zrk2 * (zgam3 - za2*cosz)
2395  zr4 = zrpp * zrkg1
2396  zr5 = zrpp * (zrk - zgam1)
2397 
2398  zt1 = zrp1 * (za1 + zrkg4)
2399  zt2 = zrm1 * (za1 - zrkg4)
2400  zt3 = zrk2 * (zgam4 + za1*cosz)
2401 
2402 ! --- ... use exponential lookup table for transmittance, or expansion
2403 ! of exponential for low optical depth
2404 
2405  zb1 = min( zrk*ztau1, 500.0 )
2406  if ( zb1 <= od_lo ) then
2407  zexm1 = f_one - zb1 + 0.5*zb1*zb1
2408  else
2409  ftind = zb1 / (bpade + zb1)
2410  itind = ftind*ntbmx + 0.5
2411  zexm1 = exp_tbl(itind)
2412  endif
2413  zexp1 = f_one / zexm1
2414 
2415  zb2 = min( sntz*ztau1, 500.0 )
2416  if ( zb2 <= od_lo ) then
2417  zexm2 = f_one - zb2 + 0.5*zb2*zb2
2418  else
2419  ftind = zb2 / (bpade + zb2)
2420  itind = ftind*ntbmx + 0.5
2421  zexm2 = exp_tbl(itind)
2422  endif
2423  zexp2 = f_one / zexm2
2424  ze1r45 = zr4*zexp1 + zr5*zexm1
2425 
2426 ! ... collimated beam
2427  if (ze1r45>=-eps1 .and. ze1r45<=eps1) then
2428  zrefb(kp) = eps1
2429  ztrab(kp) = zexm2
2430  else
2431  zden1 = zssa1 / ze1r45
2432  zrefb(kp) = max(f_zero, min(f_one, &
2433  & (zr1*zexp1 - zr2*zexm1 - zr3*zexm2)*zden1 ))
2434  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one &
2435  & - (zt1*zexp1 - zt2*zexm1 - zt3*zexp2)*zden1) ))
2436  endif
2437 
2438 ! ... diffuse beam
2439  zden1 = zr4 / (ze1r45 * zrkg1)
2440  zrefd(kp) = max(f_zero, min(f_one, &
2441  & zgam2*(zexp1 - zexm1)*zden1 ))
2442  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
2443  endif ! end if_zssaw_block
2444 
2445 ! --- ... direct beam transmittance. use exponential lookup table
2446 ! for transmittance, or expansion of exponential for low
2447 ! optical depth
2448 
2449  zr1 = ztau1 * sntz
2450  if ( zr1 <= od_lo ) then
2451  zexp3 = f_one - zr1 + 0.5*zr1*zr1
2452  else
2453  ftind = zr1 / (bpade + zr1)
2454  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2455  zexp3 = exp_tbl(itind)
2456  endif
2457 
2458  ztdbt(k) = zexp3 * ztdbt(kp)
2459  zldbt(kp) = zexp3
2460 
2461 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2462 ! (must use 'orig', unscaled cloud optical depth)
2463 
2464  zr1 = ztau0 * sntz
2465  if ( zr1 <= od_lo ) then
2466  zexp4 = f_one - zr1 + 0.5*zr1*zr1
2467  else
2468  ftind = zr1 / (bpade + zr1)
2469  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2470  zexp4 = exp_tbl(itind)
2471  endif
2472 
2473  zldbt0(k) = zexp4
2474  ztdbt0 = zexp4 * ztdbt0
2475  enddo ! end do_k_loop
2476 
2477  call swflux &
2478 ! --- inputs:
2479  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
2480  & nlay, nlp1, &
2481 ! --- outputs:
2482  & zfu, zfd &
2483  & )
2484 
2485 ! --- ... compute upward and downward fluxes at levels
2486  do k = 1, nlp1
2487  fxup0(k,ib) = fxup0(k,ib) + zsolar*zfu(k)
2488  fxdn0(k,ib) = fxdn0(k,ib) + zsolar*zfd(k)
2489  enddo
2490 
2491 !! --- ... surface downward beam/diffused flux components
2492  zb1 = zsolar*ztdbt0
2493  zb2 = zsolar*(zfd(1) - ztdbt0)
2494 
2495  if (ibd /= 0) then
2496  sfbm0(ibd) = sfbm0(ibd) + zb1
2497  sfdf0(ibd) = sfdf0(ibd) + zb2
2498  else
2499  zf1 = 0.5 * zb1
2500  zf2 = 0.5 * zb2
2501  sfbm0(1) = sfbm0(1) + zf1
2502  sfdf0(1) = sfdf0(1) + zf2
2503  sfbm0(2) = sfbm0(2) + zf1
2504  sfdf0(2) = sfdf0(2) + zf2
2505  endif
2506 ! sfbm0(ibd) = sfbm0(ibd) + zsolar*ztdbt0
2507 ! sfdf0(ibd) = sfdf0(ibd) + zsolar*(zfd(1) - ztdbt0)
2508 
2509 ! --- ... compute total sky optical parameters, layer reflectance and transmittance
2510 
2511  if ( cf1 > eps ) then
2512 
2513 ! --- ... set up toa direct beam and surface values (beam and diff)
2514  ztdbt0 = f_one
2515  zldbt(1) = f_zero
2516 
2517  do k = nlay, 1, -1
2518  kp = k + 1
2519  zc0 = f_one - cldfrc(k)
2520  zc1 = cldfrc(k)
2521  if ( zc1 > ftiny ) then ! it is a cloudy-layer
2522 
2523  ztau0 = ztaus(k) + taucw(k,ib)
2524  zssa0 = zssas(k) + ssacw(k,ib)
2525  zasy0 = zasys(k) + asycw(k,ib)
2526  zssaw = min(oneminus, zssa0 / ztau0)
2527  zasyw = zasy0 / max(ftiny, zssa0)
2528 
2529 ! --- ... delta scaling for total-sky condition
2530  za1 = zasyw * zasyw
2531  za2 = zssaw * za1
2532 
2533  ztau1 = (f_one - za2) * ztau0
2534  zssa1 = (zssaw - za2) / (f_one - za2)
2535 !org zasy1 = (zasyw - za1) / (f_one - za1)
2536  zasy1 = zasyw / (f_one + zasyw)
2537  zasy3 = 0.75 * zasy1
2538 
2539 ! --- ... general two-stream expressions
2540  if ( iswmode == 1 ) then
2541  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
2542  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
2543  zgam3 = 0.5 - zasy3 * cosz
2544  elseif ( iswmode == 2 ) then ! pifm
2545  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
2546  zgam2 = 0.75* zssa1 * (f_one- zasy1)
2547  zgam3 = 0.5 - zasy3 * cosz
2548  elseif ( iswmode == 3 ) then ! discrete ordinates
2549  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
2550  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
2551  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
2552  endif
2553  zgam4 = f_one - zgam3
2554 
2555  zrefb1 = zrefb(kp)
2556  zrefd1 = zrefd(kp)
2557  ztrab1 = ztrab(kp)
2558  ztrad1 = ztrad(kp)
2559 
2560 ! --- ... compute homogeneous reflectance and transmittance
2561 
2562  if ( zssaw >= zcrit ) then ! for conservative scattering
2563  za1 = zgam1 * cosz - zgam3
2564  za2 = zgam1 * ztau1
2565 
2566 ! --- ... use exponential lookup table for transmittance, or expansion
2567 ! of exponential for low optical depth
2568 
2569  zb1 = min( ztau1*sntz , 500.0 )
2570  if ( zb1 <= od_lo ) then
2571  zb2 = f_one - zb1 + 0.5*zb1*zb1
2572  else
2573  ftind = zb1 / (bpade + zb1)
2574  itind = ftind*ntbmx + 0.5
2575  zb2 = exp_tbl(itind)
2576  endif
2577 
2578 ! ... collimated beam
2579  zrefb(kp) = max(f_zero, min(f_one, &
2580  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
2581  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp)))
2582 
2583 ! ... isotropic incidence
2584  zrefd(kp) = max(f_zero, min(f_one, za2 / (f_one+za2) ))
2585  ztrad(kp) = max(f_zero, min(f_one, f_one - zrefd(kp) ))
2586 
2587  else ! for non-conservative scattering
2588  za1 = zgam1*zgam4 + zgam2*zgam3
2589  za2 = zgam1*zgam3 + zgam2*zgam4
2590  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
2591  zrk2= 2.0 * zrk
2592 
2593  zrp = zrk * cosz
2594  zrp1 = f_one + zrp
2595  zrm1 = f_one - zrp
2596  zrpp = f_one - zrp*zrp
2597  zrkg1= zrk + zgam1
2598  zrkg3= zrk * zgam3
2599  zrkg4= zrk * zgam4
2600 
2601  zr1 = zrm1 * (za2 + zrkg3)
2602  zr2 = zrp1 * (za2 - zrkg3)
2603  zr3 = zrk2 * (zgam3 - za2*cosz)
2604  zr4 = zrpp * zrkg1
2605  zr5 = zrpp * (zrk - zgam1)
2606 
2607  zt1 = zrp1 * (za1 + zrkg4)
2608  zt2 = zrm1 * (za1 - zrkg4)
2609  zt3 = zrk2 * (zgam4 + za1*cosz)
2610 
2611 ! --- ... use exponential lookup table for transmittance, or expansion
2612 ! of exponential for low optical depth
2613 
2614  zb1 = min( zrk*ztau1, 500.0 )
2615  if ( zb1 <= od_lo ) then
2616  zexm1 = f_one - zb1 + 0.5*zb1*zb1
2617  else
2618  ftind = zb1 / (bpade + zb1)
2619  itind = ftind*ntbmx + 0.5
2620  zexm1 = exp_tbl(itind)
2621  endif
2622  zexp1 = f_one / zexm1
2623 
2624  zb2 = min( ztau1*sntz, 500.0 )
2625  if ( zb2 <= od_lo ) then
2626  zexm2 = f_one - zb2 + 0.5*zb2*zb2
2627  else
2628  ftind = zb2 / (bpade + zb2)
2629  itind = ftind*ntbmx + 0.5
2630  zexm2 = exp_tbl(itind)
2631  endif
2632  zexp2 = f_one / zexm2
2633  ze1r45 = zr4*zexp1 + zr5*zexm1
2634 
2635 ! ... collimated beam
2636  if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then
2637  zrefb(kp) = eps1
2638  ztrab(kp) = zexm2
2639  else
2640  zden1 = zssa1 / ze1r45
2641  zrefb(kp) = max(f_zero, min(f_one, &
2642  & (zr1*zexp1-zr2*zexm1-zr3*zexm2)*zden1 ))
2643  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one - &
2644  & (zt1*zexp1-zt2*zexm1-zt3*zexp2)*zden1) ))
2645  endif
2646 
2647 ! ... diffuse beam
2648  zden1 = zr4 / (ze1r45 * zrkg1)
2649  zrefd(kp) = max(f_zero, min(f_one, &
2650  & zgam2*(zexp1 - zexm1)*zden1 ))
2651  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
2652  endif ! end if_zssaw_block
2653 
2654 ! --- ... combine clear and cloudy contributions for total sky
2655 ! and calculate direct beam transmittances
2656 
2657  zrefb(kp) = zc0*zrefb1 + zc1*zrefb(kp)
2658  zrefd(kp) = zc0*zrefd1 + zc1*zrefd(kp)
2659  ztrab(kp) = zc0*ztrab1 + zc1*ztrab(kp)
2660  ztrad(kp) = zc0*ztrad1 + zc1*ztrad(kp)
2661 
2662 ! --- ... direct beam transmittance. use exponential lookup table
2663 ! for transmittance, or expansion of exponential for low
2664 ! optical depth
2665 
2666  zr1 = ztau1 * sntz
2667  if ( zr1 <= od_lo ) then
2668  zexp3 = f_one - zr1 + 0.5*zr1*zr1
2669  else
2670  ftind = zr1 / (bpade + zr1)
2671  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2672  zexp3 = exp_tbl(itind)
2673  endif
2674 
2675  zldbt(kp) = zc0*zldbt(kp) + zc1*zexp3
2676  ztdbt(k) = zldbt(kp) * ztdbt(kp)
2677 
2678 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2679 ! (must use 'orig', unscaled cloud optical depth)
2680 
2681  zr1 = ztau0 * sntz
2682  if ( zr1 <= od_lo ) then
2683  zexp4 = f_one - zr1 + 0.5*zr1*zr1
2684  else
2685  ftind = zr1 / (bpade + zr1)
2686  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2687  zexp4 = exp_tbl(itind)
2688  endif
2689 
2690  ztdbt0 = (zc0*zldbt0(k) + zc1*zexp4) * ztdbt0
2691 
2692  else ! if_zc1_block --- it is a clear layer
2693 
2694 ! --- ... direct beam transmittance
2695  ztdbt(k) = zldbt(kp) * ztdbt(kp)
2696 
2697 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2698  ztdbt0 = zldbt0(k) * ztdbt0
2699 
2700  endif ! end if_zc1_block
2701  enddo ! end do_k_loop
2702 
2703 ! --- ... perform vertical quadrature
2704 
2705  call swflux &
2706 ! --- inputs:
2707  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
2708  & nlay, nlp1, &
2709 ! --- outputs:
2710  & zfu, zfd &
2711  & )
2712 
2713 ! --- ... compute upward and downward fluxes at levels
2714  do k = 1, nlp1
2715  fxupc(k,ib) = fxupc(k,ib) + zsolar*zfu(k)
2716  fxdnc(k,ib) = fxdnc(k,ib) + zsolar*zfd(k)
2717  enddo
2718 
2719 !! --- ... surface downward beam/diffused flux components
2720  zb1 = zsolar*ztdbt0
2721  zb2 = zsolar*(zfd(1) - ztdbt0)
2722 
2723  if (ibd /= 0) then
2724  sfbmc(ibd) = sfbmc(ibd) + zb1
2725  sfdfc(ibd) = sfdfc(ibd) + zb2
2726  else
2727  zf1 = 0.5 * zb1
2728  zf2 = 0.5 * zb2
2729  sfbmc(1) = sfbmc(1) + zf1
2730  sfdfc(1) = sfdfc(1) + zf2
2731  sfbmc(2) = sfbmc(2) + zf1
2732  sfdfc(2) = sfdfc(2) + zf2
2733  endif
2734 ! sfbmc(ibd) = sfbmc(ibd) + zsolar*ztdbt0
2735 ! sfdfc(ibd) = sfdfc(ibd) + zsolar*(zfd(1) - ztdbt0)
2736 
2737  endif ! end if_cf1_block
2738 
2739  enddo lab_do_jg
2740 
2741 ! --- ... end of g-point loop
2742 
2743  do ib = 1, nbdsw
2744  ftoadc = ftoadc + fxdn0(nlp1,ib)
2745  ftoau0 = ftoau0 + fxup0(nlp1,ib)
2746  fsfcu0 = fsfcu0 + fxup0(1,ib)
2747  fsfcd0 = fsfcd0 + fxdn0(1,ib)
2748  enddo
2749 
2750 !! --- ... uv-b surface downward flux
2751  ibd = nuvb - nblow + 1
2752  suvbf0 = fxdn0(1,ibd)
2753 
2754  if ( cf1 <= eps ) then ! clear column, set total-sky=clear-sky fluxes
2755  do ib = 1, nbdsw
2756  do k = 1, nlp1
2757  fxupc(k,ib) = fxup0(k,ib)
2758  fxdnc(k,ib) = fxdn0(k,ib)
2759  enddo
2760  enddo
2761 
2762  ftoauc = ftoau0
2763  fsfcuc = fsfcu0
2764  fsfcdc = fsfcd0
2765 
2766 !! --- ... surface downward beam/diffused flux components
2767  sfbmc(1) = sfbm0(1)
2768  sfdfc(1) = sfdf0(1)
2769  sfbmc(2) = sfbm0(2)
2770  sfdfc(2) = sfdf0(2)
2771 
2772 !! --- ... uv-b surface downward flux
2773  suvbfc = suvbf0
2774  else ! cloudy column, compute total-sky fluxes
2775  do ib = 1, nbdsw
2776  do k = 1, nlp1
2777  fxupc(k,ib) = cf1*fxupc(k,ib) + cf0*fxup0(k,ib)
2778  fxdnc(k,ib) = cf1*fxdnc(k,ib) + cf0*fxdn0(k,ib)
2779  enddo
2780  enddo
2781 
2782  do ib = 1, nbdsw
2783  ftoauc = ftoauc + fxupc(nlp1,ib)
2784  fsfcuc = fsfcuc + fxupc(1,ib)
2785  fsfcdc = fsfcdc + fxdnc(1,ib)
2786  enddo
2787 
2788 !! --- ... uv-b surface downward flux
2789  suvbfc = fxdnc(1,ibd)
2790 
2791 !! --- ... surface downward beam/diffused flux components
2792  sfbmc(1) = cf1*sfbmc(1) + cf0*sfbm0(1)
2793  sfbmc(2) = cf1*sfbmc(2) + cf0*sfbm0(2)
2794  sfdfc(1) = cf1*sfdfc(1) + cf0*sfdf0(1)
2795  sfdfc(2) = cf1*sfdfc(2) + cf0*sfdf0(2)
2796  endif ! end if_cf1_block
2797 
2798  return
2799 !...................................
2800  end subroutine spcvrtc
2801 !-----------------------------------
2802 
2806 !-----------------------------------
2807  subroutine spcvrtm &
2808 !...................................
2809 ! --- inputs:
2810  & ( ssolar,cosz,sntz,albbm,albdf,sfluxzen,cldfmc, &
2811  & cf1,cf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
2812  & nlay, nlp1, &
2813 ! --- outputs:
2814  & fxupc,fxdnc,fxup0,fxdn0, &
2815  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
2816  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
2817  & )
2818 
2819 ! =================== program usage description =================== !
2820 ! !
2821 ! purpose: computes the shortwave radiative fluxes using two-stream !
2822 ! method of h. barker and mcica, the monte-carlo independent!
2823 ! column approximation, for the representation of sub-grid !
2824 ! cloud variability (i.e. cloud overlap). !
2825 ! !
2826 ! subprograms called: swflux !
2827 ! !
2828 ! ==================== defination of variables ==================== !
2829 ! !
2830 ! inputs: size !
2831 ! ssolar - real, incoming solar flux at top 1 !
2832 ! cosz - real, cosine solar zenith angle 1 !
2833 ! sntz - real, secant solar zenith angle 1 !
2834 ! albbm - real, surface albedo for direct beam radiation 2 !
2835 ! albdf - real, surface albedo for diffused radiation 2 !
2836 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
2837 ! cldfmc - real, layer cloud fraction for g-point nlay*ngptsw!
2838 ! cf1 - real, >0: cloudy sky, otherwise: clear sky 1 !
2839 ! cf0 - real, =1-cf1 1 !
2840 ! taug - real, spectral optical depth for gases nlay*ngptsw!
2841 ! taur - real, optical depth for rayleigh scattering nlay*ngptsw!
2842 ! tauae - real, aerosols optical depth nlay*nbdsw !
2843 ! ssaae - real, aerosols single scattering albedo nlay*nbdsw !
2844 ! asyae - real, aerosols asymmetry factor nlay*nbdsw !
2845 ! taucw - real, weighted cloud optical depth nlay*nbdsw !
2846 ! ssacw - real, weighted cloud single scat albedo nlay*nbdsw !
2847 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
2848 ! nlay,nlp1 - integer, number of layers/levels 1 !
2849 ! !
2850 ! output variables: !
2851 ! fxupc - real, tot sky upward flux nlp1*nbdsw !
2852 ! fxdnc - real, tot sky downward flux nlp1*nbdsw !
2853 ! fxup0 - real, clr sky upward flux nlp1*nbdsw !
2854 ! fxdn0 - real, clr sky downward flux nlp1*nbdsw !
2855 ! ftoauc - real, tot sky toa upwd flux 1 !
2856 ! ftoau0 - real, clr sky toa upwd flux 1 !
2857 ! ftoadc - real, toa downward (incoming) solar flux 1 !
2858 ! fsfcuc - real, tot sky sfc upwd flux 1 !
2859 ! fsfcu0 - real, clr sky sfc upwd flux 1 !
2860 ! fsfcdc - real, tot sky sfc dnwd flux 1 !
2861 ! fsfcd0 - real, clr sky sfc dnwd flux 1 !
2862 ! sfbmc - real, tot sky sfc dnwd beam flux (nir/uv+vis) 2 !
2863 ! sfdfc - real, tot sky sfc dnwd diff flux (nir/uv+vis) 2 !
2864 ! sfbm0 - real, clr sky sfc dnwd beam flux (nir/uv+vis) 2 !
2865 ! sfdf0 - real, clr sky sfc dnwd diff flux (nir/uv+vis) 2 !
2866 ! suvbfc - real, tot sky sfc dnwd uv-b flux 1 !
2867 ! suvbf0 - real, clr sky sfc dnwd uv-b flux 1 !
2868 ! !
2869 ! internal variables: !
2870 ! zrefb - real, direct beam reflectivity for clear/cloudy nlp1 !
2871 ! zrefd - real, diffuse reflectivity for clear/cloudy nlp1 !
2872 ! ztrab - real, direct beam transmissivity for clear/cloudy nlp1 !
2873 ! ztrad - real, diffuse transmissivity for clear/cloudy nlp1 !
2874 ! zldbt - real, layer beam transmittance for clear/cloudy nlp1 !
2875 ! ztdbt - real, lev total beam transmittance for clr/cld nlp1 !
2876 ! !
2877 ! control parameters in module "physparam" !
2878 ! iswmode - control flag for 2-stream transfer schemes !
2879 ! = 1 delta-eddington (joseph et al., 1976) !
2880 ! = 2 pifm (zdunkowski et al., 1980) !
2881 ! = 3 discrete ordinates (liou, 1973) !
2882 ! !
2883 ! ******************************************************************* !
2884 ! original code description !
2885 ! !
2886 ! method: !
2887 ! ------- !
2888 ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations. !
2889 ! kmodts = 1 eddington (joseph et al., 1976) !
2890 ! = 2 pifm (zdunkowski et al., 1980) !
2891 ! = 3 discrete ordinates (liou, 1973) !
2892 ! !
2893 ! modifications: !
2894 ! -------------- !
2895 ! original: h. barker !
2896 ! revision: merge with rrtmg_sw: j.-j.morcrette, ecmwf, feb 2003 !
2897 ! revision: add adjustment for earth/sun distance:mjiacono,aer,oct2003!
2898 ! revision: bug fix for use of palbp and palbd: mjiacono, aer, nov2003!
2899 ! revision: bug fix to apply delta scaling to clear sky: aer, dec2004 !
2900 ! revision: code modified so that delta scaling is not done in cloudy !
2901 ! profiles if routine cldprop is used; delta scaling can be !
2902 ! applied by swithcing code below if cldprop is not used to !
2903 ! get cloud properties. aer, jan 2005 !
2904 ! revision: uniform formatting for rrtmg: mjiacono, aer, jul 2006 !
2905 ! revision: use exponential lookup table for transmittance: mjiacono, !
2906 ! aer, aug 2007 !
2907 ! !
2908 ! ******************************************************************* !
2909 ! ====================== end of description block ================= !
2910 
2911 ! --- constant parameters:
2912  real (kind=kind_phys), parameter :: zcrit = 0.9999995 ! thresold for conservative scattering
2913  real (kind=kind_phys), parameter :: zsr3 = sqrt(3.0)
2914  real (kind=kind_phys), parameter :: od_lo = 0.06
2915  real (kind=kind_phys), parameter :: eps1 = 1.0e-8
2916 
2917 ! --- inputs:
2918  integer, intent(in) :: nlay, nlp1
2919 
2920  real (kind=kind_phys), dimension(nlay,ngptsw), intent(in) :: &
2921  & taug, taur, cldfmc
2922  real (kind=kind_phys), dimension(nlay,nbdsw), intent(in) :: &
2923  & taucw, ssacw, asycw, tauae, ssaae, asyae
2924 
2925  real (kind=kind_phys), dimension(ngptsw), intent(in) :: sfluxzen
2926 
2927  real (kind=kind_phys), dimension(2), intent(in) :: albbm, albdf
2928 
2929  real (kind=kind_phys), intent(in) :: cosz, sntz, cf1, cf0, ssolar
2930 
2931 ! --- outputs:
2932  real (kind=kind_phys), dimension(nlp1,nbdsw), intent(out) :: &
2933  & fxupc, fxdnc, fxup0, fxdn0
2934 
2935  real (kind=kind_phys), dimension(2), intent(out) :: sfbmc, sfdfc, &
2936  & sfbm0, sfdf0
2937 
2938  real (kind=kind_phys), intent(out) :: suvbfc, suvbf0, ftoadc, &
2939  & ftoauc, ftoau0, fsfcuc, fsfcu0, fsfcdc, fsfcd0
2940 
2941 ! --- locals:
2942  real (kind=kind_phys), dimension(nlay) :: ztaus, zssas, zasys, &
2943  & zldbt0
2944 
2945  real (kind=kind_phys), dimension(nlp1) :: zrefb, zrefd, ztrab, &
2946  & ztrad, ztdbt, zldbt, zfu, zfd
2947 
2948  real (kind=kind_phys) :: ztau1, zssa1, zasy1, ztau0, zssa0, &
2949  & zasy0, zasy3, zssaw, zasyw, zgam1, zgam2, zgam3, zgam4, &
2950  & za1, za2, zb1, zb2, zrk, zrk2, zrp, zrp1, zrm1, zrpp, &
2951  & zrkg1, zrkg3, zrkg4, zexp1, zexm1, zexp2, zexm2, zden1, &
2952  & zexp3, zexp4, ze1r45, ftind, zsolar, ztdbt0, zr1, zr2, &
2953  & zr3, zr4, zr5, zt1, zt2, zt3, zf1, zf2
2954 
2955  integer :: ib, ibd, jb, jg, k, kp, itind
2956 !
2957 !===> ... begin here
2958 !
2959 ! --- ... initialization of output fluxes
2960 
2961  do ib = 1, nbdsw
2962  do k = 1, nlp1
2963  fxdnc(k,ib) = f_zero
2964  fxupc(k,ib) = f_zero
2965  fxdn0(k,ib) = f_zero
2966  fxup0(k,ib) = f_zero
2967  enddo
2968  enddo
2969 
2970  ftoadc = f_zero
2971  ftoauc = f_zero
2972  ftoau0 = f_zero
2973  fsfcuc = f_zero
2974  fsfcu0 = f_zero
2975  fsfcdc = f_zero
2976  fsfcd0 = f_zero
2977 
2978 !! --- ... uv-b surface downward fluxes
2979  suvbfc = f_zero
2980  suvbf0 = f_zero
2981 
2982 !! --- ... output surface flux components
2983  sfbmc(1) = f_zero
2984  sfbmc(2) = f_zero
2985  sfdfc(1) = f_zero
2986  sfdfc(2) = f_zero
2987  sfbm0(1) = f_zero
2988  sfbm0(2) = f_zero
2989  sfdf0(1) = f_zero
2990  sfdf0(2) = f_zero
2991 
2992 ! --- ... loop over all g-points in each band
2993 
2994  lab_do_jg : do jg = 1, ngptsw
2995 
2996  jb = ngb(jg)
2997  ib = jb + 1 - nblow
2998  ibd = idxsfc(jb) ! spectral band index
2999 
3000  zsolar = ssolar * sfluxzen(jg)
3001 
3002 ! --- ... set up toa direct beam and surface values (beam and diff)
3003 
3004  ztdbt(nlp1) = f_one
3005  ztdbt0 = f_one
3006 
3007  zldbt(1) = f_zero
3008  if (ibd /= 0) then
3009  zrefb(1) = albbm(ibd)
3010  zrefd(1) = albdf(ibd)
3011  else
3012  zrefb(1) = 0.5 * (albbm(1) + albbm(2))
3013  zrefd(1) = 0.5 * (albdf(1) + albdf(2))
3014  endif
3015  ztrab(1) = f_zero
3016  ztrad(1) = f_zero
3017 
3018 ! --- ... compute clear-sky optical parameters, layer reflectance and transmittance
3019 
3020  do k = nlay, 1, -1
3021  kp = k + 1
3022 
3023  ztau0 = max( ftiny, taur(k,jg)+taug(k,jg)+tauae(k,ib) )
3024  zssa0 = taur(k,jg) + tauae(k,ib)*ssaae(k,ib)
3025  zasy0 = asyae(k,ib)*ssaae(k,ib)*tauae(k,ib)
3026  zssaw = min( oneminus, zssa0 / ztau0 )
3027  zasyw = zasy0 / max( ftiny, zssa0 )
3028 
3029 ! --- ... saving clear-sky quantities for later total-sky usage
3030  ztaus(k) = ztau0
3031  zssas(k) = zssa0
3032  zasys(k) = zasy0
3033 
3034 ! --- ... delta scaling for clear-sky condition
3035  za1 = zasyw * zasyw
3036  za2 = zssaw * za1
3037 
3038  ztau1 = (f_one - za2) * ztau0
3039  zssa1 = (zssaw - za2) / (f_one - za2)
3040 !org zasy1 = (zasyw - za1) / (f_one - za1) ! this line is replaced by the next
3041  zasy1 = zasyw / (f_one + zasyw) ! to reduce truncation error
3042  zasy3 = 0.75 * zasy1
3043 
3044 ! --- ... general two-stream expressions
3045  if ( iswmode == 1 ) then
3046  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
3047  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
3048  zgam3 = 0.5 - zasy3 * cosz
3049  elseif ( iswmode == 2 ) then ! pifm
3050  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
3051  zgam2 = 0.75* zssa1 * (f_one- zasy1)
3052  zgam3 = 0.5 - zasy3 * cosz
3053  elseif ( iswmode == 3 ) then ! discrete ordinates
3054  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
3055  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
3056  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
3057  endif
3058  zgam4 = f_one - zgam3
3059 
3060 ! --- ... compute homogeneous reflectance and transmittance
3061 
3062  if ( zssaw >= zcrit ) then ! for conservative scattering
3063  za1 = zgam1 * cosz - zgam3
3064  za2 = zgam1 * ztau1
3065 
3066 ! --- ... use exponential lookup table for transmittance, or expansion
3067 ! of exponential for low optical depth
3068 
3069  zb1 = min( ztau1*sntz , 500.0 )
3070  if ( zb1 <= od_lo ) then
3071  zb2 = f_one - zb1 + 0.5*zb1*zb1
3072  else
3073  ftind = zb1 / (bpade + zb1)
3074  itind = ftind*ntbmx + 0.5
3075  zb2 = exp_tbl(itind)
3076  endif
3077 
3078 ! ... collimated beam
3079  zrefb(kp) = max(f_zero, min(f_one, &
3080  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
3081  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp) ))
3082 
3083 ! ... isotropic incidence
3084  zrefd(kp) = max(f_zero, min(f_one, za2/(f_one + za2) ))
3085  ztrad(kp) = max(f_zero, min(f_one, f_one-zrefd(kp) ))
3086 
3087  else ! for non-conservative scattering
3088  za1 = zgam1*zgam4 + zgam2*zgam3
3089  za2 = zgam1*zgam3 + zgam2*zgam4
3090  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
3091  zrk2= 2.0 * zrk
3092 
3093  zrp = zrk * cosz
3094  zrp1 = f_one + zrp
3095  zrm1 = f_one - zrp
3096  zrpp = f_one - zrp*zrp
3097  zrkg1= zrk + zgam1
3098  zrkg3= zrk * zgam3
3099  zrkg4= zrk * zgam4
3100 
3101  zr1 = zrm1 * (za2 + zrkg3)
3102  zr2 = zrp1 * (za2 - zrkg3)
3103  zr3 = zrk2 * (zgam3 - za2*cosz)
3104  zr4 = zrpp * zrkg1
3105  zr5 = zrpp * (zrk - zgam1)
3106 
3107  zt1 = zrp1 * (za1 + zrkg4)
3108  zt2 = zrm1 * (za1 - zrkg4)
3109  zt3 = zrk2 * (zgam4 + za1*cosz)
3110 
3111 ! --- ... use exponential lookup table for transmittance, or expansion
3112 ! of exponential for low optical depth
3113 
3114  zb1 = min( zrk*ztau1, 500.0 )
3115  if ( zb1 <= od_lo ) then
3116  zexm1 = f_one - zb1 + 0.5*zb1*zb1
3117  else
3118  ftind = zb1 / (bpade + zb1)
3119  itind = ftind*ntbmx + 0.5
3120  zexm1 = exp_tbl(itind)
3121  endif
3122  zexp1 = f_one / zexm1
3123 
3124  zb2 = min( sntz*ztau1, 500.0 )
3125  if ( zb2 <= od_lo ) then
3126  zexm2 = f_one - zb2 + 0.5*zb2*zb2
3127  else
3128  ftind = zb2 / (bpade + zb2)
3129  itind = ftind*ntbmx + 0.5
3130  zexm2 = exp_tbl(itind)
3131  endif
3132  zexp2 = f_one / zexm2
3133  ze1r45 = zr4*zexp1 + zr5*zexm1
3134 
3135 ! ... collimated beam
3136  if (ze1r45>=-eps1 .and. ze1r45<=eps1) then
3137  zrefb(kp) = eps1
3138  ztrab(kp) = zexm2
3139  else
3140  zden1 = zssa1 / ze1r45
3141  zrefb(kp) = max(f_zero, min(f_one, &
3142  & (zr1*zexp1 - zr2*zexm1 - zr3*zexm2)*zden1 ))
3143  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one &
3144  & - (zt1*zexp1 - zt2*zexm1 - zt3*zexp2)*zden1) ))
3145  endif
3146 
3147 ! ... diffuse beam
3148  zden1 = zr4 / (ze1r45 * zrkg1)
3149  zrefd(kp) = max(f_zero, min(f_one, &
3150  & zgam2*(zexp1 - zexm1)*zden1 ))
3151  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
3152  endif ! end if_zssaw_block
3153 
3154 ! --- ... direct beam transmittance. use exponential lookup table
3155 ! for transmittance, or expansion of exponential for low
3156 ! optical depth
3157 
3158  zr1 = ztau1 * sntz
3159  if ( zr1 <= od_lo ) then
3160  zexp3 = f_one - zr1 + 0.5*zr1*zr1
3161  else
3162  ftind = zr1 / (bpade + zr1)
3163  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3164  zexp3 = exp_tbl(itind)
3165  endif
3166 
3167  ztdbt(k) = zexp3 * ztdbt(kp)
3168  zldbt(kp) = zexp3
3169 
3170 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3171 ! (must use 'orig', unscaled cloud optical depth)
3172 
3173  zr1 = ztau0 * sntz
3174  if ( zr1 <= od_lo ) then
3175  zexp4 = f_one - zr1 + 0.5*zr1*zr1
3176  else
3177  ftind = zr1 / (bpade + zr1)
3178  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3179  zexp4 = exp_tbl(itind)
3180  endif
3181 
3182  zldbt0(k) = zexp4
3183  ztdbt0 = zexp4 * ztdbt0
3184  enddo ! end do_k_loop
3185 
3186  call swflux &
3187 ! --- inputs:
3188  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3189  & nlay, nlp1, &
3190 ! --- outputs:
3191  & zfu, zfd &
3192  & )
3193 
3194 ! --- ... compute upward and downward fluxes at levels
3195  do k = 1, nlp1
3196  fxup0(k,ib) = fxup0(k,ib) + zsolar*zfu(k)
3197  fxdn0(k,ib) = fxdn0(k,ib) + zsolar*zfd(k)
3198  enddo
3199 
3200 !! --- ... surface downward beam/diffuse flux components
3201  zb1 = zsolar*ztdbt0
3202  zb2 = zsolar*(zfd(1) - ztdbt0)
3203 
3204  if (ibd /= 0) then
3205  sfbm0(ibd) = sfbm0(ibd) + zb1
3206  sfdf0(ibd) = sfdf0(ibd) + zb2
3207  else
3208  zf1 = 0.5 * zb1
3209  zf2 = 0.5 * zb2
3210  sfbm0(1) = sfbm0(1) + zf1
3211  sfdf0(1) = sfdf0(1) + zf2
3212  sfbm0(2) = sfbm0(2) + zf1
3213  sfdf0(2) = sfdf0(2) + zf2
3214  endif
3215 ! sfbm0(ibd) = sfbm0(ibd) + zsolar*ztdbt0
3216 ! sfdf0(ibd) = sfdf0(ibd) + zsolar*(zfd(1) - ztdbt0)
3217 
3218 ! --- ... compute total sky optical parameters, layer reflectance and transmittance
3219 
3220  if ( cf1 > eps ) then
3221 
3222 ! --- ... set up toa direct beam and surface values (beam and diff)
3223  ztdbt0 = f_one
3224  zldbt(1) = f_zero
3225 
3226  do k = nlay, 1, -1
3227  kp = k + 1
3228  if ( cldfmc(k,jg) > ftiny ) then ! it is a cloudy-layer
3229 
3230  ztau0 = ztaus(k) + taucw(k,ib)
3231  zssa0 = zssas(k) + ssacw(k,ib)
3232  zasy0 = zasys(k) + asycw(k,ib)
3233  zssaw = min(oneminus, zssa0 / ztau0)
3234  zasyw = zasy0 / max(ftiny, zssa0)
3235 
3236 ! --- ... delta scaling for total-sky condition
3237  za1 = zasyw * zasyw
3238  za2 = zssaw * za1
3239 
3240  ztau1 = (f_one - za2) * ztau0
3241  zssa1 = (zssaw - za2) / (f_one - za2)
3242 !org zasy1 = (zasyw - za1) / (f_one - za1)
3243  zasy1 = zasyw / (f_one + zasyw)
3244  zasy3 = 0.75 * zasy1
3245 
3246 ! --- ... general two-stream expressions
3247  if ( iswmode == 1 ) then
3248  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
3249  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
3250  zgam3 = 0.5 - zasy3 * cosz
3251  elseif ( iswmode == 2 ) then ! pifm
3252  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
3253  zgam2 = 0.75* zssa1 * (f_one- zasy1)
3254  zgam3 = 0.5 - zasy3 * cosz
3255  elseif ( iswmode == 3 ) then ! discrete ordinates
3256  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
3257  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
3258  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
3259  endif
3260  zgam4 = f_one - zgam3
3261 
3262 ! --- ... compute homogeneous reflectance and transmittance
3263 
3264  if ( zssaw >= zcrit ) then ! for conservative scattering
3265  za1 = zgam1 * cosz - zgam3
3266  za2 = zgam1 * ztau1
3267 
3268 ! --- ... use exponential lookup table for transmittance, or expansion
3269 ! of exponential for low optical depth
3270 
3271  zb1 = min( ztau1*sntz , 500.0 )
3272  if ( zb1 <= od_lo ) then
3273  zb2 = f_one - zb1 + 0.5*zb1*zb1
3274  else
3275  ftind = zb1 / (bpade + zb1)
3276  itind = ftind*ntbmx + 0.5
3277  zb2 = exp_tbl(itind)
3278  endif
3279 
3280 ! ... collimated beam
3281  zrefb(kp) = max(f_zero, min(f_one, &
3282  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
3283  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp)))
3284 
3285 ! ... isotropic incidence
3286  zrefd(kp) = max(f_zero, min(f_one, za2 / (f_one+za2) ))
3287  ztrad(kp) = max(f_zero, min(f_one, f_one - zrefd(kp) ))
3288 
3289  else ! for non-conservative scattering
3290  za1 = zgam1*zgam4 + zgam2*zgam3
3291  za2 = zgam1*zgam3 + zgam2*zgam4
3292  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
3293  zrk2= 2.0 * zrk
3294 
3295  zrp = zrk * cosz
3296  zrp1 = f_one + zrp
3297  zrm1 = f_one - zrp
3298  zrpp = f_one - zrp*zrp
3299  zrkg1= zrk + zgam1
3300  zrkg3= zrk * zgam3
3301  zrkg4= zrk * zgam4
3302 
3303  zr1 = zrm1 * (za2 + zrkg3)
3304  zr2 = zrp1 * (za2 - zrkg3)
3305  zr3 = zrk2 * (zgam3 - za2*cosz)
3306  zr4 = zrpp * zrkg1
3307  zr5 = zrpp * (zrk - zgam1)
3308 
3309  zt1 = zrp1 * (za1 + zrkg4)
3310  zt2 = zrm1 * (za1 - zrkg4)
3311  zt3 = zrk2 * (zgam4 + za1*cosz)
3312 
3313 ! --- ... use exponential lookup table for transmittance, or expansion
3314 ! of exponential for low optical depth
3315 
3316  zb1 = min( zrk*ztau1, 500.0 )
3317  if ( zb1 <= od_lo ) then
3318  zexm1 = f_one - zb1 + 0.5*zb1*zb1
3319  else
3320  ftind = zb1 / (bpade + zb1)
3321  itind = ftind*ntbmx + 0.5
3322  zexm1 = exp_tbl(itind)
3323  endif
3324  zexp1 = f_one / zexm1
3325 
3326  zb2 = min( ztau1*sntz, 500.0 )
3327  if ( zb2 <= od_lo ) then
3328  zexm2 = f_one - zb2 + 0.5*zb2*zb2
3329  else
3330  ftind = zb2 / (bpade + zb2)
3331  itind = ftind*ntbmx + 0.5
3332  zexm2 = exp_tbl(itind)
3333  endif
3334  zexp2 = f_one / zexm2
3335  ze1r45 = zr4*zexp1 + zr5*zexm1
3336 
3337 ! ... collimated beam
3338  if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then
3339  zrefb(kp) = eps1
3340  ztrab(kp) = zexm2
3341  else
3342  zden1 = zssa1 / ze1r45
3343  zrefb(kp) = max(f_zero, min(f_one, &
3344  & (zr1*zexp1-zr2*zexm1-zr3*zexm2)*zden1 ))
3345  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one - &
3346  & (zt1*zexp1-zt2*zexm1-zt3*zexp2)*zden1) ))
3347  endif
3348 
3349 ! ... diffuse beam
3350  zden1 = zr4 / (ze1r45 * zrkg1)
3351  zrefd(kp) = max(f_zero, min(f_one, &
3352  & zgam2*(zexp1 - zexm1)*zden1 ))
3353  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
3354  endif ! end if_zssaw_block
3355 
3356 ! --- ... direct beam transmittance. use exponential lookup table
3357 ! for transmittance, or expansion of exponential for low
3358 ! optical depth
3359 
3360  zr1 = ztau1 * sntz
3361  if ( zr1 <= od_lo ) then
3362  zexp3 = f_one - zr1 + 0.5*zr1*zr1
3363  else
3364  ftind = zr1 / (bpade + zr1)
3365  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3366  zexp3 = exp_tbl(itind)
3367  endif
3368 
3369  zldbt(kp) = zexp3
3370  ztdbt(k) = zexp3 * ztdbt(kp)
3371 
3372 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3373 ! (must use 'orig', unscaled cloud optical depth)
3374 
3375  zr1 = ztau0 * sntz
3376  if ( zr1 <= od_lo ) then
3377  zexp4 = f_one - zr1 + 0.5*zr1*zr1
3378  else
3379  ftind = zr1 / (bpade + zr1)
3380  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3381  zexp4 = exp_tbl(itind)
3382  endif
3383 
3384  ztdbt0 = zexp4 * ztdbt0
3385 
3386  else ! if_cldfmc_block --- it is a clear layer
3387 
3388 ! --- ... direct beam transmittance
3389  ztdbt(k) = zldbt(kp) * ztdbt(kp)
3390 
3391 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3392  ztdbt0 = zldbt0(k) * ztdbt0
3393 
3394  endif ! end if_cldfmc_block
3395  enddo ! end do_k_loop
3396 
3397 ! --- ... perform vertical quadrature
3398 
3399  call swflux &
3400 ! --- inputs:
3401  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3402  & nlay, nlp1, &
3403 ! --- outputs:
3404  & zfu, zfd &
3405  & )
3406 
3407 ! --- ... compute upward and downward fluxes at levels
3408  do k = 1, nlp1
3409  fxupc(k,ib) = fxupc(k,ib) + zsolar*zfu(k)
3410  fxdnc(k,ib) = fxdnc(k,ib) + zsolar*zfd(k)
3411  enddo
3412 
3413 !! --- ... surface downward beam/diffused flux components
3414  zb1 = zsolar*ztdbt0
3415  zb2 = zsolar*(zfd(1) - ztdbt0)
3416 
3417  if (ibd /= 0) then
3418  sfbmc(ibd) = sfbmc(ibd) + zb1
3419  sfdfc(ibd) = sfdfc(ibd) + zb2
3420  else
3421  zf1 = 0.5 * zb1
3422  zf2 = 0.5 * zb2
3423  sfbmc(1) = sfbmc(1) + zf1
3424  sfdfc(1) = sfdfc(1) + zf2
3425  sfbmc(2) = sfbmc(2) + zf1
3426  sfdfc(2) = sfdfc(2) + zf2
3427  endif
3428 ! sfbmc(ibd) = sfbmc(ibd) + zsolar*ztdbt0
3429 ! sfdfc(ibd) = sfdfc(ibd) + zsolar*(zfd(1) - ztdbt0)
3430 
3431  endif ! end if_cf1_block
3432 
3433  enddo lab_do_jg
3434 
3435 ! --- ... end of g-point loop
3436 
3437  do ib = 1, nbdsw
3438  ftoadc = ftoadc + fxdn0(nlp1,ib)
3439  ftoau0 = ftoau0 + fxup0(nlp1,ib)
3440  fsfcu0 = fsfcu0 + fxup0(1,ib)
3441  fsfcd0 = fsfcd0 + fxdn0(1,ib)
3442  enddo
3443 
3444 !! --- ... uv-b surface downward flux
3445  ibd = nuvb - nblow + 1
3446  suvbf0 = fxdn0(1,ibd)
3447 
3448  if ( cf1 <= eps ) then ! clear column, set total-sky=clear-sky fluxes
3449  do ib = 1, nbdsw
3450  do k = 1, nlp1
3451  fxupc(k,ib) = fxup0(k,ib)
3452  fxdnc(k,ib) = fxdn0(k,ib)
3453  enddo
3454  enddo
3455 
3456  ftoauc = ftoau0
3457  fsfcuc = fsfcu0
3458  fsfcdc = fsfcd0
3459 
3460 !! --- ... surface downward beam/diffused flux components
3461  sfbmc(1) = sfbm0(1)
3462  sfdfc(1) = sfdf0(1)
3463  sfbmc(2) = sfbm0(2)
3464  sfdfc(2) = sfdf0(2)
3465 
3466 !! --- ... uv-b surface downward flux
3467  suvbfc = suvbf0
3468  else ! cloudy column, compute total-sky fluxes
3469  do ib = 1, nbdsw
3470  ftoauc = ftoauc + fxupc(nlp1,ib)
3471  fsfcuc = fsfcuc + fxupc(1,ib)
3472  fsfcdc = fsfcdc + fxdnc(1,ib)
3473  enddo
3474 
3475 !! --- ... uv-b surface downward flux
3476  suvbfc = fxdnc(1,ibd)
3477  endif ! end if_cf1_block
3478 
3479  return
3480 !...................................
3481  end subroutine spcvrtm
3482 !-----------------------------------
3483 
3485 !-----------------------------------
3486  subroutine swflux &
3487 !...................................
3488 ! --- inputs:
3489  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3490  & nlay, nlp1, &
3491 ! --- outputs:
3492  & zfu, zfd &
3493  & )
3494 
3495 ! =================== program usage description =================== !
3496 ! !
3497 ! purpose: computes the upward and downward radiation fluxes !
3498 ! !
3499 ! interface: "swflux" is called by "spcvrc" and "spcvrm" !
3500 ! !
3501 ! subroutines called : none !
3502 ! !
3503 ! ==================== defination of variables ==================== !
3504 ! !
3505 ! input variables: !
3506 ! zrefb(NLP1) - layer direct beam reflectivity !
3507 ! zrefd(NLP1) - layer diffuse reflectivity !
3508 ! ztrab(NLP1) - layer direct beam transmissivity !
3509 ! ztrad(NLP1) - layer diffuse transmissivity !
3510 ! zldbt(NLP1) - layer mean beam transmittance !
3511 ! ztdbt(NLP1) - total beam transmittance at levels !
3512 ! NLAY, NLP1 - number of layers/levels !
3513 ! !
3514 ! output variables: !
3515 ! zfu (NLP1) - upward flux at layer interface !
3516 ! zfd (NLP1) - downward flux at layer interface !
3517 ! !
3518 ! ******************************************************************* !
3519 ! ====================== end of description block ================= !
3520 
3521 ! --- inputs:
3522  integer, intent(in) :: nlay, nlp1
3523 
3524  real (kind=kind_phys), dimension(nlp1), intent(in) :: zrefb, &
3525  & zrefd, ztrab, ztrad, ztdbt, zldbt
3526 
3527 ! --- outputs:
3528  real (kind=kind_phys), dimension(nlp1), intent(out) :: zfu, zfd
3529 
3530 ! --- locals:
3531  real (kind=kind_phys), dimension(nlp1) :: zrupb,zrupd,zrdnd,ztdn
3532 
3533  real (kind=kind_phys) :: zden1
3534 
3535  integer :: k, kp
3536 !
3537 !===> ... begin here
3538 !
3539 
3540 ! --- ... link lowest layer with surface
3541 
3542  zrupb(1) = zrefb(1) ! direct beam
3543  zrupd(1) = zrefd(1) ! diffused
3544 
3545 ! --- ... pass from bottom to top
3546 
3547  do k = 1, nlay
3548  kp = k + 1
3549 
3550  zden1 = f_one / ( f_one - zrupd(k)*zrefd(kp) )
3551  zrupb(kp) = zrefb(kp) + ( ztrad(kp) * &
3552  & ( (ztrab(kp) - zldbt(kp))*zrupd(k) + &
3553  & zldbt(kp)*zrupb(k)) ) * zden1
3554  zrupd(kp) = zrefd(kp) + ztrad(kp)*ztrad(kp)*zrupd(k)*zden1
3555  enddo
3556 
3557 ! --- ... upper boundary conditions
3558 
3559  ztdn(nlp1) = f_one
3560  zrdnd(nlp1) = f_zero
3561  ztdn(nlay) = ztrab(nlp1)
3562  zrdnd(nlay) = zrefd(nlp1)
3563 
3564 ! --- ... pass from top to bottom
3565 
3566  do k = nlay, 2, -1
3567  zden1 = f_one / (f_one - zrefd(k)*zrdnd(k))
3568  ztdn(k-1) = ztdbt(k)*ztrab(k) + ( ztrad(k) * &
3569  & ( (ztdn(k) - ztdbt(k)) + ztdbt(k) * &
3570  & zrefb(k)*zrdnd(k) )) * zden1
3571  zrdnd(k-1) = zrefd(k) + ztrad(k)*ztrad(k)*zrdnd(k)*zden1
3572  enddo
3573 
3574 ! --- ... up and down-welling fluxes at levels
3575 
3576  do k = 1, nlp1
3577  zden1 = f_one / (f_one - zrdnd(k)*zrupd(k))
3578  zfu(k) = ( ztdbt(k)*zrupb(k) + &
3579  & (ztdn(k) - ztdbt(k))*zrupd(k) ) * zden1
3580  zfd(k) = ztdbt(k) + ( ztdn(k) - ztdbt(k) + &
3581  & ztdbt(k)*zrupb(k)*zrdnd(k) ) * zden1
3582  enddo
3583 
3584  return
3585 !...................................
3586  end subroutine swflux
3587 !-----------------------------------
3588 
3589 
3590 !-----------------------------------
3591  subroutine taumol &
3592 !...................................
3593 ! --- inputs:
3594  & ( colamt,colmol,fac00,fac01,fac10,fac11,jp,jt,jt1,laytrop, &
3595  & forfac,forfrac,indfor,selffac,selffrac,indself, nlay, &
3596 ! --- outputs:
3597  & sfluxzen, taug, taur &
3598  & )
3599 
3600 ! ================== program usage description ================== !
3601 ! !
3602 ! description: !
3603 ! calculate optical depths for gaseous absorption and rayleigh !
3604 ! scattering. !
3605 ! !
3606 ! subroutines called: taugb## (## = 16 - 29) !
3607 ! !
3608 ! ==================== defination of variables ==================== !
3609 ! !
3610 ! inputs: size !
3611 ! colamt - real, column amounts of absorbing gases the index !
3612 ! are for h2o, co2, o3, n2o, ch4, and o2, !
3613 ! respectively (molecules/cm**2) nlay*maxgas!
3614 ! colmol - real, total column amount (dry air+water vapor) nlay !
3615 ! facij - real, for each layer, these are factors that are !
3616 ! needed to compute the interpolation factors !
3617 ! that multiply the appropriate reference k- !
3618 ! values. a value of 0/1 for i,j indicates !
3619 ! that the corresponding factor multiplies !
3620 ! reference k-value for the lower/higher of the !
3621 ! two appropriate temperatures, and altitudes, !
3622 ! respectively. naly !
3623 ! jp - real, the index of the lower (in altitude) of the !
3624 ! two appropriate ref pressure levels needed !
3625 ! for interpolation. nlay !
3626 ! jt, jt1 - integer, the indices of the lower of the two approp !
3627 ! ref temperatures needed for interpolation (for !
3628 ! pressure levels jp and jp+1, respectively) nlay !
3629 ! laytrop - integer, tropopause layer index 1 !
3630 ! forfac - real, scale factor needed to foreign-continuum. nlay !
3631 ! forfrac - real, factor needed for temperature interpolation nlay !
3632 ! indfor - integer, index of the lower of the two appropriate !
3633 ! reference temperatures needed for foreign- !
3634 ! continuum interpolation nlay !
3635 ! selffac - real, scale factor needed to h2o self-continuum. nlay !
3636 ! selffrac- real, factor needed for temperature interpolation !
3637 ! of reference h2o self-continuum data nlay !
3638 ! indself - integer, index of the lower of the two appropriate !
3639 ! reference temperatures needed for the self- !
3640 ! continuum interpolation nlay !
3641 ! nlay - integer, number of vertical layers 1 !
3642 ! !
3643 ! output: !
3644 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
3645 ! taug - real, spectral optical depth for gases nlay*ngptsw!
3646 ! taur - real, opt depth for rayleigh scattering nlay*ngptsw!
3647 ! !
3648 ! =================================================================== !
3649 ! ************ original subprogram description *************** !
3650 ! !
3651 ! optical depths developed for the !
3652 ! !
3653 ! rapid radiative transfer model (rrtm) !
3654 ! !
3655 ! atmospheric and environmental research, inc. !
3656 ! 131 hartwell avenue !
3657 ! lexington, ma 02421 !
3658 ! !
3659 ! !
3660 ! eli j. mlawer !
3661 ! jennifer delamere !
3662 ! steven j. taubman !
3663 ! shepard a. clough !
3664 ! !
3665 ! !
3666 ! !
3667 ! email: mlawer@aer.com !
3668 ! email: jdelamer@aer.com !
3669 ! !
3670 ! the authors wish to acknowledge the contributions of the !
3671 ! following people: patrick d. brown, michael j. iacono, !
3672 ! ronald e. farren, luke chen, robert bergstrom. !
3673 ! !
3674 ! ******************************************************************* !
3675 ! !
3676 ! taumol !
3677 ! !
3678 ! this file contains the subroutines taugbn (where n goes from !
3679 ! 16 to 29). taugbn calculates the optical depths and Planck !
3680 ! fractions per g-value and layer for band n. !
3681 ! !
3682 ! output: optical depths (unitless) !
3683 ! fractions needed to compute planck functions at every layer !
3684 ! and g-value !
3685 ! !
3686 ! modifications: !
3687 ! !
3688 ! revised: adapted to f90 coding, j.-j.morcrette, ecmwf, feb 2003 !
3689 ! revised: modified for g-point reduction, mjiacono, aer, dec 2003 !
3690 ! revised: reformatted for consistency with rrtmg_lw, mjiacono, aer, !
3691 ! jul 2006 !
3692 ! !
3693 ! ******************************************************************* !
3694 ! ====================== end of description block ================= !
3695 
3696 ! --- inputs:
3697  integer, intent(in) :: nlay, laytrop
3698 
3699  integer, dimension(nlay), intent(in) :: indfor, indself, &
3700  & jp, jt, jt1
3701 
3702  real (kind=kind_phys), dimension(nlay), intent(in) :: colmol, &
3703  & fac00, fac01, fac10, fac11, forfac, forfrac, selffac, &
3704  & selffrac
3705 
3706  real (kind=kind_phys), dimension(nlay,maxgas),intent(in) :: colamt
3707 
3708 ! --- outputs:
3709  real (kind=kind_phys), dimension(ngptsw), intent(out) :: sfluxzen
3710 
3711  real (kind=kind_phys), dimension(nlay,ngptsw), intent(out) :: &
3712  & taug, taur
3713 
3714 ! --- locals:
3715  real (kind=kind_phys) :: fs, speccomb, specmult, colm1, colm2
3716 
3717  integer, dimension(nlay,nblow:nbhgh) :: id0, id1
3718 
3719  integer :: ibd, j, jb, js, k, klow, khgh, klim, ks, njb, ns
3720 !
3721 !===> ... begin here
3722 !
3723 ! --- ... loop over each spectral band
3724 
3725  do jb = nblow, nbhgh
3726 
3727 ! --- ... indices for layer optical depth
3728 
3729  do k = 1, laytrop
3730  id0(k,jb) = ((jp(k)-1)*5 + (jt(k)-1)) * nspa(jb)
3731  id1(k,jb) = ( jp(k) *5 + (jt1(k)-1)) * nspa(jb)
3732  enddo
3733 
3734  do k = laytrop+1, nlay
3735  id0(k,jb) = ((jp(k)-13)*5 + (jt(k)-1)) * nspb(jb)
3736  id1(k,jb) = ((jp(k)-12)*5 + (jt1(k)-1)) * nspb(jb)
3737  enddo
3738 
3739 ! --- ... calculate spectral flux at toa
3740 
3741  ibd = ibx(jb)
3742  njb = ng(jb)
3743  ns = ngs(jb)
3744 
3745  select case (jb)
3746 
3747  case (16, 20, 23, 25, 26, 29)
3748 
3749  do j = 1, njb
3750  sfluxzen(ns+j) = sfluxref01(j,1,ibd)
3751  enddo
3752 
3753  case (27)
3754 
3755  do j = 1, njb
3756  sfluxzen(ns+j) = scalekur * sfluxref01(j,1,ibd)
3757  enddo
3758 
3759  case default
3760 
3761  if (jb==17 .or. jb==28) then
3762 
3763  ks = nlay
3764  lab_do_k1 : do k = laytrop, nlay-1
3765  if (jp(k)<layreffr(jb) .and. jp(k+1)>=layreffr(jb)) then
3766  ks = k + 1
3767  exit lab_do_k1
3768  endif
3769  enddo lab_do_k1
3770 
3771  colm1 = colamt(ks,ix1(jb))
3772  colm2 = colamt(ks,ix2(jb))
3773  speccomb = colm1 + strrat(jb)*colm2
3774  specmult = specwt(jb) * min( oneminus, colm1/speccomb )
3775  js = 1 + int( specmult )
3776  fs = mod(specmult, f_one)
3777 
3778  do j = 1, njb
3779  sfluxzen(ns+j) = sfluxref02(j,js,ibd) &
3780  & + fs * (sfluxref02(j,js+1,ibd) - sfluxref02(j,js,ibd))
3781  enddo
3782 
3783  else
3784 
3785  ks = laytrop
3786  lab_do_k2 : do k = 1, laytrop-1
3787  if (jp(k)<layreffr(jb) .and. jp(k+1)>=layreffr(jb)) then
3788  ks = k + 1
3789  exit lab_do_k2
3790  endif
3791  enddo lab_do_k2
3792 
3793  colm1 = colamt(ks,ix1(jb))
3794  colm2 = colamt(ks,ix2(jb))
3795  speccomb = colm1 + strrat(jb)*colm2
3796  specmult = specwt(jb) * min( oneminus, colm1/speccomb )
3797  js = 1 + int( specmult )
3798  fs = mod(specmult, f_one)
3799 
3800  do j = 1, njb
3801  sfluxzen(ns+j) = sfluxref03(j,js,ibd) &
3802  & + fs * (sfluxref03(j,js+1,ibd) - sfluxref03(j,js,ibd))
3803  enddo
3804 
3805  endif
3806 
3807  end select
3808 
3809  enddo
3810 
3811 ! --- ... call taumol## to calculate layer optical depth
3812 
3813  call taumol16
3814  call taumol17
3815  call taumol18
3816  call taumol19
3817  call taumol20
3818  call taumol21
3819  call taumol22
3820  call taumol23
3821  call taumol24
3822  call taumol25
3823  call taumol26
3824  call taumol27
3825  call taumol28
3826  call taumol29
3827 
3828 
3829 ! =================
3830  contains
3831 ! =================
3832 
3834 !-----------------------------------
3835  subroutine taumol16
3836 !...................................
3837 
3838 ! ------------------------------------------------------------------ !
3839 ! band 16: 2600-3250 cm-1 (low - h2o,ch4; high - ch4) !
3840 ! ------------------------------------------------------------------ !
3841 !
3842  use module_radsw_kgb16
3843 
3844 ! --- locals:
3845 
3846  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
3847  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
3848 
3849  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
3850  integer :: inds, indf, indsp, indfp, j, js, k
3851 
3852 !
3853 !===> ... begin here
3854 !
3855 
3856 ! --- ... compute the optical depth by interpolating in ln(pressure),
3857 ! temperature, and appropriate species. below laytrop, the water
3858 ! vapor self-continuum is interpolated (in temperature) separately.
3859 
3860  do k = 1, nlay
3861  tauray = colmol(k) * rayl
3862 
3863  do j = 1, ng16
3864  taur(k,ns16+j) = tauray
3865  enddo
3866  enddo
3867 
3868  do k = 1, laytrop
3869  speccomb = colamt(k,1) + strrat(16)*colamt(k,5)
3870  specmult = 8.0 * min( oneminus, colamt(k,1)/speccomb )
3871 
3872  js = 1 + int( specmult )
3873  fs = mod( specmult, f_one )
3874  fs1= f_one - fs
3875  fac000 = fs1 * fac00(k)
3876  fac010 = fs1 * fac10(k)
3877  fac100 = fs * fac00(k)
3878  fac110 = fs * fac10(k)
3879  fac001 = fs1 * fac01(k)
3880  fac011 = fs1 * fac11(k)
3881  fac101 = fs * fac01(k)
3882  fac111 = fs * fac11(k)
3883 
3884  ind01 = id0(k,16) + js
3885  ind02 = ind01 + 1
3886  ind03 = ind01 + 9
3887  ind04 = ind01 + 10
3888  ind11 = id1(k,16) + js
3889  ind12 = ind11 + 1
3890  ind13 = ind11 + 9
3891  ind14 = ind11 + 10
3892  inds = indself(k)
3893  indf = indfor(k)
3894  indsp= inds + 1
3895  indfp= indf + 1
3896 
3897  do j = 1, ng16
3898  taug(k,ns16+j) = speccomb &
3899  & *( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
3900  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
3901  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
3902  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
3903  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
3904  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
3905  & + forfac(k) * (forref(indf,j) + forfrac(k) &
3906  & * (forref(indfp,j) - forref(indf,j))))
3907  enddo
3908  enddo
3909 
3910  do k = laytrop+1, nlay
3911  ind01 = id0(k,16) + 1
3912  ind02 = ind01 + 1
3913  ind11 = id1(k,16) + 1
3914  ind12 = ind11 + 1
3915 
3916  do j = 1, ng16
3917  taug(k,ns16+j) = colamt(k,5) &
3918  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
3919  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
3920  enddo
3921  enddo
3922 
3923  return
3924 !...................................
3925  end subroutine taumol16
3926 !-----------------------------------
3927 
3929 !-----------------------------------
3930  subroutine taumol17
3931 !...................................
3932 
3933 ! ------------------------------------------------------------------ !
3934 ! band 17: 3250-4000 cm-1 (low - h2o,co2; high - h2o,co2) !
3935 ! ------------------------------------------------------------------ !
3936 !
3937  use module_radsw_kgb17
3938 
3939 ! --- locals:
3940  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
3941  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
3942 
3943  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
3944  integer :: inds, indf, indsp, indfp, j, js, k
3945 
3946 !
3947 !===> ... begin here
3948 !
3949 
3950 ! --- ... compute the optical depth by interpolating in ln(pressure),
3951 ! temperature, and appropriate species. below laytrop, the water
3952 ! vapor self-continuum is interpolated (in temperature) separately.
3953 
3954  do k = 1, nlay
3955  tauray = colmol(k) * rayl
3956 
3957  do j = 1, ng17
3958  taur(k,ns17+j) = tauray
3959  enddo
3960  enddo
3961 
3962  do k = 1, laytrop
3963  speccomb = colamt(k,1) + strrat(17)*colamt(k,2)
3964  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
3965 
3966  js = 1 + int(specmult)
3967  fs = mod(specmult, f_one)
3968  fs1= f_one - fs
3969  fac000 = fs1 * fac00(k)
3970  fac010 = fs1 * fac10(k)
3971  fac100 = fs * fac00(k)
3972  fac110 = fs * fac10(k)
3973  fac001 = fs1 * fac01(k)
3974  fac011 = fs1 * fac11(k)
3975  fac101 = fs * fac01(k)
3976  fac111 = fs * fac11(k)
3977 
3978  ind01 = id0(k,17) + js
3979  ind02 = ind01 + 1
3980  ind03 = ind01 + 9
3981  ind04 = ind01 + 10
3982  ind11 = id1(k,17) + js
3983  ind12 = ind11 + 1
3984  ind13 = ind11 + 9
3985  ind14 = ind11 + 10
3986 
3987  inds = indself(k)
3988  indf = indfor(k)
3989  indsp= inds + 1
3990  indfp= indf + 1
3991 
3992  do j = 1, ng17
3993  taug(k,ns17+j) = speccomb &
3994  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
3995  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
3996  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
3997  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
3998  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
3999  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4000  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4001  & * (forref(indfp,j) - forref(indf,j))))
4002  enddo
4003  enddo
4004 
4005  do k = laytrop+1, nlay
4006  speccomb = colamt(k,1) + strrat(17)*colamt(k,2)
4007  specmult = 4.0 * min(oneminus, colamt(k,1) / speccomb)
4008 
4009  js = 1 + int(specmult)
4010  fs = mod(specmult, f_one)
4011  fs1= f_one - fs
4012  fac000 = fs1 * fac00(k)
4013  fac010 = fs1 * fac10(k)
4014  fac100 = fs * fac00(k)
4015  fac110 = fs * fac10(k)
4016  fac001 = fs1 * fac01(k)
4017  fac011 = fs1 * fac11(k)
4018  fac101 = fs * fac01(k)
4019  fac111 = fs * fac11(k)
4020 
4021  ind01 = id0(k,17) + js
4022  ind02 = ind01 + 1
4023  ind03 = ind01 + 5
4024  ind04 = ind01 + 6
4025  ind11 = id1(k,17) + js
4026  ind12 = ind11 + 1
4027  ind13 = ind11 + 5
4028  ind14 = ind11 + 6
4029 
4030  indf = indfor(k)
4031  indfp= indf + 1
4032 
4033  do j = 1, ng17
4034  taug(k,ns17+j) = speccomb &
4035  & * ( fac000 * absb(ind01,j) + fac100 * absb(ind02,j) &
4036  & + fac010 * absb(ind03,j) + fac110 * absb(ind04,j) &
4037  & + fac001 * absb(ind11,j) + fac101 * absb(ind12,j) &
4038  & + fac011 * absb(ind13,j) + fac111 * absb(ind14,j) ) &
4039  & + colamt(k,1) * forfac(k) * (forref(indf,j) &
4040  & + forfrac(k) * (forref(indfp,j) - forref(indf,j)))
4041  enddo
4042  enddo
4043 
4044  return
4045 !...................................
4046  end subroutine taumol17
4047 !-----------------------------------
4048 
4050 !-----------------------------------
4051  subroutine taumol18
4052 !...................................
4053 
4054 ! ------------------------------------------------------------------ !
4055 ! band 18: 4000-4650 cm-1 (low - h2o,ch4; high - ch4) !
4056 ! ------------------------------------------------------------------ !
4057 !
4058  use module_radsw_kgb18
4059 
4060 ! --- locals:
4061  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4062  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4063 
4064  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4065  integer :: inds, indf, indsp, indfp, j, js, k
4066 
4067 !
4068 !===> ... begin here
4069 !
4070 
4071 ! --- ... compute the optical depth by interpolating in ln(pressure),
4072 ! temperature, and appropriate species. below laytrop, the water
4073 ! vapor self-continuum is interpolated (in temperature) separately.
4074 
4075  do k = 1, nlay
4076  tauray = colmol(k) * rayl
4077 
4078  do j = 1, ng18
4079  taur(k,ns18+j) = tauray
4080  enddo
4081  enddo
4082 
4083  do k = 1, laytrop
4084  speccomb = colamt(k,1) + strrat(18)*colamt(k,5)
4085  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4086 
4087  js = 1 + int(specmult)
4088  fs = mod(specmult, f_one)
4089  fs1= f_one - fs
4090  fac000 = fs1 * fac00(k)
4091  fac010 = fs1 * fac10(k)
4092  fac100 = fs * fac00(k)
4093  fac110 = fs * fac10(k)
4094  fac001 = fs1 * fac01(k)
4095  fac011 = fs1 * fac11(k)
4096  fac101 = fs * fac01(k)
4097  fac111 = fs * fac11(k)
4098 
4099  ind01 = id0(k,18) + js
4100  ind02 = ind01 + 1
4101  ind03 = ind01 + 9
4102  ind04 = ind01 + 10
4103  ind11 = id1(k,18) + js
4104  ind12 = ind11 + 1
4105  ind13 = ind11 + 9
4106  ind14 = ind11 + 10
4107 
4108  inds = indself(k)
4109  indf = indfor(k)
4110  indsp= inds + 1
4111  indfp= indf + 1
4112 
4113  do j = 1, ng18
4114  taug(k,ns18+j) = speccomb &
4115  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4116  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4117  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4118  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4119  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4120  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4121  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4122  & * (forref(indfp,j) - forref(indf,j))))
4123  enddo
4124  enddo
4125 
4126  do k = laytrop+1, nlay
4127  ind01 = id0(k,18) + 1
4128  ind02 = ind01 + 1
4129  ind11 = id1(k,18) + 1
4130  ind12 = ind11 + 1
4131 
4132  do j = 1, ng18
4133  taug(k,ns18+j) = colamt(k,5) &
4134  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4135  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
4136  enddo
4137  enddo
4138 
4139  return
4140 !...................................
4141  end subroutine taumol18
4142 !-----------------------------------
4143 
4145 !-----------------------------------
4146  subroutine taumol19
4147 !...................................
4148 
4149 ! ------------------------------------------------------------------ !
4150 ! band 19: 4650-5150 cm-1 (low - h2o,co2; high - co2) !
4151 ! ------------------------------------------------------------------ !
4152 !
4153  use module_radsw_kgb19
4154 
4155 ! --- locals:
4156  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4157  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4158 
4159  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4160  integer :: inds, indf, indsp, indfp, j, js, k
4161 
4162 !
4163 !===> ... begin here
4164 !
4165 
4166 ! --- ... compute the optical depth by interpolating in ln(pressure),
4167 ! temperature, and appropriate species. below laytrop, the water
4168 ! vapor self-continuum is interpolated (in temperature) separately.
4169 
4170  do k = 1, nlay
4171  tauray = colmol(k) * rayl
4172 
4173  do j = 1, ng19
4174  taur(k,ns19+j) = tauray
4175  enddo
4176  enddo
4177 
4178  do k = 1, laytrop
4179  speccomb = colamt(k,1) + strrat(19)*colamt(k,2)
4180  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4181 
4182  js = 1 + int(specmult)
4183  fs = mod(specmult, f_one)
4184  fs1= f_one - fs
4185  fac000 = fs1 * fac00(k)
4186  fac010 = fs1 * fac10(k)
4187  fac100 = fs * fac00(k)
4188  fac110 = fs * fac10(k)
4189  fac001 = fs1 * fac01(k)
4190  fac011 = fs1 * fac11(k)
4191  fac101 = fs * fac01(k)
4192  fac111 = fs * fac11(k)
4193 
4194  ind01 = id0(k,19) + js
4195  ind02 = ind01 + 1
4196  ind03 = ind01 + 9
4197  ind04 = ind01 + 10
4198  ind11 = id1(k,19) + js
4199  ind12 = ind11 + 1
4200  ind13 = ind11 + 9
4201  ind14 = ind11 + 10
4202 
4203  inds = indself(k)
4204  indf = indfor(k)
4205  indsp= inds + 1
4206  indfp= indf + 1
4207 
4208  do j = 1, ng19
4209  taug(k,ns19+j) = speccomb &
4210  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4211  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4212  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4213  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4214  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4215  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4216  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4217  & * (forref(indfp,j) - forref(indf,j))))
4218  enddo
4219  enddo
4220 
4221  do k = laytrop+1, nlay
4222  ind01 = id0(k,19) + 1
4223  ind02 = ind01 + 1
4224  ind11 = id1(k,19) + 1
4225  ind12 = ind11 + 1
4226 
4227  do j = 1, ng19
4228  taug(k,ns19+j) = colamt(k,2) &
4229  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4230  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
4231  enddo
4232  enddo
4233 
4234 !...................................
4235  end subroutine taumol19
4236 !-----------------------------------
4237 
4239 !-----------------------------------
4240  subroutine taumol20
4241 !...................................
4242 
4243 ! ------------------------------------------------------------------ !
4244 ! band 20: 5150-6150 cm-1 (low - h2o; high - h2o) !
4245 ! ------------------------------------------------------------------ !
4246 !
4247  use module_radsw_kgb20
4248 
4249 ! --- locals:
4250  real (kind=kind_phys) :: tauray
4251 
4252  integer :: ind01, ind02, ind11, ind12
4253  integer :: inds, indf, indsp, indfp, j, k
4254 
4255 !
4256 !===> ... begin here
4257 !
4258 
4259 ! --- ... compute the optical depth by interpolating in ln(pressure),
4260 ! temperature, and appropriate species. below laytrop, the water
4261 ! vapor self-continuum is interpolated (in temperature) separately.
4262 
4263  do k = 1, nlay
4264  tauray = colmol(k) * rayl
4265 
4266  do j = 1, ng20
4267  taur(k,ns20+j) = tauray
4268  enddo
4269  enddo
4270 
4271  do k = 1, laytrop
4272  ind01 = id0(k,20) + 1
4273  ind02 = ind01 + 1
4274  ind11 = id1(k,20) + 1
4275  ind12 = ind11 + 1
4276 
4277  inds = indself(k)
4278  indf = indfor(k)
4279  indsp= inds + 1
4280  indfp= indf + 1
4281 
4282  do j = 1, ng20
4283  taug(k,ns20+j) = colamt(k,1) &
4284  & * ( (fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
4285  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j)) &
4286  & + selffac(k) * (selfref(inds,j) + selffrac(k) &
4287  & * (selfref(indsp,j) - selfref(inds,j))) &
4288  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4289  & * (forref(indfp,j) - forref(indf,j))) ) &
4290  & + colamt(k,5) * absch4(j)
4291  enddo
4292  enddo
4293 
4294  do k = laytrop+1, nlay
4295  ind01 = id0(k,20) + 1
4296  ind02 = ind01 + 1
4297  ind11 = id1(k,20) + 1
4298  ind12 = ind11 + 1
4299 
4300  indf = indfor(k)
4301  indfp= indf + 1
4302 
4303  do j = 1, ng20
4304  taug(k,ns20+j) = colamt(k,1) &
4305  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4306  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) &
4307  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4308  & * (forref(indfp,j) - forref(indf,j))) ) &
4309  & + colamt(k,5) * absch4(j)
4310  enddo
4311  enddo
4312 
4313  return
4314 !...................................
4315  end subroutine taumol20
4316 !-----------------------------------
4317 
4319 !-----------------------------------
4320  subroutine taumol21
4321 !...................................
4322 
4323 ! ------------------------------------------------------------------ !
4324 ! band 21: 6150-7700 cm-1 (low - h2o,co2; high - h2o,co2) !
4325 ! ------------------------------------------------------------------ !
4326 !
4327  use module_radsw_kgb21
4328 
4329 ! --- locals:
4330  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4331  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4332 
4333  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4334  integer :: inds, indf, indsp, indfp, j, js, k
4335 
4336 !
4337 !===> ... begin here
4338 !
4339 
4340 ! --- ... compute the optical depth by interpolating in ln(pressure),
4341 ! temperature, and appropriate species. below laytrop, the water
4342 ! vapor self-continuum is interpolated (in temperature) separately.
4343 
4344  do k = 1, nlay
4345  tauray = colmol(k) * rayl
4346 
4347  do j = 1, ng21
4348  taur(k,ns21+j) = tauray
4349  enddo
4350  enddo
4351 
4352  do k = 1, laytrop
4353  speccomb = colamt(k,1) + strrat(21)*colamt(k,2)
4354  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4355 
4356  js = 1 + int(specmult)
4357  fs = mod(specmult, f_one)
4358  fs1= f_one - fs
4359  fac000 = fs1 * fac00(k)
4360  fac010 = fs1 * fac10(k)
4361  fac100 = fs * fac00(k)
4362  fac110 = fs * fac10(k)
4363  fac001 = fs1 * fac01(k)
4364  fac011 = fs1 * fac11(k)
4365  fac101 = fs * fac01(k)
4366  fac111 = fs * fac11(k)
4367 
4368  ind01 = id0(k,21) + js
4369  ind02 = ind01 + 1
4370  ind03 = ind01 + 9
4371  ind04 = ind01 + 10
4372  ind11 = id1(k,21) + js
4373  ind12 = ind11 + 1
4374  ind13 = ind11 + 9
4375  ind14 = ind11 + 10
4376 
4377  inds = indself(k)
4378  indf = indfor(k)
4379  indsp= inds + 1
4380  indfp= indf + 1
4381 
4382  do j = 1, ng21
4383  taug(k,ns21+j) = speccomb &
4384  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4385  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4386  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4387  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4388  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4389  & + selffrac(k) * (selfref(indsp,j) - selfref(inds,j))) &
4390  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4391  & * (forref(indfp,j) - forref(indf,j))))
4392  enddo
4393  enddo
4394 
4395  do k = laytrop+1, nlay
4396  speccomb = colamt(k,1) + strrat(21)*colamt(k,2)
4397  specmult = 4.0 * min(oneminus, colamt(k,1) / speccomb)
4398 
4399  js = 1 + int(specmult)
4400  fs = mod(specmult, f_one)
4401  fs1= f_one - fs
4402  fac000 = fs1 * fac00(k)
4403  fac010 = fs1 * fac10(k)
4404  fac100 = fs * fac00(k)
4405  fac110 = fs * fac10(k)
4406  fac001 = fs1 * fac01(k)
4407  fac011 = fs1 * fac11(k)
4408  fac101 = fs * fac01(k)
4409  fac111 = fs * fac11(k)
4410 
4411  ind01 = id0(k,21) + js
4412  ind02 = ind01 + 1
4413  ind03 = ind01 + 5
4414  ind04 = ind01 + 6
4415  ind11 = id1(k,21) + js
4416  ind12 = ind11 + 1
4417  ind13 = ind11 + 5
4418  ind14 = ind11 + 6
4419 
4420  indf = indfor(k)
4421  indfp= indf + 1
4422 
4423  do j = 1, ng21
4424  taug(k,ns21+j) = speccomb &
4425  & * ( fac000 * absb(ind01,j) + fac100 * absb(ind02,j) &
4426  & + fac010 * absb(ind03,j) + fac110 * absb(ind04,j) &
4427  & + fac001 * absb(ind11,j) + fac101 * absb(ind12,j) &
4428  & + fac011 * absb(ind13,j) + fac111 * absb(ind14,j) ) &
4429  & + colamt(k,1) * forfac(k) * (forref(indf,j) &
4430  & + forfrac(k) * (forref(indfp,j) - forref(indf,j)))
4431  enddo
4432  enddo
4433 
4434 !...................................
4435  end subroutine taumol21
4436 !-----------------------------------
4437 
4439 !-----------------------------------
4440  subroutine taumol22
4441 !...................................
4442 
4443 ! ------------------------------------------------------------------ !
4444 ! band 22: 7700-8050 cm-1 (low - h2o,o2; high - o2) !
4445 ! ------------------------------------------------------------------ !
4446 !
4447  use module_radsw_kgb22
4448 
4449 ! --- locals:
4450  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4451  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111, &
4452  & o2adj, o2cont, o2tem
4453 
4454  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4455  integer :: inds, indf, indsp, indfp, j, js, k
4456 
4457 !
4458 !===> ... begin here
4459 !
4460 ! --- ... the following factor is the ratio of total o2 band intensity (lines
4461 ! and mate continuum) to o2 band intensity (line only). it is needed
4462 ! to adjust the optical depths since the k's include only lines.
4463 
4464  o2adj = 1.6
4465  o2tem = 4.35e-4 / (350.0*2.0)
4466 
4467 ! --- ... compute the optical depth by interpolating in ln(pressure),
4468 ! temperature, and appropriate species. below laytrop, the water
4469 ! vapor self-continuum is interpolated (in temperature) separately.
4470 
4471  do k = 1, nlay
4472  tauray = colmol(k) * rayl
4473 
4474  do j = 1, ng22
4475  taur(k,ns22+j) = tauray
4476  enddo
4477  enddo
4478 
4479  do k = 1, laytrop
4480  o2cont = o2tem * colamt(k,6)
4481  speccomb = colamt(k,1) + strrat(22)*colamt(k,6)
4482  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4483 
4484  js = 1 + int(specmult)
4485  fs = mod(specmult, f_one)
4486  fs1= f_one - fs
4487  fac000 = fs1 * fac00(k)
4488  fac010 = fs1 * fac10(k)
4489  fac100 = fs * fac00(k)
4490  fac110 = fs * fac10(k)
4491  fac001 = fs1 * fac01(k)
4492  fac011 = fs1 * fac11(k)
4493  fac101 = fs * fac01(k)
4494  fac111 = fs * fac11(k)
4495 
4496  ind01 = id0(k,22) + js
4497  ind02 = ind01 + 1
4498  ind03 = ind01 + 9
4499  ind04 = ind01 + 10
4500  ind11 = id1(k,22) + js
4501  ind12 = ind11 + 1
4502  ind13 = ind11 + 9
4503  ind14 = ind11 + 10
4504 
4505  inds = indself(k)
4506  indf = indfor(k)
4507  indsp= inds + 1
4508  indfp= indf + 1
4509 
4510  do j = 1, ng22
4511  taug(k,ns22+j) = speccomb &
4512  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4513  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4514  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4515  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4516  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4517  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4518  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4519  & * (forref(indfp,j) - forref(indf,j)))) + o2cont
4520  enddo
4521  enddo
4522 
4523  do k = laytrop+1, nlay
4524  o2cont = o2tem * colamt(k,6)
4525 
4526  ind01 = id0(k,22) + 1
4527  ind02 = ind01 + 1
4528  ind11 = id1(k,22) + 1
4529  ind12 = ind11 + 1
4530 
4531  do j = 1, ng22
4532  taug(k,ns22+j) = colamt(k,6) * o2adj &
4533  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4534  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) &
4535  & + o2cont
4536  enddo
4537  enddo
4538 
4539  return
4540 !...................................
4541  end subroutine taumol22
4542 !-----------------------------------
4543 
4545 !-----------------------------------
4546  subroutine taumol23
4547 !...................................
4548 
4549 ! ------------------------------------------------------------------ !
4550 ! band 23: 8050-12850 cm-1 (low - h2o; high - nothing) !
4551 ! ------------------------------------------------------------------ !
4552 !
4553  use module_radsw_kgb23
4554 
4555 ! --- locals:
4556  integer :: ind01, ind02, ind11, ind12
4557  integer :: inds, indf, indsp, indfp, j, k
4558 
4559 !
4560 !===> ... begin here
4561 !
4562 
4563 ! --- ... compute the optical depth by interpolating in ln(pressure),
4564 ! temperature, and appropriate species. below laytrop, the water
4565 ! vapor self-continuum is interpolated (in temperature) separately.
4566 
4567  do k = 1, nlay
4568  do j = 1, ng23
4569  taur(k,ns23+j) = colmol(k) * rayl(j)
4570  enddo
4571  enddo
4572 
4573  do k = 1, laytrop
4574  ind01 = id0(k,23) + 1
4575  ind02 = ind01 + 1
4576  ind11 = id1(k,23) + 1
4577  ind12 = ind11 + 1
4578 
4579  inds = indself(k)
4580  indf = indfor(k)
4581  indsp= inds + 1
4582  indfp= indf + 1
4583 
4584  do j = 1, ng23
4585  taug(k,ns23+j) = colamt(k,1) * (givfac &
4586  & * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
4587  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) &
4588  & + selffac(k) * (selfref(inds,j) + selffrac(k) &
4589  & * (selfref(indsp,j) - selfref(inds,j))) &
4590  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4591  & * (forref(indfp,j) - forref(indf,j))))
4592  enddo
4593  enddo
4594 
4595  do k = laytrop+1, nlay
4596  do j = 1, ng23
4597  taug(k,ns23+j) = f_zero
4598  enddo
4599  enddo
4600 
4601 !...................................
4602  end subroutine taumol23
4603 !-----------------------------------
4604 
4606 !-----------------------------------
4607  subroutine taumol24
4608 !...................................
4609 
4610 ! ------------------------------------------------------------------ !
4611 ! band 24: 12850-16000 cm-1 (low - h2o,o2; high - o2) !
4612 ! ------------------------------------------------------------------ !
4613 !
4614  use module_radsw_kgb24
4615 
4616 ! --- locals:
4617  real (kind=kind_phys) :: speccomb, specmult, fs, fs1, &
4618  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4619 
4620  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4621  integer :: inds, indf, indsp, indfp, j, js, k
4622 
4623 !
4624 !===> ... begin here
4625 !
4626 
4627 ! --- ... compute the optical depth by interpolating in ln(pressure),
4628 ! temperature, and appropriate species. below laytrop, the water
4629 ! vapor self-continuum is interpolated (in temperature) separately.
4630 
4631  do k = 1, laytrop
4632  speccomb = colamt(k,1) + strrat(24)*colamt(k,6)
4633  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4634 
4635  js = 1 + int(specmult)
4636  fs = mod(specmult, f_one)
4637  fs1= f_one - fs
4638  fac000 = fs1 * fac00(k)
4639  fac010 = fs1 * fac10(k)
4640  fac100 = fs * fac00(k)
4641  fac110 = fs * fac10(k)
4642  fac001 = fs1 * fac01(k)
4643  fac011 = fs1 * fac11(k)
4644  fac101 = fs * fac01(k)
4645  fac111 = fs * fac11(k)
4646 
4647  ind01 = id0(k,24) + js
4648  ind02 = ind01 + 1
4649  ind03 = ind01 + 9
4650  ind04 = ind01 + 10
4651  ind11 = id1(k,24) + js
4652  ind12 = ind11 + 1
4653  ind13 = ind11 + 9
4654  ind14 = ind11 + 10
4655 
4656  inds = indself(k)
4657  indf = indfor(k)
4658  indsp= inds + 1
4659  indfp= indf + 1
4660 
4661  do j = 1, ng24
4662  taug(k,ns24+j) = speccomb &
4663  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4664  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4665  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4666  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4667  & + colamt(k,3) * abso3a(j) + colamt(k,1) &
4668  & * (selffac(k) * (selfref(inds,j) + selffrac(k) &
4669  & * (selfref(indsp,j) - selfref(inds,j))) &
4670  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4671  & * (forref(indfp,j) - forref(indf,j))))
4672 
4673  taur(k,ns24+j) = colmol(k) &
4674  & * (rayla(j,js) + fs*(rayla(j,js+1) - rayla(j,js)))
4675  enddo
4676  enddo
4677 
4678  do k = laytrop+1, nlay
4679  ind01 = id0(k,24) + 1
4680  ind02 = ind01 + 1
4681  ind11 = id1(k,24) + 1
4682  ind12 = ind11 + 1
4683 
4684  do j = 1, ng24
4685  taug(k,ns24+j) = colamt(k,6) &
4686  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4687  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) &
4688  & + colamt(k,3) * abso3b(j)
4689 
4690  taur(k,ns24+j) = colmol(k) * raylb(j)
4691  enddo
4692  enddo
4693 
4694  return
4695 !...................................
4696  end subroutine taumol24
4697 !-----------------------------------
4698 
4700 !-----------------------------------
4701  subroutine taumol25
4702 !...................................
4703 
4704 ! ------------------------------------------------------------------ !
4705 ! band 25: 16000-22650 cm-1 (low - h2o; high - nothing) !
4706 ! ------------------------------------------------------------------ !
4707 !
4708  use module_radsw_kgb25
4709 
4710 ! --- locals:
4711  integer :: ind01, ind02, ind11, ind12
4712  integer :: j, k
4713 
4714 !
4715 !===> ... begin here
4716 !
4717 
4718 ! --- ... compute the optical depth by interpolating in ln(pressure),
4719 ! temperature, and appropriate species. below laytrop, the water
4720 ! vapor self-continuum is interpolated (in temperature) separately.
4721 
4722  do k = 1, nlay
4723  do j = 1, ng25
4724  taur(k,ns25+j) = colmol(k) * rayl(j)
4725  enddo
4726  enddo
4727 
4728  do k = 1, laytrop
4729  ind01 = id0(k,25) + 1
4730  ind02 = ind01 + 1
4731  ind11 = id1(k,25) + 1
4732  ind12 = ind11 + 1
4733 
4734  do j = 1, ng25
4735  taug(k,ns25+j) = colamt(k,1) &
4736  & * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
4737  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) &
4738  & + colamt(k,3) * abso3a(j)
4739  enddo
4740  enddo
4741 
4742  do k = laytrop+1, nlay
4743  do j = 1, ng25
4744  taug(k,ns25+j) = colamt(k,3) * abso3b(j)
4745  enddo
4746  enddo
4747 
4748  return
4749 !...................................
4750  end subroutine taumol25
4751 !-----------------------------------
4752 
4754 !-----------------------------------
4755  subroutine taumol26
4756 !...................................
4757 
4758 ! ------------------------------------------------------------------ !
4759 ! band 26: 22650-29000 cm-1 (low - nothing; high - nothing) !
4760 ! ------------------------------------------------------------------ !
4761 !
4762  use module_radsw_kgb26
4763 
4764 ! --- locals:
4765  integer :: j, k
4766 
4767 !
4768 !===> ... begin here
4769 !
4770 
4771 ! --- ... compute the optical depth by interpolating in ln(pressure),
4772 ! temperature, and appropriate species. below laytrop, the water
4773 ! vapor self-continuum is interpolated (in temperature) separately.
4774 
4775  do k = 1, nlay
4776  do j = 1, ng26
4777  taug(k,ns26+j) = f_zero
4778  taur(k,ns26+j) = colmol(k) * rayl(j)
4779  enddo
4780  enddo
4781 
4782  return
4783 !...................................
4784  end subroutine taumol26
4785 !-----------------------------------
4786 
4788 !-----------------------------------
4789  subroutine taumol27
4790 !...................................
4791 
4792 ! ------------------------------------------------------------------ !
4793 ! band 27: 29000-38000 cm-1 (low - o3; high - o3) !
4794 ! ------------------------------------------------------------------ !
4795 !
4796  use module_radsw_kgb27
4797 !
4798 ! --- locals:
4799  integer :: ind01, ind02, ind11, ind12
4800  integer :: j, k
4801 
4802 !
4803 !===> ... begin here
4804 !
4805 
4806 ! --- ... compute the optical depth by interpolating in ln(pressure),
4807 ! temperature, and appropriate species. below laytrop, the water
4808 ! vapor self-continuum is interpolated (in temperature) separately.
4809 
4810  do k = 1, nlay
4811  do j = 1, ng27
4812  taur(k,ns27+j) = colmol(k) * rayl(j)
4813  enddo
4814  enddo
4815 
4816  do k = 1, laytrop
4817  ind01 = id0(k,27) + 1
4818  ind02 = ind01 + 1
4819  ind11 = id1(k,27) + 1
4820  ind12 = ind11 + 1
4821 
4822  do j = 1, ng27
4823  taug(k,ns27+j) = colamt(k,3) &
4824  & * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
4825  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) )
4826  enddo
4827  enddo
4828 
4829  do k = laytrop+1, nlay
4830  ind01 = id0(k,27) + 1
4831  ind02 = ind01 + 1
4832  ind11 = id1(k,27) + 1
4833  ind12 = ind11 + 1
4834 
4835  do j = 1, ng27
4836  taug(k,ns27+j) = colamt(k,3) &
4837  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4838  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
4839  enddo
4840  enddo
4841 
4842  return
4843 !...................................
4844  end subroutine taumol27
4845 !-----------------------------------
4846 
4848 !-----------------------------------
4849  subroutine taumol28
4850 !...................................
4851 
4852 ! ------------------------------------------------------------------ !
4853 ! band 28: 38000-50000 cm-1 (low - o3,o2; high - o3,o2) !
4854 ! ------------------------------------------------------------------ !
4855 !
4856  use module_radsw_kgb28
4857 
4858 ! --- locals:
4859  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4860  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4861 
4862  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4863  integer :: j, js, k
4864 
4865 !
4866 !===> ... begin here
4867 !
4868 
4869 ! --- ... compute the optical depth by interpolating in ln(pressure),
4870 ! temperature, and appropriate species. below laytrop, the water
4871 ! vapor self-continuum is interpolated (in temperature) separately.
4872 
4873  do k = 1, nlay
4874  tauray = colmol(k) * rayl
4875 
4876  do j = 1, ng28
4877  taur(k,ns28+j) = tauray
4878  enddo
4879  enddo
4880 
4881  do k = 1, laytrop
4882  speccomb = colamt(k,3) + strrat(28)*colamt(k,6)
4883  specmult = 8.0 * min(oneminus, colamt(k,3) / speccomb)
4884 
4885  js = 1 + int(specmult)
4886  fs = mod(specmult, f_one)
4887  fs1= f_one - fs
4888  fac000 = fs1 * fac00(k)
4889  fac010 = fs1 * fac10(k)
4890  fac100 = fs * fac00(k)
4891  fac110 = fs * fac10(k)
4892  fac001 = fs1 * fac01(k)
4893  fac011 = fs1 * fac11(k)
4894  fac101 = fs * fac01(k)
4895  fac111 = fs * fac11(k)
4896 
4897  ind01 = id0(k,28) + js
4898  ind02 = ind01 + 1
4899  ind03 = ind01 + 9
4900  ind04 = ind01 + 10
4901  ind11 = id1(k,28) + js
4902  ind12 = ind11 + 1
4903  ind13 = ind11 + 9
4904  ind14 = ind11 + 10
4905 
4906  do j = 1, ng28
4907  taug(k,ns28+j) = speccomb &
4908  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4909  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4910  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4911  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) )
4912  enddo
4913  enddo
4914 
4915  do k = laytrop+1, nlay
4916  speccomb = colamt(k,3) + strrat(28)*colamt(k,6)
4917  specmult = 4.0 * min(oneminus, colamt(k,3) / speccomb)
4918 
4919  js = 1 + int(specmult)
4920  fs = mod(specmult, f_one)
4921  fs1= f_one - fs
4922  fac000 = fs1 * fac00(k)
4923  fac010 = fs1 * fac10(k)
4924  fac100 = fs * fac00(k)
4925  fac110 = fs * fac10(k)
4926  fac001 = fs1 * fac01(k)
4927  fac011 = fs1 * fac11(k)
4928  fac101 = fs * fac01(k)
4929  fac111 = fs * fac11(k)
4930 
4931  ind01 = id0(k,28) + js
4932  ind02 = ind01 + 1
4933  ind03 = ind01 + 5
4934  ind04 = ind01 + 6
4935  ind11 = id1(k,28) + js
4936  ind12 = ind11 + 1
4937  ind13 = ind11 + 5
4938  ind14 = ind11 + 6
4939 
4940  do j = 1, ng28
4941  taug(k,ns28+j) = speccomb &
4942  & * ( fac000 * absb(ind01,j) + fac100 * absb(ind02,j) &
4943  & + fac010 * absb(ind03,j) + fac110 * absb(ind04,j) &
4944  & + fac001 * absb(ind11,j) + fac101 * absb(ind12,j) &
4945  & + fac011 * absb(ind13,j) + fac111 * absb(ind14,j) )
4946  enddo
4947  enddo
4948 
4949  return
4950 !...................................
4951  end subroutine taumol28
4952 !-----------------------------------
4953 
4955 !-----------------------------------
4956  subroutine taumol29
4957 !...................................
4958 
4959 ! ------------------------------------------------------------------ !
4960 ! band 29: 820-2600 cm-1 (low - h2o; high - co2) !
4961 ! ------------------------------------------------------------------ !
4962 !
4963  use module_radsw_kgb29
4964 
4965 ! --- locals:
4966  real (kind=kind_phys) :: tauray
4967 
4968  integer :: ind01, ind02, ind11, ind12
4969  integer :: inds, indf, indsp, indfp, j, k
4970 
4971 !
4972 !===> ... begin here
4973 !
4974 
4975 ! --- ... compute the optical depth by interpolating in ln(pressure),
4976 ! temperature, and appropriate species. below laytrop, the water
4977 ! vapor self-continuum is interpolated (in temperature) separately.
4978 
4979  do k = 1, nlay
4980  tauray = colmol(k) * rayl
4981 
4982  do j = 1, ng29
4983  taur(k,ns29+j) = tauray
4984  enddo
4985  enddo
4986 
4987  do k = 1, laytrop
4988  ind01 = id0(k,29) + 1
4989  ind02 = ind01 + 1
4990  ind11 = id1(k,29) + 1
4991  ind12 = ind11 + 1
4992 
4993  inds = indself(k)
4994  indf = indfor(k)
4995  indsp= inds + 1
4996  indfp= indf + 1
4997 
4998  do j = 1, ng29
4999  taug(k,ns29+j) = colamt(k,1) &
5000  & * ( (fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
5001  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) &
5002  & + selffac(k) * (selfref(inds,j) + selffrac(k) &
5003  & * (selfref(indsp,j) - selfref(inds,j))) &
5004  & + forfac(k) * (forref(indf,j) + forfrac(k) &
5005  & * (forref(indfp,j) - forref(indf,j)))) &
5006  & + colamt(k,2) * absco2(j)
5007  enddo
5008  enddo
5009 
5010  do k = laytrop+1, nlay
5011  ind01 = id0(k,29) + 1
5012  ind02 = ind01 + 1
5013  ind11 = id1(k,29) + 1
5014  ind12 = ind11 + 1
5015 
5016  do j = 1, ng29
5017  taug(k,ns29+j) = colamt(k,2) &
5018  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
5019  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) &
5020  & + colamt(k,1) * absh2o(j)
5021  enddo
5022  enddo
5023 
5024  return
5025 !...................................
5026  end subroutine taumol29
5027 !-----------------------------------
5028 
5029 !...................................
5030  end subroutine taumol
5031 !-----------------------------------
5032 
5033 !
5034 !........................................!
5035  end module module_radsw_main !
5036 !========================================!
5037 
real(kind=kind_phys), parameter con_amo3
molecular wght of o3 (g/mol)
Definition: physcons.f:134
real(kind=kind_phys), parameter, public rayl
real(kind=kind_phys), parameter, public rayl
subroutine mcica_subcol
This subroutine computes the sub-colum cloud profile flag array.
Definition: radsw_main.f:1783
This module sets up absorption coeffients for band 19: 4650-5150 cm-1 (low - h2o,co2; high - co2) ...
real(kind=kind_phys), dimension(ng23), public rayl
real(kind=kind_phys), public a1s
real(kind=kind_phys), dimension(msa16, ng16), public absa
real(kind=kind_phys), dimension(msf29, ng29), public selfref
real(kind=kind_phys), parameter amdw
Definition: radsw_main.f:303
real(kind=kind_phys), parameter eps
Definition: radsw_main.f:292
real(kind=kind_phys), dimension(5), public ebari
Definition: radsw_datatb.f:186
This module sets up absorption coeffients for band 24: 12850-16000 cm-1 (low - h2o, o2; high - o2)
subroutine taumol17
The subroutine computes the optical depth in band 17: 3250-4000 cm-1 (low - h2o,co2; high - h2o...
Definition: radsw_main.f:3931
real(kind=kind_phys), dimension(msf19, ng19), public selfref
character(40), parameter vtagsw
Definition: radsw_main.f:280
real(kind=kind_phys), dimension(mfr29, ng29), public forref
integer, dimension(nblow:nbhgh), public ix2
integer, dimension(nblow:nbhgh), public ix1
real(kind=kind_phys), dimension(mfr22, ng22), public forref
real(kind=kind_phys), dimension(msb21, ng21), public absb
subroutine, public swrad
This subroutine is the main sw radiation routine.
Definition: radsw_main.f:440
real(kind=kind_phys), dimension(ng24), public abso3a
real(kind=kind_phys), dimension(msb18, ng18), public absb
real(kind=kind_phys), dimension(nblow:nbhgh), public c0r
real(kind=kind_phys), dimension(msa28, ng28), public absa
This module sets up absorption coeffients for band 27: 29000-38000 cm-1 (low - o3; high - o3) ...
real(kind=kind_phys), dimension(msa18, ng18), public absa
real(kind=kind_phys), dimension(mfr20, ng20), public forref
real(kind=kind_phys), dimension(43, nblow:nbhgh), public asyice2
Definition: radsw_datatb.f:182
subroutine cldprop
This subroutine computes the cloud optical properties for each cloudy layer and g-point interval...
Definition: radsw_main.f:1429
define type construct for optional component downward fluxes at surface
Definition: radsw_param.f:110
real(kind=kind_phys), dimension(43, nblow:nbhgh), public extice2
Definition: radsw_datatb.f:182
real(kind=kind_phys), dimension(msf24, ng24), public selfref
real(kind=kind_phys), dimension(ng24, mfx24), public rayla
real(kind=kind_phys), dimension(msb28, ng28), public absb
real(kind=kind_phys), dimension(msb24, ng24), public absb
integer, parameter nbhgh
band range upper limit
Definition: radsw_param.f:130
This module sets up absorption coeffients for band 20: 5150-6150 cm-1 (low - h2o; high - h2o) ...
integer, parameter ntbmx
index upper limit of trans table
Definition: radsw_param.f:140
real(kind=kind_phys), public a0s
real(kind=kind_phys), dimension(ngmax, mfs01, mfb01), target, public sfluxref01
real(kind=kind_phys), dimension(ng26), public rayl
integer, parameter ngptsw
total number of g-point in all bands
Definition: radsw_param.f:134
real(kind=kind_phys), dimension(ng24), public abso3b
integer, save isubcsw
sub-column cloud approx flag in sw radiation
Definition: physparam.f:226
real(kind=kind_phys), dimension(msf22, ng22), public selfref
real(kind=kind_phys), dimension(58, nblow:nbhgh), public extliq1
Definition: radsw_datatb.f:180
This module contains coefficients of cloud optical properties for each of the spectral bands...
Definition: radsw_datatb.f:126
real(kind=kind_phys), dimension(nblow:nbhgh), public c0s
real(kind=kind_phys), dimension(msa23, ng23), public absa
This module sets up absorption coeffients for band 23: 8050-12850 cm-1 (low - h2o; high - nothing) ...
real(kind=kind_phys), dimension(msa22, ng22), public absa
real(kind=kind_phys), dimension(ng29), public absco2
define type construct for radiation fluxes at toa
Definition: radsw_param.f:70
real(kind=kind_phys), dimension(msb27, ng27), public absb
real(kind=kind_phys), dimension(ng25), public abso3a
integer, save iovrsw
cloud overlapping control flag for sw
Definition: physparam.f:197
real(kind=kind_phys), parameter stpfac
Definition: radsw_main.f:295
real(kind=kind_phys), dimension(msb29, ng29), public absb
This module contains some the most frequently used math and physics constants for gcm models...
Definition: physcons.f:29
real(kind=kind_phys), dimension(ngmax, mfs02, mfb02), target, public sfluxref02
integer, parameter iswrgas
sw rare gases effect control flag (ch4,n2o,o2,...) =0:no; =1:yes.
Definition: physparam.f:62
define type construct for radiation fluxes at surface
Definition: radsw_param.f:82
This module sets up absorption coeffients for band 18: 4000-4650 cm-1 (low - h2o,ch4; high - ch4) ...
real(kind=kind_phys), dimension(msb16, ng16), public absb
integer, parameter nbdsw
Definition: radsw_param.f:144
subroutine taumol22
The subroutine computes the optical depth in band 22: 7700-8050 cm-1 (low - h2o,o2; high - o2) ...
Definition: radsw_main.f:4441
real(kind=kind_phys), dimension(nblow:nbhgh), public b0s
real(kind=kind_phys), dimension(msb20, ng20), public absb
real(kind=kind_phys), dimension(msf23, ng23), public selfref
real(kind=kind_phys), parameter, public rayl
subroutine spcvrtm
This subroutine computes the shortwave radiative fluxes using two-stream method of h...
Definition: radsw_main.f:2809
integer, dimension(nblow:nbhgh) nspb
Definition: radsw_main.f:307
real(kind=kind_phys), dimension(43, nblow:nbhgh), public ssaice2
Definition: radsw_datatb.f:182
subroutine setcoef
This subroutine computes various coefficients needed in radiative transfer calculation.
Definition: radsw_main.f:1945
real(kind=kind_phys), dimension(msb22, ng22), public absb
real(kind=kind_phys), dimension(ng25), public abso3b
subroutine taumol23
The subroutine computes the optical depth in band 23: 8050-12850 cm-1 (low - h2o; high - nothing) ...
Definition: radsw_main.f:4547
integer, save iswcliq
sw optical property for liquid clouds =0:input cld opt depth, ignoring iswcice setting =1:input c...
Definition: physparam.f:68
real(kind=kind_phys), dimension(59) preflog
Definition: radsw_datatb.f:69
real(kind=kind_phys), dimension(msa21, ng21), public absa
integer, parameter ipsdsw0
Definition: radsw_main.f:349
This module defines commonly used control variables/parameters in physics related programs...
Definition: physparam.f:27
subroutine taumol26
The subroutine computes the optical depth in band 26: 22650-29000 cm-1 (low - nothing; high - nothing...
Definition: radsw_main.f:4756
real(kind=kind_phys), dimension(ng24), public raylb
integer, dimension(nblow:nbhgh) idxsfc
Definition: radsw_main.f:307
This module contains reference temperature and pressure.
Definition: radsw_datatb.f:58
real(kind=kind_phys), parameter, public rayl
integer, parameter iswrate
sw heating rate unit control flag =1:k/day; =2:k/second.
Definition: physparam.f:58
This module contains SW band parameters set up.
Definition: radsw_param.f:58
real(kind=kind_phys), parameter ftiny
Definition: radsw_main.f:296
integer, parameter nblow
band range lower limit
Definition: radsw_param.f:128
subroutine taumol19
The subroutine computes the optical depth in band 19: 4650-5150 cm-1 (low - h2o,co2; high - co2) ...
Definition: radsw_main.f:4147
real(kind=kind_phys), parameter amdo3
Definition: radsw_main.f:304
This module sets up absorption coeffients for band 26: 22650-29000 cm-1 (low - nothing; high - nothin...
subroutine taumol28
The subroutine computes the optical depth in band 28: 38000-50000 cm-1 (low - o3,o2; high - o3...
Definition: radsw_main.f:4850
real(kind=kind_phys), dimension(msb19, ng19), public absb
subroutine taumol18
The subroutine computes the optical depth in band 18: 4000-4650 cm-1 (low - h2o,ch4; high - ch4) ...
Definition: radsw_main.f:4052
real(kind=kind_phys), dimension(nblow:nbhgh), public strrat
real(kind=kind_phys), dimension(46, nblow:nbhgh), public extice3
Definition: radsw_datatb.f:184
real(kind=kind_phys), dimension(mfr16, ng16), public forref
real(kind=kind_phys), dimension(msf20, ng20), public selfref
real(kind=kind_phys), parameter, public rayl
real(kind=kind_phys), dimension(5), public fbari
Definition: radsw_datatb.f:186
real(kind=kind_phys), parameter, public rayl
real(kind=kind_phys), dimension(5), public dbari
Definition: radsw_datatb.f:186
real(kind=kind_phys), dimension(ng20), public absch4
real(kind=kind_phys), dimension(msf17, ng17), public selfref
integer, save ivflip
vertical profile indexing flag
Definition: physparam.f:224
This module sets up absorption coeffients for band 29: 820-2600 cm-1 (low - h2o; high - co2) ...
define type construct for optional radiation flux profiles
Definition: radsw_param.f:96
This module sets up absorption coeffients for band 17: 3250-4000 cm-1 (low - h2o,co2; high - h2o...
real(kind=kind_phys), dimension(mfr23, ng23), public forref
subroutine taumol29
The subroutine computes the optical depth in band 29: 820-2600 cm-1 (low - h2o; high - co2) ...
Definition: radsw_main.f:4957
integer, dimension(nblow:nbhgh) idxebc
Definition: radsw_main.f:307
real(kind=kind_phys), dimension(58, nblow:nbhgh), public asyliq1
Definition: radsw_datatb.f:180
real(kind=kind_phys), dimension(46, nblow:nbhgh), public asyice3
Definition: radsw_datatb.f:184
real(kind=kind_phys), dimension(msa27, ng27), public absa
subroutine, public rswinit
This subroutine initializes non-varying module variables, conversion factors, and look-up tables...
Definition: radsw_main.f:1272
real(kind=kind_phys), dimension(mfr21, ng21), public forref
subroutine taumol27
The subroutine computes the optical depth in band 27: 29000-38000 cm-1 (low - o3; high - o3) ...
Definition: radsw_main.f:4790
This module sets up absorption coeffients for band 22: 7700-8050 cm-1 (low - h2o, o2; high - o2) ...
real(kind=kind_phys), dimension(ngmax, mfs03, mfb03), target, public sfluxref03
real(kind=kind_phys), dimension(ng29), public absh2o
real(kind=kind_phys), parameter bpade
Definition: radsw_main.f:294
real(kind=kind_phys), parameter con_amw
molecular wght of water vapor (g/mol)
Definition: physcons.f:132
real(kind=kind_phys), dimension(mfr18, ng18), public forref
subroutine taumol
Definition: radsw_main.f:3593
real(kind=kind_phys), parameter s0
Definition: radsw_main.f:297
This module contains spectral distribution of solar radiation flux used to obtain the incoming solar ...
real(kind=kind_phys), dimension(msf16, ng16), public selfref
real(kind=kind_phys), dimension(msb17, ng17), public absb
subroutine taumol20
The subroutine computes the optical depth in band 20: 5150-6150 cm-1 (low - h2o; high - h2o) ...
Definition: radsw_main.f:4241
subroutine taumol21
The subroutine computes the optical depth in band 21: 6150-7700 cm-1 (low - h2o,co2; high - h2o...
Definition: radsw_main.f:4321
real(kind=kind_phys), dimension(ng27), public rayl
real(kind=kind_phys), parameter con_amd
molecular wght of dry air (g/mol)
Definition: physcons.f:130
real(kind=kind_phys), dimension(mfr19, ng19), public forref
This module sets up absorption coeffients for band 25: 16000-22650 cm-1 (low - h2o; high - nothing) ...
This module includes ncep&#39;s modifications of the rrtm-sw radiation code from aer inc.
Definition: radsw_main.f:260
real(kind=kind_phys), parameter, public rayl
real(kind=kind_phys), dimension(msa20, ng20), public absa
real(kind=kind_phys), dimension(58, nblow:nbhgh), public ssaliq1
Definition: radsw_datatb.f:180
real(kind=kind_phys), dimension(msa17, ng17), public absa
real(kind=kind_phys), parameter con_g
gravity (m/s2)
Definition: physcons.f:52
real(kind=kind_phys), dimension(mfr24, ng24), public forref
real(kind=kind_phys), parameter f_zero
Definition: radsw_main.f:299
real(kind=kind_phys), dimension(ng25), public rayl
integer, dimension(nblow:nbhgh) ngs
Definition: radsw_param.f:166
integer, dimension(nblow:nbhgh), public ibx
subroutine taumol25
The subroutine computes the optical depth in band 25: 16000-22650 cm-1 (low - h2o; high - nothing) ...
Definition: radsw_main.f:4702
integer, dimension(nblow:nbhgh) nspa
Definition: radsw_main.f:307
real(kind=kind_phys), dimension(46, nblow:nbhgh), public ssaice3
Definition: radsw_datatb.f:184
real(kind=kind_phys), dimension(msa25, ng25), public absa
real(kind=kind_phys), parameter, public scalekur
real(kind=kind_phys), parameter, public rayl
subroutine spcvrtc
This subroutine computes the shortwave radiative fluxes using two-stream method.
Definition: radsw_main.f:2100
real(kind=kind_phys), parameter, public givfac
real(kind=kind_phys), dimension(0:ntbmx) exp_tbl
Definition: radsw_main.f:340
integer, dimension(nblow:nbhgh) ng
Definition: radsw_param.f:153
real(kind=kind_phys) heatfac
Definition: radsw_main.f:345
real(kind=kind_phys), dimension(mfr17, ng17), public forref
real(kind=kind_phys), public a0r
This module sets up absorption coeffients for band 21: 6150-7700 cm-1 (low - h2o,co2; high - h2o...
real(kind=kind_phys), dimension(5), public cbari
Definition: radsw_datatb.f:186
integer, parameter iswmode
sw control flag for 2-stream transfer scheme =1:delta-eddington (joseph et al., 1976) =2:pifm (zd...
Definition: physparam.f:86
real(kind=kind_phys), dimension(msf18, ng18), public selfref
real(kind=kind_phys), parameter, public rayl
real(kind=kind_phys), parameter con_cp
spec heat air at p (J/kg/K)
Definition: physcons.f:73
integer, dimension(ngptsw) ngb
Definition: radsw_param.f:171
real(kind=kind_phys), parameter f_one
Definition: radsw_main.f:300
real(kind=kind_phys), parameter oneminus
Definition: radsw_main.f:293
integer, parameter nuvb
Definition: radsw_main.f:329
integer, parameter maxgas
max num of absorbing gases
Definition: radsw_param.f:138
real(kind=kind_phys), dimension(nblow:nbhgh), public b1s
integer, dimension(nblow:nbhgh), public layreffr
real(kind=kind_phys), parameter con_avgd
avogadro constant (1/mol)
Definition: physcons.f:125
real(kind=kind_phys), dimension(nblow:nbhgh), public b0r
This module sets up absorption coeffients for band 28: 38000-50000 cm-1 (low - o3,o2; high - o3,o2)
real(kind=kind_phys), dimension(msa29, ng29), public absa
real(kind=kind_phys), dimension(msa19, ng19), public absa
real(kind=kind_phys), dimension(5), public abari
Definition: radsw_datatb.f:186
This module sets up absorption coefficients for band 16: 2600-3250 cm-1 (low - h2o, ch4; high - ch4)
integer, save iswcice
sw optical property for ice clouds (only iswcliq>0) =0:not defined yet =1:input cip...
Definition: physparam.f:77
real(kind=kind_phys), dimension(msf21, ng21), public selfref
real(kind=kind_phys), dimension(59) tref
Definition: radsw_datatb.f:69
real(kind=kind_phys), dimension(5), public bbari
Definition: radsw_datatb.f:186
real(kind=kind_phys), dimension(nblow:nbhgh), public specwt
subroutine taumol24
The subroutine computes the optical depth in band 24: 12850-16000 cm-1 (low - h2o,o2; high - o2)
Definition: radsw_main.f:4608
subroutine taumol16
The subroutine computes the optical depth in band 16: 2600-3250 cm-1 (low - h2o,ch4; high - ch4) ...
Definition: radsw_main.f:3836
integer, save icldflg
cloud optical property scheme control flag
Definition: physparam.f:193
subroutine swflux
This subroutine computes the upward and downward radiation fluxes.
Definition: radsw_main.f:3488
real(kind=kind_phys), dimension(msa24, ng24), public absa