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