GFS Physics Documentation
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 
405 ! FEB 2017 A.Cheng - add odpth output, effective radius input !
406 !========================================!
407  module module_radsw_main !
408 !........................................!
409 !
410  use physparam, only : iswrate, iswrgas, iswcliq, iswcice, &
411  & isubcsw, icldflg, iovrsw, ivflip, &
412  & iswmode, kind_phys
413  use physcons, only : con_g, con_cp, con_avgd, con_amd, &
414  & con_amw, con_amo3
415 
417  use mersenne_twister, only : random_setseed, random_number, &
418  & random_stat
419  use module_radsw_ref, only : preflog, tref
421 !
422  implicit none
423 !
424  private
425 !
426 ! --- version tag and last revision date
427  character(40), parameter :: &
428  & VTAGSW='NCEP SW v5.1 Nov 2012 -RRTMG-SW v3.8 '
429 ! & VTAGSW='NCEP SW v5.0 Aug 2012 -RRTMG-SW v3.8 '
430 ! & VTAGSW='RRTMG-SW v3.8 Nov 2009'
431 ! & VTAGSW='RRTMG-SW v3.7 Nov 2009'
432 ! & VTAGSW='RRTMG-SW v3.61 Oct 2008'
433 ! & VTAGSW='RRTMG-SW v3.5 Oct 2008'
434 ! & VTAGSW='RRTM-SW 112v2.3 Apr 2007'
435 ! & VTAGSW='RRTM-SW 112v2.3 Mar 2005'
436 ! & VTAGSW='RRTM-SW 112v2.0 Jul 2004'
437 
439 
440  real (kind=kind_phys), parameter :: eps = 1.0e-6
441  real (kind=kind_phys), parameter :: oneminus= 1.0 - eps
443  real (kind=kind_phys), parameter :: bpade = 1.0/0.278
444  real (kind=kind_phys), parameter :: stpfac = 296.0/1013.0
445  real (kind=kind_phys), parameter :: ftiny = 1.0e-12
446  real (kind=kind_phys), parameter :: flimit = 1.0e-20
448  real (kind=kind_phys), parameter :: s0 = 1368.22
449 
450  real (kind=kind_phys), parameter :: f_zero = 0.0
451  real (kind=kind_phys), parameter :: f_one = 1.0
452 
454  real (kind=kind_phys), parameter :: amdw = con_amd/con_amw
455  real (kind=kind_phys), parameter :: amdo3 = con_amd/con_amo3
456 
458  integer, dimension(nblow:nbhgh) :: nspa, nspb
460  integer, dimension(nblow:nbhgh) :: idxsfc
462  integer, dimension(nblow:nbhgh) :: idxebc
463 
464  data nspa(:) / 9, 9, 9, 9, 1, 9, 9, 1, 9, 1, 0, 1, 9, 1 /
465  data nspb(:) / 1, 5, 1, 1, 1, 5, 1, 0, 1, 0, 0, 1, 5, 1 /
466 
467 ! data idxsfc(:) / 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1 / ! band index for sfc flux
468  data idxsfc(:) / 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1 / ! band index for sfc flux
469  data idxebc(:) / 5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 1, 1, 1, 5 / ! band index for cld prop
470 
471 ! --- band wavenumber intervals
472 ! real (kind=kind_phys), dimension(nblow:nbhgh):: wavenum1,wavenum2
473 ! data wavenum1(:) / &
474 ! & 2600.0, 3250.0, 4000.0, 4650.0, 5150.0, 6150.0, 7700.0, &
475 ! & 8050.0,12850.0,16000.0,22650.0,29000.0,38000.0, 820.0 /
476 ! data wavenum2(:) / &
477 ! 3250.0, 4000.0, 4650.0, 5150.0, 6150.0, 7700.0, 8050.0, &
478 ! & 12850.0,16000.0,22650.0,29000.0,38000.0,50000.0, 2600.0 /
479 ! real (kind=kind_phys), dimension(nblow:nbhgh) :: delwave
480 ! data delwave(:) / &
481 ! & 650.0, 750.0, 650.0, 500.0, 1000.0, 1550.0, 350.0, &
482 ! & 4800.0, 3150.0, 6650.0, 6350.0, 9000.0,12000.0, 1780.0 /
483 
485  integer, parameter :: nuvb = 27
486 
488  logical :: lhswb = .false.
489  logical :: lhsw0 = .false.
490  logical :: lflxprf= .false.
491  logical :: lfdncmp= .false.
492 
493 
495  real (kind=kind_phys) :: exp_tbl(0:ntbmx)
496 
497 
500  real (kind=kind_phys) :: heatfac
501 
502 
504  integer, parameter :: ipsdsw0 = 1
505 
506 ! --- public accessable subprograms
507 
508  public swrad, rswinit
509 
510 
511 ! =================
512  contains
513 ! =================
514 
585 !-----------------------------------
586  subroutine swrad &
587  & ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, & ! --- inputs
588  & clouds,icseed,aerosols,sfcalb, &
589  & cosz,solcon,nday,idxday, &
590  & npts, nlay, nlp1, lprnt, &
591  & hswc,topflx,sfcflx, & ! --- outputs
592  & hsw0,hswb,flxprf,fdncmp & ! --- optional
593  & )
595 ! ==================== defination of variables ==================== !
596 ! !
597 ! input variables: !
598 ! plyr (npts,nlay) : model layer mean pressure in mb !
599 ! plvl (npts,nlp1) : model level pressure in mb !
600 ! tlyr (npts,nlay) : model layer mean temperature in k !
601 ! tlvl (npts,nlp1) : model level temperature in k (not in use) !
602 ! qlyr (npts,nlay) : layer specific humidity in gm/gm *see inside !
603 ! olyr (npts,nlay) : layer ozone concentration in gm/gm !
604 ! gasvmr(npts,nlay,:): atmospheric constent gases: !
605 ! (check module_radiation_gases for definition) !
606 ! gasvmr(:,:,1) - co2 volume mixing ratio !
607 ! gasvmr(:,:,2) - n2o volume mixing ratio !
608 ! gasvmr(:,:,3) - ch4 volume mixing ratio !
609 ! gasvmr(:,:,4) - o2 volume mixing ratio !
610 ! gasvmr(:,:,5) - co volume mixing ratio (not used) !
611 ! gasvmr(:,:,6) - cfc11 volume mixing ratio (not used) !
612 ! gasvmr(:,:,7) - cfc12 volume mixing ratio (not used) !
613 ! gasvmr(:,:,8) - cfc22 volume mixing ratio (not used) !
614 ! gasvmr(:,:,9) - ccl4 volume mixing ratio (not used) !
615 ! clouds(npts,nlay,:): cloud profile !
616 ! (check module_radiation_clouds for definition) !
617 ! --- for iswcliq > 0 --- !
618 ! clouds(:,:,1) - layer total cloud fraction !
619 ! clouds(:,:,2) - layer in-cloud liq water path (g/m**2) !
620 ! clouds(:,:,3) - mean eff radius for liq cloud (micron) !
621 ! clouds(:,:,4) - layer in-cloud ice water path (g/m**2) !
622 ! clouds(:,:,5) - mean eff radius for ice cloud (micron) !
623 ! clouds(:,:,6) - layer rain drop water path (g/m**2) !
624 ! clouds(:,:,7) - mean eff radius for rain drop (micron) !
625 ! clouds(:,:,8) - layer snow flake water path (g/m**2) !
626 ! clouds(:,:,9) - mean eff radius for snow flake (micron) !
627 ! --- for iswcliq = 0 --- !
628 ! clouds(:,:,1) - layer total cloud fraction !
629 ! clouds(:,:,2) - layer cloud optical depth !
630 ! clouds(:,:,3) - layer cloud single scattering albedo !
631 ! clouds(:,:,4) - layer cloud asymmetry factor !
632 ! icseed(npts) : auxiliary special cloud related array !
633 ! when module variable isubcsw=2, it provides !
634 ! permutation seed for each column profile that !
635 ! are used for generating random numbers. !
636 ! when isubcsw /=2, it will not be used. !
637 ! aerosols(npts,nlay,nbdsw,:) : aerosol optical properties !
638 ! (check module_radiation_aerosols for definition) !
639 ! (:,:,:,1) - optical depth !
640 ! (:,:,:,2) - single scattering albedo !
641 ! (:,:,:,3) - asymmetry parameter !
642 ! sfcalb(npts, : ) : surface albedo in fraction !
643 ! (check module_radiation_surface for definition) !
644 ! ( :, 1 ) - near ir direct beam albedo !
645 ! ( :, 2 ) - near ir diffused albedo !
646 ! ( :, 3 ) - uv+vis direct beam albedo !
647 ! ( :, 4 ) - uv+vis diffused albedo !
648 ! cosz (npts) : cosine of solar zenith angle !
649 ! solcon : solar constant (w/m**2) !
650 ! NDAY : num of daytime points !
651 ! idxday(npts) : index array for daytime points !
652 ! npts : number of horizontal points !
653 ! nlay,nlp1 : vertical layer/lavel numbers !
654 ! lprnt : logical check print flag !
655 ! !
656 ! output variables: !
657 ! hswc (npts,nlay): total sky heating rates (k/sec or k/day) !
658 ! topflx(npts) : radiation fluxes at toa (w/m**2), components: !
659 ! (check module_radsw_parameters for definition) !
660 ! upfxc - total sky upward flux at toa !
661 ! dnflx - total sky downward flux at toa !
662 ! upfx0 - clear sky upward flux at toa !
663 ! sfcflx(npts) : radiation fluxes at sfc (w/m**2), components: !
664 ! (check module_radsw_parameters for definition) !
665 ! upfxc - total sky upward flux at sfc !
666 ! dnfxc - total sky downward flux at sfc !
667 ! upfx0 - clear sky upward flux at sfc !
668 ! dnfx0 - clear sky downward flux at sfc !
669 ! !
670 !!optional outputs variables: !
671 ! hswb(npts,nlay,nbdsw): spectral band total sky heating rates !
672 ! hsw0 (npts,nlay): clear sky heating rates (k/sec or k/day) !
673 ! flxprf(npts,nlp1): level radiation fluxes (w/m**2), components: !
674 ! (check module_radsw_parameters for definition) !
675 ! dnfxc - total sky downward flux at interface !
676 ! upfxc - total sky upward flux at interface !
677 ! dnfx0 - clear sky downward flux at interface !
678 ! upfx0 - clear sky upward flux at interface !
679 ! fdncmp(npts) : component surface downward fluxes (w/m**2): !
680 ! (check module_radsw_parameters for definition) !
681 ! uvbfc - total sky downward uv-b flux at sfc !
682 ! uvbf0 - clear sky downward uv-b flux at sfc !
683 ! nirbm - downward surface nir direct beam flux !
684 ! nirdf - downward surface nir diffused flux !
685 ! visbm - downward surface uv+vis direct beam flux !
686 ! visdf - downward surface uv+vis diffused flux !
687 ! !
688 ! external module variables: (in physparam) !
689 ! iswrgas - control flag for rare gases (ch4,n2o,o2, etc.) !
690 ! =0: do not include rare gases !
691 ! >0: include all rare gases !
692 ! iswcliq - control flag for liq-cloud optical properties !
693 ! =0: input cloud optical depth, fixed ssa, asy !
694 ! =1: use hu and stamnes(1993) method for liq cld !
695 ! =2: not used !
696 ! iswcice - control flag for ice-cloud optical properties !
697 ! *** if iswcliq==0, iswcice is ignored !
698 ! =1: use ebert and curry (1992) scheme for ice clouds !
699 ! =2: use streamer v3.0 (2001) method for ice clouds !
700 ! =3: use fu's method (1996) for ice clouds !
701 ! iswmode - control flag for 2-stream transfer scheme !
702 ! =1; delta-eddington (joseph et al., 1976) !
703 ! =2: pifm (zdunkowski et al., 1980) !
704 ! =3: discrete ordinates (liou, 1973) !
705 ! isubcsw - sub-column cloud approximation control flag !
706 ! =0: no sub-col cld treatment, use grid-mean cld quantities !
707 ! =1: mcica sub-col, prescribed seeds to get random numbers !
708 ! =2: mcica sub-col, providing array icseed for random numbers!
709 ! iovrsw - cloud overlapping control flag !
710 ! =0: random overlapping clouds !
711 ! =1: maximum/random overlapping clouds !
712 ! =2: maximum overlap cloud !
713 ! ivflip - control flg for direction of vertical index !
714 ! =0: index from toa to surface !
715 ! =1: index from surface to toa !
716 ! !
717 ! module parameters, control variables: !
718 ! nblow,nbhgh - lower and upper limits of spectral bands !
719 ! maxgas - maximum number of absorbing gaseous !
720 ! ngptsw - total number of g-point subintervals !
721 ! ng## - number of g-points in band (##=16-29) !
722 ! ngb(ngptsw) - band indices for each g-point !
723 ! bpade - pade approximation constant (1/0.278) !
724 ! nspa,nspb(nblow:nbhgh) !
725 ! - number of lower/upper ref atm's per band !
726 ! ipsdsw0 - permutation seed for mcica sub-col clds !
727 ! !
728 ! major local variables: !
729 ! pavel (nlay) - layer pressures (mb) !
730 ! delp (nlay) - layer pressure thickness (mb) !
731 ! tavel (nlay) - layer temperatures (k) !
732 ! coldry (nlay) - dry air column amount !
733 ! (1.e-20*molecules/cm**2) !
734 ! cldfrc (nlay) - layer cloud fraction (norm by tot cld) !
735 ! cldfmc (nlay,ngptsw) - layer cloud fraction for g-point !
736 ! taucw (nlay,nbdsw) - cloud optical depth !
737 ! ssacw (nlay,nbdsw) - cloud single scattering albedo (weighted) !
738 ! asycw (nlay,nbdsw) - cloud asymmetry factor (weighted) !
739 ! tauaer (nlay,nbdsw) - aerosol optical depths !
740 ! ssaaer (nlay,nbdsw) - aerosol single scattering albedo !
741 ! asyaer (nlay,nbdsw) - aerosol asymmetry factor !
742 ! colamt (nlay,maxgas) - column amounts of absorbing gases !
743 ! 1 to maxgas are for h2o, co2, o3, n2o, !
744 ! ch4, o2, co, respectively (mol/cm**2) !
745 ! facij (nlay) - indicator of interpolation factors !
746 ! =0/1: indicate lower/higher temp & height !
747 ! selffac(nlay) - scale factor for self-continuum, equals !
748 ! (w.v. density)/(atm density at 296K,1013 mb) !
749 ! selffrac(nlay) - factor for temp interpolation of ref !
750 ! self-continuum data !
751 ! indself(nlay) - index of the lower two appropriate ref !
752 ! temp for the self-continuum interpolation !
753 ! forfac (nlay) - scale factor for w.v. foreign-continuum !
754 ! forfrac(nlay) - factor for temp interpolation of ref !
755 ! w.v. foreign-continuum data !
756 ! indfor (nlay) - index of the lower two appropriate ref !
757 ! temp for the foreign-continuum interp !
758 ! laytrop - layer at which switch is made from one !
759 ! combination of key species to another !
760 ! jp(nlay),jt(nlay),jt1(nlay) !
761 ! - lookup table indexes !
762 ! flxucb(nlp1,nbdsw) - spectral bnd total-sky upward flx (w/m2) !
763 ! flxdcb(nlp1,nbdsw) - spectral bnd total-sky downward flx (w/m2)!
764 ! flxu0b(nlp1,nbdsw) - spectral bnd clear-sky upward flx (w/m2) !
765 ! flxd0b(nlp1,nbdsw) - spectral b d clear-sky downward flx (w/m2)!
766 ! !
767 ! !
768 ! ===================== end of definitions ==================== !
769 
770 ! --- inputs:
771  integer, intent(in) :: npts, nlay, nlp1, NDAY
772 
773  integer, dimension(:), intent(in) :: idxday, icseed
774 
775  logical, intent(in) :: lprnt
776 
777  real (kind=kind_phys), dimension(npts,nlp1), intent(in) :: &
778  & plvl, tlvl
779  real (kind=kind_phys), dimension(npts,nlay), intent(in) :: &
780  & plyr, tlyr, qlyr, olyr
781  real (kind=kind_phys), dimension(npts,4), intent(in) :: sfcalb
782 
783  real (kind=kind_phys), dimension(npts,nlay,9),intent(in):: gasvmr
784  real (kind=kind_phys), dimension(npts,nlay,11):: clouds
785  real (kind=kind_phys), dimension(npts,nlay,nbdsw,3),intent(in):: &
786  & aerosols
787 
788  real (kind=kind_phys), intent(in) :: cosz(npts), solcon
789 
790 ! --- outputs:
791  real (kind=kind_phys), dimension(npts,nlay), intent(out) :: hswc
792 
793  type(topfsw_type), dimension(npts), intent(out) :: topflx
794  type(sfcfsw_type), dimension(npts), intent(out) :: sfcflx
795 
796 !! --- optional outputs:
797  real (kind=kind_phys), dimension(npts,nlay,nbdsw), optional, &
798  & intent(out) :: hswb
799 
800  real (kind=kind_phys), dimension(npts,nlay), optional, &
801  & intent(out) :: hsw0
802  type (profsw_type), dimension(npts,nlp1), optional, &
803  & intent(out) :: flxprf
804  type (cmpfsw_type), dimension(npts), optional, &
805  & intent(out) :: fdncmp
806 
807 ! --- locals:
808  real (kind=kind_phys), dimension(nlay,ngptsw) :: cldfmc, &
809  & taug, taur
810  real (kind=kind_phys), dimension(nlp1,nbdsw):: fxupc, fxdnc, &
811  & fxup0, fxdn0
812 
813  real (kind=kind_phys), dimension(nlay,nbdsw) :: &
814  & tauae, ssaae, asyae, taucw, ssacw, asycw
815 
816  real (kind=kind_phys), dimension(ngptsw) :: sfluxzen
817 
818  real (kind=kind_phys), dimension(nlay) :: cldfrc, delp, &
819  & pavel, tavel, coldry, colmol, h2ovmr, o3vmr, temcol, &
820  & cliqp, reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, &
821  & cfrac, fac00, fac01, fac10, fac11, forfac, forfrac, &
822  & selffac, selffrac, rfdelp
823 
824  real (kind=kind_phys), dimension(nlp1) :: fnet, flxdc, flxuc, &
825  & flxd0, flxu0
826 
827  real (kind=kind_phys), dimension(2) :: albbm, albdf, sfbmc, &
828  & sfbm0, sfdfc, sfdf0
829 
830  real (kind=kind_phys) :: cosz1, sntz1, tem0, tem1, tem2, s0fac, &
831  & ssolar, zcf0, zcf1, ftoau0, ftoauc, ftoadc, &
832  & fsfcu0, fsfcuc, fsfcd0, fsfcdc, suvbfc, suvbf0
833 
834 ! --- column amount of absorbing gases:
835 ! (:,m) m = 1-h2o, 2-co2, 3-o3, 4-n2o, 5-ch4, 6-o2, 7-co
836  real (kind=kind_phys) :: colamt(nlay,maxgas)
837 
838  integer, dimension(npts) :: ipseed
839  integer, dimension(nlay) :: indfor, indself, jp, jt, jt1
840 
841  integer :: i, ib, ipt, j1, k, kk, laytrop, mb
842 !
843 !===> ... begin here
844 !
845 
846  lhswb = present ( hswb )
847  lhsw0 = present ( hsw0 )
848  lflxprf= present ( flxprf )
849  lfdncmp= present ( fdncmp )
850 
852 ! *** s0, the solar constant at toa in w/m**2, is hard-coded with
853 ! each spectra band, the total flux is about 1368.22 w/m**2.
854 
855  s0fac = solcon / s0
856 
858 
859  hswc(:,:) = f_zero
860  topflx = topfsw_type( f_zero, f_zero, f_zero )
861  sfcflx = sfcfsw_type( f_zero, f_zero, f_zero, f_zero )
862 
863 !! --- ... initial optional outputs
864  if ( lflxprf ) then
865  flxprf = profsw_type( f_zero, f_zero, f_zero, f_zero )
866  endif
867 
868  if ( lfdncmp ) then
869  fdncmp = cmpfsw_type(f_zero,f_zero,f_zero,f_zero,f_zero,f_zero)
870  endif
871 
872  if ( lhsw0 ) then
873  hsw0(:,:) = f_zero
874  endif
875 
876  if ( lhswb ) then
877  hswb(:,:,:) = f_zero
878  endif
879 
882 
883  if ( isubcsw == 1 ) then ! advance prescribed permutation seed
884  do i = 1, npts
885  ipseed(i) = ipsdsw0 + i
886  enddo
887  elseif ( isubcsw == 2 ) then ! use input array of permutaion seeds
888  do i = 1, npts
889  ipseed(i) = icseed(i)
890  enddo
891  endif
892 
893  if ( lprnt ) then
894  write(0,*)' In radsw, isubcsw, ipsdsw0,ipseed =', &
895  & isubcsw, ipsdsw0, ipseed
896  endif
897 
898 ! --- ... loop over each daytime grid point
899 
900  lab_do_ipt : do ipt = 1, nday
901 
902  j1 = idxday(ipt)
903 
904  cosz1 = cosz(j1)
905  sntz1 = f_one / cosz(j1)
906  ssolar = s0fac * cosz(j1)
907 
909  albbm(1) = sfcalb(j1,1)
910  albdf(1) = sfcalb(j1,2)
911  albbm(2) = sfcalb(j1,3)
912  albdf(2) = sfcalb(j1,4)
913 
915 ! the vertical index of internal array is from surface to top
916 
917  if (ivflip == 0) then ! input from toa to sfc
918 
919  tem1 = 100.0 * con_g
920  tem2 = 1.0e-20 * 1.0e3 * con_avgd
921 
922  do k = 1, nlay
923  kk = nlp1 - k
924  pavel(k) = plyr(j1,kk)
925  tavel(k) = tlyr(j1,kk)
926  delp(k) = plvl(j1,kk+1) - plvl(j1,kk)
932 
933 !test use
934 ! h2ovmr(k)= max(f_zero,qlyr(j1,kk)*amdw) ! input mass mixing ratio
935 ! h2ovmr(k)= max(f_zero,qlyr(j1,kk)) ! input vol mixing ratio
936 ! o3vmr (k)= max(f_zero,olyr(j1,kk)) ! input vol mixing ratio
937 !ncep model use
938  h2ovmr(k)= max(f_zero,qlyr(j1,kk)*amdw/(f_one-qlyr(j1,kk))) ! input specific humidity
939  o3vmr(k)= max(f_zero,olyr(j1,kk)*amdo3) ! input mass mixing ratio
940 
941  tem0 = (f_one - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
942  coldry(k) = tem2 * delp(k) / (tem1*tem0*(f_one + h2ovmr(k)))
943  temcol(k) = 1.0e-12 * coldry(k)
944 
945  colamt(k,1) = max(f_zero, coldry(k)*h2ovmr(k)) ! h2o
946  colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(j1,kk,1)) ! co2
947  colamt(k,3) = max(f_zero, coldry(k)*o3vmr(k)) ! o3
948  colmol(k) = coldry(k) + colamt(k,1)
949  enddo
950 
951 ! --- ... set up gas column amount, convert from volume mixing ratio
952 ! to molec/cm2 based on coldry (scaled to 1.0e-20)
953 
954  if (iswrgas > 0) then
955  do k = 1, nlay
956  kk = nlp1 - k
957  colamt(k,4) = max(temcol(k), coldry(k)*gasvmr(j1,kk,2)) ! n2o
958  colamt(k,5) = max(temcol(k), coldry(k)*gasvmr(j1,kk,3)) ! ch4
959  colamt(k,6) = max(temcol(k), coldry(k)*gasvmr(j1,kk,4)) ! o2
960 ! colamt(k,7) = max(temcol(k), coldry(k)*gasvmr(j1,kk,5)) ! co - notused
961  enddo
962  else
963  do k = 1, nlay
964  colamt(k,4) = temcol(k) ! n2o
965  colamt(k,5) = temcol(k) ! ch4
966  colamt(k,6) = temcol(k) ! o2
967 ! colamt(k,7) = temcol(k) ! co - notused
968  enddo
969  endif
970 
972 
973  do k = 1, nlay
974  kk = nlp1 - k
975  do ib = 1, nbdsw
976  tauae(k,ib) = aerosols(j1,kk,ib,1)
977  ssaae(k,ib) = aerosols(j1,kk,ib,2)
978  asyae(k,ib) = aerosols(j1,kk,ib,3)
979  enddo
980  enddo
981 
983  if (iswcliq > 0) then ! use prognostic cloud method
984  do k = 1, nlay
985  kk = nlp1 - k
986  cfrac(k) = clouds(j1,kk,1) ! cloud fraction
987  cliqp(k) = clouds(j1,kk,2) ! cloud liq path
988  reliq(k) = clouds(j1,kk,3) ! liq partical effctive radius
989  cicep(k) = clouds(j1,kk,4) ! cloud ice path
990  reice(k) = clouds(j1,kk,5) ! ice partical effctive radius
991  cdat1(k) = clouds(j1,kk,6) ! cloud rain drop path
992  cdat2(k) = clouds(j1,kk,7) ! rain partical effctive radius
993  cdat3(k) = clouds(j1,kk,8) ! cloud snow path
994  cdat4(k) = clouds(j1,kk,9) ! snow partical effctive radius
995  enddo
996  else ! use diagnostic cloud method
997  do k = 1, nlay
998  kk = nlp1 - k
999  cfrac(k) = clouds(j1,kk,1) ! cloud fraction
1000  cdat1(k) = clouds(j1,kk,2) ! cloud optical depth
1001  cdat2(k) = clouds(j1,kk,3) ! cloud single scattering albedo
1002  cdat3(k) = clouds(j1,kk,4) ! cloud asymmetry factor
1003  enddo
1004  endif ! end if_iswcliq
1005 
1006  else ! input from sfc to toa
1007 
1008  tem1 = 100.0 * con_g
1009  tem2 = 1.0e-20 * 1.0e3 * con_avgd
1010 
1011  do k = 1, nlay
1012  pavel(k) = plyr(j1,k)
1013  tavel(k) = tlyr(j1,k)
1014  delp(k) = plvl(j1,k) - plvl(j1,k+1)
1015 
1016 ! --- ... set absorber amount
1017 !test use
1018 ! h2ovmr(k)= max(f_zero,qlyr(j1,k)*amdw) ! input mass mixing ratio
1019 ! h2ovmr(k)= max(f_zero,qlyr(j1,k)) ! input vol mixing ratio
1020 ! o3vmr (k)= max(f_zero,olyr(j1,k)) ! input vol mixing ratio
1021 !ncep model use
1022  h2ovmr(k)= max(f_zero,qlyr(j1,k)*amdw/(f_one-qlyr(j1,k))) ! input specific humidity
1023  o3vmr(k)= max(f_zero,olyr(j1,k)*amdo3) ! input mass mixing ratio
1024 
1025  tem0 = (f_one - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
1026  coldry(k) = tem2 * delp(k) / (tem1*tem0*(f_one + h2ovmr(k)))
1027  temcol(k) = 1.0e-12 * coldry(k)
1028 
1029  colamt(k,1) = max(f_zero, coldry(k)*h2ovmr(k)) ! h2o
1030  colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(j1,k,1)) ! co2
1031  colamt(k,3) = max(f_zero, coldry(k)*o3vmr(k)) ! o3
1032  colmol(k) = coldry(k) + colamt(k,1)
1033  enddo
1034 
1035 
1036  if (lprnt) then
1037  if (ipt == 1) then
1038  write(0,*)' pavel=',pavel
1039  write(0,*)' tavel=',tavel
1040  write(0,*)' delp=',delp
1041  write(0,*)' h2ovmr=',h2ovmr*1000
1042  write(0,*)' o3vmr=',o3vmr*1000000
1043  endif
1044  endif
1045 
1046 ! --- ... set up gas column amount, convert from volume mixing ratio
1047 ! to molec/cm2 based on coldry (scaled to 1.0e-20)
1048 
1049  if (iswrgas > 0) then
1050  do k = 1, nlay
1051  colamt(k,4) = max(temcol(k), coldry(k)*gasvmr(j1,k,2)) ! n2o
1052  colamt(k,5) = max(temcol(k), coldry(k)*gasvmr(j1,k,3)) ! ch4
1053  colamt(k,6) = max(temcol(k), coldry(k)*gasvmr(j1,k,4)) ! o2
1054 ! colamt(k,7) = max(temcol(k), coldry(k)*gasvmr(j1,k,5)) ! co - notused
1055  enddo
1056  else
1057  do k = 1, nlay
1058  colamt(k,4) = temcol(k) ! n2o
1059  colamt(k,5) = temcol(k) ! ch4
1060  colamt(k,6) = temcol(k) ! o2
1061 ! colamt(k,7) = temcol(k) ! co - notused
1062  enddo
1063  endif
1064 
1065 ! --- ... set aerosol optical properties
1066 
1067  do ib = 1, nbdsw
1068  do k = 1, nlay
1069  tauae(k,ib) = aerosols(j1,k,ib,1)
1070  ssaae(k,ib) = aerosols(j1,k,ib,2)
1071  asyae(k,ib) = aerosols(j1,k,ib,3)
1072  enddo
1073  enddo
1074 
1075  if (iswcliq > 0) then ! use prognostic cloud method
1076  do k = 1, nlay
1077  cfrac(k) = clouds(j1,k,1) ! cloud fraction
1078  cliqp(k) = clouds(j1,k,2) ! cloud liq path
1079  reliq(k) = clouds(j1,k,3) ! liq partical effctive radius
1080  cicep(k) = clouds(j1,k,4) ! cloud ice path
1081  reice(k) = clouds(j1,k,5) ! ice partical effctive radius
1082  cdat1(k) = clouds(j1,k,6) ! cloud rain drop path
1083  cdat2(k) = clouds(j1,k,7) ! rain partical effctive radius
1084  cdat3(k) = clouds(j1,k,8) ! cloud snow path
1085  cdat4(k) = clouds(j1,k,9) ! snow partical effctive radius
1086  enddo
1087  else ! use diagnostic cloud method
1088  do k = 1, nlay
1089  cfrac(k) = clouds(j1,k,1) ! cloud fraction
1090  cdat1(k) = clouds(j1,k,2) ! cloud optical depth
1091  cdat2(k) = clouds(j1,k,3) ! cloud single scattering albedo
1092  cdat3(k) = clouds(j1,k,4) ! cloud asymmetry factor
1093  enddo
1094  endif ! end if_iswcliq
1095 
1096  endif ! if_ivflip
1097 
1102 
1103  zcf0 = f_one
1104  zcf1 = f_one
1105  if (iovrsw == 0) then ! random overlapping
1106  do k = 1, nlay
1107  zcf0 = zcf0 * (f_one - cfrac(k))
1108  enddo
1109  else if (iovrsw == 1) then ! max/ran overlapping
1110  do k = 1, nlay
1111  if (cfrac(k) > ftiny) then ! cloudy layer
1112  zcf1 = min( zcf1, f_one-cfrac(k) )
1113  elseif (zcf1 < f_one) then ! clear layer
1114  zcf0 = zcf0 * zcf1
1115  zcf1 = f_one
1116  endif
1117  enddo
1118  zcf0 = zcf0 * zcf1
1119  else if (iovrsw == 2) then ! maximum overlapping
1120  do k = 1, nlay
1121  zcf0 = min( zcf0, f_one-cfrac(k) )
1122  enddo
1123  endif
1124 
1125  if (zcf0 <= ftiny) zcf0 = f_zero
1126  if (zcf0 > oneminus) zcf0 = f_one
1127  zcf1 = f_one - zcf0
1128 
1131 
1132  if (zcf1 > f_zero) then ! cloudy sky column
1133 
1134  call cldprop &
1135 ! --- inputs:
1136  & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, &
1137  & zcf1, nlay, ipseed(j1), &
1138 ! --- outputs:
1139  & taucw, ssacw, asycw, cldfrc, cldfmc &
1140  & )
1141 
1142  else ! clear sky column
1143  cldfrc(:) = f_zero
1144  cldfmc(:,:)= f_zero
1145  do i = 1, nbdsw
1146  do k = 1, nlay
1147  taucw(k,i) = f_zero
1148  ssacw(k,i) = f_zero
1149  asycw(k,i) = f_zero
1150  enddo
1151  enddo
1152  endif ! end if_zcf1_block
1153  do k = 1, nlay
1154  clouds(j1,k,10) = taucw(k,10)
1155  end do
1158  call setcoef &
1159 ! --- inputs:
1160  & ( pavel,tavel,h2ovmr, nlay,nlp1, &
1161 ! --- outputs:
1162  & laytrop,jp,jt,jt1,fac00,fac01,fac10,fac11, &
1163  & selffac,selffrac,indself,forfac,forfrac,indfor &
1164  & )
1165 
1168  call taumol &
1169 ! --- inputs:
1170  & ( colamt,colmol,fac00,fac01,fac10,fac11,jp,jt,jt1,laytrop, &
1171  & forfac,forfrac,indfor,selffac,selffrac,indself, nlay, &
1172 ! --- outputs:
1173  & sfluxzen, taug, taur &
1174  & )
1175 
1181 
1182  if ( isubcsw <= 0 ) then ! use standard cloud scheme
1183 
1184  call spcvrtc &
1185 ! --- inputs:
1186  & ( ssolar,cosz1,sntz1,albbm,albdf,sfluxzen,cldfrc, &
1187  & zcf1,zcf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
1188  & nlay, nlp1, &
1189 ! --- outputs:
1190  & fxupc,fxdnc,fxup0,fxdn0, &
1191  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
1192  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
1193  & )
1194 
1195  else ! use mcica cloud scheme
1196 
1197  call spcvrtm &
1198 ! --- inputs:
1199  & ( ssolar,cosz1,sntz1,albbm,albdf,sfluxzen,cldfmc, &
1200  & zcf1,zcf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
1201  & nlay, nlp1, &
1202 ! --- outputs:
1203  & fxupc,fxdnc,fxup0,fxdn0, &
1204  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
1205  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
1206  & )
1207 
1208  endif
1209 
1211 ! --- ... sum up total spectral fluxes for total-sky
1212 
1213  do k = 1, nlp1
1214  flxuc(k) = f_zero
1215  flxdc(k) = f_zero
1216 
1217  do ib = 1, nbdsw
1218  flxuc(k) = flxuc(k) + fxupc(k,ib)
1219  flxdc(k) = flxdc(k) + fxdnc(k,ib)
1220  enddo
1221  enddo
1222 
1223 !! --- ... optional clear sky fluxes
1224 
1225  if ( lhsw0 .or. lflxprf ) then
1226  do k = 1, nlp1
1227  flxu0(k) = f_zero
1228  flxd0(k) = f_zero
1229 
1230  do ib = 1, nbdsw
1231  flxu0(k) = flxu0(k) + fxup0(k,ib)
1232  flxd0(k) = flxd0(k) + fxdn0(k,ib)
1233  enddo
1234  enddo
1235  endif
1236 
1237 ! --- ... prepare for final outputs
1238 
1239  do k = 1, nlay
1240  rfdelp(k) = heatfac / delp(k)
1241  enddo
1242 
1243  if ( lfdncmp ) then
1244 !! --- ... optional uv-b surface downward flux
1245  fdncmp(j1)%uvbf0 = suvbf0
1246  fdncmp(j1)%uvbfc = suvbfc
1247 
1248 !! --- ... optional beam and diffuse sfc fluxes
1249  fdncmp(j1)%nirbm = sfbmc(1)
1250  fdncmp(j1)%nirdf = sfdfc(1)
1251  fdncmp(j1)%visbm = sfbmc(2)
1252  fdncmp(j1)%visdf = sfdfc(2)
1253  endif ! end if_lfdncmp
1254 
1255 ! --- ... toa and sfc fluxes
1256 
1257  topflx(j1)%upfxc = ftoauc
1258  topflx(j1)%dnfxc = ftoadc
1259  topflx(j1)%upfx0 = ftoau0
1260 
1261  sfcflx(j1)%upfxc = fsfcuc
1262  sfcflx(j1)%dnfxc = fsfcdc
1263  sfcflx(j1)%upfx0 = fsfcu0
1264  sfcflx(j1)%dnfx0 = fsfcd0
1265 
1266  if (ivflip == 0) then ! output from toa to sfc
1267 
1268 ! --- ... compute heating rates
1269 
1270  fnet(1) = flxdc(1) - flxuc(1)
1271 
1272  do k = 2, nlp1
1273  kk = nlp1 - k + 1
1274  fnet(k) = flxdc(k) - flxuc(k)
1275  hswc(j1,kk) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1276  enddo
1277 
1278 !! --- ... optional flux profiles
1279 
1280  if ( lflxprf ) then
1281  do k = 1, nlp1
1282  kk = nlp1 - k + 1
1283  flxprf(j1,kk)%upfxc = flxuc(k)
1284  flxprf(j1,kk)%dnfxc = flxdc(k)
1285  flxprf(j1,kk)%upfx0 = flxu0(k)
1286  flxprf(j1,kk)%dnfx0 = flxd0(k)
1287  enddo
1288  endif
1289 
1290 !! --- ... optional clear sky heating rates
1291 
1292  if ( lhsw0 ) then
1293  fnet(1) = flxd0(1) - flxu0(1)
1294 
1295  do k = 2, nlp1
1296  kk = nlp1 - k + 1
1297  fnet(k) = flxd0(k) - flxu0(k)
1298  hsw0(j1,kk) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1299  enddo
1300  endif
1301 
1302 !! --- ... optional spectral band heating rates
1303 
1304  if ( lhswb ) then
1305  do mb = 1, nbdsw
1306  fnet(1) = fxdnc(1,mb) - fxupc(1,mb)
1307 
1308  do k = 2, nlp1
1309  kk = nlp1 - k + 1
1310  fnet(k) = fxdnc(k,mb) - fxupc(k,mb)
1311  hswb(j1,kk,mb) = (fnet(k) - fnet(k-1)) * rfdelp(k-1)
1312  enddo
1313  enddo
1314  endif
1315 
1316  else ! output from sfc to toa
1317 
1318 ! --- ... compute heating rates
1319 
1320  fnet(1) = flxdc(1) - flxuc(1)
1321 
1322  do k = 2, nlp1
1323  fnet(k) = flxdc(k) - flxuc(k)
1324  hswc(j1,k-1) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1325  enddo
1326 
1327 !! --- ... optional flux profiles
1328 
1329  if ( lflxprf ) then
1330  do k = 1, nlp1
1331  flxprf(j1,k)%upfxc = flxuc(k)
1332  flxprf(j1,k)%dnfxc = flxdc(k)
1333  flxprf(j1,k)%upfx0 = flxu0(k)
1334  flxprf(j1,k)%dnfx0 = flxd0(k)
1335  enddo
1336  endif
1337 
1338 !! --- ... optional clear sky heating rates
1339 
1340  if ( lhsw0 ) then
1341  fnet(1) = flxd0(1) - flxu0(1)
1342 
1343  do k = 2, nlp1
1344  fnet(k) = flxd0(k) - flxu0(k)
1345  hsw0(j1,k-1) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1346  enddo
1347  endif
1348 
1349 !! --- ... optional spectral band heating rates
1350 
1351  if ( lhswb ) then
1352  do mb = 1, nbdsw
1353  fnet(1) = fxdnc(1,mb) - fxupc(1,mb)
1354 
1355  do k = 1, nlay
1356  fnet(k+1) = fxdnc(k+1,mb) - fxupc(k+1,mb)
1357  hswb(j1,k,mb) = (fnet(k+1) - fnet(k)) * rfdelp(k)
1358  enddo
1359  enddo
1360  endif
1361 
1362  endif ! if_ivflip
1363 
1364  enddo lab_do_ipt
1365 
1366  return
1367 !...................................
1368  end subroutine swrad
1369 !-----------------------------------
1371 
1372 
1376 !-----------------------------------
1377  subroutine rswinit &
1378  & ( me ) ! --- inputs:
1379 ! --- outputs: (none)
1380 
1381 ! =================== program usage description =================== !
1382 ! !
1383 ! purpose: initialize non-varying module variables, conversion factors,!
1384 ! and look-up tables. !
1385 ! !
1386 ! subprograms called: none !
1387 ! !
1388 ! ==================== defination of variables ==================== !
1389 ! !
1390 ! inputs: !
1391 ! me - print control for parallel process !
1392 ! !
1393 ! outputs: (none) !
1394 ! !
1395 ! external module variables: (in physparam) !
1396 ! iswrate - heating rate unit selections !
1397 ! =1: output in k/day !
1398 ! =2: output in k/second !
1399 ! iswrgas - control flag for rare gases (ch4,n2o,o2, etc.) !
1400 ! =0: do not include rare gases !
1401 ! >0: include all rare gases !
1402 ! iswcliq - liquid cloud optical properties contrl flag !
1403 ! =0: input cloud opt depth from diagnostic scheme !
1404 ! >0: input cwp,rew, and other cloud content parameters !
1405 ! isubcsw - sub-column cloud approximation control flag !
1406 ! =0: no sub-col cld treatment, use grid-mean cld quantities !
1407 ! =1: mcica sub-col, prescribed seeds to get random numbers !
1408 ! =2: mcica sub-col, providing array icseed for random numbers!
1409 ! icldflg - cloud scheme control flag !
1410 ! =0: diagnostic scheme gives cloud tau, omiga, and g. !
1411 ! =1: prognostic scheme gives cloud liq/ice path, etc. !
1412 ! iovrsw - clouds vertical overlapping control flag !
1413 ! =0: random overlapping clouds !
1414 ! =1: maximum/random overlapping clouds !
1415 ! =2: maximum overlap cloud !
1416 ! iswmode - control flag for 2-stream transfer scheme !
1417 ! =1; delta-eddington (joseph et al., 1976) !
1418 ! =2: pifm (zdunkowski et al., 1980) !
1419 ! =3: discrete ordinates (liou, 1973) !
1420 ! !
1421 ! ******************************************************************* !
1422 ! !
1423 ! definitions: !
1424 ! arrays for 10000-point look-up tables: !
1425 ! tau_tbl clear-sky optical depth !
1426 ! exp_tbl exponential lookup table for transmittance !
1427 ! !
1428 ! ******************************************************************* !
1429 ! !
1430 ! ====================== end of description block ================= !
1431 
1432 ! --- inputs:
1433  integer, intent(in) :: me
1434 
1435 ! --- outputs: none
1436 
1437 ! --- locals:
1438  real (kind=kind_phys), parameter :: expeps = 1.e-20
1439 
1440  integer :: i
1441 
1442  real (kind=kind_phys) :: tfn, tau
1443 
1444 !
1445 !===> ... begin here
1446 !
1447  if ( iovrsw<0 .or. iovrsw>2 ) then
1448  print *,' *** Error in specification of cloud overlap flag', &
1449  & ' IOVRSW=',iovrsw,' in RSWINIT !!'
1450  stop
1451  endif
1452 
1453  if (me == 0) then
1454  print *,' - Using AER Shortwave Radiation, Version: ',vtagsw
1455 
1456  if (iswmode == 1) then
1457  print *,' --- Delta-eddington 2-stream transfer scheme'
1458  else if (iswmode == 2) then
1459  print *,' --- PIFM 2-stream transfer scheme'
1460  else if (iswmode == 3) then
1461  print *,' --- Discrete ordinates 2-stream transfer scheme'
1462  endif
1463 
1464  if (iswrgas <= 0) then
1465  print *,' --- Rare gases absorption is NOT included in SW'
1466  else
1467  print *,' --- Include rare gases N2O, CH4, O2, absorptions',&
1468  & ' in SW'
1469  endif
1470 
1471  if ( isubcsw == 0 ) then
1472  print *,' --- Using standard grid average clouds, no ', &
1473  & 'sub-column clouds approximation applied'
1474  elseif ( isubcsw == 1 ) then
1475  print *,' --- Using MCICA sub-colum clouds approximation ', &
1476  & 'with a prescribed sequence of permutation seeds'
1477  elseif ( isubcsw == 2 ) then
1478  print *,' --- Using MCICA sub-colum clouds approximation ', &
1479  & 'with provided input array of permutation seeds'
1480  else
1481  print *,' *** Error in specification of sub-column cloud ', &
1482  & ' control flag isubcsw =',isubcsw,' !!'
1483  stop
1484  endif
1485  endif
1486 
1487 ! --- ... check cloud flags for consistency
1488 
1489  if ((icldflg == 0 .and. iswcliq /= 0) .or. &
1490  & (icldflg == 1 .and. iswcliq == 0)) then
1491  print *,' *** Model cloud scheme inconsistent with SW', &
1492  & ' radiation cloud radiative property setup !!'
1493  stop
1494  endif
1495 
1496 ! --- ... setup constant factors for heating rate
1497 ! the 1.0e-2 is to convert pressure from mb to N/m**2
1498 
1499  if (iswrate == 1) then
1500 ! heatfac = 8.4391
1501 ! heatfac = con_g * 86400. * 1.0e-2 / con_cp ! (in k/day)
1502  heatfac = con_g * 864.0 / con_cp ! (in k/day)
1503  else
1504  heatfac = con_g * 1.0e-2 / con_cp ! (in k/second)
1505  endif
1506 
1507 ! --- ... define exponential lookup tables for transmittance. tau is
1508 ! computed as a function of the tau transition function, and
1509 ! transmittance is calculated as a function of tau. all tables
1510 ! are computed at intervals of 0.0001. the inverse of the
1511 ! constant used in the Pade approximation to the tau transition
1512 ! function is set to bpade.
1513 
1514  exp_tbl(0) = 1.0
1515  exp_tbl(ntbmx) = expeps
1516 
1517  do i = 1, ntbmx-1
1518  tfn = float(i) / float(ntbmx-i)
1519  tau = bpade * tfn
1520  exp_tbl(i) = exp( -tau )
1521  enddo
1522 
1523  return
1524 !...................................
1525  end subroutine rswinit
1526 !-----------------------------------
1527 
1562 !-----------------------------------
1563  subroutine cldprop &
1564  & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, & ! --- inputs
1565  & cf1, nlay, ipseed, &
1566  & taucw, ssacw, asycw, cldfrc, cldfmc & ! --- output
1567  & )
1569 ! =================== program usage description =================== !
1570 ! !
1571 ! Purpose: Compute the cloud optical properties for each cloudy layer !
1572 ! and g-point interval. !
1573 ! !
1574 ! subprograms called: none !
1575 ! !
1576 ! ==================== defination of variables ==================== !
1577 ! !
1578 ! inputs: size !
1579 ! cfrac - real, layer cloud fraction nlay !
1580 ! ..... for iswcliq > 0 (prognostic cloud sckeme) - - - !
1581 ! cliqp - real, layer in-cloud liq water path (g/m**2) nlay !
1582 ! reliq - real, mean eff radius for liq cloud (micron) nlay !
1583 ! cicep - real, layer in-cloud ice water path (g/m**2) nlay !
1584 ! reice - real, mean eff radius for ice cloud (micron) nlay !
1585 ! cdat1 - real, layer rain drop water path (g/m**2) nlay !
1586 ! cdat2 - real, effective radius for rain drop (micron) nlay !
1587 ! cdat3 - real, layer snow flake water path(g/m**2) nlay !
1588 ! cdat4 - real, mean eff radius for snow flake(micron) nlay !
1589 ! ..... for iswcliq = 0 (diagnostic cloud sckeme) - - - !
1590 ! cdat1 - real, layer cloud optical depth nlay !
1591 ! cdat2 - real, layer cloud single scattering albedo nlay !
1592 ! cdat3 - real, layer cloud asymmetry factor nlay !
1593 ! cdat4 - real, optional use nlay !
1594 ! cliqp - real, not used nlay !
1595 ! cicep - real, not used nlay !
1596 ! reliq - real, not used nlay !
1597 ! reice - real, not used nlay !
1598 ! !
1599 ! cf1 - real, effective total cloud cover at surface 1 !
1600 ! nlay - integer, vertical layer number 1 !
1601 ! ipseed- permutation seed for generating random numbers (isubcsw>0) !
1602 ! !
1603 ! outputs: !
1604 ! taucw - real, cloud optical depth, w/o delta scaled nlay*nbdsw !
1605 ! ssacw - real, weighted cloud single scattering albedo nlay*nbdsw !
1606 ! (ssa = ssacw / taucw) !
1607 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
1608 ! (asy = asycw / ssacw) !
1609 ! cldfrc - real, cloud fraction of grid mean value nlay !
1610 ! cldfmc - real, cloud fraction for each sub-column nlay*ngptsw!
1611 ! !
1612 ! !
1613 ! explanation of the method for each value of iswcliq, and iswcice. !
1614 ! set up in module "physparam" !
1615 ! !
1616 ! iswcliq=0 : input cloud optical property (tau, ssa, asy). !
1617 ! (used for diagnostic cloud method) !
1618 ! iswcliq>0 : input cloud liq/ice path and effective radius, also !
1619 ! require the user of 'iswcice' to specify the method !
1620 ! used to compute aborption due to water/ice parts. !
1621 ! ................................................................... !
1622 ! !
1623 ! iswcliq=1 : liquid water cloud optical properties are computed !
1624 ! as in hu and stamnes (1993), j. clim., 6, 728-742. !
1625 ! !
1626 ! iswcice used only when iswcliq > 0 !
1627 ! the cloud ice path (g/m2) and ice effective radius !
1628 ! (microns) are inputs. !
1629 ! iswcice=1 : ice cloud optical properties are computed as in !
1630 ! ebert and curry (1992), jgr, 97, 3831-3836. !
1631 ! iswcice=2 : ice cloud optical properties are computed as in !
1632 ! streamer v3.0 (2001), key, streamer user's guide, !
1633 ! cooperative institude for meteorological studies,95pp!
1634 ! iswcice=3 : ice cloud optical properties are computed as in !
1635 ! fu (1996), j. clim., 9. !
1636 ! !
1637 ! other cloud control module variables: !
1638 ! isubcsw =0: standard cloud scheme, no sub-col cloud approximation !
1639 ! >0: mcica sub-col cloud scheme using ipseed as permutation!
1640 ! seed for generating rundom numbers !
1641 ! !
1642 ! ====================== end of description block ================= !
1643 !
1645 
1646 ! --- inputs:
1647  integer, intent(in) :: nlay, ipseed
1648  real (kind=kind_phys), intent(in) :: cf1
1649 
1650  real (kind=kind_phys), dimension(nlay), intent(in) :: cliqp, &
1651  & reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, cfrac
1652 
1653 ! --- outputs:
1654  real (kind=kind_phys), dimension(nlay,ngptsw), intent(out) :: &
1655  & cldfmc
1656  real (kind=kind_phys), dimension(nlay,nbdsw), intent(out) :: &
1657  & taucw, ssacw, asycw
1658  real (kind=kind_phys), dimension(nlay), intent(out) :: cldfrc
1659 
1660 ! --- locals:
1661  real (kind=kind_phys), dimension(nblow:nbhgh) :: tauliq, tauice, &
1662  & ssaliq, ssaice, ssaran, ssasnw, asyliq, asyice, &
1663  & asyran, asysnw
1664  real (kind=kind_phys), dimension(nlay) :: cldf
1665 
1666  real (kind=kind_phys) :: dgeice, factor, fint, tauran, tausnw, &
1667  & cldliq, refliq, cldice, refice, cldran, cldsnw, refsnw, &
1668  & extcoliq, ssacoliq, asycoliq, extcoice, ssacoice, asycoice,&
1669  & dgesnw
1670 
1671  logical :: lcloudy(nlay,ngptsw)
1672  integer :: ia, ib, ig, jb, k, index
1673 
1674 !
1675 !===> ... begin here
1676 !
1677  do ib = 1, nbdsw
1678  do k = 1, nlay
1679  taucw(k,ib) = f_zero
1680  ssacw(k,ib) = f_one
1681  asycw(k,ib) = f_zero
1682  enddo
1683  enddo
1684 
1686 
1687  lab_if_iswcliq : if (iswcliq > 0) then
1688 
1689  lab_do_k : do k = 1, nlay
1690  lab_if_cld : if (cfrac(k) > ftiny) then
1691 
1703 
1704  cldran = cdat1(k)
1705 ! refran = cdat2(k)
1706  cldsnw = cdat3(k)
1707  refsnw = cdat4(k)
1708  dgesnw = 1.0315 * refsnw ! for fu's snow formula
1709 
1710  tauran = cldran * a0r
1711 
1712 ! --- if use fu's formula it needs to be normalized by snow/ice density
1713 ! !not use snow density = 0.1 g/cm**3 = 0.1 g/(mu * m**2)
1714 ! use ice density = 0.9167 g/cm**3 = 0.9167 g/(mu * m**2)
1715 ! 1/0.9167 = 1.09087
1716 ! factor 1.5396=8/(3*sqrt(3)) converts reff to generalized ice particle size
1717 ! use newer factor value 1.0315
1718  if (cldsnw>f_zero .and. refsnw>10.0_kind_phys) then
1719 ! tausnw = cldsnw * (a0s + a1s/refsnw)
1720  tausnw = cldsnw*1.09087*(a0s + a1s/dgesnw) ! fu's formula
1721  else
1722  tausnw = f_zero
1723  endif
1724 
1725  do ib = nblow, nbhgh
1726  ssaran(ib) = tauran * (f_one - b0r(ib))
1727  ssasnw(ib) = tausnw * (f_one - (b0s(ib)+b1s(ib)*dgesnw))
1728  asyran(ib) = ssaran(ib) * c0r(ib)
1729  asysnw(ib) = ssasnw(ib) * c0s(ib)
1730  enddo
1731 
1732  cldliq = cliqp(k)
1733  cldice = cicep(k)
1734  refliq = reliq(k)
1735  refice = reice(k)
1736 
1737 ! --- ... calculation of absorption coefficients due to water clouds.
1738 
1739  if ( cldliq <= f_zero ) then
1740  do ib = nblow, nbhgh
1741  tauliq(ib) = f_zero
1742  ssaliq(ib) = f_zero
1743  asyliq(ib) = f_zero
1744  enddo
1745  else
1746  if ( iswcliq == 1 ) then
1747  factor = refliq - 1.5
1748  index = max( 1, min( 57, int( factor ) ))
1749  fint = factor - float(index)
1750 
1751  do ib = nblow, nbhgh
1752  extcoliq = max(f_zero, extliq1(index,ib) &
1753  & + fint*(extliq1(index+1,ib)-extliq1(index,ib)) )
1754  ssacoliq = max(f_zero, min(f_one, ssaliq1(index,ib) &
1755  & + fint*(ssaliq1(index+1,ib)-ssaliq1(index,ib)) ))
1756 
1757  asycoliq = max(f_zero, min(f_one, asyliq1(index,ib) &
1758  & + fint*(asyliq1(index+1,ib)-asyliq1(index,ib)) ))
1759 ! forcoliq = asycoliq * asycoliq
1760 
1761  tauliq(ib) = cldliq * extcoliq
1762  ssaliq(ib) = tauliq(ib) * ssacoliq
1763  asyliq(ib) = ssaliq(ib) * asycoliq
1764  enddo
1765  endif ! end if_iswcliq_block
1766  endif ! end if_cldliq_block
1767 
1768 ! --- ... calculation of absorption coefficients due to ice clouds.
1769 
1770  if ( cldice <= f_zero ) then
1771  do ib = nblow, nbhgh
1772  tauice(ib) = f_zero
1773  ssaice(ib) = f_zero
1774  asyice(ib) = f_zero
1775  enddo
1776  else
1777 
1778 ! --- ... ebert and curry approach for all particle sizes though somewhat
1779 ! unjustified for large ice particles
1780 
1781  if ( iswcice == 1 ) then
1782  refice = min(130.0_kind_phys,max(13.0_kind_phys,refice))
1783 
1784  do ib = nblow, nbhgh
1785  ia = idxebc(ib) ! eb_&_c band index for ice cloud coeff
1786 
1787  extcoice = max(f_zero, abari(ia)+bbari(ia)/refice )
1788  ssacoice = max(f_zero, min(f_one, &
1789  & f_one-cbari(ia)-dbari(ia)*refice ))
1790  asycoice = max(f_zero, min(f_one, &
1791  & ebari(ia)+fbari(ia)*refice ))
1792 ! forcoice = asycoice * asycoice
1793 
1794  tauice(ib) = cldice * extcoice
1795  ssaice(ib) = tauice(ib) * ssacoice
1796  asyice(ib) = ssaice(ib) * asycoice
1797  enddo
1798 
1799 ! --- ... streamer approach for ice effective radius between 5.0 and 131.0 microns
1800 
1801  elseif ( iswcice == 2 ) then
1802  refice = min(131.0_kind_phys,max(5.0_kind_phys,refice))
1803 
1804  factor = (refice - 2.0) / 3.0
1805  index = max( 1, min( 42, int( factor ) ))
1806  fint = factor - float(index)
1807 
1808  do ib = nblow, nbhgh
1809  extcoice = max(f_zero, extice2(index,ib) &
1810  & + fint*(extice2(index+1,ib)-extice2(index,ib)) )
1811  ssacoice = max(f_zero, min(f_one, ssaice2(index,ib) &
1812  & + fint*(ssaice2(index+1,ib)-ssaice2(index,ib)) ))
1813  asycoice = max(f_zero, min(f_one, asyice2(index,ib) &
1814  & + fint*(asyice2(index+1,ib)-asyice2(index,ib)) ))
1815 ! forcoice = asycoice * asycoice
1816 
1817  tauice(ib) = cldice * extcoice
1818  ssaice(ib) = tauice(ib) * ssacoice
1819  asyice(ib) = ssaice(ib) * asycoice
1820  enddo
1821 
1822 ! --- ... fu's approach for ice effective radius between 4.8 and 135 microns
1823 ! (generalized effective size from 5 to 140 microns)
1824 
1825  elseif ( iswcice == 3 ) then
1826  dgeice = max( 5.0, min( 140.0, 1.0315*refice ))
1827 
1828  factor = (dgeice - 2.0) / 3.0
1829  index = max( 1, min( 45, int( factor ) ))
1830  fint = factor - float(index)
1831 
1832  do ib = nblow, nbhgh
1833  extcoice = max(f_zero, extice3(index,ib) &
1834  & + fint*(extice3(index+1,ib)-extice3(index,ib)) )
1835  ssacoice = max(f_zero, min(f_one, ssaice3(index,ib) &
1836  & + fint*(ssaice3(index+1,ib)-ssaice3(index,ib)) ))
1837  asycoice = max(f_zero, min(f_one, asyice3(index,ib) &
1838  & + fint*(asyice3(index+1,ib)-asyice3(index,ib)) ))
1839 ! fdelta = max(f_zero, min(f_one, fdlice3(index,ib) &
1840 ! & + fint*(fdlice3(index+1,ib)-fdlice3(index,ib)) ))
1841 ! forcoice = min( asycoice, fdelta+0.5/ssacoice ) ! see fu 1996 p. 2067
1842 
1843  tauice(ib) = cldice * extcoice
1844  ssaice(ib) = tauice(ib) * ssacoice
1845  asyice(ib) = ssaice(ib) * asycoice
1846  enddo
1847 
1848  endif ! end if_iswcice_block
1849  endif ! end if_cldice_block
1850 
1851  do ib = 1, nbdsw
1852  jb = nblow + ib - 1
1853  taucw(k,ib) = tauliq(jb)+tauice(jb)+tauran+tausnw
1854  ssacw(k,ib) = ssaliq(jb)+ssaice(jb)+ssaran(jb)+ssasnw(jb)
1855  asycw(k,ib) = asyliq(jb)+asyice(jb)+asyran(jb)+asysnw(jb)
1856  enddo
1857 
1858  endif lab_if_cld
1859  enddo lab_do_k
1860 
1861  else lab_if_iswcliq
1862 
1863  do k = 1, nlay
1864  if (cfrac(k) > ftiny) then
1865  do ib = 1, nbdsw
1866  taucw(k,ib) = cdat1(k)
1867  ssacw(k,ib) = cdat1(k) * cdat2(k)
1868  asycw(k,ib) = ssacw(k,ib) * cdat3(k)
1869  enddo
1870  endif
1871  enddo
1872 
1873  endif lab_if_iswcliq
1874 
1877 
1878  if ( isubcsw > 0 ) then ! mcica sub-col clouds approx
1879 
1880  cldf(:) = cfrac(:)
1881  where (cldf(:) < ftiny)
1882  cldf(:) = f_zero
1883  end where
1884 
1885 ! --- ... call sub-column cloud generator
1886 
1887  call mcica_subcol &
1888 ! --- inputs:
1889  & ( cldf, nlay, ipseed, &
1890 ! --- outputs:
1891  & lcloudy &
1892  & )
1893 
1894  do ig = 1, ngptsw
1895  do k = 1, nlay
1896  if ( lcloudy(k,ig) ) then
1897  cldfmc(k,ig) = f_one
1898  else
1899  cldfmc(k,ig) = f_zero
1900  endif
1901  enddo
1902  enddo
1903 
1904  else ! non-mcica, normalize cloud
1905 
1906  do k = 1, nlay
1907  cldfrc(k) = cfrac(k) / cf1
1908  enddo
1909  endif ! end if_isubcsw_block
1910 
1911  return
1912 !...................................
1913  end subroutine cldprop
1914 !-----------------------------------
1916 
1917 
1923 ! ----------------------------------
1924  subroutine mcica_subcol &
1925  & ( cldf, nlay, ipseed, & ! --- inputs
1926  & lcloudy & ! --- outputs
1927  & )
1929 ! ==================== defination of variables ==================== !
1930 ! !
1931 ! input variables: size !
1932 ! cldf - real, layer cloud fraction nlay !
1933 ! nlay - integer, number of model vertical layers 1 !
1934 ! ipseed - integer, permute seed for random num generator 1 !
1935 ! ** note : if the cloud generator is called multiple times, need !
1936 ! to permute the seed between each call; if between calls !
1937 ! for lw and sw, use values differ by the number of g-pts. !
1938 ! !
1939 ! output variables: !
1940 ! lcloudy - logical, sub-colum cloud profile flag array nlay*ngptsw!
1941 ! !
1942 ! other control flags from module variables: !
1943 ! iovrsw : control flag for cloud overlapping method !
1944 ! =0:random; =1:maximum/random; =2:maximum !
1945 ! !
1946 ! !
1947 ! ===================== end of definitions ==================== !
1948 
1949  implicit none
1950 
1951 ! --- inputs:
1952  integer, intent(in) :: nlay, ipseed
1953 
1954  real (kind=kind_phys), dimension(nlay), intent(in) :: cldf
1955 
1956 ! --- outputs:
1957  logical, dimension(nlay,ngptsw), intent(out):: lcloudy
1958 
1959 ! --- locals:
1960  real (kind=kind_phys) :: cdfunc(nlay,ngptsw), tem1, &
1961  & rand2d(nlay*ngptsw), rand1d(ngptsw)
1962 
1963  type(random_stat) :: stat ! for thread safe random generator
1964 
1965  integer :: k, n, k1
1966 !
1967 !===> ... begin here
1968 !
1969 ! --- ... advance randum number generator by ipseed values
1970 
1971  call random_setseed &
1972 ! --- inputs:
1973  & ( ipseed, &
1974 ! --- outputs:
1975  & stat &
1976  & )
1977 
1978 ! --- ... sub-column set up according to overlapping assumption
1979 
1980  select case ( iovrsw )
1981 
1982  case( 0 ) ! random overlap, pick a random value at every level
1983 
1984  call random_number &
1985 ! --- inputs: ( none )
1986 ! --- outputs:
1987  & ( rand2d, stat )
1988 
1989  k1 = 0
1990  do n = 1, ngptsw
1991  do k = 1, nlay
1992  k1 = k1 + 1
1993  cdfunc(k,n) = rand2d(k1)
1994  enddo
1995  enddo
1996 
1997  case( 1 ) ! max-ran overlap
1998 
1999  call random_number &
2000 ! --- inputs: ( none )
2001 ! --- outputs:
2002  & ( rand2d, stat )
2003 
2004  k1 = 0
2005  do n = 1, ngptsw
2006  do k = 1, nlay
2007  k1 = k1 + 1
2008  cdfunc(k,n) = rand2d(k1)
2009  enddo
2010  enddo
2011 
2012 ! --- first pick a random number for bottom/top layer.
2013 ! then walk up the column: (aer's code)
2014 ! if layer below is cloudy, use the same rand num in the layer below
2015 ! if layer below is clear, use a new random number
2016 
2017 ! --- from bottom up
2018  do k = 2, nlay
2019  k1 = k - 1
2020  tem1 = f_one - cldf(k1)
2021 
2022  do n = 1, ngptsw
2023  if ( cdfunc(k1,n) > tem1 ) then
2024  cdfunc(k,n) = cdfunc(k1,n)
2025  else
2026  cdfunc(k,n) = cdfunc(k,n) * tem1
2027  endif
2028  enddo
2029  enddo
2030 
2031 ! --- then walk down the column: (if use original author's method)
2032 ! if layer above is cloudy, use the same rand num in the layer above
2033 ! if layer above is clear, use a new random number
2034 
2035 ! --- from top down
2036 ! do k = nlay-1, 1, -1
2037 ! k1 = k + 1
2038 ! tem1 = f_one - cldf(k1)
2039 
2040 ! do n = 1, ngptsw
2041 ! if ( cdfunc(k1,n) > tem1 ) then
2042 ! cdfunc(k,n) = cdfunc(k1,n)
2043 ! else
2044 ! cdfunc(k,n) = cdfunc(k,n) * tem1
2045 ! endif
2046 ! enddo
2047 ! enddo
2048 
2049  case( 2 ) ! maximum overlap, pick same random numebr at every level
2050 
2051  call random_number &
2052 ! --- inputs: ( none )
2053 ! --- outputs:
2054  & ( rand1d, stat )
2055 
2056  do n = 1, ngptsw
2057  tem1 = rand1d(n)
2058 
2059  do k = 1, nlay
2060  cdfunc(k,n) = tem1
2061  enddo
2062  enddo
2063 
2064  end select
2065 
2066 ! --- ... generate subcolumns for homogeneous clouds
2067 
2068  do k = 1, nlay
2069  tem1 = f_one - cldf(k)
2070 
2071  do n = 1, ngptsw
2072  lcloudy(k,n) = cdfunc(k,n) >= tem1
2073  enddo
2074  enddo
2075 
2076  return
2077 ! ..................................
2078  end subroutine mcica_subcol
2079 ! ----------------------------------
2080 
2105 ! ----------------------------------
2106  subroutine setcoef &
2107  & ( pavel,tavel,h2ovmr, nlay,nlp1, & ! --- inputs
2108  & laytrop,jp,jt,jt1,fac00,fac01,fac10,fac11, & ! --- outputs
2109  & selffac,selffrac,indself,forfac,forfrac,indfor &
2110  & )
2112 ! =================== program usage description =================== !
2113 ! !
2114 ! purpose: compute various coefficients needed in radiative transfer !
2115 ! calculations. !
2116 ! !
2117 ! subprograms called: none !
2118 ! !
2119 ! ==================== defination of variables ==================== !
2120 ! !
2121 ! inputs: -size- !
2122 ! pavel - real, layer pressures (mb) nlay !
2123 ! tavel - real, layer temperatures (k) nlay !
2124 ! h2ovmr - real, layer w.v. volum mixing ratio (kg/kg) nlay !
2125 ! nlay/nlp1 - integer, total number of vertical layers, levels 1 !
2126 ! !
2127 ! outputs: !
2128 ! laytrop - integer, tropopause layer index (unitless) 1 !
2129 ! jp - real, indices of lower reference pressure nlay !
2130 ! jt, jt1 - real, indices of lower reference temperatures nlay !
2131 ! at levels of jp and jp+1 !
2132 ! facij - real, factors multiply the reference ks, nlay !
2133 ! i,j=0/1 for lower/higher of the 2 appropriate !
2134 ! temperatures and altitudes. !
2135 ! selffac - real, scale factor for w. v. self-continuum nlay !
2136 ! equals (w. v. density)/(atmospheric density !
2137 ! at 296k and 1013 mb) !
2138 ! selffrac - real, factor for temperature interpolation of nlay !
2139 ! reference w. v. self-continuum data !
2140 ! indself - integer, index of lower ref temp for selffac nlay !
2141 ! forfac - real, scale factor for w. v. foreign-continuum nlay !
2142 ! forfrac - real, factor for temperature interpolation of nlay !
2143 ! reference w.v. foreign-continuum data !
2144 ! indfor - integer, index of lower ref temp for forfac nlay !
2145 ! !
2146 ! ====================== end of definitions =================== !
2147 
2148 ! --- inputs:
2149  integer, intent(in) :: nlay, nlp1
2150 
2151  real (kind=kind_phys), dimension(:), intent(in) :: pavel, tavel, &
2152  & h2ovmr
2153 
2154 ! --- outputs:
2155  integer, dimension(nlay), intent(out) :: indself, indfor, &
2156  & jp, jt, jt1
2157  integer, intent(out) :: laytrop
2158 
2159  real (kind=kind_phys), dimension(nlay), intent(out) :: fac00, &
2160  & fac01, fac10, fac11, selffac, selffrac, forfac, forfrac
2161 
2162 ! --- locals:
2163  real (kind=kind_phys) :: plog, fp, fp1, ft, ft1, tem1, tem2
2164 
2165  integer :: i, k, jp1
2166 !
2167 !===> ... begin here
2168 !
2169  laytrop= nlay
2170 
2171  do k = 1, nlay
2172 
2173  forfac(k) = pavel(k)*stpfac / (tavel(k)*(f_one + h2ovmr(k)))
2174 
2175 ! --- ... find the two reference pressures on either side of the
2176 ! layer pressure. store them in jp and jp1. store in fp the
2177 ! fraction of the difference (in ln(pressure)) between these
2178 ! two values that the layer pressure lies.
2179 
2180  plog = log(pavel(k))
2181  jp(k) = max(1, min(58, int(36.0 - 5.0*(plog+0.04)) ))
2182  jp1 = jp(k) + 1
2183  fp = 5.0 * (preflog(jp(k)) - plog)
2184 
2185 ! --- ... determine, for each reference pressure (jp and jp1), which
2186 ! reference temperature (these are different for each reference
2187 ! pressure) is nearest the layer temperature but does not exceed it.
2188 ! store these indices in jt and jt1, resp. store in ft (resp. ft1)
2189 ! the fraction of the way between jt (jt1) and the next highest
2190 ! reference temperature that the layer temperature falls.
2191 
2192  tem1 = (tavel(k) - tref(jp(k))) / 15.0
2193  tem2 = (tavel(k) - tref(jp1 )) / 15.0
2194  jt(k) = max(1, min(4, int(3.0 + tem1) ))
2195  jt1(k) = max(1, min(4, int(3.0 + tem2) ))
2196  ft = tem1 - float(jt(k) - 3)
2197  ft1 = tem2 - float(jt1(k) - 3)
2198 
2199 ! --- ... we have now isolated the layer ln pressure and temperature,
2200 ! between two reference pressures and two reference temperatures
2201 ! (for each reference pressure). we multiply the pressure
2202 ! fraction fp with the appropriate temperature fractions to get
2203 ! the factors that will be needed for the interpolation that yields
2204 ! the optical depths (performed in routines taugbn for band n).
2205 
2206  fp1 = f_one - fp
2207  fac10(k) = fp1 * ft
2208  fac00(k) = fp1 * (f_one - ft)
2209  fac11(k) = fp * ft1
2210  fac01(k) = fp * (f_one - ft1)
2211 
2212 ! --- ... if the pressure is less than ~100mb, perform a different
2213 ! set of species interpolations.
2214 
2215  if ( plog > 4.56 ) then
2216 
2217  laytrop = k
2218 
2219 ! --- ... set up factors needed to separately include the water vapor
2220 ! foreign-continuum in the calculation of absorption coefficient.
2221 
2222  tem1 = (332.0 - tavel(k)) / 36.0
2223  indfor(k) = min(2, max(1, int(tem1)))
2224  forfrac(k) = tem1 - float(indfor(k))
2225 
2226 ! --- ... set up factors needed to separately include the water vapor
2227 ! self-continuum in the calculation of absorption coefficient.
2228 
2229  tem2 = (tavel(k) - 188.0) / 7.2
2230  indself(k) = min(9, max(1, int(tem2)-7))
2231  selffrac(k) = tem2 - float(indself(k) + 7)
2232  selffac(k) = h2ovmr(k) * forfac(k)
2233 
2234  else
2235 
2236 ! --- ... set up factors needed to separately include the water vapor
2237 ! foreign-continuum in the calculation of absorption coefficient.
2238 
2239  tem1 = (tavel(k) - 188.0) / 36.0
2240  indfor(k) = 3
2241  forfrac(k) = tem1 - f_one
2242 
2243  indself(k) = 0
2244  selffrac(k) = f_zero
2245  selffac(k) = f_zero
2246 
2247  endif
2248 
2249  enddo ! end_do_k_loop
2250 
2251  return
2252 ! ..................................
2253  end subroutine setcoef
2254 ! ----------------------------------
2255 
2295 !-----------------------------------
2296  subroutine spcvrtc &
2297  & ( ssolar,cosz,sntz,albbm,albdf,sfluxzen,cldfrc, & ! --- inputs
2298  & cf1,cf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
2299  & nlay, nlp1, &
2300  & fxupc,fxdnc,fxup0,fxdn0, & ! --- outputs
2301  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
2302  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
2303  & )
2305 ! =================== program usage description =================== !
2306 ! !
2307 ! purpose: computes the shortwave radiative fluxes using two-stream !
2308 ! method !
2309 ! !
2310 ! subprograms called: vrtqdr !
2311 ! !
2312 ! ==================== defination of variables ==================== !
2313 ! !
2314 ! inputs: size !
2315 ! ssolar - real, incoming solar flux at top 1 !
2316 ! cosz - real, cosine solar zenith angle 1 !
2317 ! sntz - real, secant solar zenith angle 1 !
2318 ! albbm - real, surface albedo for direct beam radiation 2 !
2319 ! albdf - real, surface albedo for diffused radiation 2 !
2320 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
2321 ! cldfrc - real, layer cloud fraction nlay !
2322 ! cf1 - real, >0: cloudy sky, otherwise: clear sky 1 !
2323 ! cf0 - real, =1-cf1 1 !
2324 ! taug - real, spectral optical depth for gases nlay*ngptsw!
2325 ! taur - real, optical depth for rayleigh scattering nlay*ngptsw!
2326 ! tauae - real, aerosols optical depth nlay*nbdsw !
2327 ! ssaae - real, aerosols single scattering albedo nlay*nbdsw !
2328 ! asyae - real, aerosols asymmetry factor nlay*nbdsw !
2329 ! taucw - real, weighted cloud optical depth nlay*nbdsw !
2330 ! ssacw - real, weighted cloud single scat albedo nlay*nbdsw !
2331 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
2332 ! nlay,nlp1 - integer, number of layers/levels 1 !
2333 ! !
2334 ! output variables: !
2335 ! fxupc - real, tot sky upward flux nlp1*nbdsw !
2336 ! fxdnc - real, tot sky downward flux nlp1*nbdsw !
2337 ! fxup0 - real, clr sky upward flux nlp1*nbdsw !
2338 ! fxdn0 - real, clr sky downward flux nlp1*nbdsw !
2339 ! ftoauc - real, tot sky toa upwd flux 1 !
2340 ! ftoau0 - real, clr sky toa upwd flux 1 !
2341 ! ftoadc - real, toa downward (incoming) solar flux 1 !
2342 ! fsfcuc - real, tot sky sfc upwd flux 1 !
2343 ! fsfcu0 - real, clr sky sfc upwd flux 1 !
2344 ! fsfcdc - real, tot sky sfc dnwd flux 1 !
2345 ! fsfcd0 - real, clr sky sfc dnwd flux 1 !
2346 ! sfbmc - real, tot sky sfc dnwd beam flux (nir/uv+vis) 2 !
2347 ! sfdfc - real, tot sky sfc dnwd diff flux (nir/uv+vis) 2 !
2348 ! sfbm0 - real, clr sky sfc dnwd beam flux (nir/uv+vis) 2 !
2349 ! sfdf0 - real, clr sky sfc dnwd diff flux (nir/uv+vis) 2 !
2350 ! suvbfc - real, tot sky sfc dnwd uv-b flux 1 !
2351 ! suvbf0 - real, clr sky sfc dnwd uv-b flux 1 !
2352 ! !
2353 ! internal variables: !
2354 ! zrefb - real, direct beam reflectivity for clear/cloudy nlp1 !
2355 ! zrefd - real, diffuse reflectivity for clear/cloudy nlp1 !
2356 ! ztrab - real, direct beam transmissivity for clear/cloudy nlp1 !
2357 ! ztrad - real, diffuse transmissivity for clear/cloudy nlp1 !
2358 ! zldbt - real, layer beam transmittance for clear/cloudy nlp1 !
2359 ! ztdbt - real, lev total beam transmittance for clr/cld nlp1 !
2360 ! !
2361 ! control parameters in module "physparam" !
2362 ! iswmode - control flag for 2-stream transfer schemes !
2363 ! = 1 delta-eddington (joseph et al., 1976) !
2364 ! = 2 pifm (zdunkowski et al., 1980) !
2365 ! = 3 discrete ordinates (liou, 1973) !
2366 ! !
2367 ! ******************************************************************* !
2368 ! original code description !
2369 ! !
2370 ! method: !
2371 ! ------- !
2372 ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations. !
2373 ! kmodts = 1 eddington (joseph et al., 1976) !
2374 ! = 2 pifm (zdunkowski et al., 1980) !
2375 ! = 3 discrete ordinates (liou, 1973) !
2376 ! !
2377 ! modifications: !
2378 ! -------------- !
2379 ! original: h. barker !
2380 ! revision: merge with rrtmg_sw: j.-j.morcrette, ecmwf, feb 2003 !
2381 ! revision: add adjustment for earth/sun distance:mjiacono,aer,oct2003!
2382 ! revision: bug fix for use of palbp and palbd: mjiacono, aer, nov2003!
2383 ! revision: bug fix to apply delta scaling to clear sky: aer, dec2004 !
2384 ! revision: code modified so that delta scaling is not done in cloudy !
2385 ! profiles if routine cldprop is used; delta scaling can be !
2386 ! applied by swithcing code below if cldprop is not used to !
2387 ! get cloud properties. aer, jan 2005 !
2388 ! revision: uniform formatting for rrtmg: mjiacono, aer, jul 2006 !
2389 ! revision: use exponential lookup table for transmittance: mjiacono, !
2390 ! aer, aug 2007 !
2391 ! !
2392 ! ******************************************************************* !
2393 ! ====================== end of description block ================= !
2394 
2395 ! --- constant parameters:
2396  real (kind=kind_phys), parameter :: zcrit = 0.9999995 ! thresold for conservative scattering
2397  real (kind=kind_phys), parameter :: zsr3 = sqrt(3.0)
2398  real (kind=kind_phys), parameter :: od_lo = 0.06
2399  real (kind=kind_phys), parameter :: eps1 = 1.0e-8
2400 
2401 ! --- inputs:
2402  integer, intent(in) :: nlay, nlp1
2403 
2404  real (kind=kind_phys), dimension(nlay,ngptsw), intent(in) :: &
2405  & taug, taur
2406  real (kind=kind_phys), dimension(nlay,nbdsw), intent(in) :: &
2407  & taucw, ssacw, asycw, tauae, ssaae, asyae
2408 
2409  real (kind=kind_phys), dimension(ngptsw), intent(in) :: sfluxzen
2410  real (kind=kind_phys), dimension(nlay), intent(in) :: cldfrc
2411 
2412  real (kind=kind_phys), dimension(2), intent(in) :: albbm, albdf
2413 
2414  real (kind=kind_phys), intent(in) :: cosz, sntz, cf1, cf0, ssolar
2415 
2416 ! --- outputs:
2417  real (kind=kind_phys), dimension(nlp1,nbdsw), intent(out) :: &
2418  & fxupc, fxdnc, fxup0, fxdn0
2419 
2420  real (kind=kind_phys), dimension(2), intent(out) :: sfbmc, sfdfc, &
2421  & sfbm0, sfdf0
2422 
2423  real (kind=kind_phys), intent(out) :: suvbfc, suvbf0, ftoadc, &
2424  & ftoauc, ftoau0, fsfcuc, fsfcu0, fsfcdc, fsfcd0
2425 
2426 ! --- locals:
2427  real (kind=kind_phys), dimension(nlay) :: ztaus, zssas, zasys, &
2428  & zldbt0
2429 
2430  real (kind=kind_phys), dimension(nlp1) :: zrefb, zrefd, ztrab, &
2431  & ztrad, ztdbt, zldbt, zfu, zfd
2432 
2433  real (kind=kind_phys) :: ztau1, zssa1, zasy1, ztau0, zssa0, &
2434  & zasy0, zasy3, zssaw, zasyw, zgam1, zgam2, zgam3, zgam4, &
2435  & zc0, zc1, za1, za2, zb1, zb2, zrk, zrk2, zrp, zrp1, zrm1, &
2436  & zrpp, zrkg1, zrkg3, zrkg4, zexp1, zexm1, zexp2, zexm2, &
2437  & zexp3, zexp4, zden1, ze1r45, ftind, zsolar, zrefb1, &
2438  & zrefd1, ztrab1, ztrad1, ztdbt0, zr1, zr2, zr3, zr4, zr5, &
2439  & zt1, zt2, zt3, zf1, zf2, zrpp1
2440 
2441  integer :: ib, ibd, jb, jg, k, kp, itind
2442 !
2443 !===> ... begin here
2444 
2446  do ib = 1, nbdsw
2447  do k = 1, nlp1
2448  fxdnc(k,ib) = f_zero
2449  fxupc(k,ib) = f_zero
2450  fxdn0(k,ib) = f_zero
2451  fxup0(k,ib) = f_zero
2452  enddo
2453  enddo
2454 
2455  ftoadc = f_zero
2456  ftoauc = f_zero
2457  ftoau0 = f_zero
2458  fsfcuc = f_zero
2459  fsfcu0 = f_zero
2460  fsfcdc = f_zero
2461  fsfcd0 = f_zero
2462 
2463 !! --- ... uv-b surface downward fluxes
2464  suvbfc = f_zero
2465  suvbf0 = f_zero
2466 
2467 !! --- ... output surface flux components
2468  sfbmc(1) = f_zero
2469  sfbmc(2) = f_zero
2470  sfdfc(1) = f_zero
2471  sfdfc(2) = f_zero
2472  sfbm0(1) = f_zero
2473  sfbm0(2) = f_zero
2474  sfdf0(1) = f_zero
2475  sfdf0(2) = f_zero
2476 
2477 ! --- ... loop over all g-points in each band
2478 
2479  lab_do_jg : do jg = 1, ngptsw
2480 
2481  jb = ngb(jg)
2482  ib = jb + 1 - nblow
2483  ibd = idxsfc(jb)
2484 
2485  zsolar = ssolar * sfluxzen(jg)
2486 
2487 ! --- ... set up toa direct beam and surface values (beam and diff)
2488 
2489  ztdbt(nlp1) = f_one
2490  ztdbt0 = f_one
2491 
2492  zldbt(1) = f_zero
2493  if (ibd /= 0) then
2494  zrefb(1) = albbm(ibd)
2495  zrefd(1) = albdf(ibd)
2496  else
2497  zrefb(1) = 0.5 * (albbm(1) + albbm(2))
2498  zrefd(1) = 0.5 * (albdf(1) + albdf(2))
2499  endif
2500  ztrab(1) = f_zero
2501  ztrad(1) = f_zero
2502 
2513 
2514  do k = nlay, 1, -1
2515  kp = k + 1
2516 
2517  ztau0 = max( ftiny, taur(k,jg)+taug(k,jg)+tauae(k,ib) )
2518  zssa0 = taur(k,jg) + tauae(k,ib)*ssaae(k,ib)
2519  zasy0 = asyae(k,ib)*ssaae(k,ib)*tauae(k,ib)
2520  zssaw = min( oneminus, zssa0 / ztau0 )
2521  zasyw = zasy0 / max( ftiny, zssa0 )
2522 
2523 ! --- ... saving clear-sky quantities for later total-sky usage
2524  ztaus(k) = ztau0
2525  zssas(k) = zssa0
2526  zasys(k) = zasy0
2527 
2528 ! --- ... delta scaling for clear-sky condition
2529  za1 = zasyw * zasyw
2530  za2 = zssaw * za1
2531 
2532  ztau1 = (f_one - za2) * ztau0
2533  zssa1 = (zssaw - za2) / (f_one - za2)
2534 !org zasy1 = (zasyw - za1) / (f_one - za1) ! this line is replaced by the next
2535  zasy1 = zasyw / (f_one + zasyw) ! to reduce truncation error
2536  zasy3 = 0.75 * zasy1
2537 
2538 ! --- ... general two-stream expressions
2539  if ( iswmode == 1 ) then
2540  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
2541  zgam2 =-0.25 + zssa1 * (f_one - zasy3)
2542  zgam3 = 0.5 - zasy3 * cosz
2543  elseif ( iswmode == 2 ) then ! pifm
2544  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
2545  zgam2 = 0.75* zssa1 * (f_one- zasy1)
2546  zgam3 = 0.5 - zasy3 * cosz
2547  elseif ( iswmode == 3 ) then ! discrete ordinates
2548  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
2549  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
2550  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
2551  endif
2552  zgam4 = f_one - zgam3
2553 
2554 ! --- ... compute homogeneous reflectance and transmittance
2555 
2556  if ( zssaw >= zcrit ) then ! for conservative scattering
2557  za1 = zgam1 * cosz - zgam3
2558  za2 = zgam1 * ztau1
2559 
2560 ! --- ... use exponential lookup table for transmittance, or expansion
2561 ! of exponential for low optical depth
2562 
2563  zb1 = min( ztau1*sntz , 500.0 )
2564  if ( zb1 <= od_lo ) then
2565  zb2 = f_one - zb1 + 0.5*zb1*zb1
2566  else
2567  ftind = zb1 / (bpade + zb1)
2568  itind = ftind*ntbmx + 0.5
2569  zb2 = exp_tbl(itind)
2570  endif
2571 
2572 ! ... collimated beam
2573  zrefb(kp) = max(f_zero, min(f_one, &
2574  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
2575  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp) ))
2576 
2577 ! ... isotropic incidence
2578  zrefd(kp) = max(f_zero, min(f_one, za2/(f_one + za2) ))
2579  ztrad(kp) = max(f_zero, min(f_one, f_one-zrefd(kp) ))
2580 
2581  else ! for non-conservative scattering
2582  za1 = zgam1*zgam4 + zgam2*zgam3
2583  za2 = zgam1*zgam3 + zgam2*zgam4
2584  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
2585  zrk2= 2.0 * zrk
2586 
2587  zrp = zrk * cosz
2588  zrp1 = f_one + zrp
2589  zrm1 = f_one - zrp
2590  zrpp1= f_one - zrp*zrp
2591  zrpp = sign( max(flimit, abs(zrpp1)), zrpp1 ) ! avoid numerical singularity
2592  zrkg1= zrk + zgam1
2593  zrkg3= zrk * zgam3
2594  zrkg4= zrk * zgam4
2595 
2596  zr1 = zrm1 * (za2 + zrkg3)
2597  zr2 = zrp1 * (za2 - zrkg3)
2598  zr3 = zrk2 * (zgam3 - za2*cosz)
2599  zr4 = zrpp * zrkg1
2600  zr5 = zrpp * (zrk - zgam1)
2601 
2602  zt1 = zrp1 * (za1 + zrkg4)
2603  zt2 = zrm1 * (za1 - zrkg4)
2604  zt3 = zrk2 * (zgam4 + za1*cosz)
2605 
2606 ! --- ... use exponential lookup table for transmittance, or expansion
2607 ! of exponential for low optical depth
2608 
2609  zb1 = min( zrk*ztau1, 500.0 )
2610  if ( zb1 <= od_lo ) then
2611  zexm1 = f_one - zb1 + 0.5*zb1*zb1
2612  else
2613  ftind = zb1 / (bpade + zb1)
2614  itind = ftind*ntbmx + 0.5
2615  zexm1 = exp_tbl(itind)
2616  endif
2617  zexp1 = f_one / zexm1
2618 
2619  zb2 = min( sntz*ztau1, 500.0 )
2620  if ( zb2 <= od_lo ) then
2621  zexm2 = f_one - zb2 + 0.5*zb2*zb2
2622  else
2623  ftind = zb2 / (bpade + zb2)
2624  itind = ftind*ntbmx + 0.5
2625  zexm2 = exp_tbl(itind)
2626  endif
2627  zexp2 = f_one / zexm2
2628  ze1r45 = zr4*zexp1 + zr5*zexm1
2629 
2630 ! ... collimated beam
2631  if (ze1r45>=-eps1 .and. ze1r45<=eps1) then
2632  zrefb(kp) = eps1
2633  ztrab(kp) = zexm2
2634  else
2635  zden1 = zssa1 / ze1r45
2636  zrefb(kp) = max(f_zero, min(f_one, &
2637  & (zr1*zexp1 - zr2*zexm1 - zr3*zexm2)*zden1 ))
2638  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one &
2639  & - (zt1*zexp1 - zt2*zexm1 - zt3*zexp2)*zden1) ))
2640  endif
2641 
2642 ! ... diffuse beam
2643  zden1 = zr4 / (ze1r45 * zrkg1)
2644  zrefd(kp) = max(f_zero, min(f_one, &
2645  & zgam2*(zexp1 - zexm1)*zden1 ))
2646  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
2647  endif ! end if_zssaw_block
2648 
2649 ! --- ... direct beam transmittance. use exponential lookup table
2650 ! for transmittance, or expansion of exponential for low
2651 ! optical depth
2652 
2653  zr1 = ztau1 * sntz
2654  if ( zr1 <= od_lo ) then
2655  zexp3 = f_one - zr1 + 0.5*zr1*zr1
2656  else
2657  ftind = zr1 / (bpade + zr1)
2658  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2659  zexp3 = exp_tbl(itind)
2660  endif
2661 
2662  ztdbt(k) = zexp3 * ztdbt(kp)
2663  zldbt(kp) = zexp3
2664 
2665 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2666 ! (must use 'orig', unscaled cloud optical depth)
2667 
2668  zr1 = ztau0 * sntz
2669  if ( zr1 <= od_lo ) then
2670  zexp4 = f_one - zr1 + 0.5*zr1*zr1
2671  else
2672  ftind = zr1 / (bpade + zr1)
2673  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2674  zexp4 = exp_tbl(itind)
2675  endif
2676 
2677  zldbt0(k) = zexp4
2678  ztdbt0 = zexp4 * ztdbt0
2679  enddo ! end do_k_loop
2680 
2681  call vrtqdr &
2682 ! --- inputs:
2683  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
2684  & nlay, nlp1, &
2685 ! --- outputs:
2686  & zfu, zfd &
2687  & )
2688 
2689 ! --- ... compute upward and downward fluxes at levels
2690  do k = 1, nlp1
2691  fxup0(k,ib) = fxup0(k,ib) + zsolar*zfu(k)
2692  fxdn0(k,ib) = fxdn0(k,ib) + zsolar*zfd(k)
2693  enddo
2694 
2695 !! --- ... surface downward beam/diffused flux components
2696  zb1 = zsolar*ztdbt0
2697  zb2 = zsolar*(zfd(1) - ztdbt0)
2698 
2699  if (ibd /= 0) then
2700  sfbm0(ibd) = sfbm0(ibd) + zb1
2701  sfdf0(ibd) = sfdf0(ibd) + zb2
2702  else
2703  zf1 = 0.5 * zb1
2704  zf2 = 0.5 * zb2
2705  sfbm0(1) = sfbm0(1) + zf1
2706  sfdf0(1) = sfdf0(1) + zf2
2707  sfbm0(2) = sfbm0(2) + zf1
2708  sfdf0(2) = sfdf0(2) + zf2
2709  endif
2710 ! sfbm0(ibd) = sfbm0(ibd) + zsolar*ztdbt0
2711 ! sfdf0(ibd) = sfdf0(ibd) + zsolar*(zfd(1) - ztdbt0)
2712 
2722 
2723  if ( cf1 > eps ) then
2724 
2725 ! --- ... set up toa direct beam and surface values (beam and diff)
2726  ztdbt0 = f_one
2727  zldbt(1) = f_zero
2728 
2729  do k = nlay, 1, -1
2730  kp = k + 1
2731  zc0 = f_one - cldfrc(k)
2732  zc1 = cldfrc(k)
2733  if ( zc1 > ftiny ) then ! it is a cloudy-layer
2734 
2735  ztau0 = ztaus(k) + taucw(k,ib)
2736  zssa0 = zssas(k) + ssacw(k,ib)
2737  zasy0 = zasys(k) + asycw(k,ib)
2738  zssaw = min(oneminus, zssa0 / ztau0)
2739  zasyw = zasy0 / max(ftiny, zssa0)
2740 
2741 ! --- ... delta scaling for total-sky condition
2742  za1 = zasyw * zasyw
2743  za2 = zssaw * za1
2744 
2745  ztau1 = (f_one - za2) * ztau0
2746  zssa1 = (zssaw - za2) / (f_one - za2)
2747 !org zasy1 = (zasyw - za1) / (f_one - za1)
2748  zasy1 = zasyw / (f_one + zasyw)
2749  zasy3 = 0.75 * zasy1
2750 
2751 ! --- ... general two-stream expressions
2752  if ( iswmode == 1 ) then
2753  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
2754  zgam2 =-0.25 + zssa1 * (f_one - zasy3)
2755  zgam3 = 0.5 - zasy3 * cosz
2756  elseif ( iswmode == 2 ) then ! pifm
2757  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
2758  zgam2 = 0.75* zssa1 * (f_one- zasy1)
2759  zgam3 = 0.5 - zasy3 * cosz
2760  elseif ( iswmode == 3 ) then ! discrete ordinates
2761  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
2762  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
2763  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
2764  endif
2765  zgam4 = f_one - zgam3
2766 
2767  zrefb1 = zrefb(kp)
2768  zrefd1 = zrefd(kp)
2769  ztrab1 = ztrab(kp)
2770  ztrad1 = ztrad(kp)
2771 
2772 ! --- ... compute homogeneous reflectance and transmittance
2773 
2774  if ( zssaw >= zcrit ) then ! for conservative scattering
2775  za1 = zgam1 * cosz - zgam3
2776  za2 = zgam1 * ztau1
2777 
2778 ! --- ... use exponential lookup table for transmittance, or expansion
2779 ! of exponential for low optical depth
2780 
2781  zb1 = min( ztau1*sntz , 500.0 )
2782  if ( zb1 <= od_lo ) then
2783  zb2 = f_one - zb1 + 0.5*zb1*zb1
2784  else
2785  ftind = zb1 / (bpade + zb1)
2786  itind = ftind*ntbmx + 0.5
2787  zb2 = exp_tbl(itind)
2788  endif
2789 
2790 ! ... collimated beam
2791  zrefb(kp) = max(f_zero, min(f_one, &
2792  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
2793  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp)))
2794 
2795 ! ... isotropic incidence
2796  zrefd(kp) = max(f_zero, min(f_one, za2 / (f_one+za2) ))
2797  ztrad(kp) = max(f_zero, min(f_one, f_one - zrefd(kp) ))
2798 
2799  else ! for non-conservative scattering
2800  za1 = zgam1*zgam4 + zgam2*zgam3
2801  za2 = zgam1*zgam3 + zgam2*zgam4
2802  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
2803  zrk2= 2.0 * zrk
2804 
2805  zrp = zrk * cosz
2806  zrp1 = f_one + zrp
2807  zrm1 = f_one - zrp
2808  zrpp1= f_one - zrp*zrp
2809  zrpp = sign( max(flimit, abs(zrpp1)), zrpp1 ) ! avoid numerical singularity
2810  zrkg1= zrk + zgam1
2811  zrkg3= zrk * zgam3
2812  zrkg4= zrk * zgam4
2813 
2814  zr1 = zrm1 * (za2 + zrkg3)
2815  zr2 = zrp1 * (za2 - zrkg3)
2816  zr3 = zrk2 * (zgam3 - za2*cosz)
2817  zr4 = zrpp * zrkg1
2818  zr5 = zrpp * (zrk - zgam1)
2819 
2820  zt1 = zrp1 * (za1 + zrkg4)
2821  zt2 = zrm1 * (za1 - zrkg4)
2822  zt3 = zrk2 * (zgam4 + za1*cosz)
2823 
2824 ! --- ... use exponential lookup table for transmittance, or expansion
2825 ! of exponential for low optical depth
2826 
2827  zb1 = min( zrk*ztau1, 500.0 )
2828  if ( zb1 <= od_lo ) then
2829  zexm1 = f_one - zb1 + 0.5*zb1*zb1
2830  else
2831  ftind = zb1 / (bpade + zb1)
2832  itind = ftind*ntbmx + 0.5
2833  zexm1 = exp_tbl(itind)
2834  endif
2835  zexp1 = f_one / zexm1
2836 
2837  zb2 = min( ztau1*sntz, 500.0 )
2838  if ( zb2 <= od_lo ) then
2839  zexm2 = f_one - zb2 + 0.5*zb2*zb2
2840  else
2841  ftind = zb2 / (bpade + zb2)
2842  itind = ftind*ntbmx + 0.5
2843  zexm2 = exp_tbl(itind)
2844  endif
2845  zexp2 = f_one / zexm2
2846  ze1r45 = zr4*zexp1 + zr5*zexm1
2847 
2848 ! ... collimated beam
2849  if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then
2850  zrefb(kp) = eps1
2851  ztrab(kp) = zexm2
2852  else
2853  zden1 = zssa1 / ze1r45
2854  zrefb(kp) = max(f_zero, min(f_one, &
2855  & (zr1*zexp1-zr2*zexm1-zr3*zexm2)*zden1 ))
2856  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one - &
2857  & (zt1*zexp1-zt2*zexm1-zt3*zexp2)*zden1) ))
2858  endif
2859 
2860 ! ... diffuse beam
2861  zden1 = zr4 / (ze1r45 * zrkg1)
2862  zrefd(kp) = max(f_zero, min(f_one, &
2863  & zgam2*(zexp1 - zexm1)*zden1 ))
2864  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
2865  endif ! end if_zssaw_block
2866 
2867 ! --- ... combine clear and cloudy contributions for total sky
2868 ! and calculate direct beam transmittances
2869 
2870  zrefb(kp) = zc0*zrefb1 + zc1*zrefb(kp)
2871  zrefd(kp) = zc0*zrefd1 + zc1*zrefd(kp)
2872  ztrab(kp) = zc0*ztrab1 + zc1*ztrab(kp)
2873  ztrad(kp) = zc0*ztrad1 + zc1*ztrad(kp)
2874 
2875 ! --- ... direct beam transmittance. use exponential lookup table
2876 ! for transmittance, or expansion of exponential for low
2877 ! optical depth
2878 
2879  zr1 = ztau1 * sntz
2880  if ( zr1 <= od_lo ) then
2881  zexp3 = f_one - zr1 + 0.5*zr1*zr1
2882  else
2883  ftind = zr1 / (bpade + zr1)
2884  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2885  zexp3 = exp_tbl(itind)
2886  endif
2887 
2888  zldbt(kp) = zc0*zldbt(kp) + zc1*zexp3
2889  ztdbt(k) = zldbt(kp) * ztdbt(kp)
2890 
2891 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2892 ! (must use 'orig', unscaled cloud optical depth)
2893 
2894  zr1 = ztau0 * sntz
2895  if ( zr1 <= od_lo ) then
2896  zexp4 = f_one - zr1 + 0.5*zr1*zr1
2897  else
2898  ftind = zr1 / (bpade + zr1)
2899  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2900  zexp4 = exp_tbl(itind)
2901  endif
2902 
2903  ztdbt0 = (zc0*zldbt0(k) + zc1*zexp4) * ztdbt0
2904 
2905  else ! if_zc1_block --- it is a clear layer
2906 
2907 ! --- ... direct beam transmittance
2908  ztdbt(k) = zldbt(kp) * ztdbt(kp)
2909 
2910 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2911  ztdbt0 = zldbt0(k) * ztdbt0
2912 
2913  endif ! end if_zc1_block
2914  enddo ! end do_k_loop
2915 
2916 ! --- ... perform vertical quadrature
2917 
2918  call vrtqdr &
2919 ! --- inputs:
2920  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
2921  & nlay, nlp1, &
2922 ! --- outputs:
2923  & zfu, zfd &
2924  & )
2925 
2926 ! --- ... compute upward and downward fluxes at levels
2927  do k = 1, nlp1
2928  fxupc(k,ib) = fxupc(k,ib) + zsolar*zfu(k)
2929  fxdnc(k,ib) = fxdnc(k,ib) + zsolar*zfd(k)
2930  enddo
2931 
2933 ! --- ... surface downward beam/diffused flux components
2934  zb1 = zsolar*ztdbt0
2935  zb2 = zsolar*(zfd(1) - ztdbt0)
2936 
2937  if (ibd /= 0) then
2938  sfbmc(ibd) = sfbmc(ibd) + zb1
2939  sfdfc(ibd) = sfdfc(ibd) + zb2
2940  else
2941  zf1 = 0.5 * zb1
2942  zf2 = 0.5 * zb2
2943  sfbmc(1) = sfbmc(1) + zf1
2944  sfdfc(1) = sfdfc(1) + zf2
2945  sfbmc(2) = sfbmc(2) + zf1
2946  sfdfc(2) = sfdfc(2) + zf2
2947  endif
2948 ! sfbmc(ibd) = sfbmc(ibd) + zsolar*ztdbt0
2949 ! sfdfc(ibd) = sfdfc(ibd) + zsolar*(zfd(1) - ztdbt0)
2950 
2951  endif ! end if_cf1_block
2952 
2953  enddo lab_do_jg
2954 
2955 ! --- ... end of g-point loop
2956 
2957  do ib = 1, nbdsw
2958  ftoadc = ftoadc + fxdn0(nlp1,ib)
2959  ftoau0 = ftoau0 + fxup0(nlp1,ib)
2960  fsfcu0 = fsfcu0 + fxup0(1,ib)
2961  fsfcd0 = fsfcd0 + fxdn0(1,ib)
2962  enddo
2963 
2964 !! --- ... uv-b surface downward flux
2965  ibd = nuvb - nblow + 1
2966  suvbf0 = fxdn0(1,ibd)
2967 
2968  if ( cf1 <= eps ) then ! clear column, set total-sky=clear-sky fluxes
2969  do ib = 1, nbdsw
2970  do k = 1, nlp1
2971  fxupc(k,ib) = fxup0(k,ib)
2972  fxdnc(k,ib) = fxdn0(k,ib)
2973  enddo
2974  enddo
2975 
2976  ftoauc = ftoau0
2977  fsfcuc = fsfcu0
2978  fsfcdc = fsfcd0
2979 
2980 !! --- ... surface downward beam/diffused flux components
2981  sfbmc(1) = sfbm0(1)
2982  sfdfc(1) = sfdf0(1)
2983  sfbmc(2) = sfbm0(2)
2984  sfdfc(2) = sfdf0(2)
2985 
2986 !! --- ... uv-b surface downward flux
2987  suvbfc = suvbf0
2988  else ! cloudy column, compute total-sky fluxes
2989  do ib = 1, nbdsw
2990  do k = 1, nlp1
2991  fxupc(k,ib) = cf1*fxupc(k,ib) + cf0*fxup0(k,ib)
2992  fxdnc(k,ib) = cf1*fxdnc(k,ib) + cf0*fxdn0(k,ib)
2993  enddo
2994  enddo
2995 
2996  do ib = 1, nbdsw
2997  ftoauc = ftoauc + fxupc(nlp1,ib)
2998  fsfcuc = fsfcuc + fxupc(1,ib)
2999  fsfcdc = fsfcdc + fxdnc(1,ib)
3000  enddo
3001 
3002 !! --- ... uv-b surface downward flux
3003  suvbfc = fxdnc(1,ibd)
3004 
3005 !! --- ... surface downward beam/diffused flux components
3006  sfbmc(1) = cf1*sfbmc(1) + cf0*sfbm0(1)
3007  sfbmc(2) = cf1*sfbmc(2) + cf0*sfbm0(2)
3008  sfdfc(1) = cf1*sfdfc(1) + cf0*sfdf0(1)
3009  sfdfc(2) = cf1*sfdfc(2) + cf0*sfdf0(2)
3010  endif ! end if_cf1_block
3011 
3012  return
3013 !...................................
3014  end subroutine spcvrtc
3015 !-----------------------------------
3017 
3057 !-----------------------------------
3058  subroutine spcvrtm &
3059  & ( ssolar,cosz,sntz,albbm,albdf,sfluxzen,cldfmc, & ! --- inputs
3060  & cf1,cf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
3061  & nlay, nlp1, &
3062  & fxupc,fxdnc,fxup0,fxdn0, & ! --- outputs
3063  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
3064  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
3065  & )
3067 ! =================== program usage description =================== !
3068 ! !
3069 ! purpose: computes the shortwave radiative fluxes using two-stream !
3070 ! method of h. barker and mcica, the monte-carlo independent!
3071 ! column approximation, for the representation of sub-grid !
3072 ! cloud variability (i.e. cloud overlap). !
3073 ! !
3074 ! subprograms called: vrtqdr !
3075 ! !
3076 ! ==================== defination of variables ==================== !
3077 ! !
3078 ! inputs: size !
3079 ! ssolar - real, incoming solar flux at top 1 !
3080 ! cosz - real, cosine solar zenith angle 1 !
3081 ! sntz - real, secant solar zenith angle 1 !
3082 ! albbm - real, surface albedo for direct beam radiation 2 !
3083 ! albdf - real, surface albedo for diffused radiation 2 !
3084 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
3085 ! cldfmc - real, layer cloud fraction for g-point nlay*ngptsw!
3086 ! cf1 - real, >0: cloudy sky, otherwise: clear sky 1 !
3087 ! cf0 - real, =1-cf1 1 !
3088 ! taug - real, spectral optical depth for gases nlay*ngptsw!
3089 ! taur - real, optical depth for rayleigh scattering nlay*ngptsw!
3090 ! tauae - real, aerosols optical depth nlay*nbdsw !
3091 ! ssaae - real, aerosols single scattering albedo nlay*nbdsw !
3092 ! asyae - real, aerosols asymmetry factor nlay*nbdsw !
3093 ! taucw - real, weighted cloud optical depth nlay*nbdsw !
3094 ! ssacw - real, weighted cloud single scat albedo nlay*nbdsw !
3095 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
3096 ! nlay,nlp1 - integer, number of layers/levels 1 !
3097 ! !
3098 ! output variables: !
3099 ! fxupc - real, tot sky upward flux nlp1*nbdsw !
3100 ! fxdnc - real, tot sky downward flux nlp1*nbdsw !
3101 ! fxup0 - real, clr sky upward flux nlp1*nbdsw !
3102 ! fxdn0 - real, clr sky downward flux nlp1*nbdsw !
3103 ! ftoauc - real, tot sky toa upwd flux 1 !
3104 ! ftoau0 - real, clr sky toa upwd flux 1 !
3105 ! ftoadc - real, toa downward (incoming) solar flux 1 !
3106 ! fsfcuc - real, tot sky sfc upwd flux 1 !
3107 ! fsfcu0 - real, clr sky sfc upwd flux 1 !
3108 ! fsfcdc - real, tot sky sfc dnwd flux 1 !
3109 ! fsfcd0 - real, clr sky sfc dnwd flux 1 !
3110 ! sfbmc - real, tot sky sfc dnwd beam flux (nir/uv+vis) 2 !
3111 ! sfdfc - real, tot sky sfc dnwd diff flux (nir/uv+vis) 2 !
3112 ! sfbm0 - real, clr sky sfc dnwd beam flux (nir/uv+vis) 2 !
3113 ! sfdf0 - real, clr sky sfc dnwd diff flux (nir/uv+vis) 2 !
3114 ! suvbfc - real, tot sky sfc dnwd uv-b flux 1 !
3115 ! suvbf0 - real, clr sky sfc dnwd uv-b flux 1 !
3116 ! !
3117 ! internal variables: !
3118 ! zrefb - real, direct beam reflectivity for clear/cloudy nlp1 !
3119 ! zrefd - real, diffuse reflectivity for clear/cloudy nlp1 !
3120 ! ztrab - real, direct beam transmissivity for clear/cloudy nlp1 !
3121 ! ztrad - real, diffuse transmissivity for clear/cloudy nlp1 !
3122 ! zldbt - real, layer beam transmittance for clear/cloudy nlp1 !
3123 ! ztdbt - real, lev total beam transmittance for clr/cld nlp1 !
3124 ! !
3125 ! control parameters in module "physparam" !
3126 ! iswmode - control flag for 2-stream transfer schemes !
3127 ! = 1 delta-eddington (joseph et al., 1976) !
3128 ! = 2 pifm (zdunkowski et al., 1980) !
3129 ! = 3 discrete ordinates (liou, 1973) !
3130 ! !
3131 ! ******************************************************************* !
3132 ! original code description !
3133 ! !
3134 ! method: !
3135 ! ------- !
3136 ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations. !
3137 ! kmodts = 1 eddington (joseph et al., 1976) !
3138 ! = 2 pifm (zdunkowski et al., 1980) !
3139 ! = 3 discrete ordinates (liou, 1973) !
3140 ! !
3141 ! modifications: !
3142 ! -------------- !
3143 ! original: h. barker !
3144 ! revision: merge with rrtmg_sw: j.-j.morcrette, ecmwf, feb 2003 !
3145 ! revision: add adjustment for earth/sun distance:mjiacono,aer,oct2003!
3146 ! revision: bug fix for use of palbp and palbd: mjiacono, aer, nov2003!
3147 ! revision: bug fix to apply delta scaling to clear sky: aer, dec2004 !
3148 ! revision: code modified so that delta scaling is not done in cloudy !
3149 ! profiles if routine cldprop is used; delta scaling can be !
3150 ! applied by swithcing code below if cldprop is not used to !
3151 ! get cloud properties. aer, jan 2005 !
3152 ! revision: uniform formatting for rrtmg: mjiacono, aer, jul 2006 !
3153 ! revision: use exponential lookup table for transmittance: mjiacono, !
3154 ! aer, aug 2007 !
3155 ! !
3156 ! ******************************************************************* !
3157 ! ====================== end of description block ================= !
3158 
3159 ! --- constant parameters:
3160  real (kind=kind_phys), parameter :: zcrit = 0.9999995 ! thresold for conservative scattering
3161  real (kind=kind_phys), parameter :: zsr3 = sqrt(3.0)
3162  real (kind=kind_phys), parameter :: od_lo = 0.06
3163  real (kind=kind_phys), parameter :: eps1 = 1.0e-8
3164 
3165 ! --- inputs:
3166  integer, intent(in) :: nlay, nlp1
3167 
3168  real (kind=kind_phys), dimension(nlay,ngptsw), intent(in) :: &
3169  & taug, taur, cldfmc
3170  real (kind=kind_phys), dimension(nlay,nbdsw), intent(in) :: &
3171  & taucw, ssacw, asycw, tauae, ssaae, asyae
3172 
3173  real (kind=kind_phys), dimension(ngptsw), intent(in) :: sfluxzen
3174 
3175  real (kind=kind_phys), dimension(2), intent(in) :: albbm, albdf
3176 
3177  real (kind=kind_phys), intent(in) :: cosz, sntz, cf1, cf0, ssolar
3178 
3179 ! --- outputs:
3180  real (kind=kind_phys), dimension(nlp1,nbdsw), intent(out) :: &
3181  & fxupc, fxdnc, fxup0, fxdn0
3182 
3183  real (kind=kind_phys), dimension(2), intent(out) :: sfbmc, sfdfc, &
3184  & sfbm0, sfdf0
3185 
3186  real (kind=kind_phys), intent(out) :: suvbfc, suvbf0, ftoadc, &
3187  & ftoauc, ftoau0, fsfcuc, fsfcu0, fsfcdc, fsfcd0
3188 
3189 ! --- locals:
3190  real (kind=kind_phys), dimension(nlay) :: ztaus, zssas, zasys, &
3191  & zldbt0
3192 
3193  real (kind=kind_phys), dimension(nlp1) :: zrefb, zrefd, ztrab, &
3194  & ztrad, ztdbt, zldbt, zfu, zfd
3195 
3196  real (kind=kind_phys) :: ztau1, zssa1, zasy1, ztau0, zssa0, &
3197  & zasy0, zasy3, zssaw, zasyw, zgam1, zgam2, zgam3, zgam4, &
3198  & za1, za2, zb1, zb2, zrk, zrk2, zrp, zrp1, zrm1, zrpp, &
3199  & zrkg1, zrkg3, zrkg4, zexp1, zexm1, zexp2, zexm2, zden1, &
3200  & zexp3, zexp4, ze1r45, ftind, zsolar, ztdbt0, zr1, zr2, &
3201  & zr3, zr4, zr5, zt1, zt2, zt3, zf1, zf2, zrpp1
3202 
3203  integer :: ib, ibd, jb, jg, k, kp, itind
3204 !
3205 !===> ... begin here
3206 !
3208 
3209  do ib = 1, nbdsw
3210  do k = 1, nlp1
3211  fxdnc(k,ib) = f_zero
3212  fxupc(k,ib) = f_zero
3213  fxdn0(k,ib) = f_zero
3214  fxup0(k,ib) = f_zero
3215  enddo
3216  enddo
3217 
3218  ftoadc = f_zero
3219  ftoauc = f_zero
3220  ftoau0 = f_zero
3221  fsfcuc = f_zero
3222  fsfcu0 = f_zero
3223  fsfcdc = f_zero
3224  fsfcd0 = f_zero
3225 
3226 !! --- ... uv-b surface downward fluxes
3227  suvbfc = f_zero
3228  suvbf0 = f_zero
3229 
3230 !! --- ... output surface flux components
3231  sfbmc(1) = f_zero
3232  sfbmc(2) = f_zero
3233  sfdfc(1) = f_zero
3234  sfdfc(2) = f_zero
3235  sfbm0(1) = f_zero
3236  sfbm0(2) = f_zero
3237  sfdf0(1) = f_zero
3238  sfdf0(2) = f_zero
3239 
3240 ! --- ... loop over all g-points in each band
3241 
3242  lab_do_jg : do jg = 1, ngptsw
3243 
3244  jb = ngb(jg)
3245  ib = jb + 1 - nblow
3246  ibd = idxsfc(jb) ! spectral band index
3247 
3248  zsolar = ssolar * sfluxzen(jg)
3249 
3250 ! --- ... set up toa direct beam and surface values (beam and diff)
3251 
3252  ztdbt(nlp1) = f_one
3253  ztdbt0 = f_one
3254 
3255  zldbt(1) = f_zero
3256  if (ibd /= 0) then
3257  zrefb(1) = albbm(ibd)
3258  zrefd(1) = albdf(ibd)
3259  else
3260  zrefb(1) = 0.5 * (albbm(1) + albbm(2))
3261  zrefd(1) = 0.5 * (albdf(1) + albdf(2))
3262  endif
3263  ztrab(1) = f_zero
3264  ztrad(1) = f_zero
3265 
3275 
3276  do k = nlay, 1, -1
3277  kp = k + 1
3278 
3279  ztau0 = max( ftiny, taur(k,jg)+taug(k,jg)+tauae(k,ib) )
3280  zssa0 = taur(k,jg) + tauae(k,ib)*ssaae(k,ib)
3281  zasy0 = asyae(k,ib)*ssaae(k,ib)*tauae(k,ib)
3282  zssaw = min( oneminus, zssa0 / ztau0 )
3283  zasyw = zasy0 / max( ftiny, zssa0 )
3284 
3285 ! --- ... saving clear-sky quantities for later total-sky usage
3286  ztaus(k) = ztau0
3287  zssas(k) = zssa0
3288  zasys(k) = zasy0
3289 
3290 ! --- ... delta scaling for clear-sky condition
3291  za1 = zasyw * zasyw
3292  za2 = zssaw * za1
3293 
3294  ztau1 = (f_one - za2) * ztau0
3295  zssa1 = (zssaw - za2) / (f_one - za2)
3296 !org zasy1 = (zasyw - za1) / (f_one - za1) ! this line is replaced by the next
3297  zasy1 = zasyw / (f_one + zasyw) ! to reduce truncation error
3298  zasy3 = 0.75 * zasy1
3299 
3300 ! --- ... general two-stream expressions
3301  if ( iswmode == 1 ) then
3302  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
3303  zgam2 =-0.25 + zssa1 * (f_one - zasy3)
3304  zgam3 = 0.5 - zasy3 * cosz
3305  elseif ( iswmode == 2 ) then ! pifm
3306  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
3307  zgam2 = 0.75* zssa1 * (f_one- zasy1)
3308  zgam3 = 0.5 - zasy3 * cosz
3309  elseif ( iswmode == 3 ) then ! discrete ordinates
3310  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
3311  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
3312  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
3313  endif
3314  zgam4 = f_one - zgam3
3315 
3316 ! --- ... compute homogeneous reflectance and transmittance
3317 
3318  if ( zssaw >= zcrit ) then ! for conservative scattering
3319  za1 = zgam1 * cosz - zgam3
3320  za2 = zgam1 * ztau1
3321 
3322 ! --- ... use exponential lookup table for transmittance, or expansion
3323 ! of exponential for low optical depth
3324 
3325  zb1 = min( ztau1*sntz , 500.0 )
3326  if ( zb1 <= od_lo ) then
3327  zb2 = f_one - zb1 + 0.5*zb1*zb1
3328  else
3329  ftind = zb1 / (bpade + zb1)
3330  itind = ftind*ntbmx + 0.5
3331  zb2 = exp_tbl(itind)
3332  endif
3333 
3334 ! ... collimated beam
3335  zrefb(kp) = max(f_zero, min(f_one, &
3336  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
3337  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp) ))
3338 
3339 ! ... isotropic incidence
3340  zrefd(kp) = max(f_zero, min(f_one, za2/(f_one + za2) ))
3341  ztrad(kp) = max(f_zero, min(f_one, f_one-zrefd(kp) ))
3342 
3343  else ! for non-conservative scattering
3344  za1 = zgam1*zgam4 + zgam2*zgam3
3345  za2 = zgam1*zgam3 + zgam2*zgam4
3346  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
3347  zrk2= 2.0 * zrk
3348 
3349  zrp = zrk * cosz
3350  zrp1 = f_one + zrp
3351  zrm1 = f_one - zrp
3352  zrpp1= f_one - zrp*zrp
3353  zrpp = sign( max(flimit, abs(zrpp1)), zrpp1 ) ! avoid numerical singularity
3354  zrkg1= zrk + zgam1
3355  zrkg3= zrk * zgam3
3356  zrkg4= zrk * zgam4
3357 
3358  zr1 = zrm1 * (za2 + zrkg3)
3359  zr2 = zrp1 * (za2 - zrkg3)
3360  zr3 = zrk2 * (zgam3 - za2*cosz)
3361  zr4 = zrpp * zrkg1
3362  zr5 = zrpp * (zrk - zgam1)
3363 
3364  zt1 = zrp1 * (za1 + zrkg4)
3365  zt2 = zrm1 * (za1 - zrkg4)
3366  zt3 = zrk2 * (zgam4 + za1*cosz)
3367 
3368 ! --- ... use exponential lookup table for transmittance, or expansion
3369 ! of exponential for low optical depth
3370 
3371  zb1 = min( zrk*ztau1, 500.0 )
3372  if ( zb1 <= od_lo ) then
3373  zexm1 = f_one - zb1 + 0.5*zb1*zb1
3374  else
3375  ftind = zb1 / (bpade + zb1)
3376  itind = ftind*ntbmx + 0.5
3377  zexm1 = exp_tbl(itind)
3378  endif
3379  zexp1 = f_one / zexm1
3380 
3381  zb2 = min( sntz*ztau1, 500.0 )
3382  if ( zb2 <= od_lo ) then
3383  zexm2 = f_one - zb2 + 0.5*zb2*zb2
3384  else
3385  ftind = zb2 / (bpade + zb2)
3386  itind = ftind*ntbmx + 0.5
3387  zexm2 = exp_tbl(itind)
3388  endif
3389  zexp2 = f_one / zexm2
3390  ze1r45 = zr4*zexp1 + zr5*zexm1
3391 
3392 ! ... collimated beam
3393  if (ze1r45>=-eps1 .and. ze1r45<=eps1) then
3394  zrefb(kp) = eps1
3395  ztrab(kp) = zexm2
3396  else
3397  zden1 = zssa1 / ze1r45
3398  zrefb(kp) = max(f_zero, min(f_one, &
3399  & (zr1*zexp1 - zr2*zexm1 - zr3*zexm2)*zden1 ))
3400  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one &
3401  & - (zt1*zexp1 - zt2*zexm1 - zt3*zexp2)*zden1) ))
3402  endif
3403 
3404 ! ... diffuse beam
3405  zden1 = zr4 / (ze1r45 * zrkg1)
3406  zrefd(kp) = max(f_zero, min(f_one, &
3407  & zgam2*(zexp1 - zexm1)*zden1 ))
3408  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
3409  endif ! end if_zssaw_block
3410 
3411 ! --- ... direct beam transmittance. use exponential lookup table
3412 ! for transmittance, or expansion of exponential for low
3413 ! optical depth
3414 
3415  zr1 = ztau1 * sntz
3416  if ( zr1 <= od_lo ) then
3417  zexp3 = f_one - zr1 + 0.5*zr1*zr1
3418  else
3419  ftind = zr1 / (bpade + zr1)
3420  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3421  zexp3 = exp_tbl(itind)
3422  endif
3423 
3424  ztdbt(k) = zexp3 * ztdbt(kp)
3425  zldbt(kp) = zexp3
3426 
3427 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3428 ! (must use 'orig', unscaled cloud optical depth)
3429 
3430  zr1 = ztau0 * sntz
3431  if ( zr1 <= od_lo ) then
3432  zexp4 = f_one - zr1 + 0.5*zr1*zr1
3433  else
3434  ftind = zr1 / (bpade + zr1)
3435  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3436  zexp4 = exp_tbl(itind)
3437  endif
3438 
3439  zldbt0(k) = zexp4
3440  ztdbt0 = zexp4 * ztdbt0
3441  enddo ! end do_k_loop
3442 
3443  call vrtqdr &
3444 ! --- inputs:
3445  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3446  & nlay, nlp1, &
3447 ! --- outputs:
3448  & zfu, zfd &
3449  & )
3450 
3451 ! --- ... compute upward and downward fluxes at levels
3452  do k = 1, nlp1
3453  fxup0(k,ib) = fxup0(k,ib) + zsolar*zfu(k)
3454  fxdn0(k,ib) = fxdn0(k,ib) + zsolar*zfd(k)
3455  enddo
3456 
3457 !! --- ... surface downward beam/diffuse flux components
3458  zb1 = zsolar*ztdbt0
3459  zb2 = zsolar*(zfd(1) - ztdbt0)
3460 
3461  if (ibd /= 0) then
3462  sfbm0(ibd) = sfbm0(ibd) + zb1
3463  sfdf0(ibd) = sfdf0(ibd) + zb2
3464  else
3465  zf1 = 0.5 * zb1
3466  zf2 = 0.5 * zb2
3467  sfbm0(1) = sfbm0(1) + zf1
3468  sfdf0(1) = sfdf0(1) + zf2
3469  sfbm0(2) = sfbm0(2) + zf1
3470  sfdf0(2) = sfdf0(2) + zf2
3471  endif
3472 ! sfbm0(ibd) = sfbm0(ibd) + zsolar*ztdbt0
3473 ! sfdf0(ibd) = sfdf0(ibd) + zsolar*(zfd(1) - ztdbt0)
3474 
3484 
3485  if ( cf1 > eps ) then
3486 
3487 ! --- ... set up toa direct beam and surface values (beam and diff)
3488  ztdbt0 = f_one
3489  zldbt(1) = f_zero
3490 
3491  do k = nlay, 1, -1
3492  kp = k + 1
3493  if ( cldfmc(k,jg) > ftiny ) then ! it is a cloudy-layer
3494 
3495  ztau0 = ztaus(k) + taucw(k,ib)
3496  zssa0 = zssas(k) + ssacw(k,ib)
3497  zasy0 = zasys(k) + asycw(k,ib)
3498  zssaw = min(oneminus, zssa0 / ztau0)
3499  zasyw = zasy0 / max(ftiny, zssa0)
3500 
3501 ! --- ... delta scaling for total-sky condition
3502  za1 = zasyw * zasyw
3503  za2 = zssaw * za1
3504 
3505  ztau1 = (f_one - za2) * ztau0
3506  zssa1 = (zssaw - za2) / (f_one - za2)
3507 !org zasy1 = (zasyw - za1) / (f_one - za1)
3508  zasy1 = zasyw / (f_one + zasyw)
3509  zasy3 = 0.75 * zasy1
3510 
3511 ! --- ... general two-stream expressions
3512  if ( iswmode == 1 ) then
3513  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
3514  zgam2 =-0.25 + zssa1 * (f_one - zasy3)
3515  zgam3 = 0.5 - zasy3 * cosz
3516  elseif ( iswmode == 2 ) then ! pifm
3517  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
3518  zgam2 = 0.75* zssa1 * (f_one- zasy1)
3519  zgam3 = 0.5 - zasy3 * cosz
3520  elseif ( iswmode == 3 ) then ! discrete ordinates
3521  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
3522  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
3523  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
3524  endif
3525  zgam4 = f_one - zgam3
3526 
3527 ! --- ... compute homogeneous reflectance and transmittance
3528 
3529  if ( zssaw >= zcrit ) then ! for conservative scattering
3530  za1 = zgam1 * cosz - zgam3
3531  za2 = zgam1 * ztau1
3532 
3533 ! --- ... use exponential lookup table for transmittance, or expansion
3534 ! of exponential for low optical depth
3535 
3536  zb1 = min( ztau1*sntz , 500.0 )
3537  if ( zb1 <= od_lo ) then
3538  zb2 = f_one - zb1 + 0.5*zb1*zb1
3539  else
3540  ftind = zb1 / (bpade + zb1)
3541  itind = ftind*ntbmx + 0.5
3542  zb2 = exp_tbl(itind)
3543  endif
3544 
3545 ! ... collimated beam
3546  zrefb(kp) = max(f_zero, min(f_one, &
3547  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
3548  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp)))
3549 
3550 ! ... isotropic incidence
3551  zrefd(kp) = max(f_zero, min(f_one, za2 / (f_one+za2) ))
3552  ztrad(kp) = max(f_zero, min(f_one, f_one - zrefd(kp) ))
3553 
3554  else ! for non-conservative scattering
3555  za1 = zgam1*zgam4 + zgam2*zgam3
3556  za2 = zgam1*zgam3 + zgam2*zgam4
3557  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
3558  zrk2= 2.0 * zrk
3559 
3560  zrp = zrk * cosz
3561  zrp1 = f_one + zrp
3562  zrm1 = f_one - zrp
3563  zrpp1= f_one - zrp*zrp
3564  zrpp = sign( max(flimit, abs(zrpp1)), zrpp1 ) ! avoid numerical singularity
3565  zrkg1= zrk + zgam1
3566  zrkg3= zrk * zgam3
3567  zrkg4= zrk * zgam4
3568 
3569  zr1 = zrm1 * (za2 + zrkg3)
3570  zr2 = zrp1 * (za2 - zrkg3)
3571  zr3 = zrk2 * (zgam3 - za2*cosz)
3572  zr4 = zrpp * zrkg1
3573  zr5 = zrpp * (zrk - zgam1)
3574 
3575  zt1 = zrp1 * (za1 + zrkg4)
3576  zt2 = zrm1 * (za1 - zrkg4)
3577  zt3 = zrk2 * (zgam4 + za1*cosz)
3578 
3579 ! --- ... use exponential lookup table for transmittance, or expansion
3580 ! of exponential for low optical depth
3581 
3582  zb1 = min( zrk*ztau1, 500.0 )
3583  if ( zb1 <= od_lo ) then
3584  zexm1 = f_one - zb1 + 0.5*zb1*zb1
3585  else
3586  ftind = zb1 / (bpade + zb1)
3587  itind = ftind*ntbmx + 0.5
3588  zexm1 = exp_tbl(itind)
3589  endif
3590  zexp1 = f_one / zexm1
3591 
3592  zb2 = min( ztau1*sntz, 500.0 )
3593  if ( zb2 <= od_lo ) then
3594  zexm2 = f_one - zb2 + 0.5*zb2*zb2
3595  else
3596  ftind = zb2 / (bpade + zb2)
3597  itind = ftind*ntbmx + 0.5
3598  zexm2 = exp_tbl(itind)
3599  endif
3600  zexp2 = f_one / zexm2
3601  ze1r45 = zr4*zexp1 + zr5*zexm1
3602 
3603 ! ... collimated beam
3604  if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then
3605  zrefb(kp) = eps1
3606  ztrab(kp) = zexm2
3607  else
3608  zden1 = zssa1 / ze1r45
3609  zrefb(kp) = max(f_zero, min(f_one, &
3610  & (zr1*zexp1-zr2*zexm1-zr3*zexm2)*zden1 ))
3611  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one - &
3612  & (zt1*zexp1-zt2*zexm1-zt3*zexp2)*zden1) ))
3613  endif
3614 
3615 ! ... diffuse beam
3616  zden1 = zr4 / (ze1r45 * zrkg1)
3617  zrefd(kp) = max(f_zero, min(f_one, &
3618  & zgam2*(zexp1 - zexm1)*zden1 ))
3619  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
3620  endif ! end if_zssaw_block
3621 
3622 ! --- ... direct beam transmittance. use exponential lookup table
3623 ! for transmittance, or expansion of exponential for low
3624 ! optical depth
3625 
3626  zr1 = ztau1 * sntz
3627  if ( zr1 <= od_lo ) then
3628  zexp3 = f_one - zr1 + 0.5*zr1*zr1
3629  else
3630  ftind = zr1 / (bpade + zr1)
3631  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3632  zexp3 = exp_tbl(itind)
3633  endif
3634 
3635  zldbt(kp) = zexp3
3636  ztdbt(k) = zexp3 * ztdbt(kp)
3637 
3638 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3639 ! (must use 'orig', unscaled cloud optical depth)
3640 
3641  zr1 = ztau0 * sntz
3642  if ( zr1 <= od_lo ) then
3643  zexp4 = f_one - zr1 + 0.5*zr1*zr1
3644  else
3645  ftind = zr1 / (bpade + zr1)
3646  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3647  zexp4 = exp_tbl(itind)
3648  endif
3649 
3650  ztdbt0 = zexp4 * ztdbt0
3651 
3652  else ! if_cldfmc_block --- it is a clear layer
3653 
3654 ! --- ... direct beam transmittance
3655  ztdbt(k) = zldbt(kp) * ztdbt(kp)
3656 
3657 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3658  ztdbt0 = zldbt0(k) * ztdbt0
3659 
3660  endif ! end if_cldfmc_block
3661  enddo ! end do_k_loop
3662 
3663 ! --- ... perform vertical quadrature
3664 
3665  call vrtqdr &
3666 ! --- inputs:
3667  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3668  & nlay, nlp1, &
3669 ! --- outputs:
3670  & zfu, zfd &
3671  & )
3672 
3673 ! --- ... compute upward and downward fluxes at levels
3674  do k = 1, nlp1
3675  fxupc(k,ib) = fxupc(k,ib) + zsolar*zfu(k)
3676  fxdnc(k,ib) = fxdnc(k,ib) + zsolar*zfd(k)
3677  enddo
3678 
3680 ! --- ... surface downward beam/diffused flux components
3681  zb1 = zsolar*ztdbt0
3682  zb2 = zsolar*(zfd(1) - ztdbt0)
3683 
3684  if (ibd /= 0) then
3685  sfbmc(ibd) = sfbmc(ibd) + zb1
3686  sfdfc(ibd) = sfdfc(ibd) + zb2
3687  else
3688  zf1 = 0.5 * zb1
3689  zf2 = 0.5 * zb2
3690  sfbmc(1) = sfbmc(1) + zf1
3691  sfdfc(1) = sfdfc(1) + zf2
3692  sfbmc(2) = sfbmc(2) + zf1
3693  sfdfc(2) = sfdfc(2) + zf2
3694  endif
3695 ! sfbmc(ibd) = sfbmc(ibd) + zsolar*ztdbt0
3696 ! sfdfc(ibd) = sfdfc(ibd) + zsolar*(zfd(1) - ztdbt0)
3697 
3698  endif ! end if_cf1_block
3699 
3700  enddo lab_do_jg
3701 
3702 ! --- ... end of g-point loop
3703 
3704  do ib = 1, nbdsw
3705  ftoadc = ftoadc + fxdn0(nlp1,ib)
3706  ftoau0 = ftoau0 + fxup0(nlp1,ib)
3707  fsfcu0 = fsfcu0 + fxup0(1,ib)
3708  fsfcd0 = fsfcd0 + fxdn0(1,ib)
3709  enddo
3710 
3711 !! --- ... uv-b surface downward flux
3712  ibd = nuvb - nblow + 1
3713  suvbf0 = fxdn0(1,ibd)
3714 
3715  if ( cf1 <= eps ) then ! clear column, set total-sky=clear-sky fluxes
3716  do ib = 1, nbdsw
3717  do k = 1, nlp1
3718  fxupc(k,ib) = fxup0(k,ib)
3719  fxdnc(k,ib) = fxdn0(k,ib)
3720  enddo
3721  enddo
3722 
3723  ftoauc = ftoau0
3724  fsfcuc = fsfcu0
3725  fsfcdc = fsfcd0
3726 
3727 !! --- ... surface downward beam/diffused flux components
3728  sfbmc(1) = sfbm0(1)
3729  sfdfc(1) = sfdf0(1)
3730  sfbmc(2) = sfbm0(2)
3731  sfdfc(2) = sfdf0(2)
3732 
3733 !! --- ... uv-b surface downward flux
3734  suvbfc = suvbf0
3735  else ! cloudy column, compute total-sky fluxes
3736  do ib = 1, nbdsw
3737  ftoauc = ftoauc + fxupc(nlp1,ib)
3738  fsfcuc = fsfcuc + fxupc(1,ib)
3739  fsfcdc = fsfcdc + fxdnc(1,ib)
3740  enddo
3741 
3742 !! --- ... uv-b surface downward flux
3743  suvbfc = fxdnc(1,ibd)
3744  endif ! end if_cf1_block
3745 
3746  return
3747 !...................................
3748  end subroutine spcvrtm
3749 !-----------------------------------
3750 
3764 !-----------------------------------
3765  subroutine vrtqdr &
3766  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, & ! inputs
3767  & nlay, nlp1, &
3768  & zfu, zfd & ! outputs:
3769  & )
3771 ! =================== program usage description =================== !
3772 ! !
3773 ! purpose: computes the upward and downward radiation fluxes !
3774 ! !
3775 ! interface: "vrtqdr" is called by "spcvrc" and "spcvrm" !
3776 ! !
3777 ! subroutines called : none !
3778 ! !
3779 ! ==================== defination of variables ==================== !
3780 ! !
3781 ! input variables: !
3782 ! zrefb(NLP1) - layer direct beam reflectivity !
3783 ! zrefd(NLP1) - layer diffuse reflectivity !
3784 ! ztrab(NLP1) - layer direct beam transmissivity !
3785 ! ztrad(NLP1) - layer diffuse transmissivity !
3786 ! zldbt(NLP1) - layer mean beam transmittance !
3787 ! ztdbt(NLP1) - total beam transmittance at levels !
3788 ! NLAY, NLP1 - number of layers/levels !
3789 ! !
3790 ! output variables: !
3791 ! zfu (NLP1) - upward flux at layer interface !
3792 ! zfd (NLP1) - downward flux at layer interface !
3793 ! !
3794 ! ******************************************************************* !
3795 ! ====================== end of description block ================= !
3796 
3797 ! --- inputs:
3798  integer, intent(in) :: nlay, nlp1
3799 
3800  real (kind=kind_phys), dimension(nlp1), intent(in) :: zrefb, &
3801  & zrefd, ztrab, ztrad, ztdbt, zldbt
3802 
3803 ! --- outputs:
3804  real (kind=kind_phys), dimension(nlp1), intent(out) :: zfu, zfd
3805 
3806 ! --- locals:
3807  real (kind=kind_phys), dimension(nlp1) :: zrupb,zrupd,zrdnd,ztdn
3808 
3809  real (kind=kind_phys) :: zden1
3810 
3811  integer :: k, kp
3812 !
3813 !===> ... begin here
3814 !
3815 
3817  zrupb(1) = zrefb(1) ! direct beam
3818  zrupd(1) = zrefd(1) ! diffused
3819 
3821  do k = 1, nlay
3822  kp = k + 1
3823 
3824  zden1 = f_one / ( f_one - zrupd(k)*zrefd(kp) )
3825  zrupb(kp) = zrefb(kp) + ( ztrad(kp) * &
3826  & ( (ztrab(kp) - zldbt(kp))*zrupd(k) + &
3827  & zldbt(kp)*zrupb(k)) ) * zden1
3828  zrupd(kp) = zrefd(kp) + ztrad(kp)*ztrad(kp)*zrupd(k)*zden1
3829  enddo
3830 
3832  ztdn(nlp1) = f_one
3833  zrdnd(nlp1) = f_zero
3834  ztdn(nlay) = ztrab(nlp1)
3835  zrdnd(nlay) = zrefd(nlp1)
3836 
3838  do k = nlay, 2, -1
3839  zden1 = f_one / (f_one - zrefd(k)*zrdnd(k))
3840  ztdn(k-1) = ztdbt(k)*ztrab(k) + ( ztrad(k) * &
3841  & ( (ztdn(k) - ztdbt(k)) + ztdbt(k) * &
3842  & zrefb(k)*zrdnd(k) )) * zden1
3843  zrdnd(k-1) = zrefd(k) + ztrad(k)*ztrad(k)*zrdnd(k)*zden1
3844  enddo
3845 
3847  do k = 1, nlp1
3848  zden1 = f_one / (f_one - zrdnd(k)*zrupd(k))
3849  zfu(k) = ( ztdbt(k)*zrupb(k) + &
3850  & (ztdn(k) - ztdbt(k))*zrupd(k) ) * zden1
3851  zfd(k) = ztdbt(k) + ( ztdn(k) - ztdbt(k) + &
3852  & ztdbt(k)*zrupb(k)*zrdnd(k) ) * zden1
3853  enddo
3854 
3855  return
3856 !...................................
3857  end subroutine vrtqdr
3858 !-----------------------------------
3860 
3900 !-----------------------------------
3901  subroutine taumol &
3902  & ( colamt,colmol,fac00,fac01,fac10,fac11,jp,jt,jt1,laytrop, & ! --- inputs
3903  & forfac,forfrac,indfor,selffac,selffrac,indself, nlay, &
3904  & sfluxzen, taug, taur & ! --- outputs
3905  & )
3907 ! ================== program usage description ================== !
3908 ! !
3909 ! description: !
3910 ! calculate optical depths for gaseous absorption and rayleigh !
3911 ! scattering. !
3912 ! !
3913 ! subroutines called: taugb## (## = 16 - 29) !
3914 ! !
3915 ! ==================== defination of variables ==================== !
3916 ! !
3917 ! inputs: size !
3918 ! colamt - real, column amounts of absorbing gases the index !
3919 ! are for h2o, co2, o3, n2o, ch4, and o2, !
3920 ! respectively (molecules/cm**2) nlay*maxgas!
3921 ! colmol - real, total column amount (dry air+water vapor) nlay !
3922 ! facij - real, for each layer, these are factors that are !
3923 ! needed to compute the interpolation factors !
3924 ! that multiply the appropriate reference k- !
3925 ! values. a value of 0/1 for i,j indicates !
3926 ! that the corresponding factor multiplies !
3927 ! reference k-value for the lower/higher of the !
3928 ! two appropriate temperatures, and altitudes, !
3929 ! respectively. naly !
3930 ! jp - real, the index of the lower (in altitude) of the !
3931 ! two appropriate ref pressure levels needed !
3932 ! for interpolation. nlay !
3933 ! jt, jt1 - integer, the indices of the lower of the two approp !
3934 ! ref temperatures needed for interpolation (for !
3935 ! pressure levels jp and jp+1, respectively) nlay !
3936 ! laytrop - integer, tropopause layer index 1 !
3937 ! forfac - real, scale factor needed to foreign-continuum. nlay !
3938 ! forfrac - real, factor needed for temperature interpolation nlay !
3939 ! indfor - integer, index of the lower of the two appropriate !
3940 ! reference temperatures needed for foreign- !
3941 ! continuum interpolation nlay !
3942 ! selffac - real, scale factor needed to h2o self-continuum. nlay !
3943 ! selffrac- real, factor needed for temperature interpolation !
3944 ! of reference h2o self-continuum data nlay !
3945 ! indself - integer, index of the lower of the two appropriate !
3946 ! reference temperatures needed for the self- !
3947 ! continuum interpolation nlay !
3948 ! nlay - integer, number of vertical layers 1 !
3949 ! !
3950 ! output: !
3951 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
3952 ! taug - real, spectral optical depth for gases nlay*ngptsw!
3953 ! taur - real, opt depth for rayleigh scattering nlay*ngptsw!
3954 ! !
3955 ! =================================================================== !
3956 ! ************ original subprogram description *************** !
3957 ! !
3958 ! optical depths developed for the !
3959 ! !
3960 ! rapid radiative transfer model (rrtm) !
3961 ! !
3962 ! atmospheric and environmental research, inc. !
3963 ! 131 hartwell avenue !
3964 ! lexington, ma 02421 !
3965 ! !
3966 ! !
3967 ! eli j. mlawer !
3968 ! jennifer delamere !
3969 ! steven j. taubman !
3970 ! shepard a. clough !
3971 ! !
3972 ! !
3973 ! !
3974 ! email: mlawer@aer.com !
3975 ! email: jdelamer@aer.com !
3976 ! !
3977 ! the authors wish to acknowledge the contributions of the !
3978 ! following people: patrick d. brown, michael j. iacono, !
3979 ! ronald e. farren, luke chen, robert bergstrom. !
3980 ! !
3981 ! ******************************************************************* !
3982 ! !
3983 ! taumol !
3984 ! !
3985 ! this file contains the subroutines taugbn (where n goes from !
3986 ! 16 to 29). taugbn calculates the optical depths and Planck !
3987 ! fractions per g-value and layer for band n. !
3988 ! !
3989 ! output: optical depths (unitless) !
3990 ! fractions needed to compute planck functions at every layer !
3991 ! and g-value !
3992 ! !
3993 ! modifications: !
3994 ! !
3995 ! revised: adapted to f90 coding, j.-j.morcrette, ecmwf, feb 2003 !
3996 ! revised: modified for g-point reduction, mjiacono, aer, dec 2003 !
3997 ! revised: reformatted for consistency with rrtmg_lw, mjiacono, aer, !
3998 ! jul 2006 !
3999 ! !
4000 ! ******************************************************************* !
4001 ! ====================== end of description block ================= !
4002 
4003 ! --- inputs:
4004  integer, intent(in) :: nlay, laytrop
4005 
4006  integer, dimension(nlay), intent(in) :: indfor, indself, &
4007  & jp, jt, jt1
4008 
4009  real (kind=kind_phys), dimension(nlay), intent(in) :: colmol, &
4010  & fac00, fac01, fac10, fac11, forfac, forfrac, selffac, &
4011  & selffrac
4012 
4013  real (kind=kind_phys), dimension(nlay,maxgas),intent(in) :: colamt
4014 
4015 ! --- outputs:
4016  real (kind=kind_phys), dimension(ngptsw), intent(out) :: sfluxzen
4017 
4018  real (kind=kind_phys), dimension(nlay,ngptsw), intent(out) :: &
4019  & taug, taur
4020 
4021 ! --- locals:
4022  real (kind=kind_phys) :: fs, speccomb, specmult, colm1, colm2
4023 
4024  integer, dimension(nlay,nblow:nbhgh) :: id0, id1
4025 
4026  integer :: ibd, j, jb, js, k, klow, khgh, klim, ks, njb, ns
4027 !
4028 !===> ... begin here
4029 !
4030 ! --- ... loop over each spectral band
4031 
4032  do jb = nblow, nbhgh
4033 
4034 ! --- ... indices for layer optical depth
4035 
4036  do k = 1, laytrop
4037  id0(k,jb) = ((jp(k)-1)*5 + (jt(k)-1)) * nspa(jb)
4038  id1(k,jb) = ( jp(k) *5 + (jt1(k)-1)) * nspa(jb)
4039  enddo
4040 
4041  do k = laytrop+1, nlay
4042  id0(k,jb) = ((jp(k)-13)*5 + (jt(k)-1)) * nspb(jb)
4043  id1(k,jb) = ((jp(k)-12)*5 + (jt1(k)-1)) * nspb(jb)
4044  enddo
4045 
4046 ! --- ... calculate spectral flux at toa
4047 
4048  ibd = ibx(jb)
4049  njb = ng(jb)
4050  ns = ngs(jb)
4051 
4052  select case (jb)
4053 
4054  case (16, 20, 23, 25, 26, 29)
4055 
4056  do j = 1, njb
4057  sfluxzen(ns+j) = sfluxref01(j,1,ibd)
4058  enddo
4059 
4060  case (27)
4061 
4062  do j = 1, njb
4063  sfluxzen(ns+j) = scalekur * sfluxref01(j,1,ibd)
4064  enddo
4065 
4066  case default
4067 
4068  if (jb==17 .or. jb==28) then
4069 
4070  ks = nlay
4071  lab_do_k1 : do k = laytrop, nlay-1
4072  if (jp(k)<layreffr(jb) .and. jp(k+1)>=layreffr(jb)) then
4073  ks = k + 1
4074  exit lab_do_k1
4075  endif
4076  enddo lab_do_k1
4077 
4078  colm1 = colamt(ks,ix1(jb))
4079  colm2 = colamt(ks,ix2(jb))
4080  speccomb = colm1 + strrat(jb)*colm2
4081  specmult = specwt(jb) * min( oneminus, colm1/speccomb )
4082  js = 1 + int( specmult )
4083  fs = mod(specmult, f_one)
4084 
4085  do j = 1, njb
4086  sfluxzen(ns+j) = sfluxref02(j,js,ibd) &
4087  & + fs * (sfluxref02(j,js+1,ibd) - sfluxref02(j,js,ibd))
4088  enddo
4089 
4090  else
4091 
4092  ks = laytrop
4093  lab_do_k2 : do k = 1, laytrop-1
4094  if (jp(k)<layreffr(jb) .and. jp(k+1)>=layreffr(jb)) then
4095  ks = k + 1
4096  exit lab_do_k2
4097  endif
4098  enddo lab_do_k2
4099 
4100  colm1 = colamt(ks,ix1(jb))
4101  colm2 = colamt(ks,ix2(jb))
4102  speccomb = colm1 + strrat(jb)*colm2
4103  specmult = specwt(jb) * min( oneminus, colm1/speccomb )
4104  js = 1 + int( specmult )
4105  fs = mod(specmult, f_one)
4106 
4107  do j = 1, njb
4108  sfluxzen(ns+j) = sfluxref03(j,js,ibd) &
4109  & + fs * (sfluxref03(j,js+1,ibd) - sfluxref03(j,js,ibd))
4110  enddo
4111 
4112  endif
4113 
4114  end select
4115 
4116  enddo
4117 
4119 
4121  call taumol16
4123  call taumol17
4125  call taumol18
4127  call taumol19
4129  call taumol20
4131  call taumol21
4133  call taumol22
4135  call taumol23
4137  call taumol24
4139  call taumol25
4141  call taumol26
4143  call taumol27
4145  call taumol28
4147  call taumol29
4148 
4149 
4150 ! =================
4151  contains
4152 ! =================
4153 
4156 !-----------------------------------
4157  subroutine taumol16
4158 !...................................
4159 
4160 ! ------------------------------------------------------------------ !
4161 ! band 16: 2600-3250 cm-1 (low - h2o,ch4; high - ch4) !
4162 ! ------------------------------------------------------------------ !
4163 !
4164  use module_radsw_kgb16
4165 
4166 ! --- locals:
4167 
4168  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4169  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4170 
4171  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4172  integer :: inds, indf, indsp, indfp, j, js, k
4173 
4174 !
4175 !===> ... begin here
4176 !
4177 
4178 ! --- ... compute the optical depth by interpolating in ln(pressure),
4179 ! temperature, and appropriate species. below laytrop, the water
4180 ! vapor self-continuum is interpolated (in temperature) separately.
4181 
4182  do k = 1, nlay
4183  tauray = colmol(k) * rayl
4184 
4185  do j = 1, ng16
4186  taur(k,ns16+j) = tauray
4187  enddo
4188  enddo
4189 
4190  do k = 1, laytrop
4191  speccomb = colamt(k,1) + strrat(16)*colamt(k,5)
4192  specmult = 8.0 * min( oneminus, colamt(k,1)/speccomb )
4193 
4194  js = 1 + int( specmult )
4195  fs = mod( specmult, f_one )
4196  fs1= f_one - fs
4197  fac000 = fs1 * fac00(k)
4198  fac010 = fs1 * fac10(k)
4199  fac100 = fs * fac00(k)
4200  fac110 = fs * fac10(k)
4201  fac001 = fs1 * fac01(k)
4202  fac011 = fs1 * fac11(k)
4203  fac101 = fs * fac01(k)
4204  fac111 = fs * fac11(k)
4205 
4206  ind01 = id0(k,16) + js
4207  ind02 = ind01 + 1
4208  ind03 = ind01 + 9
4209  ind04 = ind01 + 10
4210  ind11 = id1(k,16) + js
4211  ind12 = ind11 + 1
4212  ind13 = ind11 + 9
4213  ind14 = ind11 + 10
4214  inds = indself(k)
4215  indf = indfor(k)
4216  indsp= inds + 1
4217  indfp= indf + 1
4218 
4219  do j = 1, ng16
4220  taug(k,ns16+j) = speccomb &
4221  & *( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4222  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4223  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4224  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4225  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4226  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4227  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4228  & * (forref(indfp,j) - forref(indf,j))))
4229  enddo
4230  enddo
4231 
4232  do k = laytrop+1, nlay
4233  ind01 = id0(k,16) + 1
4234  ind02 = ind01 + 1
4235  ind11 = id1(k,16) + 1
4236  ind12 = ind11 + 1
4237 
4238  do j = 1, ng16
4239  taug(k,ns16+j) = colamt(k,5) &
4240  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4241  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
4242  enddo
4243  enddo
4244 
4245  return
4246 !...................................
4247  end subroutine taumol16
4248 !-----------------------------------
4249 
4250 
4253 !-----------------------------------
4254  subroutine taumol17
4255 !...................................
4256 
4257 ! ------------------------------------------------------------------ !
4258 ! band 17: 3250-4000 cm-1 (low - h2o,co2; high - h2o,co2) !
4259 ! ------------------------------------------------------------------ !
4260 !
4261  use module_radsw_kgb17
4262 
4263 ! --- locals:
4264  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4265  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4266 
4267  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4268  integer :: inds, indf, indsp, indfp, j, js, k
4269 
4270 !
4271 !===> ... begin here
4272 !
4273 
4274 ! --- ... compute the optical depth by interpolating in ln(pressure),
4275 ! temperature, and appropriate species. below laytrop, the water
4276 ! vapor self-continuum is interpolated (in temperature) separately.
4277 
4278  do k = 1, nlay
4279  tauray = colmol(k) * rayl
4280 
4281  do j = 1, ng17
4282  taur(k,ns17+j) = tauray
4283  enddo
4284  enddo
4285 
4286  do k = 1, laytrop
4287  speccomb = colamt(k,1) + strrat(17)*colamt(k,2)
4288  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4289 
4290  js = 1 + int(specmult)
4291  fs = mod(specmult, f_one)
4292  fs1= f_one - fs
4293  fac000 = fs1 * fac00(k)
4294  fac010 = fs1 * fac10(k)
4295  fac100 = fs * fac00(k)
4296  fac110 = fs * fac10(k)
4297  fac001 = fs1 * fac01(k)
4298  fac011 = fs1 * fac11(k)
4299  fac101 = fs * fac01(k)
4300  fac111 = fs * fac11(k)
4301 
4302  ind01 = id0(k,17) + js
4303  ind02 = ind01 + 1
4304  ind03 = ind01 + 9
4305  ind04 = ind01 + 10
4306  ind11 = id1(k,17) + js
4307  ind12 = ind11 + 1
4308  ind13 = ind11 + 9
4309  ind14 = ind11 + 10
4310 
4311  inds = indself(k)
4312  indf = indfor(k)
4313  indsp= inds + 1
4314  indfp= indf + 1
4315 
4316  do j = 1, ng17
4317  taug(k,ns17+j) = speccomb &
4318  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4319  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4320  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4321  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4322  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4323  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4324  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4325  & * (forref(indfp,j) - forref(indf,j))))
4326  enddo
4327  enddo
4328 
4329  do k = laytrop+1, nlay
4330  speccomb = colamt(k,1) + strrat(17)*colamt(k,2)
4331  specmult = 4.0 * min(oneminus, colamt(k,1) / speccomb)
4332 
4333  js = 1 + int(specmult)
4334  fs = mod(specmult, f_one)
4335  fs1= f_one - fs
4336  fac000 = fs1 * fac00(k)
4337  fac010 = fs1 * fac10(k)
4338  fac100 = fs * fac00(k)
4339  fac110 = fs * fac10(k)
4340  fac001 = fs1 * fac01(k)
4341  fac011 = fs1 * fac11(k)
4342  fac101 = fs * fac01(k)
4343  fac111 = fs * fac11(k)
4344 
4345  ind01 = id0(k,17) + js
4346  ind02 = ind01 + 1
4347  ind03 = ind01 + 5
4348  ind04 = ind01 + 6
4349  ind11 = id1(k,17) + js
4350  ind12 = ind11 + 1
4351  ind13 = ind11 + 5
4352  ind14 = ind11 + 6
4353 
4354  indf = indfor(k)
4355  indfp= indf + 1
4356 
4357  do j = 1, ng17
4358  taug(k,ns17+j) = speccomb &
4359  & * ( fac000 * absb(ind01,j) + fac100 * absb(ind02,j) &
4360  & + fac010 * absb(ind03,j) + fac110 * absb(ind04,j) &
4361  & + fac001 * absb(ind11,j) + fac101 * absb(ind12,j) &
4362  & + fac011 * absb(ind13,j) + fac111 * absb(ind14,j) ) &
4363  & + colamt(k,1) * forfac(k) * (forref(indf,j) &
4364  & + forfrac(k) * (forref(indfp,j) - forref(indf,j)))
4365  enddo
4366  enddo
4367 
4368  return
4369 !...................................
4370  end subroutine taumol17
4371 !-----------------------------------
4372 
4373 
4376 !-----------------------------------
4377  subroutine taumol18
4378 !...................................
4379 
4380 ! ------------------------------------------------------------------ !
4381 ! band 18: 4000-4650 cm-1 (low - h2o,ch4; high - ch4) !
4382 ! ------------------------------------------------------------------ !
4383 !
4384  use module_radsw_kgb18
4385 
4386 ! --- locals:
4387  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4388  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4389 
4390  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4391  integer :: inds, indf, indsp, indfp, j, js, k
4392 
4393 !
4394 !===> ... begin here
4395 !
4396 
4397 ! --- ... compute the optical depth by interpolating in ln(pressure),
4398 ! temperature, and appropriate species. below laytrop, the water
4399 ! vapor self-continuum is interpolated (in temperature) separately.
4400 
4401  do k = 1, nlay
4402  tauray = colmol(k) * rayl
4403 
4404  do j = 1, ng18
4405  taur(k,ns18+j) = tauray
4406  enddo
4407  enddo
4408 
4409  do k = 1, laytrop
4410  speccomb = colamt(k,1) + strrat(18)*colamt(k,5)
4411  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4412 
4413  js = 1 + int(specmult)
4414  fs = mod(specmult, f_one)
4415  fs1= f_one - fs
4416  fac000 = fs1 * fac00(k)
4417  fac010 = fs1 * fac10(k)
4418  fac100 = fs * fac00(k)
4419  fac110 = fs * fac10(k)
4420  fac001 = fs1 * fac01(k)
4421  fac011 = fs1 * fac11(k)
4422  fac101 = fs * fac01(k)
4423  fac111 = fs * fac11(k)
4424 
4425  ind01 = id0(k,18) + js
4426  ind02 = ind01 + 1
4427  ind03 = ind01 + 9
4428  ind04 = ind01 + 10
4429  ind11 = id1(k,18) + js
4430  ind12 = ind11 + 1
4431  ind13 = ind11 + 9
4432  ind14 = ind11 + 10
4433 
4434  inds = indself(k)
4435  indf = indfor(k)
4436  indsp= inds + 1
4437  indfp= indf + 1
4438 
4439  do j = 1, ng18
4440  taug(k,ns18+j) = speccomb &
4441  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4442  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4443  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4444  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4445  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4446  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4447  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4448  & * (forref(indfp,j) - forref(indf,j))))
4449  enddo
4450  enddo
4451 
4452  do k = laytrop+1, nlay
4453  ind01 = id0(k,18) + 1
4454  ind02 = ind01 + 1
4455  ind11 = id1(k,18) + 1
4456  ind12 = ind11 + 1
4457 
4458  do j = 1, ng18
4459  taug(k,ns18+j) = colamt(k,5) &
4460  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4461  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
4462  enddo
4463  enddo
4464 
4465  return
4466 !...................................
4467  end subroutine taumol18
4468 !-----------------------------------
4469 
4472 !-----------------------------------
4473  subroutine taumol19
4474 !...................................
4475 
4476 ! ------------------------------------------------------------------ !
4477 ! band 19: 4650-5150 cm-1 (low - h2o,co2; high - co2) !
4478 ! ------------------------------------------------------------------ !
4479 !
4480  use module_radsw_kgb19
4481 
4482 ! --- locals:
4483  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4484  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4485 
4486  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4487  integer :: inds, indf, indsp, indfp, j, js, k
4488 
4489 !
4490 !===> ... begin here
4491 !
4492 
4493 ! --- ... compute the optical depth by interpolating in ln(pressure),
4494 ! temperature, and appropriate species. below laytrop, the water
4495 ! vapor self-continuum is interpolated (in temperature) separately.
4496 
4497  do k = 1, nlay
4498  tauray = colmol(k) * rayl
4499 
4500  do j = 1, ng19
4501  taur(k,ns19+j) = tauray
4502  enddo
4503  enddo
4504 
4505  do k = 1, laytrop
4506  speccomb = colamt(k,1) + strrat(19)*colamt(k,2)
4507  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4508 
4509  js = 1 + int(specmult)
4510  fs = mod(specmult, f_one)
4511  fs1= f_one - fs
4512  fac000 = fs1 * fac00(k)
4513  fac010 = fs1 * fac10(k)
4514  fac100 = fs * fac00(k)
4515  fac110 = fs * fac10(k)
4516  fac001 = fs1 * fac01(k)
4517  fac011 = fs1 * fac11(k)
4518  fac101 = fs * fac01(k)
4519  fac111 = fs * fac11(k)
4520 
4521  ind01 = id0(k,19) + js
4522  ind02 = ind01 + 1
4523  ind03 = ind01 + 9
4524  ind04 = ind01 + 10
4525  ind11 = id1(k,19) + js
4526  ind12 = ind11 + 1
4527  ind13 = ind11 + 9
4528  ind14 = ind11 + 10
4529 
4530  inds = indself(k)
4531  indf = indfor(k)
4532  indsp= inds + 1
4533  indfp= indf + 1
4534 
4535  do j = 1, ng19
4536  taug(k,ns19+j) = speccomb &
4537  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4538  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4539  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4540  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4541  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4542  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4543  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4544  & * (forref(indfp,j) - forref(indf,j))))
4545  enddo
4546  enddo
4547 
4548  do k = laytrop+1, nlay
4549  ind01 = id0(k,19) + 1
4550  ind02 = ind01 + 1
4551  ind11 = id1(k,19) + 1
4552  ind12 = ind11 + 1
4553 
4554  do j = 1, ng19
4555  taug(k,ns19+j) = colamt(k,2) &
4556  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4557  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
4558  enddo
4559  enddo
4560 
4561 !...................................
4562  end subroutine taumol19
4563 !-----------------------------------
4564 
4565 
4568 !-----------------------------------
4569  subroutine taumol20
4570 !...................................
4571 
4572 ! ------------------------------------------------------------------ !
4573 ! band 20: 5150-6150 cm-1 (low - h2o; high - h2o) !
4574 ! ------------------------------------------------------------------ !
4575 !
4576  use module_radsw_kgb20
4577 
4578 ! --- locals:
4579  real (kind=kind_phys) :: tauray
4580 
4581  integer :: ind01, ind02, ind11, ind12
4582  integer :: inds, indf, indsp, indfp, j, k
4583 
4584 !
4585 !===> ... begin here
4586 !
4587 
4588 ! --- ... compute the optical depth by interpolating in ln(pressure),
4589 ! temperature, and appropriate species. below laytrop, the water
4590 ! vapor self-continuum is interpolated (in temperature) separately.
4591 
4592  do k = 1, nlay
4593  tauray = colmol(k) * rayl
4594 
4595  do j = 1, ng20
4596  taur(k,ns20+j) = tauray
4597  enddo
4598  enddo
4599 
4600  do k = 1, laytrop
4601  ind01 = id0(k,20) + 1
4602  ind02 = ind01 + 1
4603  ind11 = id1(k,20) + 1
4604  ind12 = ind11 + 1
4605 
4606  inds = indself(k)
4607  indf = indfor(k)
4608  indsp= inds + 1
4609  indfp= indf + 1
4610 
4611  do j = 1, ng20
4612  taug(k,ns20+j) = colamt(k,1) &
4613  & * ( (fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
4614  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j)) &
4615  & + selffac(k) * (selfref(inds,j) + selffrac(k) &
4616  & * (selfref(indsp,j) - selfref(inds,j))) &
4617  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4618  & * (forref(indfp,j) - forref(indf,j))) ) &
4619  & + colamt(k,5) * absch4(j)
4620  enddo
4621  enddo
4622 
4623  do k = laytrop+1, nlay
4624  ind01 = id0(k,20) + 1
4625  ind02 = ind01 + 1
4626  ind11 = id1(k,20) + 1
4627  ind12 = ind11 + 1
4628 
4629  indf = indfor(k)
4630  indfp= indf + 1
4631 
4632  do j = 1, ng20
4633  taug(k,ns20+j) = colamt(k,1) &
4634  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4635  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) &
4636  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4637  & * (forref(indfp,j) - forref(indf,j))) ) &
4638  & + colamt(k,5) * absch4(j)
4639  enddo
4640  enddo
4641 
4642  return
4643 !...................................
4644  end subroutine taumol20
4645 !-----------------------------------
4646 
4647 
4650 !-----------------------------------
4651  subroutine taumol21
4652 !...................................
4653 
4654 ! ------------------------------------------------------------------ !
4655 ! band 21: 6150-7700 cm-1 (low - h2o,co2; high - h2o,co2) !
4656 ! ------------------------------------------------------------------ !
4657 !
4658  use module_radsw_kgb21
4659 
4660 ! --- locals:
4661  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4662  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4663 
4664  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4665  integer :: inds, indf, indsp, indfp, j, js, k
4666 
4667 !
4668 !===> ... begin here
4669 !
4670 
4671 ! --- ... compute the optical depth by interpolating in ln(pressure),
4672 ! temperature, and appropriate species. below laytrop, the water
4673 ! vapor self-continuum is interpolated (in temperature) separately.
4674 
4675  do k = 1, nlay
4676  tauray = colmol(k) * rayl
4677 
4678  do j = 1, ng21
4679  taur(k,ns21+j) = tauray
4680  enddo
4681  enddo
4682 
4683  do k = 1, laytrop
4684  speccomb = colamt(k,1) + strrat(21)*colamt(k,2)
4685  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4686 
4687  js = 1 + int(specmult)
4688  fs = mod(specmult, f_one)
4689  fs1= f_one - fs
4690  fac000 = fs1 * fac00(k)
4691  fac010 = fs1 * fac10(k)
4692  fac100 = fs * fac00(k)
4693  fac110 = fs * fac10(k)
4694  fac001 = fs1 * fac01(k)
4695  fac011 = fs1 * fac11(k)
4696  fac101 = fs * fac01(k)
4697  fac111 = fs * fac11(k)
4698 
4699  ind01 = id0(k,21) + js
4700  ind02 = ind01 + 1
4701  ind03 = ind01 + 9
4702  ind04 = ind01 + 10
4703  ind11 = id1(k,21) + js
4704  ind12 = ind11 + 1
4705  ind13 = ind11 + 9
4706  ind14 = ind11 + 10
4707 
4708  inds = indself(k)
4709  indf = indfor(k)
4710  indsp= inds + 1
4711  indfp= indf + 1
4712 
4713  do j = 1, ng21
4714  taug(k,ns21+j) = speccomb &
4715  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4716  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4717  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4718  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4719  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4720  & + selffrac(k) * (selfref(indsp,j) - selfref(inds,j))) &
4721  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4722  & * (forref(indfp,j) - forref(indf,j))))
4723  enddo
4724  enddo
4725 
4726  do k = laytrop+1, nlay
4727  speccomb = colamt(k,1) + strrat(21)*colamt(k,2)
4728  specmult = 4.0 * min(oneminus, colamt(k,1) / speccomb)
4729 
4730  js = 1 + int(specmult)
4731  fs = mod(specmult, f_one)
4732  fs1= f_one - fs
4733  fac000 = fs1 * fac00(k)
4734  fac010 = fs1 * fac10(k)
4735  fac100 = fs * fac00(k)
4736  fac110 = fs * fac10(k)
4737  fac001 = fs1 * fac01(k)
4738  fac011 = fs1 * fac11(k)
4739  fac101 = fs * fac01(k)
4740  fac111 = fs * fac11(k)
4741 
4742  ind01 = id0(k,21) + js
4743  ind02 = ind01 + 1
4744  ind03 = ind01 + 5
4745  ind04 = ind01 + 6
4746  ind11 = id1(k,21) + js
4747  ind12 = ind11 + 1
4748  ind13 = ind11 + 5
4749  ind14 = ind11 + 6
4750 
4751  indf = indfor(k)
4752  indfp= indf + 1
4753 
4754  do j = 1, ng21
4755  taug(k,ns21+j) = speccomb &
4756  & * ( fac000 * absb(ind01,j) + fac100 * absb(ind02,j) &
4757  & + fac010 * absb(ind03,j) + fac110 * absb(ind04,j) &
4758  & + fac001 * absb(ind11,j) + fac101 * absb(ind12,j) &
4759  & + fac011 * absb(ind13,j) + fac111 * absb(ind14,j) ) &
4760  & + colamt(k,1) * forfac(k) * (forref(indf,j) &
4761  & + forfrac(k) * (forref(indfp,j) - forref(indf,j)))
4762  enddo
4763  enddo
4764 
4765 !...................................
4766  end subroutine taumol21
4767 !-----------------------------------
4768 
4769 
4772 !-----------------------------------
4773  subroutine taumol22
4774 !...................................
4775 
4776 ! ------------------------------------------------------------------ !
4777 ! band 22: 7700-8050 cm-1 (low - h2o,o2; high - o2) !
4778 ! ------------------------------------------------------------------ !
4779 !
4780  use module_radsw_kgb22
4781 
4782 ! --- locals:
4783  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
4784  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111, &
4785  & o2adj, o2cont, o2tem
4786 
4787  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4788  integer :: inds, indf, indsp, indfp, j, js, k
4789 
4790 !
4791 !===> ... begin here
4792 !
4793 ! --- ... the following factor is the ratio of total o2 band intensity (lines
4794 ! and mate continuum) to o2 band intensity (line only). it is needed
4795 ! to adjust the optical depths since the k's include only lines.
4796 
4797  o2adj = 1.6
4798  o2tem = 4.35e-4 / (350.0*2.0)
4799 
4800 
4801 ! --- ... compute the optical depth by interpolating in ln(pressure),
4802 ! temperature, and appropriate species. below laytrop, the water
4803 ! vapor self-continuum is interpolated (in temperature) separately.
4804 
4805  do k = 1, nlay
4806  tauray = colmol(k) * rayl
4807 
4808  do j = 1, ng22
4809  taur(k,ns22+j) = tauray
4810  enddo
4811  enddo
4812 
4813  do k = 1, laytrop
4814  o2cont = o2tem * colamt(k,6)
4815  speccomb = colamt(k,1) + strrat(22)*colamt(k,6)
4816  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4817 
4818  js = 1 + int(specmult)
4819  fs = mod(specmult, f_one)
4820  fs1= f_one - fs
4821  fac000 = fs1 * fac00(k)
4822  fac010 = fs1 * fac10(k)
4823  fac100 = fs * fac00(k)
4824  fac110 = fs * fac10(k)
4825  fac001 = fs1 * fac01(k)
4826  fac011 = fs1 * fac11(k)
4827  fac101 = fs * fac01(k)
4828  fac111 = fs * fac11(k)
4829 
4830  ind01 = id0(k,22) + js
4831  ind02 = ind01 + 1
4832  ind03 = ind01 + 9
4833  ind04 = ind01 + 10
4834  ind11 = id1(k,22) + js
4835  ind12 = ind11 + 1
4836  ind13 = ind11 + 9
4837  ind14 = ind11 + 10
4838 
4839  inds = indself(k)
4840  indf = indfor(k)
4841  indsp= inds + 1
4842  indfp= indf + 1
4843 
4844  do j = 1, ng22
4845  taug(k,ns22+j) = speccomb &
4846  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
4847  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
4848  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
4849  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
4850  & + colamt(k,1) * (selffac(k) * (selfref(inds,j) &
4851  & + selffrac(k) * (selfref(indsp,j)-selfref(inds,j))) &
4852  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4853  & * (forref(indfp,j) - forref(indf,j)))) + o2cont
4854  enddo
4855  enddo
4856 
4857  do k = laytrop+1, nlay
4858  o2cont = o2tem * colamt(k,6)
4859 
4860  ind01 = id0(k,22) + 1
4861  ind02 = ind01 + 1
4862  ind11 = id1(k,22) + 1
4863  ind12 = ind11 + 1
4864 
4865  do j = 1, ng22
4866  taug(k,ns22+j) = colamt(k,6) * o2adj &
4867  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
4868  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) &
4869  & + o2cont
4870  enddo
4871  enddo
4872 
4873  return
4874 !...................................
4875  end subroutine taumol22
4876 !-----------------------------------
4877 
4878 
4881 !-----------------------------------
4882  subroutine taumol23
4883 !...................................
4884 
4885 ! ------------------------------------------------------------------ !
4886 ! band 23: 8050-12850 cm-1 (low - h2o; high - nothing) !
4887 ! ------------------------------------------------------------------ !
4888 !
4889  use module_radsw_kgb23
4890 
4891 ! --- locals:
4892  integer :: ind01, ind02, ind11, ind12
4893  integer :: inds, indf, indsp, indfp, j, k
4894 
4895 !
4896 !===> ... begin here
4897 !
4898 
4899 ! --- ... compute the optical depth by interpolating in ln(pressure),
4900 ! temperature, and appropriate species. below laytrop, the water
4901 ! vapor self-continuum is interpolated (in temperature) separately.
4902 
4903  do k = 1, nlay
4904  do j = 1, ng23
4905  taur(k,ns23+j) = colmol(k) * rayl(j)
4906  enddo
4907  enddo
4908 
4909  do k = 1, laytrop
4910  ind01 = id0(k,23) + 1
4911  ind02 = ind01 + 1
4912  ind11 = id1(k,23) + 1
4913  ind12 = ind11 + 1
4914 
4915  inds = indself(k)
4916  indf = indfor(k)
4917  indsp= inds + 1
4918  indfp= indf + 1
4919 
4920  do j = 1, ng23
4921  taug(k,ns23+j) = colamt(k,1) * (givfac &
4922  & * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
4923  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) &
4924  & + selffac(k) * (selfref(inds,j) + selffrac(k) &
4925  & * (selfref(indsp,j) - selfref(inds,j))) &
4926  & + forfac(k) * (forref(indf,j) + forfrac(k) &
4927  & * (forref(indfp,j) - forref(indf,j))))
4928  enddo
4929  enddo
4930 
4931  do k = laytrop+1, nlay
4932  do j = 1, ng23
4933  taug(k,ns23+j) = f_zero
4934  enddo
4935  enddo
4936 
4937 !...................................
4938  end subroutine taumol23
4939 !-----------------------------------
4940 
4941 
4944 !-----------------------------------
4945  subroutine taumol24
4946 !...................................
4947 
4948 ! ------------------------------------------------------------------ !
4949 ! band 24: 12850-16000 cm-1 (low - h2o,o2; high - o2) !
4950 ! ------------------------------------------------------------------ !
4951 !
4952  use module_radsw_kgb24
4953 
4954 ! --- locals:
4955  real (kind=kind_phys) :: speccomb, specmult, fs, fs1, &
4956  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
4957 
4958  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
4959  integer :: inds, indf, indsp, indfp, j, js, k
4960 
4961 !
4962 !===> ... begin here
4963 !
4964 
4965 ! --- ... compute the optical depth by interpolating in ln(pressure),
4966 ! temperature, and appropriate species. below laytrop, the water
4967 ! vapor self-continuum is interpolated (in temperature) separately.
4968 
4969  do k = 1, laytrop
4970  speccomb = colamt(k,1) + strrat(24)*colamt(k,6)
4971  specmult = 8.0 * min(oneminus, colamt(k,1) / speccomb)
4972 
4973  js = 1 + int(specmult)
4974  fs = mod(specmult, f_one)
4975  fs1= f_one - fs
4976  fac000 = fs1 * fac00(k)
4977  fac010 = fs1 * fac10(k)
4978  fac100 = fs * fac00(k)
4979  fac110 = fs * fac10(k)
4980  fac001 = fs1 * fac01(k)
4981  fac011 = fs1 * fac11(k)
4982  fac101 = fs * fac01(k)
4983  fac111 = fs * fac11(k)
4984 
4985  ind01 = id0(k,24) + js
4986  ind02 = ind01 + 1
4987  ind03 = ind01 + 9
4988  ind04 = ind01 + 10
4989  ind11 = id1(k,24) + js
4990  ind12 = ind11 + 1
4991  ind13 = ind11 + 9
4992  ind14 = ind11 + 10
4993 
4994  inds = indself(k)
4995  indf = indfor(k)
4996  indsp= inds + 1
4997  indfp= indf + 1
4998 
4999  do j = 1, ng24
5000  taug(k,ns24+j) = speccomb &
5001  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
5002  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
5003  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
5004  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) ) &
5005  & + colamt(k,3) * abso3a(j) + colamt(k,1) &
5006  & * (selffac(k) * (selfref(inds,j) + selffrac(k) &
5007  & * (selfref(indsp,j) - selfref(inds,j))) &
5008  & + forfac(k) * (forref(indf,j) + forfrac(k) &
5009  & * (forref(indfp,j) - forref(indf,j))))
5010 
5011  taur(k,ns24+j) = colmol(k) &
5012  & * (rayla(j,js) + fs*(rayla(j,js+1) - rayla(j,js)))
5013  enddo
5014  enddo
5015 
5016  do k = laytrop+1, nlay
5017  ind01 = id0(k,24) + 1
5018  ind02 = ind01 + 1
5019  ind11 = id1(k,24) + 1
5020  ind12 = ind11 + 1
5021 
5022  do j = 1, ng24
5023  taug(k,ns24+j) = colamt(k,6) &
5024  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
5025  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) &
5026  & + colamt(k,3) * abso3b(j)
5027 
5028  taur(k,ns24+j) = colmol(k) * raylb(j)
5029  enddo
5030  enddo
5031 
5032  return
5033 !...................................
5034  end subroutine taumol24
5035 !-----------------------------------
5036 
5037 
5040 !-----------------------------------
5041  subroutine taumol25
5042 !...................................
5043 
5044 ! ------------------------------------------------------------------ !
5045 ! band 25: 16000-22650 cm-1 (low - h2o; high - nothing) !
5046 ! ------------------------------------------------------------------ !
5047 !
5048  use module_radsw_kgb25
5049 
5050 ! --- locals:
5051  integer :: ind01, ind02, ind11, ind12
5052  integer :: j, k
5053 
5054 !
5055 !===> ... begin here
5056 !
5057 
5058 ! --- ... compute the optical depth by interpolating in ln(pressure),
5059 ! temperature, and appropriate species. below laytrop, the water
5060 ! vapor self-continuum is interpolated (in temperature) separately.
5061 
5062  do k = 1, nlay
5063  do j = 1, ng25
5064  taur(k,ns25+j) = colmol(k) * rayl(j)
5065  enddo
5066  enddo
5067 
5068  do k = 1, laytrop
5069  ind01 = id0(k,25) + 1
5070  ind02 = ind01 + 1
5071  ind11 = id1(k,25) + 1
5072  ind12 = ind11 + 1
5073 
5074  do j = 1, ng25
5075  taug(k,ns25+j) = colamt(k,1) &
5076  & * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
5077  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) &
5078  & + colamt(k,3) * abso3a(j)
5079  enddo
5080  enddo
5081 
5082  do k = laytrop+1, nlay
5083  do j = 1, ng25
5084  taug(k,ns25+j) = colamt(k,3) * abso3b(j)
5085  enddo
5086  enddo
5087 
5088  return
5089 !...................................
5090  end subroutine taumol25
5091 !-----------------------------------
5092 
5093 
5096 !-----------------------------------
5097  subroutine taumol26
5098 !...................................
5099 
5100 ! ------------------------------------------------------------------ !
5101 ! band 26: 22650-29000 cm-1 (low - nothing; high - nothing) !
5102 ! ------------------------------------------------------------------ !
5103 !
5104  use module_radsw_kgb26
5105 
5106 ! --- locals:
5107  integer :: j, k
5108 
5109 !
5110 !===> ... begin here
5111 !
5112 
5113 ! --- ... compute the optical depth by interpolating in ln(pressure),
5114 ! temperature, and appropriate species. below laytrop, the water
5115 ! vapor self-continuum is interpolated (in temperature) separately.
5116 
5117  do k = 1, nlay
5118  do j = 1, ng26
5119  taug(k,ns26+j) = f_zero
5120  taur(k,ns26+j) = colmol(k) * rayl(j)
5121  enddo
5122  enddo
5123 
5124  return
5125 !...................................
5126  end subroutine taumol26
5127 !-----------------------------------
5128 
5129 
5132 !-----------------------------------
5133  subroutine taumol27
5134 !...................................
5135 
5136 ! ------------------------------------------------------------------ !
5137 ! band 27: 29000-38000 cm-1 (low - o3; high - o3) !
5138 ! ------------------------------------------------------------------ !
5139 !
5140  use module_radsw_kgb27
5141 !
5142 ! --- locals:
5143  integer :: ind01, ind02, ind11, ind12
5144  integer :: j, k
5145 
5146 !
5147 !===> ... begin here
5148 !
5149 
5150 ! --- ... compute the optical depth by interpolating in ln(pressure),
5151 ! temperature, and appropriate species. below laytrop, the water
5152 ! vapor self-continuum is interpolated (in temperature) separately.
5153 
5154  do k = 1, nlay
5155  do j = 1, ng27
5156  taur(k,ns27+j) = colmol(k) * rayl(j)
5157  enddo
5158  enddo
5159 
5160  do k = 1, laytrop
5161  ind01 = id0(k,27) + 1
5162  ind02 = ind01 + 1
5163  ind11 = id1(k,27) + 1
5164  ind12 = ind11 + 1
5165 
5166  do j = 1, ng27
5167  taug(k,ns27+j) = colamt(k,3) &
5168  & * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
5169  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) )
5170  enddo
5171  enddo
5172 
5173  do k = laytrop+1, nlay
5174  ind01 = id0(k,27) + 1
5175  ind02 = ind01 + 1
5176  ind11 = id1(k,27) + 1
5177  ind12 = ind11 + 1
5178 
5179  do j = 1, ng27
5180  taug(k,ns27+j) = colamt(k,3) &
5181  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
5182  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
5183  enddo
5184  enddo
5185 
5186  return
5187 !...................................
5188  end subroutine taumol27
5189 !-----------------------------------
5190 
5191 
5194 !-----------------------------------
5195  subroutine taumol28
5196 !...................................
5197 
5198 ! ------------------------------------------------------------------ !
5199 ! band 28: 38000-50000 cm-1 (low - o3,o2; high - o3,o2) !
5200 ! ------------------------------------------------------------------ !
5201 !
5202  use module_radsw_kgb28
5203 
5204 ! --- locals:
5205  real (kind=kind_phys) :: speccomb, specmult, tauray, fs, fs1, &
5206  & fac000,fac001,fac010,fac011, fac100,fac101,fac110,fac111
5207 
5208  integer :: ind01, ind02, ind03, ind04, ind11, ind12, ind13, ind14
5209  integer :: j, js, k
5210 
5211 !
5212 !===> ... begin here
5213 !
5214 
5215 ! --- ... compute the optical depth by interpolating in ln(pressure),
5216 ! temperature, and appropriate species. below laytrop, the water
5217 ! vapor self-continuum is interpolated (in temperature) separately.
5218 
5219  do k = 1, nlay
5220  tauray = colmol(k) * rayl
5221 
5222  do j = 1, ng28
5223  taur(k,ns28+j) = tauray
5224  enddo
5225  enddo
5226 
5227  do k = 1, laytrop
5228  speccomb = colamt(k,3) + strrat(28)*colamt(k,6)
5229  specmult = 8.0 * min(oneminus, colamt(k,3) / speccomb)
5230 
5231  js = 1 + int(specmult)
5232  fs = mod(specmult, f_one)
5233  fs1= f_one - fs
5234  fac000 = fs1 * fac00(k)
5235  fac010 = fs1 * fac10(k)
5236  fac100 = fs * fac00(k)
5237  fac110 = fs * fac10(k)
5238  fac001 = fs1 * fac01(k)
5239  fac011 = fs1 * fac11(k)
5240  fac101 = fs * fac01(k)
5241  fac111 = fs * fac11(k)
5242 
5243  ind01 = id0(k,28) + js
5244  ind02 = ind01 + 1
5245  ind03 = ind01 + 9
5246  ind04 = ind01 + 10
5247  ind11 = id1(k,28) + js
5248  ind12 = ind11 + 1
5249  ind13 = ind11 + 9
5250  ind14 = ind11 + 10
5251 
5252  do j = 1, ng28
5253  taug(k,ns28+j) = speccomb &
5254  & * ( fac000 * absa(ind01,j) + fac100 * absa(ind02,j) &
5255  & + fac010 * absa(ind03,j) + fac110 * absa(ind04,j) &
5256  & + fac001 * absa(ind11,j) + fac101 * absa(ind12,j) &
5257  & + fac011 * absa(ind13,j) + fac111 * absa(ind14,j) )
5258  enddo
5259  enddo
5260 
5261  do k = laytrop+1, nlay
5262  speccomb = colamt(k,3) + strrat(28)*colamt(k,6)
5263  specmult = 4.0 * min(oneminus, colamt(k,3) / speccomb)
5264 
5265  js = 1 + int(specmult)
5266  fs = mod(specmult, f_one)
5267  fs1= f_one - fs
5268  fac000 = fs1 * fac00(k)
5269  fac010 = fs1 * fac10(k)
5270  fac100 = fs * fac00(k)
5271  fac110 = fs * fac10(k)
5272  fac001 = fs1 * fac01(k)
5273  fac011 = fs1 * fac11(k)
5274  fac101 = fs * fac01(k)
5275  fac111 = fs * fac11(k)
5276 
5277  ind01 = id0(k,28) + js
5278  ind02 = ind01 + 1
5279  ind03 = ind01 + 5
5280  ind04 = ind01 + 6
5281  ind11 = id1(k,28) + js
5282  ind12 = ind11 + 1
5283  ind13 = ind11 + 5
5284  ind14 = ind11 + 6
5285 
5286  do j = 1, ng28
5287  taug(k,ns28+j) = speccomb &
5288  & * ( fac000 * absb(ind01,j) + fac100 * absb(ind02,j) &
5289  & + fac010 * absb(ind03,j) + fac110 * absb(ind04,j) &
5290  & + fac001 * absb(ind11,j) + fac101 * absb(ind12,j) &
5291  & + fac011 * absb(ind13,j) + fac111 * absb(ind14,j) )
5292  enddo
5293  enddo
5294 
5295  return
5296 !...................................
5297  end subroutine taumol28
5298 !-----------------------------------
5299 
5300 
5303 !-----------------------------------
5304  subroutine taumol29
5305 !...................................
5306 
5307 ! ------------------------------------------------------------------ !
5308 ! band 29: 820-2600 cm-1 (low - h2o; high - co2) !
5309 ! ------------------------------------------------------------------ !
5310 !
5311  use module_radsw_kgb29
5312 
5313 ! --- locals:
5314  real (kind=kind_phys) :: tauray
5315 
5316  integer :: ind01, ind02, ind11, ind12
5317  integer :: inds, indf, indsp, indfp, j, k
5318 
5319 !
5320 !===> ... begin here
5321 !
5322 
5323 ! --- ... compute the optical depth by interpolating in ln(pressure),
5324 ! temperature, and appropriate species. below laytrop, the water
5325 ! vapor self-continuum is interpolated (in temperature) separately.
5326 
5327  do k = 1, nlay
5328  tauray = colmol(k) * rayl
5329 
5330  do j = 1, ng29
5331  taur(k,ns29+j) = tauray
5332  enddo
5333  enddo
5334 
5335  do k = 1, laytrop
5336  ind01 = id0(k,29) + 1
5337  ind02 = ind01 + 1
5338  ind11 = id1(k,29) + 1
5339  ind12 = ind11 + 1
5340 
5341  inds = indself(k)
5342  indf = indfor(k)
5343  indsp= inds + 1
5344  indfp= indf + 1
5345 
5346  do j = 1, ng29
5347  taug(k,ns29+j) = colamt(k,1) &
5348  & * ( (fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) &
5349  & + fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) &
5350  & + selffac(k) * (selfref(inds,j) + selffrac(k) &
5351  & * (selfref(indsp,j) - selfref(inds,j))) &
5352  & + forfac(k) * (forref(indf,j) + forfrac(k) &
5353  & * (forref(indfp,j) - forref(indf,j)))) &
5354  & + colamt(k,2) * absco2(j)
5355  enddo
5356  enddo
5357 
5358  do k = laytrop+1, nlay
5359  ind01 = id0(k,29) + 1
5360  ind02 = ind01 + 1
5361  ind11 = id1(k,29) + 1
5362  ind12 = ind11 + 1
5363 
5364  do j = 1, ng29
5365  taug(k,ns29+j) = colamt(k,2) &
5366  & * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) &
5367  & + fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) &
5368  & + colamt(k,1) * absh2o(j)
5369  enddo
5370  enddo
5371 
5372  return
5373 !...................................
5374  end subroutine taumol29
5375 !-----------------------------------
5376 
5377 !...................................
5378  end subroutine taumol
5379 !-----------------------------------
5380 !! @}
5381 
5382 !
5383 !........................................!
5384  end module module_radsw_main !
5385 !========================================!
5386 !! @}
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...
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) ...
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(5), public ebari
asymmetry coefficients
Definition: radsw_datatb.f:278
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:4255
real(kind=kind_phys) heatfac
the factor for heating rates (in k/day, or k/sec set by subroutine &#39;rswinit&#39;)
Definition: radsw_main.f:500
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...
integer, dimension(nblow:nbhgh), public ix2
indexes for 2nd entries of the two key species for each of the spectral bands
real(kind=kind_phys), dimension(0:ntbmx) exp_tbl
those data will be set up only once by "rswinit"
Definition: radsw_main.f:495
integer, dimension(nblow:nbhgh), public ix1
indexes for 1st entries of the two key species for each of the spectral bands
derived type for special components of surface SW fluxes
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...
real(kind=kind_phys), dimension(nblow:nbhgh), public c0r
asymmetry coefficients
subroutine vrtqdr(zrefb, zrefd, ztrab, ztrad, zldbt, ztdbt, NLAY, NLP1, zfu, zfd)
This subroutine is called by spcvrtc() and spcvrtm(), and computes the upward and downward radiation ...
Definition: radsw_main.f:3770
This module sets up absorption coeffients for band 27: 29000-38000 cm-1 (low - o3; high - o3) ...
real(kind=kind_phys), dimension(43, nblow:nbhgh), public asyice2
asymmetry coefficients
Definition: radsw_datatb.f:249
real(kind=kind_phys), dimension(43, nblow:nbhgh), public extice2
extinction coefficients
Definition: radsw_datatb.f:243
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 index
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 optical depth and transmittance tables
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...
real(kind=kind_phys), public a0s
optical depth coefficients
real(kind=kind_phys), dimension(ngmax, mfs01, mfb01), target, public sfluxref01
spectral solar fluxes, j=1,2,...,7 for SW band number of 16,20,23,25,26,27,29
integer, parameter nuvb
uv-b band index
Definition: radsw_main.f:485
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
real(kind=kind_phys), dimension(58, nblow:nbhgh), public extliq1
extinction coefficients
Definition: radsw_datatb.f:232
This module contains cloud radiative property coefficients.
Definition: radsw_datatb.f:176
real(kind=kind_phys), dimension(nblow:nbhgh), public c0s
asymmetry coefficients
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...
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 ...
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(ngmax, mfs02, mfb02), target, public sfluxref02
spectral solar fluxes, j=1,2 for SW band number of 17 and 28
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 ...
subroutine, public swrad(plyr, plvl, tlyr, tlvl, qlyr, olyr, gasvmr, clouds, icseed, aerosols, sfcalb, cosz, solcon, NDAY, idxday, npts, nlay, nlp1, lprnt, hswc, topflx, sfcflx, HSW0, HSWB, FLXPRF, FDNCMP)
This subroutine is the main SW radiation routine.
Definition: radsw_main.f:594
subroutine spcvrtc(ssolar, cosz, sntz, albbm, albdf, sfluxzen, cldfrc, cf1, cf0, taug, taur, tauae, ssaae, asyae, taucw, ssacw, asycw, nlay, nlp1, fxupc, fxdnc, fxup0, fxdn0, ftoauc, ftoau0, ftoadc, fsfcuc, fsfcu0, fsfcdc, fsfcd0, sfbmc, sfdfc, sfbm0, sfdf0, suvbfc, suvbf0)
This subroutine computes the shortwave radiative fluxes using two-stream method.
Definition: radsw_main.f:2304
This module sets up absorption coeffients for band 18: 4000-4650 cm-1 (low - h2o,ch4; high - ch4) ...
subroutine taumol22
The subroutine computes the optical depth in band 22: 7700-8050 cm-1 (low - h2o,o2; high - o2) ...
Definition: radsw_main.f:4774
real(kind=kind_phys), dimension(nblow:nbhgh), public b0s
single scattering albedo coefficients
derived type for SW fluxes&#39; column profiles (at layer interfaces)
Definition: radsw_param.f:98
integer, dimension(nblow:nbhgh) idxebc
band index for cld prop
Definition: radsw_main.f:462
subroutine cldprop(cfrac, cliqp, reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, cf1, nlay, ipseed, taucw, ssacw, asycw, cldfrc, cldfmc)
This subroutine computes the cloud optical properties for each cloudy layer and g-point interval...
Definition: radsw_main.f:1568
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), dimension(43, nblow:nbhgh), public ssaice2
single scattering albedo coefficients
Definition: radsw_datatb.f:246
subroutine taumol23
The subroutine computes the optical depth in band 23: 8050-12850 cm-1 (low - h2o; high - nothing) ...
Definition: radsw_main.f:4883
real(kind=kind_phys), dimension(59) preflog
logarithm,ln(p), of a 59-level standard pressure profile
Definition: radsw_datatb.f:83
subroutine taumol26
The subroutine computes the optical depth in band 26: 22650-29000 cm-1 (low - nothing; high - nothing...
Definition: radsw_main.f:5098
real(kind=kind_phys), dimension(ng25), public abso3b
o3
This module contains the reference pressures (in logarithm form) at 59 vertical levels (TOA is omitte...
Definition: radsw_datatb.f:73
integer, dimension(nblow:nbhgh) idxsfc
band index for sfc flux
Definition: radsw_main.f:460
subroutine taumol16
The subroutine computes the optical depth in band 16: 2600-3250 cm-1 (low - h2o,ch4; high - ch4) ...
Definition: radsw_main.f:4158
subroutine, public rswinit(me)
This subroutine initializes non-varying module variables, conversion factors, and look-up tables...
Definition: radsw_main.f:1379
This module is for specifying the band structures and program parameters used by the RRTMG-SW scheme...
Definition: radsw_param.f:66
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 index
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:4474
This module sets up absorption coeffients for band 26: 22650-29000 cm-1 (low - nothing; high - nothin...
subroutine setcoef(pavel, tavel, h2ovmr, nlay, nlp1, laytrop, jp, jt, jt1, fac00, fac01, fac10, fac11, selffac, selffrac, indself, forfac, forfrac, indfor)
This subroutine computes various coefficients needed in radiative transfer calculation.
Definition: radsw_main.f:2111
subroutine taumol28
The subroutine computes the optical depth in band 28: 38000-50000 cm-1 (low - o3,o2; high - o3...
Definition: radsw_main.f:5196
subroutine taumol18
The subroutine computes the optical depth in band 18: 4000-4650 cm-1 (low - h2o,ch4; high - ch4) ...
Definition: radsw_main.f:4378
real(kind=kind_phys), dimension(nblow:nbhgh), public strrat
coefficients for computing binary absorbing species
real(kind=kind_phys), dimension(46, nblow:nbhgh), public extice3
extinction coefficients
Definition: radsw_datatb.f:255
real(kind=kind_phys), dimension(5), public fbari
asymmetry coefficients
Definition: radsw_datatb.f:280
real(kind=kind_phys), dimension(5), public dbari
single scattering albedo coefficients
Definition: radsw_datatb.f:276
subroutine spcvrtm(ssolar, cosz, sntz, albbm, albdf, sfluxzen, cldfmc, cf1, cf0, taug, taur, tauae, ssaae, asyae, taucw, ssacw, asycw, nlay, nlp1, fxupc, fxdnc, fxup0, fxdn0, ftoauc, ftoau0, ftoadc, fsfcuc, fsfcu0, fsfcdc, fsfcd0, sfbmc, sfdfc, sfbm0, sfdf0, suvbfc, suvbf0)
This subroutine computes the shortwave radiative fluxes using two-stream method of h...
Definition: radsw_main.f:3066
This module sets up absorption coeffients for band 29: 820-2600 cm-1 (low - h2o; high - co2) ...
This module sets up absorption coeffients for band 17: 3250-4000 cm-1 (low - h2o,co2; high - h2o...
subroutine taumol29
The subroutine computes the optical depth in band 29: 820-2600 cm-1 (low - h2o; high - co2) ...
Definition: radsw_main.f:5305
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...
real(kind=kind_phys), dimension(58, nblow:nbhgh), public asyliq1
asymmetry coefficients
Definition: radsw_datatb.f:238
real(kind=kind_phys), dimension(46, nblow:nbhgh), public asyice3
asymmetry coefficients
Definition: radsw_datatb.f:261
subroutine taumol27
The subroutine computes the optical depth in band 27: 29000-38000 cm-1 (low - o3; high - o3) ...
Definition: radsw_main.f:5134
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:448
real(kind=kind_phys), dimension(ngmax, mfs03, mfb03), target, public sfluxref03
spectral solar fluxes, j=1,2,...,5 for SW band number of 18,19,21,22,24
integer, parameter ipsdsw0
initial permutation seed used for sub-column cloud scheme
Definition: radsw_main.f:504
This module contains various indexes and coefficients for SW spectral bands, as well as the spectral ...
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), 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:4570
subroutine taumol21
The subroutine computes the optical depth in band 21: 6150-7700 cm-1 (low - h2o,co2; high - h2o...
Definition: radsw_main.f:4652
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(58, nblow:nbhgh), public ssaliq1
single scattering albedo coefficients
Definition: radsw_datatb.f:235
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:443
real(kind=kind_phys), parameter, public rayl
rayleigh extinction coefficient at
integer, dimension(nblow:nbhgh) ngs
array contains values of NS16-NS29
Definition: radsw_param.f:168
integer, dimension(nblow:nbhgh), public ibx
band index (3rd index in array sfluxref described below)
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
subroutine taumol25
The subroutine computes the optical depth in band 25: 16000-22650 cm-1 (low - h2o; high - nothing) ...
Definition: radsw_main.f:5042
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...
real(kind=kind_phys), dimension(46, nblow:nbhgh), public ssaice3
single scattering albedo coefficients
Definition: radsw_datatb.f:258
subroutine mcica_subcol(cldf, nlay, ipseed, lcloudy)
This subroutine computes the sub-colum cloud profile flag array.
Definition: radsw_main.f:1928
real(kind=kind_phys), public a0r
optical depth coefficients
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), dimension(5), public cbari
single scattering albedo coefficients
Definition: radsw_datatb.f:274
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.
derived type for SW fluxes at TOA
Definition: radsw_param.f:76
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
reverse checking of band index for each g-point
Definition: radsw_param.f:173
real(kind=kind_phys), dimension(ng29), public absh2o
h2o
integer, parameter maxgas
maximum number 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...
integer, dimension(nblow:nbhgh), public layreffr
reference pressure level for each of the spectral bands
subroutine taumol(colamt, colmol, fac00, fac01, fac10, fac11, jp, jt, jt1, laytrop, forfac, forfrac, indfor, selffac, selffrac, indself, nlay, sfluxzen, taug, taur)
This subroutine calculates optical depths for gaseous absorption and rayleigh scattering subroutine...
Definition: radsw_main.f:3906
real(kind=kind_phys), dimension(nblow:nbhgh), public b0r
single scattering albedo coefficients
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...
real(kind=kind_phys), dimension(5), public abari
extinction coefficients
Definition: radsw_datatb.f:270
This module sets up absorption coefficients for band 16: 2600-3250 cm-1 (low - h2o, ch4; high - ch4)
derived type for SW fluxes at surface
Definition: radsw_param.f:86
real(kind=kind_phys), dimension(ng24), public abso3b
o3
real(kind=kind_phys), dimension(59) tref
MLS standard atmosphere temperature at the standard pressure levels.
Definition: radsw_datatb.f:85
real(kind=kind_phys), dimension(5), public bbari
extinction coefficients
Definition: radsw_datatb.f:272
real(kind=kind_phys), dimension(nblow:nbhgh), public specwt
weighting parameters for major absorbers in each band
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:4946
real(kind=kind_phys), dimension(ng29), public absco2
co2
real(kind=kind_phys), dimension(ng25), public rayl
rayleigh extinction coefficient
real(kind=kind_phys), parameter, public rayl
rayleigh extinction coefficient at