Line data Source code
1 :
2 : module radsw
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Purpose: Solar radiation calculations.
6 : !
7 : !-----------------------------------------------------------------------
8 : use shr_kind_mod, only: r8 => shr_kind_r8
9 : use ppgrid, only: pcols, pver, pverp
10 : use cam_abortutils, only: endrun
11 : use cam_history, only: outfld
12 : use scamMod, only: single_column,scm_crm_mode,have_asdir, &
13 : asdirobs, have_asdif, asdifobs, have_aldir, &
14 : aldirobs, have_aldif, aldifobs
15 : use cam_logfile, only: iulog
16 : use radconstants, only: nswbands, get_sw_spectral_boundaries, &
17 : idx_sw_diag, indxsl
18 :
19 : implicit none
20 :
21 : private
22 : save
23 :
24 : ! Public methods
25 :
26 : public ::&
27 : radsw_init, &! initialize constants
28 : radcswmx ! driver for solar radiation code
29 :
30 : ! Private module data
31 :
32 : real(r8) :: gravit ! Acceleration of gravity
33 : real(r8) :: rga ! 1./gravit
34 : real(r8) :: sslp ! Standard sea-level pressure
35 :
36 : !===============================================================================
37 : CONTAINS
38 : !===============================================================================
39 :
40 749232 : subroutine radcswmx(lchnk ,ncol , &
41 : E_pint ,E_pmid ,E_h2ommr ,E_o3mmr , &
42 : E_o2mmr ,E_cld ,E_cicewp ,E_cliqwp ,E_rel , &
43 : E_rei ,eccf ,E_coszrs ,solin , &
44 : E_asdir ,E_asdif ,E_aldir ,E_aldif ,nmxrgn , &
45 : pmxrgn ,qrs,qrsc,fsnt ,fsntc ,fsdtoa, fsntoa, &
46 : fsutoa ,fsntoac, fsnirtoa,fsnrtoac,fsnrtoaq,fsns , &
47 : fsnsc ,fsdsc ,fsds ,sols ,soll , &
48 : solsd ,solld , fns ,fcns , &
49 : Nday ,Nnite ,IdxDay ,IdxNite, E_co2mmr, &
50 : E_aer_tau, E_aer_tau_w, E_aer_tau_w_g, E_aer_tau_w_f, tauxcl_out, tauxci_out)
51 : !-----------------------------------------------------------------------
52 : !
53 :
54 : ! Purpose:
55 : ! Solar radiation code
56 : !
57 : ! Method:
58 : ! Basic method is Delta-Eddington as described in:
59 : !
60 : ! Briegleb, Bruce P., 1992: Delta-Eddington
61 : ! Approximation for Solar Radiation in the NCAR Community Climate Model,
62 : ! Journal of Geophysical Research, Vol 97, D7, pp7603-7612).
63 : !
64 : ! Five changes to the basic method described above are:
65 : ! (1) addition of sulfate aerosols (Kiehl and Briegleb, 1993)
66 : ! (2) the distinction between liquid and ice particle clouds
67 : ! (Kiehl et al, 1996);
68 : ! (3) provision for calculating TOA fluxes with spectral response to
69 : ! match Nimbus-7 visible/near-IR radiometers (Collins, 1998);
70 : ! (4) max-random overlap (Collins, 2001)
71 : ! (5) The near-IR absorption by H2O was updated in 2003 by Collins,
72 : ! Lee-Taylor, and Edwards for consistency with the new line data in
73 : ! Hitran 2000 and the H2O continuum version CKD 2.4. Modifications
74 : ! were optimized by reducing RMS errors in heating rates relative
75 : ! to a series of benchmark calculations for the 5 standard AFGL
76 : ! atmospheres. The benchmarks were performed using DISORT2 combined
77 : ! with GENLN3. The near-IR scattering optical depths for Rayleigh
78 : ! scattering were also adjusted, as well as the correction for
79 : ! stratospheric heating by H2O.
80 : !
81 : ! The treatment of maximum-random overlap is described in the
82 : ! comment block "INDEX CALCULATIONS FOR MAX OVERLAP".
83 : !
84 : ! Divides solar spectrum into 19 intervals from 0.2-5.0 micro-meters.
85 : ! solar flux fractions specified for each interval. allows for
86 : ! seasonally and diurnally varying solar input. Includes molecular,
87 : ! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud,
88 : ! and surface absorption. Computes delta-eddington reflections and
89 : ! transmissions assuming homogeneously mixed layers. Adds the layers
90 : ! assuming scattering between layers to be isotropic, and distinguishes
91 : ! direct solar beam from scattered radiation.
92 : !
93 : ! Longitude loops are broken into 1 or 2 sections, so that only daylight
94 : ! (i.e. coszrs > 0) computations are done.
95 : !
96 : ! Note that an extra layer above the model top layer is added.
97 : !
98 : ! cgs units are used.
99 : !
100 : ! Special diagnostic calculation of the clear sky surface and total column
101 : ! absorbed flux is also done for cloud forcing diagnostics.
102 : !
103 : !-----------------------------------------------------------------------
104 : !
105 : ! D. Parks (NEC) 09/11/03
106 : ! Restructuring of routine to support SX vector architecture.
107 : !
108 : ! Possible improvements:
109 : !
110 : ! 1. Look at vectorizing index calculations for maximum overlap.
111 : !
112 : ! 2. Consider making innermost loop in flux computations the number
113 : ! of spectral intervals. Given that NS is fixed at 19, the trade-off
114 : ! will be stride one memory accesses of length 19 versus indirect
115 : ! addressing (list vector - gather/scatter) with potential vector
116 : ! lenghts of the number of day light points. Vectorizing on the number
117 : ! of spectral intervals seems worthwhile for low resolution models (T42),
118 : ! but might be inefficient with higher resolutions.
119 : !
120 : ! 3. Move the linearization of daylight points (compression/expansion) out
121 : ! of radcswmx and into d_p_coupling. This would eliminate the cost of
122 : ! routines CmpDayNite and ExpDayNite.
123 : !
124 : ! 4. Look at expliciting computing all streams in upward propagation of
125 : ! radiation. There would be additional floating point operations in
126 : ! exchange for the elimination of indirect addressing.
127 : !
128 : !-----------------------------------------------------------------------
129 :
130 : use rad_solar_var, only: get_variability
131 : use cmparray_mod, only: CmpDayNite, ExpDayNite
132 : use quicksort, only: quick_sort
133 : use phys_control, only: phys_getopts
134 : use solar_irrad_data, only: sol_tsi, do_spctrl_scaling, ref_tsi
135 : use radconstants, only: frcsol, ph2o, pco2, po2
136 :
137 : !-----------------------Constants for new band (640-700 nm)-------------
138 : real(r8) v_raytau_35
139 : real(r8) v_raytau_64
140 : real(r8) v_abo3_35
141 : real(r8) v_abo3_64
142 : parameter( &
143 : v_raytau_35 = 0.155208_r8, &
144 : v_raytau_64 = 0.0392_r8, &
145 : v_abo3_35 = 2.4058030e+01_r8, &
146 : v_abo3_64 = 2.210e+01_r8 &
147 : )
148 :
149 :
150 : !-------------Parameters for accelerating max-random solution-------------
151 : !
152 : ! The solution time scales like prod(j:1->N) (1 + n_j) where
153 : ! N = number of max-overlap regions (nmxrgn)
154 : ! n_j = number of unique cloud amounts in region j
155 : !
156 : ! Therefore the solution cost can be reduced by decreasing n_j.
157 : ! cldmin reduces n_j by treating cloud amounts < cldmin as clear sky.
158 : ! cldeps reduces n_j by treating cloud amounts identical to log(1/cldeps)
159 : ! decimal places as identical
160 : !
161 : ! areamin reduces the cost by dropping configurations that occupy
162 : ! a surface area < areamin of the model grid box. The surface area
163 : ! for a configuration C(j,k_j), where j is the region number and k_j is the
164 : ! index for a unique cloud amount (in descending order from biggest to
165 : ! smallest clouds) in region j, is
166 : !
167 : ! A = prod(j:1->N) [C(j,k_j) - C(j,k_j+1)]
168 : !
169 : ! where C(j,0) = 1.0 and C(j,n_j+1) = 0.0.
170 : !
171 : ! nconfgmax reduces the cost and improves load balancing by setting an upper
172 : ! bound on the number of cloud configurations in the solution. If the number
173 : ! of configurations exceeds nconfgmax, the nconfgmax configurations with the
174 : ! largest area are retained, and the fluxes are normalized by the total area
175 : ! of these nconfgmax configurations. For the current max/random overlap
176 : ! assumption (see subroutine cldovrlap), 30 levels, and cloud-amount
177 : ! parameterization, the mean and RMS number of configurations are
178 : ! both roughly 5. nconfgmax has been set to the mean+2*RMS number, or 15.
179 : !
180 : ! Minimum cloud amount (as a fraction of the grid-box area) to
181 : ! distinguish from clear sky
182 : !
183 : real(r8) cldmin
184 : parameter (cldmin = 1.0e-80_r8)
185 : !
186 : ! Minimimum horizontal area (as a fraction of the grid-box area) to retain
187 : ! for a unique cloud configuration in the max-random solution
188 : !
189 : real(r8) areamin
190 : parameter (areamin = 0.01_r8)
191 : !
192 : ! Decimal precision of cloud amount (0 -> preserve full resolution;
193 : ! 10^-n -> preserve n digits of cloud amount)
194 : !
195 : real(r8) cldeps
196 : parameter (cldeps = 0.0_r8)
197 : !
198 : ! Maximum number of configurations to include in solution
199 : !
200 : integer nconfgmax
201 : parameter (nconfgmax = 15)
202 : !------------------------------Commons----------------------------------
203 : !
204 : ! Input arguments
205 : !
206 : integer, intent(in) :: lchnk ! chunk identifier
207 : integer, intent(in) :: ncol ! number of atmospheric columns
208 :
209 : integer,intent(in) :: Nday ! Number of daylight columns
210 : integer,intent(in) :: Nnite ! Number of night columns
211 : integer,intent(in), dimension(pcols) :: IdxDay ! Indicies of daylight coumns
212 : integer,intent(in), dimension(pcols) :: IdxNite ! Indicies of night coumns
213 :
214 :
215 : real(r8), intent(in) :: E_pmid(pcols,pver) ! Level pressure
216 : real(r8), intent(in) :: E_pint(pcols,pverp) ! Interface pressure
217 : real(r8), intent(in) :: E_h2ommr(pcols,pver) ! Specific humidity (h2o mass mix ratio)
218 : real(r8), intent(in) :: E_o3mmr(pcols,pver) ! Ozone mass mixing ratio
219 : real(r8), intent(in) :: E_o2mmr(pcols) ! oxygen mass mixing ratio
220 : !
221 : real(r8), intent(in) :: E_cld(pcols,pver) ! Fractional cloud cover
222 : real(r8), intent(in) :: E_cicewp(pcols,pver) ! in-cloud cloud ice water path
223 : real(r8), intent(in) :: E_cliqwp(pcols,pver) ! in-cloud cloud liquid water path
224 : real(r8), intent(in) :: E_rel(pcols,pver) ! Liquid effective drop size (microns)
225 : real(r8), intent(in) :: E_rei(pcols,pver) ! Ice effective drop size (microns)
226 : !
227 : real(r8), intent(in) :: eccf ! Eccentricity factor (1./earth-sun dist^2)
228 : real(r8), intent(in) :: E_coszrs(pcols) ! Cosine solar zenith angle
229 : real(r8), intent(in) :: E_asdir(pcols) ! 0.2-0.7 micro-meter srfc alb: direct rad
230 : real(r8), intent(in) :: E_aldir(pcols) ! 0.7-5.0 micro-meter srfc alb: direct rad
231 : real(r8), intent(in) :: E_asdif(pcols) ! 0.2-0.7 micro-meter srfc alb: diffuse rad
232 : real(r8), intent(in) :: E_aldif(pcols) ! 0.7-5.0 micro-meter srfc alb: diffuse rad
233 : real(r8), intent(in) :: E_co2mmr(pcols) ! co2 column mean mmr
234 :
235 : !
236 : ! Aerosol radiative property arrays
237 : !
238 : real(r8),intent(in) :: E_aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth
239 : real(r8),intent(in) :: E_aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
240 : real(r8),intent(in) :: E_aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * w * tau
241 : real(r8),intent(in) :: E_aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
242 :
243 : !
244 : ! IN/OUT arguments
245 : !
246 : real(r8), intent(inout) :: pmxrgn(pcols,pverp) ! Maximum values of pressure for each
247 : ! ! maximally overlapped region.
248 : ! ! 0->pmxrgn(i,1) is range of pressure for
249 : ! ! 1st region,pmxrgn(i,1)->pmxrgn(i,2) for
250 : ! ! 2nd region, etc
251 : integer, intent(inout) :: nmxrgn(pcols) ! Number of maximally overlapped regions
252 : !
253 : ! Output arguments
254 : !
255 :
256 : real(r8), intent(out) :: solin(pcols) ! Incident solar flux
257 : real(r8), intent(out) :: qrs (pcols,pver) ! Solar heating rate
258 : real(r8), intent(out) :: qrsc(pcols,pver) ! Clearsky solar heating rate
259 : real(r8), intent(out) :: fsns(pcols) ! Surface absorbed solar flux
260 : real(r8), intent(out) :: fsnt(pcols) ! Total column absorbed solar flux
261 : real(r8), intent(out) :: fsntoa(pcols) ! Net solar flux at TOA
262 : real(r8), intent(out) :: fsutoa(pcols) ! Upward solar flux at TOA
263 : real(r8), intent(out) :: fsds(pcols) ! Flux shortwave downwelling surface
264 : !
265 : real(r8), intent(out) :: fsnsc(pcols) ! Clear sky surface absorbed solar flux
266 : real(r8), intent(out) :: fsdsc(pcols) ! Clear sky surface downwelling solar flux
267 : real(r8), intent(out) :: fsntc(pcols) ! Clear sky total column absorbed solar flx
268 : real(r8), intent(out) :: fsdtoa(pcols) ! Downwelling solar flux at TOA
269 : real(r8), intent(out) :: fsntoac(pcols) ! Clear sky net solar flx at TOA
270 : real(r8), intent(out) :: sols(pcols) ! Direct solar rad on surface (< 0.7)
271 : real(r8), intent(out) :: soll(pcols) ! Direct solar rad on surface (>= 0.7)
272 : real(r8), intent(out) :: solsd(pcols) ! Diffuse solar rad on surface (< 0.7)
273 : real(r8), intent(out) :: solld(pcols) ! Diffuse solar rad on surface (>= 0.7)
274 : real(r8), intent(out) :: fsnirtoa(pcols) ! Near-IR flux absorbed at toa
275 : real(r8), intent(out) :: fsnrtoac(pcols) ! Clear sky near-IR flux absorbed at toa
276 : real(r8), intent(out) :: fsnrtoaq(pcols) ! Net near-IR flux at toa >= 0.7 microns
277 :
278 : real(r8), intent(out) :: fns(pcols,pverp) ! net flux at interfaces
279 : real(r8), intent(out) :: fcns(pcols,pverp) ! net clear-sky flux at interfaces
280 :
281 : real(r8), intent(out) :: tauxcl_out(pcols,pver) ! liquid cloud visible sw optical depth
282 : real(r8), intent(out) :: tauxci_out(pcols,pver) ! ice cloud visible sw optical depth
283 :
284 : !
285 : !---------------------------Local variables-----------------------------
286 : !
287 : ! Local and reordered copies of the intent(in) variables
288 : !
289 : real(r8):: aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth
290 : real(r8):: aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
291 : real(r8):: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * w * tau
292 : real(r8):: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
293 : real(r8) :: pmid(pcols,pver) ! Level pressure
294 : real(r8) :: pint(pcols,pverp) ! Interface pressure
295 : real(r8) :: h2ommr(pcols,pver) ! Specific humidity (h2o mass mix ratio)
296 : real(r8) :: o3mmr(pcols,pver) ! Ozone mass mixing ratio
297 : !
298 : real(r8) :: cld(pcols,pver) ! Fractional cloud cover
299 : real(r8) :: cicewp(pcols,pver) ! in-cloud cloud ice water path
300 : real(r8) :: cliqwp(pcols,pver) ! in-cloud cloud liquid water path
301 : real(r8) :: rel(pcols,pver) ! Liquid effective drop size (microns)
302 : real(r8) :: rei(pcols,pver) ! Ice effective drop size (microns)
303 : !
304 : real(r8) :: coszrs(pcols) ! Cosine solar zenith angle
305 : real(r8) :: asdir(pcols) ! 0.2-0.7 micro-meter srfc alb: direct rad
306 : real(r8) :: aldir(pcols) ! 0.7-5.0 micro-meter srfc alb: direct rad
307 : real(r8) :: asdif(pcols) ! 0.2-0.7 micro-meter srfc alb: diffuse rad
308 : real(r8) :: aldif(pcols) ! 0.7-5.0 micro-meter srfc alb: diffuse rad
309 : real(r8) :: co2mmr(pcols) ! co2 column mean mmr
310 : real(r8) :: o2mmr(pcols) ! o2 column mean mmr
311 :
312 : real(r8) :: tot_irrad
313 : !
314 : ! Max/random overlap variables
315 : !
316 : real(r8) asort(pverp) ! 1 - cloud amounts to be sorted for max ovrlp.
317 : real(r8) atmp ! Temporary storage for sort when nxs = 2
318 : real(r8) cld0 ! 1 - (cld amt) used to make wstr, cstr, nstr
319 : real(r8) totwgt(pcols) ! Total of xwgts = total fractional area of
320 : ! grid-box covered by cloud configurations
321 : ! included in solution to fluxes
322 :
323 : real(r8) wgtv(nconfgmax) ! Weights for fluxes
324 : ! 1st index is configuration number
325 : real(r8) wstr(pverp,pverp) ! area weighting factors for streams
326 : ! 1st index is for stream #,
327 : ! 2nd index is for region #
328 :
329 : real(r8) xexpt ! solar direct beam trans. for layer above
330 : real(r8) xrdnd ! diffuse reflectivity for layer above
331 : real(r8) xrupd ! diffuse reflectivity for layer below
332 : real(r8) xrups ! direct-beam reflectivity for layer below
333 : real(r8) xtdnt ! total trans for layers above
334 :
335 : real(r8) xwgt ! product of cloud amounts
336 :
337 : real(r8) yexpt ! solar direct beam trans. for layer above
338 : real(r8) yrdnd ! diffuse reflectivity for layer above
339 : real(r8) yrupd ! diffuse reflectivity for layer below
340 : real(r8) ytdnd ! dif-beam transmission for layers above
341 : real(r8) ytupd ! dif-beam transmission for layers below
342 :
343 : real(r8) zexpt ! solar direct beam trans. for layer above
344 : real(r8) zrdnd ! diffuse reflectivity for layer above
345 : real(r8) zrupd ! diffuse reflectivity for layer below
346 : real(r8) zrups ! direct-beam reflectivity for layer below
347 : real(r8) ztdnt ! total trans for layers above
348 :
349 : logical new_term ! Flag for configurations to include in fluxes
350 : logical region_found ! flag for identifying regions
351 :
352 : integer ccon(nconfgmax,0:pverp,pcols)
353 : ! flags for presence of clouds
354 : ! 1st index is for level # (including
355 : ! layer above top of model and at surface)
356 : ! 2nd index is for configuration #
357 : integer cstr(0:pverp,pverp)
358 : ! flags for presence of clouds
359 : ! 1st index is for level # (including
360 : ! layer above top of model and at surface)
361 : ! 2nd index is for stream #
362 : integer icond(nconfgmax,0:pverp,pcols)
363 : ! Indices for copying rad. properties from
364 : ! one identical downward cld config.
365 : ! to another in adding method (step 2)
366 : ! 1st index is for interface # (including
367 : ! layer above top of model and at surface)
368 : ! 2nd index is for configuration # range
369 : integer iconu(nconfgmax,0:pverp,pcols)
370 : ! Indices for copying rad. properties from
371 : ! one identical upward configuration
372 : ! to another in adding method (step 2)
373 : ! 1st index is for interface # (including
374 : ! layer above top of model and at surface)
375 : ! 2nd index is for configuration # range
376 : integer iconfig ! Counter for random-ovrlap configurations
377 : integer irgn ! Index for max-overlap regions
378 : integer is0 ! Lower end of stream index range
379 : integer is1 ! Upper end of stream index range
380 : integer isn ! Stream index
381 : integer istr(pverp+1) ! index for stream #s during flux calculation
382 : integer istrtd(0:nconfgmax+1,0:pverp,pcols)
383 : ! indices into icond
384 : ! 1st index is for interface # (including
385 : ! layer above top of model and at surface)
386 : ! 2nd index is for configuration # range
387 : integer istrtu(0:nconfgmax+1,0:pverp,pcols)
388 : ! indices into iconu
389 : ! 1st index is for interface # (including
390 : ! layer above top of model and at surface)
391 : ! 2nd index is for configuration # range
392 : integer j ! Configuration index
393 : integer jj ! Configuration index
394 : integer k1 ! Level index
395 : integer k2 ! Level index
396 : integer ksort(pverp) ! Level indices of cloud amounts to be sorted
397 : integer ktmp ! Temporary storage for sort when nxs = 2
398 : integer kx1(0:pverp) ! Level index for top of max-overlap region
399 : integer kx2(0:pverp) ! Level index for bottom of max-overlap region
400 : integer l ! Index
401 : integer l0 ! Index
402 : integer mrgn ! Counter for nrgn
403 : integer mstr ! Counter for nstr
404 : integer n0 ! Number of configurations with ccon(:,k,:)==0
405 : integer n1 ! Number of configurations with ccon(:,k,:)==1
406 : integer nconfig(pcols) ! Number of random-ovrlap configurations
407 : integer nconfigm ! Value of config before testing for areamin,
408 : ! nconfgmax
409 : integer npasses ! number of passes over the indexing loop
410 : integer nrgn ! Number of max overlap regions at current
411 : ! longitude
412 : integer nstr(pverp) ! Number of unique cloud configurations
413 : ! ("streams") in a max-overlapped region
414 : ! 1st index is for region #
415 : integer nuniq ! # of unique cloud configurations
416 : integer nuniqd(0:pverp,pcols) ! # of unique cloud configurations: TOA
417 : ! to level k
418 : integer nuniqu(0:pverp,pcols) ! # of unique cloud configurations: surface
419 : ! to level k
420 : integer nxs ! Number of cloudy layers between k1 and k2
421 : integer ptr0(nconfgmax) ! Indices of configurations with ccon(:,k,:)==0
422 : integer ptr1(nconfgmax) ! Indices of configurations with ccon(:,k,:)==1
423 : integer ptrc(nconfgmax) ! Pointer for configurations sorted by wgtv
424 : integer, dimension(1) :: min_idx ! required for return val of func minloc
425 :
426 : !
427 : ! Other
428 : !
429 : integer ns ! Spectral loop index
430 : integer i ! Longitude loop index
431 : integer k ! Level loop index
432 : integer km1 ! k - 1
433 : integer kp1 ! k + 1
434 : integer n ! Loop index for daylight
435 : integer ksz ! dust size bin index
436 : integer kaer ! aerosol group index
437 : !
438 : ! A. Slingo's data for cloud particle radiative properties (from 'A GCM
439 : ! Parameterization for the Shortwave Properties of Water Clouds' JAS
440 : ! vol. 46 may 1989 pp 1419-1427)
441 : !
442 :
443 : real(r8), parameter :: abarl(4) = (/ 2.817e-02_r8, 2.682e-02_r8,2.264e-02_r8,1.281e-02_r8/)
444 : real(r8), parameter :: bbarl(4) = (/ 1.305_r8 , 1.346_r8 ,1.454_r8 ,1.641_r8 /)
445 : real(r8), parameter :: cbarl(4) = (/-5.62e-08_r8 ,-6.94e-06_r8 ,4.64e-04_r8 ,0.201_r8 /)
446 : real(r8), parameter :: dbarl(4) = (/ 1.63e-07_r8 , 2.35e-05_r8 ,1.24e-03_r8 ,7.56e-03_r8 /)
447 : real(r8), parameter :: ebarl(4) = (/ 0.829_r8 , 0.794_r8 ,0.754_r8 ,0.826_r8 /)
448 : real(r8), parameter :: fbarl(4) = (/ 2.482e-03_r8, 4.226e-03_r8,6.560e-03_r8,4.353e-03_r8/)
449 :
450 : real(r8) :: abarli ! A coefficient for current spectral band
451 : real(r8) :: bbarli ! B coefficient for current spectral band
452 : real(r8) :: cbarli ! C coefficient for current spectral band
453 : real(r8) :: dbarli ! D coefficient for current spectral band
454 : real(r8) :: ebarli ! E coefficient for current spectral band
455 : real(r8) :: fbarli ! F coefficient for current spectral band
456 : !
457 : ! Caution... A. Slingo recommends no less than 4.0 micro-meters nor
458 : ! greater than 20 micro-meters
459 : !
460 : ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
461 : !
462 : real(r8), parameter :: abari(4) = (/ 3.448e-03_r8, 3.448e-03_r8,3.448e-03_r8,3.448e-03_r8/)
463 : real(r8), parameter :: bbari(4) = (/ 2.431_r8 , 2.431_r8 ,2.431_r8 ,2.431_r8 /)
464 : real(r8), parameter :: cbari(4) = (/ 1.00e-05_r8 , 1.10e-04_r8 ,1.861e-02_r8,.46658_r8 /)
465 : real(r8), parameter :: dbari(4) = (/ 0.0_r8 , 1.405e-05_r8,8.328e-04_r8,2.05e-05_r8 /)
466 : real(r8), parameter :: ebari(4) = (/ 0.7661_r8 , 0.7730_r8 ,0.794_r8 ,0.9595_r8 /)
467 : real(r8), parameter :: fbari(4) = (/ 5.851e-04_r8, 5.665e-04_r8,7.267e-04_r8,1.076e-04_r8/)
468 :
469 : real(r8) :: abarii ! A coefficient for current spectral band
470 : real(r8) :: bbarii ! B coefficient for current spectral band
471 : real(r8) :: cbarii ! C coefficient for current spectral band
472 : real(r8) :: dbarii ! D coefficient for current spectral band
473 : real(r8) :: ebarii ! E coefficient for current spectral band
474 : real(r8) :: fbarii ! F coefficient for current spectral band
475 :
476 : !
477 : ! UPDATE TO H2O NEAR-IR: Delta optimized for Hitran 2K and CKD 2.4
478 : !
479 : real(r8), parameter :: delta = 0.0014257179260883_r8
480 : !
481 : ! END UPDATE
482 : !
483 : real(r8) :: albdir(pcols,nswbands) ! Current spc intrvl srf alb to direct rad
484 : real(r8) :: albdif(pcols,nswbands) ! Current spc intrvl srf alb to diffuse rad
485 : !
486 : ! Next series depends on spectral interval
487 : !
488 : real(r8) :: wgtint ! Weight for specific spectral interval
489 :
490 : !
491 : ! weight for 0.64 - 0.7 microns appropriate to clear skies over oceans
492 : !
493 : real(r8), parameter :: nirwgt(nswbands) = &
494 : (/ 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, &
495 : 0.0_r8, 0.0_r8, 0.0_r8, 0.320518_r8, 1.0_r8, 1.0_r8, &
496 : 1.0_r8, 1.0_r8, 1.0_r8, 1.0_r8, 1.0_r8, &
497 : 1.0_r8, 1.0_r8, 1.0_r8 /)
498 :
499 : !
500 : ! UPDATE TO H2O NEAR-IR: Rayleigh scattering optimized for Hitran 2K & CKD 2.4
501 : !
502 : real(r8), parameter :: raytau(nswbands) = &
503 : (/ 4.020_r8, 2.180_r8, 1.700_r8, 1.450_r8, 1.250_r8, &
504 : 1.085_r8, 0.730_r8, v_raytau_35, v_raytau_64, &
505 : 0.02899756_r8, 0.01356763_r8, 0.00537341_r8, &
506 : 0.00228515_r8, 0.00105028_r8, 0.00046631_r8, &
507 : 0.00025734_r8, &
508 : .0001_r8, .0001_r8, .0001_r8/)
509 : !
510 : ! END UPDATE
511 : !
512 :
513 : !
514 : ! Absorption coefficients
515 : !
516 : !
517 : ! UPDATE TO H2O NEAR-IR: abh2o optimized for Hitran 2K and CKD 2.4
518 : !
519 : real(r8), parameter :: abh2o(nswbands) = &
520 : (/ .000_r8, .000_r8, .000_r8, .000_r8, .000_r8, &
521 : .000_r8, .000_r8, .000_r8, .000_r8, &
522 : 0.00256608_r8, 0.06310504_r8, 0.42287445_r8, 2.45397941_r8, &
523 : 11.20070807_r8, 47.66091389_r8, 240.19010243_r8, &
524 : .000_r8, .000_r8, .000_r8/)
525 : !
526 : ! END UPDATE
527 : !
528 :
529 : real(r8), parameter :: abo3(nswbands) = &
530 : (/5.370e+04_r8, 13.080e+04_r8, 9.292e+04_r8, 4.530e+04_r8, 1.616e+04_r8, &
531 : 4.441e+03_r8, 1.775e+02_r8, v_abo3_35, v_abo3_64, .000_r8, &
532 : .000_r8, .000_r8 , .000_r8 , .000_r8 , .000_r8, &
533 : .000_r8, .000_r8 , .000_r8 , .000_r8 /)
534 :
535 : real(r8), parameter :: abco2(nswbands) = &
536 : (/ .000_r8, .000_r8, .000_r8, .000_r8, .000_r8, &
537 : .000_r8, .000_r8, .000_r8, .000_r8, .000_r8, &
538 : .000_r8, .000_r8, .000_r8, .000_r8, .000_r8, &
539 : .000_r8, .094_r8, .196_r8, 1.963_r8/)
540 :
541 : real(r8), parameter :: abo2(nswbands) = &
542 : (/ .000_r8, .000_r8, .000_r8, .000_r8, .000_r8, &
543 : .000_r8, .000_r8, .000_r8,1.11e-05_r8,6.69e-05_r8, &
544 : .000_r8, .000_r8, .000_r8, .000_r8, .000_r8, &
545 : .000_r8, .000_r8, .000_r8, .000_r8/)
546 : !
547 : ! Diagnostic and accumulation arrays; note that fswup, and
548 : ! fswdn are not used in the computation,but are retained for future use.
549 : !
550 : real(r8) solflx(pcols) ! Solar flux in current interval
551 : real(r8) totfld (pcols,0:pver) ! Spectrally summed flux divergence
552 : real(r8) totfldc(pcols,0:pver) ! Spectrally summed flux divergence (clearsky)
553 : real(r8) fswup(pcols,0:pverp) ! Spectrally summed up flux
554 : real(r8) fswdn(pcols,0:pverp) ! Spectrally summed down flux
555 : !
556 : ! Cloud radiative property arrays
557 : !
558 : real(r8) tauxcl(pcols,0:pver) ! water cloud extinction optical depth
559 : real(r8) tauxci(pcols,0:pver) ! ice cloud extinction optical depth
560 : real(r8) wcl(pcols,0:pver) ! liquid cloud single scattering albedo
561 : real(r8) gcl(pcols,0:pver) ! liquid cloud asymmetry parameter
562 : real(r8) fcl(pcols,0:pver) ! liquid cloud forward scattered fraction
563 : real(r8) wci(pcols,0:pver) ! ice cloud single scattering albedo
564 : real(r8) gci(pcols,0:pver) ! ice cloud asymmetry parameter
565 : real(r8) fci(pcols,0:pver) ! ice cloud forward scattered fraction
566 :
567 : !
568 : ! Various arrays and other constants:
569 : !
570 : real(r8) pflx(pcols,0:pverp) ! Interface press, including extra layer
571 : real(r8) zenfac(pcols) ! Square root of cos solar zenith angle
572 : real(r8) sqrco2(pcols) ! Square root of the co2 mass mixg ratio
573 : real(r8) tmp1 ! Temporary constant array
574 : real(r8) tmp2 ! Temporary constant array
575 : real(r8) pdel ! Pressure difference across layer
576 : real(r8) path ! Mass path of layer
577 : real(r8) ptop ! Lower interface pressure of extra layer
578 : real(r8) ptho2 ! Used to compute mass path of o2
579 : real(r8) ptho3 ! Used to compute mass path of o3
580 : real(r8) pthco2 ! Used to compute mass path of co2
581 : real(r8) pthh2o ! Used to compute mass path of h2o
582 : real(r8) h2ostr ! Inverse sq. root h2o mass mixing ratio
583 :
584 : real(r8) wavmin(nswbands) ! Spectral interval minimum wavelength
585 : real(r8) wavmax(nswbands) ! Spectral interval maximum wavelength
586 : real(r8) wavmid(nswbands) ! Spectral interval middle wavelength
587 : real(r8) trayoslp ! Rayleigh optical depth/standard pressure
588 : real(r8) tmp1l ! Temporary constant array
589 : real(r8) tmp2l ! Temporary constant array
590 : real(r8) tmp3l ! Temporary constant array
591 : real(r8) tmp1i ! Temporary constant array
592 : real(r8) tmp2i ! Temporary constant array
593 : real(r8) tmp3i ! Temporary constant array
594 : real(r8) rdenom ! Multiple scattering term
595 : real(r8) rdirexp ! layer direct ref times exp transmission
596 : real(r8) tdnmexp ! total transmission - exp transmission
597 : real(r8) psf(nswbands) ! Frac of solar flux in spect interval
598 : !
599 : ! Layer absorber amounts; note that 0 refers to the extra layer added
600 : ! above the top model layer
601 : !
602 : real(r8) uh2o(pcols,0:pver) ! Layer absorber amount of h2o
603 : real(r8) uo3(pcols,0:pver) ! Layer absorber amount of o3
604 : real(r8) uco2(pcols,0:pver) ! Layer absorber amount of co2
605 : real(r8) uo2(pcols,0:pver) ! Layer absorber amount of o2
606 : !
607 : ! Total column absorber amounts:
608 : !
609 : real(r8) uth2o(pcols) ! Total column absorber amount of h2o
610 : real(r8) uto3(pcols) ! Total column absorber amount of o3
611 : real(r8) utco2(pcols) ! Total column absorber amount of co2
612 : real(r8) uto2(pcols) ! Total column absorber amount of o2
613 : !
614 : ! These arrays are defined for pver model layers; 0 refers to the extra
615 : ! layer on top:
616 : !
617 : real(r8) rdir(nswbands,pcols,0:pver) ! Layer reflectivity to direct rad
618 : real(r8) rdif(nswbands,pcols,0:pver) ! Layer reflectivity to diffuse rad
619 : real(r8) tdir(nswbands,pcols,0:pver) ! Layer transmission to direct rad
620 : real(r8) tdif(nswbands,pcols,0:pver) ! Layer transmission to diffuse rad
621 : real(r8) explay(nswbands,pcols,0:pver) ! Solar beam exp trans. for layer
622 :
623 : real(r8) rdirc(nswbands,pcols,0:pver) ! Clear Layer reflec. to direct rad
624 : real(r8) rdifc(nswbands,pcols,0:pver) ! Clear Layer reflec. to diffuse rad
625 : real(r8) tdirc(nswbands,pcols,0:pver) ! Clear Layer trans. to direct rad
626 : real(r8) tdifc(nswbands,pcols,0:pver) ! Clear Layer trans. to diffuse rad
627 : real(r8) explayc(nswbands,pcols,0:pver) ! Solar beam exp trans. clear layer
628 :
629 : real(r8) fus(pcols,pverp) ! Upward flux (added for CRM)
630 : real(r8) fds(pcols,pverp) ! Downward flux (added for CRM)
631 : real(r8) fusc(pcols,pverp) ! Upward clear-sky flux (added for CRM)
632 : real(r8) fdsc(pcols,pverp) ! Downward clear-sky flux (added for CRM)
633 :
634 : real(r8) flxdiv ! Flux divergence for layer
635 :
636 : !
637 : ! Temporary arrays for either clear or cloudy values.
638 : !
639 : real(r8), dimension(nswbands) :: Trdir
640 : real(r8), dimension(nswbands) :: Trdif
641 : real(r8), dimension(nswbands) :: Ttdir
642 : real(r8), dimension(nswbands) :: Ttdif
643 : real(r8), dimension(nswbands) :: Texplay
644 : !
645 : !
646 : ! Radiative Properties:
647 : !
648 : ! There are 1 classes of properties:
649 : ! (1. All-sky bulk properties
650 : ! (2. Clear-sky properties
651 : !
652 : ! The first set of properties are generated during step 2 of the solution.
653 : !
654 : ! These arrays are defined at model interfaces; in 1st index (for level #),
655 : ! 0 is the top of the extra layer above the model top, and
656 : ! pverp is the earth surface. 2nd index is for cloud configuration
657 : ! defined over a whole column.
658 : !
659 : real(r8) exptdn(nswbands,0:pverp,nconfgmax,pcols) ! Sol. beam trans from layers above
660 : real(r8) rdndif(nswbands,0:pverp,nconfgmax,pcols) ! Ref to dif rad for layers above
661 : real(r8) rupdif(nswbands,0:pverp,nconfgmax,pcols) ! Ref to dif rad for layers below
662 : real(r8) rupdir(nswbands,0:pverp,nconfgmax,pcols) ! Ref to dir rad for layers below
663 : real(r8) tdntot(nswbands,0:pverp,nconfgmax,pcols) ! Total trans for layers above
664 : !
665 : ! Bulk properties used during the clear-sky calculation.
666 : !
667 : real(r8) exptdnc(pcols,0:pverp) ! clr: Sol. beam trans from layers above
668 : real(r8) rdndifc(pcols,0:pverp) ! clr: Ref to dif rad for layers above
669 : real(r8) rupdifc(pcols,0:pverp) ! clr: Ref to dif rad for layers below
670 : real(r8) rupdirc(pcols,0:pverp) ! clr: Ref to dir rad for layers below
671 : real(r8) tdntotc(pcols,0:pverp) ! clr: Total trans for layers above
672 :
673 : real(r8) fluxup(nswbands,0:pverp,pcols) ! Up flux at model interface
674 : real(r8) fluxdn(nswbands,0:pverp,pcols) ! Down flux at model interface
675 : real(r8) wexptdn(nswbands,pcols) ! Direct solar beam trans. to surface
676 : !
677 : ! Scalars used in vectorization
678 : !
679 : integer :: kk
680 : !
681 : ! Arrays used in vectorization
682 : !
683 : real(r8) v_wgtv(nconfgmax,pcols) ! Weights for fluxes
684 :
685 : real(r8) :: rdiff, ro, rn
686 : rdiff(ro,rn) = abs((ro-rn)/merge(ro,1.0_r8,ro /= 0.0_r8))
687 :
688 : ! solar variability factor
689 : real(r8) :: sfac(nswbands)
690 :
691 : character(len=16) :: microp_scheme ! microphysics scheme
692 :
693 : !-----------------------------------------------------------------------
694 : ! START OF CALCULATION
695 : !-----------------------------------------------------------------------
696 :
697 749232 : call phys_getopts(microp_scheme_out=microp_scheme)
698 :
699 : !
700 : ! Initialize output fields:
701 : !
702 12510432 : fsds(1:ncol) = 0.0_r8
703 :
704 12510432 : fsnirtoa(1:ncol) = 0.0_r8
705 12510432 : fsnrtoac(1:ncol) = 0.0_r8
706 12510432 : fsnrtoaq(1:ncol) = 0.0_r8
707 :
708 12510432 : fsns(1:ncol) = 0.0_r8
709 12510432 : fsnsc(1:ncol) = 0.0_r8
710 12510432 : fsdsc(1:ncol) = 0.0_r8
711 :
712 12510432 : fsnt(1:ncol) = 0.0_r8
713 12510432 : fsntc(1:ncol) = 0.0_r8
714 12510432 : fsntoa(1:ncol) = 0.0_r8
715 12510432 : fsdtoa(1:ncol) = 0.0_r8
716 12510432 : fsutoa(1:ncol) = 0.0_r8
717 12510432 : fsntoac(1:ncol) = 0.0_r8
718 :
719 12510432 : solin(1:ncol) = 0.0_r8
720 :
721 12510432 : sols(1:ncol) = 0.0_r8
722 12510432 : soll(1:ncol) = 0.0_r8
723 12510432 : solsd(1:ncol) = 0.0_r8
724 12510432 : solld(1:ncol) = 0.0_r8
725 :
726 326020464 : qrs (1:ncol,1:pver) = 0.0_r8
727 326020464 : qrsc(1:ncol,1:pver) = 0.0_r8
728 338530896 : fns(1:ncol,1:pverp) = 0.0_r8
729 338530896 : fcns(1:ncol,1:pverp) = 0.0_r8
730 749232 : if (single_column.and.scm_crm_mode) then
731 0 : fus(1:ncol,1:pverp) = 0.0_r8
732 0 : fds(1:ncol,1:pverp) = 0.0_r8
733 0 : fusc(:ncol,:pverp) = 0.0_r8
734 0 : fdsc(:ncol,:pverp) = 0.0_r8
735 : endif
736 :
737 749232 : tauxcl_out(1:pcols,1:pver) = 0.0_r8
738 749232 : tauxci_out(1:pcols,1:pver) = 0.0_r8
739 :
740 : !
741 : ! If night everywhere, return:
742 : !
743 749232 : if ( Nday == 0 ) then
744 : return
745 : endif
746 :
747 : !
748 : ! Rearrange input arrays
749 : !
750 :
751 391783 : call CmpDayNite(E_pmid, pmid, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
752 391783 : call CmpDayNite(E_pint, pint, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
753 391783 : call CmpDayNite(E_h2ommr, h2ommr, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
754 391783 : call CmpDayNite(E_o3mmr, o3mmr, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
755 391783 : call CmpDayNite(E_cld, cld, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
756 391783 : call CmpDayNite(E_cicewp, cicewp, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
757 391783 : call CmpDayNite(E_cliqwp, cliqwp, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
758 391783 : call CmpDayNite(E_rel, rel, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
759 391783 : call CmpDayNite(E_rei, rei, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
760 391783 : call CmpDayNite(E_coszrs, coszrs, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
761 391783 : call CmpDayNite(E_asdir, asdir, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
762 391783 : call CmpDayNite(E_aldir, aldir, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
763 391783 : call CmpDayNite(E_asdif, asdif, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
764 391783 : call CmpDayNite(E_aldif, aldif, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
765 391783 : call CmpDayNite(E_co2mmr, co2mmr, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
766 391783 : call CmpDayNite(E_o2mmr, o2mmr, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
767 :
768 391783 : call CmpDayNite(pmxrgn, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
769 391783 : call CmpDayNite(nmxrgn, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
770 391783 : call CmpDayNite(E_aer_tau, aer_tau, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
771 391783 : call CmpDayNite(E_aer_tau_w, aer_tau_w, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
772 391783 : call CmpDayNite(E_aer_tau_w_g, aer_tau_w_g, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
773 391783 : call CmpDayNite(E_aer_tau_w_f, aer_tau_w_f, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
774 :
775 391783 : if (scm_crm_mode) then
776 : ! overwrite albedos for CRM
777 0 : if(have_asdir) asdir = asdirobs(1)
778 0 : if(have_asdif) asdif = asdifobs(1)
779 0 : if(have_aldir) aldir = aldirobs(1)
780 0 : if(have_aldif) aldif = aldifobs(1)
781 : endif
782 : !
783 : ! Perform other initializations
784 : !
785 391783 : tmp1 = 0.5_r8/(gravit*sslp)
786 391783 : tmp2 = delta/gravit
787 :
788 6272383 : do i = 1, Nday
789 6272383 : sqrco2(i) = sqrt(co2mmr(i))
790 : end do
791 :
792 391783 : if ( do_spctrl_scaling ) then
793 0 : call get_variability(sfac)
794 0 : tot_irrad = ref_tsi
795 : else
796 391783 : tot_irrad = sol_tsi
797 : endif
798 :
799 10969924 : do k=1,pverp
800 169746124 : do i=1,Nday
801 169354341 : pflx(i,k) = pint(i,k)
802 : end do
803 : end do
804 :
805 6272383 : do i=1,Nday
806 : !
807 : ! Define solar incident radiation and interface pressures:
808 : !
809 :
810 5880600 : solin(i) = tot_irrad*1.e3_r8*eccf*coszrs(i)
811 :
812 5880600 : pflx(i,0) = 0._r8
813 : !
814 : ! Compute optical paths:
815 : !
816 5880600 : ptop = pflx(i,1)
817 5880600 : ptho2 = o2mmr(i) * ptop / gravit
818 5880600 : ptho3 = o3mmr(i,1) * ptop / gravit
819 5880600 : pthco2 = sqrco2(i) * (ptop / gravit)
820 5880600 : h2ostr = sqrt( 1._r8 / h2ommr(i,1) )
821 5880600 : zenfac(i) = sqrt(coszrs(i))
822 : pthh2o = ptop**2*tmp1 + (ptop*rga)* &
823 5880600 : (h2ostr*zenfac(i)*delta)
824 5880600 : uh2o(i,0) = h2ommr(i,1)*pthh2o
825 5880600 : uco2(i,0) = zenfac(i)*pthco2
826 5880600 : uo2 (i,0) = zenfac(i)*ptho2
827 6272383 : uo3 (i,0) = ptho3
828 : !
829 : ! End do i=1,Nday
830 : !
831 : end do
832 :
833 10578141 : do k=1,pver
834 :
835 163473741 : do i=1,Nday
836 :
837 152895600 : pdel = pflx(i,k+1) - pflx(i,k)
838 152895600 : path = pdel / gravit
839 152895600 : ptho2 = o2mmr(i) * path
840 152895600 : ptho3 = o3mmr(i,k) * path
841 152895600 : pthco2 = sqrco2(i) * path
842 152895600 : h2ostr = sqrt(1.0_r8/h2ommr(i,k))
843 152895600 : pthh2o = (pflx(i,k+1)**2 - pflx(i,k)**2)*tmp1 + pdel*h2ostr*zenfac(i)*tmp2
844 152895600 : uh2o(i,k) = h2ommr(i,k)*pthh2o
845 152895600 : uco2(i,k) = zenfac(i)*pthco2
846 152895600 : uo2 (i,k) = zenfac(i)*ptho2
847 163081958 : uo3 (i,k) = ptho3
848 :
849 : !
850 : ! End do i=1,Nday
851 : !
852 : end do
853 : !
854 : ! End k=1,pver
855 : !
856 : end do
857 : !
858 : ! Compute column absorber amounts for the clear sky computation:
859 : !
860 6272383 : do i=1,Nday
861 :
862 5880600 : uth2o(i) = 0.0_r8
863 5880600 : uto3(i) = 0.0_r8
864 5880600 : utco2(i) = 0.0_r8
865 5880600 : uto2(i) = 0.0_r8
866 :
867 159167983 : do k=1,pver
868 152895600 : uth2o(i) = uth2o(i) + uh2o(i,k)
869 152895600 : uto3(i) = uto3(i) + uo3(i,k)
870 152895600 : utco2(i) = utco2(i) + uco2(i,k)
871 158776200 : uto2(i) = uto2(i) + uo2(i,k)
872 : !
873 : ! End k=1,pver
874 : !
875 : end do
876 : !
877 : ! End do i=1,Nday
878 : !
879 : end do
880 : !
881 : ! Set cloud properties for top (0) layer; so long as tauxcl is zero,
882 : ! there is no cloud above top of model; the other cloud properties
883 : ! are arbitrary:
884 : !
885 6272383 : do i=1,Nday
886 :
887 5880600 : tauxcl(i,0) = 0._r8
888 5880600 : wcl(i,0) = 0.999999_r8
889 5880600 : gcl(i,0) = 0.85_r8
890 5880600 : fcl(i,0) = 0.725_r8
891 5880600 : tauxci(i,0) = 0._r8
892 5880600 : wci(i,0) = 0.999999_r8
893 5880600 : gci(i,0) = 0.85_r8
894 6272383 : fci(i,0) = 0.725_r8
895 : !
896 : ! End do i=1,Nday
897 : !
898 : end do
899 : !
900 : ! Begin spectral loop
901 : !
902 7835660 : do ns=1,nswbands
903 : !
904 : ! Set cloud extinction optical depth, single scatter albedo,
905 : ! asymmetry parameter, and forward scattered fraction:
906 : !
907 7443877 : abarli = abarl(indxsl(ns))
908 7443877 : bbarli = bbarl(indxsl(ns))
909 7443877 : cbarli = cbarl(indxsl(ns))
910 7443877 : dbarli = dbarl(indxsl(ns))
911 7443877 : ebarli = ebarl(indxsl(ns))
912 7443877 : fbarli = fbarl(indxsl(ns))
913 : !
914 7443877 : abarii = abari(indxsl(ns))
915 7443877 : bbarii = bbari(indxsl(ns))
916 7443877 : cbarii = cbari(indxsl(ns))
917 7443877 : dbarii = dbari(indxsl(ns))
918 7443877 : ebarii = ebari(indxsl(ns))
919 7443877 : fbarii = fbari(indxsl(ns))
920 : !
921 : ! adjustfraction within spectral interval to allow for the possibility of
922 : ! sub-divisions within a particular interval:
923 : !
924 7443877 : psf(ns) = 1.0_r8
925 7443877 : if(ph2o(ns)/=0._r8) psf(ns) = psf(ns)*ph2o(ns)
926 7443877 : if(pco2(ns)/=0._r8) psf(ns) = psf(ns)*pco2(ns)
927 7443877 : if(po2 (ns)/=0._r8) psf(ns) = psf(ns)*po2 (ns)
928 :
929 200984679 : do k=1,pver
930 :
931 3106001079 : do i=1,Nday
932 :
933 : ! liquid
934 : ! note that optical properties for liquid valid only
935 : ! in range of 4.2 > rel > 16 micron (Slingo 89)
936 2905016400 : if ( microp_scheme == 'MG' ) then
937 0 : tmp2l = 1._r8 - cbarli - dbarli*min(max(4.2_r8,rel(i,k)),16._r8)
938 0 : tmp3l = fbarli*min(max(4.2_r8,rel(i,k)),16._r8)
939 : else
940 2905016400 : tmp2l = 1._r8 - cbarli - dbarli*rel(i,k)
941 2905016400 : tmp3l = fbarli*rel(i,k)
942 : endif
943 :
944 : ! ice
945 : ! note that optical properties for ice valid only
946 : ! in range of 13 > rei > 130 micron (Ebert and Curry 92)
947 2905016400 : if ( microp_scheme == 'MG' ) then
948 0 : tmp2i = 1._r8 - cbarii - dbarii*min(max(13._r8,rei(i,k)),130._r8)
949 0 : tmp3i = fbarii*min(max(13._r8,rei(i,k)),130._r8)
950 : else
951 2905016400 : tmp2i = 1._r8 - cbarii - dbarii*rei(i,k)
952 2905016400 : tmp3i = fbarii*rei(i,k)
953 : endif
954 :
955 2905016400 : if (cld(i,k) >= cldmin .and. cld(i,k) >= cldeps) then
956 :
957 : ! liquid
958 735757748 : if ( microp_scheme == 'MG' ) then
959 0 : tmp1l = abarli + bbarli/min(max(4.2_r8,rel(i,k)),16._r8)
960 : else
961 735757748 : tmp1l = abarli + bbarli/rel(i,k)
962 : endif
963 :
964 : ! ice
965 735757748 : if ( microp_scheme == 'MG' ) then
966 0 : tmp1i = abarii + bbarii/max(13._r8,min(rei(i,k),130._r8))
967 : else
968 735757748 : tmp1i = abarii + bbarii/rei(i,k)
969 : endif
970 :
971 735757748 : tauxcl(i,k) = cliqwp(i,k)*tmp1l
972 735757748 : tauxci(i,k) = cicewp(i,k)*tmp1i
973 :
974 735757748 : if (ns .eq. idx_sw_diag) then
975 38724092 : tauxcl_out(i,k) = cliqwp(i,k)*tmp1l
976 38724092 : tauxci_out(i,k) = cicewp(i,k)*tmp1i
977 : endif
978 : else
979 2169258652 : tauxcl(i,k) = 0.0_r8
980 2169258652 : tauxci(i,k) = 0.0_r8
981 : endif
982 :
983 : ! Do not let single scatter albedo be 1. Delta-eddington solution
984 : ! for non-conservative case has different analytic form from solution
985 : ! for conservative case, and raddedmx is written for non-conservative case.
986 2905016400 : wcl(i,k) = min(tmp2l,.999999_r8)
987 2905016400 : gcl(i,k) = ebarli + tmp3l
988 2905016400 : fcl(i,k) = gcl(i,k)*gcl(i,k)
989 :
990 2905016400 : wci(i,k) = min(tmp2i,.999999_r8)
991 2905016400 : gci(i,k) = ebarii + tmp3i
992 3098557202 : fci(i,k) = gci(i,k)*gci(i,k)
993 :
994 : end do ! End do i=1,Nday
995 : end do ! End do k=1,pver
996 :
997 :
998 : !
999 : ! Set reflectivities for surface based on mid-point wavelength
1000 : !
1001 7443877 : call get_sw_spectral_boundaries(wavmin, wavmax, 'micrometer')
1002 7443877 : wavmid(ns) = 0.5_r8*(wavmin(ns) + wavmax(ns))
1003 : !
1004 : ! Wavelength less than 0.7 micro-meter
1005 : !
1006 7443877 : if (wavmid(ns) < 0.7_r8 ) then
1007 56451447 : do i=1,Nday
1008 52925400 : albdir(i,ns) = asdir(i)
1009 56451447 : albdif(i,ns) = asdif(i)
1010 : end do
1011 : !
1012 : ! Wavelength greater than 0.7 micro-meter
1013 : !
1014 : else
1015 62723830 : do i=1,Nday
1016 58806000 : albdir(i,ns) = aldir(i)
1017 62723830 : albdif(i,ns) = aldif(i)
1018 : end do
1019 : end if
1020 7443877 : trayoslp = raytau(ns)/sslp
1021 : !
1022 : ! Layer input properties now completely specified; compute the
1023 : ! delta-Eddington solution reflectivities and transmissivities
1024 : ! for each layer
1025 : !
1026 : call raddedmx(coszrs ,Nday , &
1027 : abh2o(ns),abo3(ns) ,abco2(ns),abo2(ns) , &
1028 : uh2o ,uo3 ,uco2 ,uo2 , &
1029 : trayoslp ,pflx ,ns , &
1030 : tauxcl ,wcl ,gcl ,fcl , &
1031 : tauxci ,wci ,gci ,fci , &
1032 : aer_tau(:,:,ns) ,aer_tau_w(:,:,ns) ,aer_tau_w_g(:,:,ns) ,aer_tau_w_f(:,:,ns) , &
1033 : rdir ,rdif ,tdir ,tdif ,explay , &
1034 7835660 : rdirc ,rdifc ,tdirc ,tdifc ,explayc )
1035 : !
1036 : ! End spectral loop
1037 : !
1038 : end do
1039 : !
1040 : !----------------------------------------------------------------------
1041 : !
1042 : ! Solution for max/random cloud overlap.
1043 : !
1044 : ! Steps:
1045 : ! (1. delta-Eddington solution for each layer (called above)
1046 : !
1047 : ! (2. The adding method is used to
1048 : ! compute the reflectivity and transmissivity to direct and diffuse
1049 : ! radiation from the top and bottom of the atmosphere for each
1050 : ! cloud configuration. This calculation is based upon the
1051 : ! max-random overlap assumption.
1052 : !
1053 : ! (3. to solve for the fluxes, combine the
1054 : ! bulk properties of the atmosphere above/below the region.
1055 : !
1056 : ! Index calculations for steps 2-3 are performed outside spectral
1057 : ! loop to avoid redundant calculations. Index calculations (with
1058 : ! application of areamin & nconfgmax conditions) are performed
1059 : ! first to identify the minimum subset of terms for the configurations
1060 : ! satisfying the areamin & nconfgmax conditions. This minimum set is
1061 : ! used to identify the corresponding minimum subset of terms in
1062 : ! steps 2 and 3.
1063 : !
1064 6268528 : do iconfig = 1, nconfgmax
1065 94085745 : ccon(iconfig,0,1:Nday) = 0
1066 94085745 : ccon(iconfig,pverp,1:Nday) = 0
1067 :
1068 94085745 : icond(iconfig,0,1:Nday) = iconfig
1069 94477528 : iconu(iconfig,pverp,1:Nday) = iconfig
1070 : end do
1071 : !
1072 : ! Construction of nuniqu/d, istrtu/d, iconu/d using binary tree
1073 : !
1074 6272383 : nuniqd(0,1:Nday) = 1
1075 6272383 : nuniqu(pverp,1:Nday) = 1
1076 :
1077 6272383 : istrtd(1,0,1:Nday) = 1
1078 6272383 : istrtu(1,pverp,1:Nday) = 1
1079 :
1080 :
1081 : !CSD$ PARALLEL DO PRIVATE( npasses, kx2, mrgn, region_found, k1, k2, kx1, nxs, ksort, asort ) &
1082 : !CSD$ PRIVATE ( ktmp, atmp, cstr, mstr, nstr, cld0, wstr, nrgn, nconfigm, istr, new_term, xwgt ) &
1083 : !CSD$ PRIVATE ( j, ptrc, wgtv, km1, nuniq, is0, is1, n0, n1, ptr0, ptr1, kp1, i, irgn ) &
1084 : !CSD$ PRIVATE ( k, l, iconfig, l0, isn )
1085 6272383 : do i=1,Nday
1086 :
1087 : !----------------------------------------------------------------------
1088 : ! INDEX CALCULATIONS FOR MAX OVERLAP
1089 : !
1090 : ! The column is divided into sets of adjacent layers, called regions,
1091 : ! in which the clouds are maximally overlapped. The clouds are
1092 : ! randomly overlapped between different regions. The number of
1093 : ! regions in a column is set by nmxrgn, and the range of pressures
1094 : ! included in each region is set by pmxrgn.
1095 : !
1096 : ! The following calculations determine the number of unique cloud
1097 : ! configurations (assuming maximum overlap), called "streams",
1098 : ! within each region. Each stream consists of a vector of binary
1099 : ! clouds (either 0 or 100% cloud cover). Over the depth of the region,
1100 : ! each stream requires a separate calculation of radiative properties. These
1101 : ! properties are generated using the adding method from
1102 : ! the radiative properties for each layer calculated by raddedmx.
1103 : !
1104 : ! The upward and downward-propagating streams are treated
1105 : ! separately.
1106 : !
1107 : ! We will refer to a particular configuration of binary clouds
1108 : ! within a single max-overlapped region as a "stream". We will
1109 : ! refer to a particular arrangement of binary clouds over the entire column
1110 : ! as a "configuration".
1111 : !
1112 : ! This section of the code generates the following information:
1113 : ! (1. nrgn : the true number of max-overlap regions (need not = nmxrgn)
1114 : ! (2. nstr : the number of streams in a region (>=1)
1115 : ! (3. cstr : flags for presence of clouds at each layer in each stream
1116 : ! (4. wstr : the fractional horizontal area of a grid box covered
1117 : ! by each stream
1118 : ! (5. kx1,2 : level indices for top/bottom of each region
1119 : !
1120 : ! The max-overlap calculation proceeds in 3 stages:
1121 : ! (1. compute layer radiative properties in raddedmx.
1122 : ! (2. combine these properties between layers
1123 : ! (3. combine properties to compute fluxes at each interface.
1124 : !
1125 : ! Most of the indexing information calculated here is used in steps 2-3
1126 : ! after the call to raddedmx.
1127 : !
1128 : ! Initialize indices for layers to be max-overlapped
1129 : !
1130 : ! Loop to handle fix in totwgt=0. For original overlap config
1131 : ! from npasses = 0.
1132 : !
1133 : npasses = 0
1134 0 : do
1135 20810214 : do irgn = 0, nmxrgn(i)
1136 20810214 : kx2(irgn) = 0
1137 : end do
1138 : mrgn = 0
1139 : !
1140 : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
1141 : !
1142 14929614 : do irgn = 1, nmxrgn(i)
1143 : !
1144 : ! Calculate min/max layer indices inside region.
1145 : !
1146 9049014 : region_found = .false.
1147 9049014 : if (kx2(irgn-1) < pver) then
1148 9049014 : k1 = kx2(irgn-1)+1
1149 9049014 : kx1(irgn) = k1
1150 9049014 : kx2(irgn) = k1-1
1151 44578071 : do k2 = pver, k1, -1
1152 44578071 : if (pmid(i,k2) <= pmxrgn(i,irgn)) then
1153 9049014 : kx2(irgn) = k2
1154 9049014 : mrgn = mrgn+1
1155 : region_found = .true.
1156 : exit
1157 : end if
1158 : end do
1159 : else
1160 : exit
1161 : endif
1162 :
1163 5880600 : if (region_found) then
1164 : !
1165 : ! Sort cloud areas and corresponding level indices.
1166 : !
1167 9049014 : nxs = 0
1168 : if (cldeps > 0) then
1169 : do k = k1,k2
1170 : if (cld(i,k) >= cldmin .and. cld(i,k) >= cldeps) then
1171 : nxs = nxs+1
1172 : ksort(nxs) = k
1173 : !
1174 : ! We need indices for clouds in order of largest to smallest, so
1175 : ! sort 1-cld in ascending order
1176 : !
1177 : asort(nxs) = 1.0_r8-(floor(cld(i,k)/cldeps)*cldeps)
1178 : end if
1179 : end do
1180 : else
1181 161944614 : do k = k1,k2
1182 161944614 : if (cld(i,k) >= cldmin) then
1183 38724092 : nxs = nxs+1
1184 38724092 : ksort(nxs) = k
1185 : !
1186 : ! We need indices for clouds in order of largest to smallest, so
1187 : ! sort 1-cld in ascending order
1188 : !
1189 38724092 : asort(nxs) = 1.0_r8-cld(i,k)
1190 : end if
1191 : end do
1192 : endif
1193 : !
1194 : ! If nxs eq 1, no need to sort.
1195 : ! If nxs eq 2, sort by swapping if necessary
1196 : ! If nxs ge 3, sort using local sort routine
1197 : !
1198 9049014 : if (nxs == 2) then
1199 2052291 : if (asort(2) < asort(1)) then
1200 1550125 : ktmp = ksort(1)
1201 1550125 : ksort(1) = ksort(2)
1202 1550125 : ksort(2) = ktmp
1203 :
1204 1550125 : atmp = asort(1)
1205 1550125 : asort(1) = asort(2)
1206 1550125 : asort(2) = atmp
1207 : endif
1208 6996723 : else if (nxs >= 3) then
1209 5534586 : call quick_sort(asort(1:nxs),ksort(1:nxs))
1210 : endif
1211 : !
1212 : ! Construct wstr, cstr, nstr for this region
1213 : !
1214 953873116 : cstr(k1:k2,1:nxs+1) = 0
1215 : mstr = 1
1216 : cld0 = 0.0_r8
1217 47773106 : do l = 1, nxs
1218 38724092 : if (asort(l) /= cld0) then
1219 37598493 : wstr(mstr,mrgn) = asort(l) - cld0
1220 37598493 : cld0 = asort(l)
1221 37598493 : mstr = mstr + 1
1222 : endif
1223 207609127 : cstr(ksort(l),mstr:nxs+1) = 1
1224 : end do
1225 9049014 : nstr(mrgn) = mstr
1226 9049014 : wstr(mstr,mrgn) = 1.0_r8 - cld0
1227 : !
1228 : ! End test of region_found = true
1229 : !
1230 : endif
1231 : !
1232 : ! End loop over regions irgn for max-overlap
1233 : !
1234 : end do
1235 5880600 : nrgn = mrgn
1236 : !
1237 : ! Finish construction of cstr for additional top layer
1238 : !
1239 37613302 : cstr(0,1:nstr(1)) = 0
1240 : !
1241 : ! INDEX COMPUTATIONS FOR STEP 2-3
1242 : ! This section of the code generates the following information:
1243 : ! (1. totwgt step 3 total frac. area of configurations satisfying
1244 : ! areamin & nconfgmax criteria
1245 : ! (2. wgtv step 3 frac. area of configurations
1246 : ! (3. ccon step 2 binary flag for clouds in each configuration
1247 : ! (4. nconfig steps 2-3 number of configurations
1248 : ! (5. nuniqu/d step 2 Number of unique cloud configurations for
1249 : ! up/downwelling rad. between surface/TOA
1250 : ! and level k
1251 : ! (6. istrtu/d step 2 Indices into iconu/d
1252 : ! (7. iconu/d step 2 Cloud configurations which are identical
1253 : ! for up/downwelling rad. between surface/TOA
1254 : ! and level k
1255 : !
1256 : ! Number of configurations (all permutations of streams in each region)
1257 : !
1258 14929614 : nconfigm = product(nstr(1: nrgn))
1259 : !
1260 : ! Construction of totwgt, wgtv, ccon, nconfig
1261 : !
1262 14929614 : istr(1: nrgn) = 1
1263 5880600 : nconfig(i) = 0
1264 5880600 : totwgt(i) = 0.0_r8
1265 5880600 : new_term = .true.
1266 96641319 : do iconfig = 1, nconfigm
1267 : xwgt = 1.0_r8
1268 281580814 : do mrgn = 1, nrgn
1269 281580814 : xwgt = xwgt * wstr(istr(mrgn),mrgn)
1270 : end do
1271 90760719 : if (xwgt >= areamin) then
1272 38641753 : nconfig(i) = nconfig(i) + 1
1273 38641753 : if (nconfig(i) <= nconfgmax) then
1274 37842590 : j = nconfig(i)
1275 37842590 : ptrc(nconfig(i)) = nconfig(i)
1276 : else
1277 799163 : nconfig(i) = nconfgmax
1278 799163 : if (new_term) then
1279 666497 : min_idx = minloc(wgtv)
1280 666497 : j = min_idx(1)
1281 : endif
1282 799163 : if (wgtv(j) < xwgt) then
1283 593723 : totwgt(i) = totwgt(i) - wgtv(j)
1284 : new_term = .true.
1285 : else
1286 : new_term = .false.
1287 : endif
1288 : endif
1289 37842590 : if (new_term) then
1290 38436313 : wgtv(j) = xwgt
1291 38436313 : totwgt(i) = totwgt(i) + xwgt
1292 104373107 : do mrgn = 1, nrgn
1293 1103717245 : ccon(j,kx1(mrgn):kx2(mrgn),i) = cstr(kx1(mrgn):kx2(mrgn),istr(mrgn))
1294 : end do
1295 : endif
1296 : endif
1297 :
1298 90760719 : mrgn = nrgn
1299 90760719 : istr(mrgn) = istr(mrgn) + 1
1300 113795766 : do while (istr(mrgn) > nstr(mrgn) .and. mrgn > 1)
1301 17154447 : istr(mrgn) = 1
1302 17154447 : mrgn = mrgn - 1
1303 17154447 : istr(mrgn) = istr(mrgn) + 1
1304 : end do
1305 : !
1306 : ! End do iconfig = 1, nconfigm
1307 : !
1308 : end do
1309 : !
1310 : ! If totwgt(i) = 0 implement maximum overlap and make another pass
1311 : ! if totwgt(i) = 0 on this second pass then terminate.
1312 : !
1313 5880600 : if (totwgt(i) > 0._r8) then
1314 : exit
1315 : else
1316 0 : npasses = npasses + 1
1317 0 : if (npasses >= 2 ) then
1318 0 : write(iulog,*)'RADCSWMX: Maximum overlap of column ','failed'
1319 0 : call endrun('RADCSWMX')
1320 : endif
1321 0 : nmxrgn(i)=1
1322 0 : pmxrgn(i,1)=1.0e30_r8
1323 : end if
1324 : !
1325 : ! End npasses = 0, do
1326 : !
1327 : end do
1328 : !
1329 : ! Finish construction of ccon
1330 : !
1331 :
1332 5880600 : istrtd(2,0,i) = nconfig(i)+1
1333 5880600 : istrtu(2,pverp,i) = nconfig(i)+1
1334 :
1335 164656800 : do k = 1, pverp
1336 158776200 : km1 = k-1
1337 158776200 : nuniq = 0
1338 158776200 : istrtd(1,k,i) = 1
1339 541164606 : do l0 = 1, nuniqd(km1,i)
1340 382388406 : is0 = istrtd(l0,km1,i)
1341 382388406 : is1 = istrtd(l0+1,km1,i)-1
1342 382388406 : n0 = 0
1343 382388406 : n1 = 0
1344 1404138336 : do isn = is0, is1
1345 1021749930 : j = icond(isn,km1,i)
1346 1404138336 : if (ccon(j,k,i) == 0) then
1347 887016785 : n0 = n0 + 1
1348 887016785 : ptr0(n0) = j
1349 : else ! if (ccon(j,k,i) == 1) then
1350 134733145 : n1 = n1 + 1
1351 134733145 : ptr1(n1) = j
1352 : endif
1353 : end do
1354 382388406 : if (n0 > 0) then
1355 328633987 : nuniq = nuniq + 1
1356 328633987 : istrtd(nuniq+1,k,i) = istrtd(nuniq,k,i)+n0
1357 1215650772 : icond(istrtd(nuniq,k,i):istrtd(nuniq+1,k,i)-1,k,i) = ptr0(1:n0)
1358 : endif
1359 541164606 : if (n1 > 0) then
1360 85716409 : nuniq = nuniq + 1
1361 85716409 : istrtd(nuniq+1,k,i) = istrtd(nuniq,k,i)+n1
1362 220449554 : icond(istrtd(nuniq,k,i):istrtd(nuniq+1,k,i)-1,k,i) = ptr1(1:n1)
1363 : endif
1364 : end do
1365 164656800 : nuniqd(k,i) = nuniq
1366 : end do
1367 : !
1368 : ! Find 'transition point' in downward configurations where the number
1369 : ! of 'configurations' changes from 1. This is used to optimize the
1370 : ! construction of the upward configurations.
1371 : ! Note: k1 == transition point
1372 : !
1373 :
1374 69526800 : do k = pverp,0,-1
1375 69526800 : if ( nuniqd(k,i) == 1) then
1376 : k1 = k
1377 : exit
1378 : end if
1379 : end do
1380 :
1381 63880230 : do k = pver, k1+1,-1
1382 57999630 : kp1 = k+1
1383 57999630 : nuniq = 0
1384 57999630 : istrtu(1,k,i) = 1
1385 253767407 : do l0 = 1, nuniqu(kp1,i)
1386 195767777 : is0 = istrtu(l0,kp1,i)
1387 195767777 : is1 = istrtu(l0+1,kp1,i)-1
1388 195767777 : n0 = 0
1389 195767777 : n1 = 0
1390 643879007 : do isn = is0, is1
1391 448111230 : j = iconu(isn,kp1,i)
1392 643879007 : if (ccon(j,k,i) == 0) then
1393 314719835 : n0 = n0 + 1
1394 314719835 : ptr0(n0) = j
1395 : else ! if (ccon(j,k,i) == 1) then
1396 133391395 : n1 = n1 + 1
1397 133391395 : ptr1(n1) = j
1398 : endif
1399 : end do
1400 195767777 : if (n0 > 0) then
1401 145872508 : nuniq = nuniq + 1
1402 145872508 : istrtu(nuniq+1,k,i) = istrtu(nuniq,k,i)+n0
1403 460592343 : iconu(istrtu(nuniq,k,i):istrtu(nuniq+1,k,i)-1,k,i) = ptr0(1:n0)
1404 : endif
1405 253767407 : if (n1 > 0) then
1406 81857259 : nuniq = nuniq + 1
1407 81857259 : istrtu(nuniq+1,k,i) = istrtu(nuniq,k,i)+n1
1408 215248654 : iconu(istrtu(nuniq,k,i):istrtu(nuniq+1,k,i)-1,k,i) = ptr1(1:n1)
1409 : endif
1410 : end do
1411 63880230 : nuniqu(k,i) = nuniq
1412 : end do
1413 : !
1414 : ! Copy identical configurations from 'transition point' to surface.
1415 : !
1416 5880600 : k1 = min(pverp-1,k1)
1417 5880600 : nuniq = nuniqu(k1+1,i)
1418 106657170 : do k = k1,0,-1
1419 100776570 : nuniqu(k,i) = nuniq
1420 674415270 : iconu(1:nuniq,k,i) = iconu(1:nuniq,k1+1,i)
1421 781072440 : istrtu(1:nuniq+1,k,i) = istrtu(1:nuniq+1,k1+1,i)
1422 : end do
1423 :
1424 44114973 : v_wgtv(1:nconfig(i),i) = wgtv(1:nconfig(i))
1425 :
1426 : !
1427 : !----------------------------------------------------------------------
1428 : ! End of index calculations
1429 : !----------------------------------------------------------------------
1430 : !
1431 : ! End do i=1,Nday
1432 : !
1433 : end do
1434 : !CSD$ END PARALLEL
1435 :
1436 : !----------------------------------------------------------------------
1437 : ! Start of flux calculations
1438 : !----------------------------------------------------------------------
1439 : !
1440 : ! Initialize spectrally integrated totals:
1441 : !
1442 169746124 : totfld (1:Nday,0:pver) = 0.0_r8
1443 169746124 : totfldc(1:Nday,0:pver) = 0.0_r8
1444 169746124 : fswup (1:Nday,0:pver) = 0.0_r8
1445 169746124 : fswdn (1:Nday,0:pver) = 0.0_r8
1446 :
1447 6272383 : fswup (1:Nday,pverp) = 0.0_r8
1448 6272383 : fswdn (1:Nday,pverp) = 0.0_r8
1449 : !
1450 : ! Start spectral interval
1451 : !
1452 : !old do ns = 1,nswbands
1453 : !old wgtint = nirwgt(ns)
1454 :
1455 6272383 : do i=1,Nday
1456 :
1457 : !----------------------------------------------------------------------
1458 : ! STEP 2
1459 : !
1460 : !
1461 : ! Apply adding method to solve for radiative properties
1462 : !
1463 : ! first initialize the bulk properties at toa
1464 : !
1465 :
1466 : ! nswbands, 0:pverp, nconfgmax, pcols
1467 :
1468 762732400 : rdndif(:,0,1:nconfig(i),i) = 0.0_r8
1469 762732400 : exptdn(:,0,1:nconfig(i),i) = 1.0_r8
1470 763124183 : tdntot(:,0,1:nconfig(i),i) = 1.0_r8
1471 : !
1472 : ! End do i=1,Nday
1473 : !
1474 : end do
1475 : !
1476 : ! solve for properties involving downward propagation of radiation.
1477 : ! the bulk properties are:
1478 : !
1479 : ! (1. exptdn sol. beam dwn. trans from layers above
1480 : ! (2. rdndif ref to dif rad for layers above
1481 : ! (3. tdntot total trans for layers above
1482 : !
1483 :
1484 : !CSD$ PARALLEL DO PRIVATE( km1, is0, is1, j, jj, Ttdif, Trdif, Trdir, Ttdir, Texplay ) &
1485 : !CSD$ PRIVATE( xexpt, xrdnd, tdnmexp, ytdnd, yrdnd, rdenom, rdirexp, zexpt, zrdnd, ztdnt ) &
1486 : !CSD$ PRIVATE( i, k, l0, ns, isn )
1487 6272383 : do i = 1, Nday
1488 165048583 : do k = 1, pverp
1489 158776200 : km1 = k - 1
1490 547045206 : do l0 = 1, nuniqd(km1,i)
1491 382388406 : is0 = istrtd(l0,km1,i)
1492 382388406 : is1 = istrtd(l0+1,km1,i)-1
1493 :
1494 382388406 : j = icond(is0,km1,i)
1495 :
1496 : !
1497 : ! If cloud in layer, use cloudy layer radiative properties (ccon == 1)
1498 : ! If clear layer, use clear-sky layer radiative properties (ccon /= 1)
1499 : !
1500 382388406 : if ( ccon(j,km1,i) == 1 ) then
1501 1714328180 : Ttdif(:) = tdif(:,i,km1)
1502 1714328180 : Trdif(:) = rdif(:,i,km1)
1503 1714328180 : Trdir(:) = rdir(:,i,km1)
1504 1714328180 : Ttdir(:) = tdir(:,i,km1)
1505 1714328180 : Texplay(:) = explay(:,i,km1)
1506 : else
1507 5933439940 : Ttdif(:) = tdifc(:,i,km1)
1508 5933439940 : Trdif(:) = rdifc(:,i,km1)
1509 5933439940 : Trdir(:) = rdirc(:,i,km1)
1510 5933439940 : Ttdir(:) = tdirc(:,i,km1)
1511 5933439940 : Texplay(:) = explayc(:,i,km1)
1512 : end if
1513 :
1514 7647768120 : do ns = 1, nswbands
1515 7265379714 : xexpt = exptdn(ns,km1,j,i)
1516 7265379714 : xrdnd = rdndif(ns,km1,j,i)
1517 7265379714 : tdnmexp = tdntot(ns,km1,j,i) - xexpt
1518 :
1519 7265379714 : ytdnd = Ttdif(ns)
1520 7265379714 : yrdnd = Trdif(ns)
1521 :
1522 7265379714 : rdenom = 1._r8/(1._r8-yrdnd*xrdnd)
1523 7265379714 : rdirexp = Trdir(ns)*xexpt
1524 :
1525 7265379714 : zexpt = xexpt * Texplay(ns)
1526 7265379714 : zrdnd = yrdnd + xrdnd*(ytdnd**2)*rdenom
1527 : ztdnt = xexpt*Ttdir(ns) + ytdnd* &
1528 7265379714 : (tdnmexp + xrdnd*rdirexp)*rdenom
1529 :
1530 7265379714 : exptdn(ns,k,j,i) = zexpt
1531 7265379714 : rdndif(ns,k,j,i) = zrdnd
1532 7647768120 : tdntot(ns,k,j,i) = ztdnt
1533 : end do ! ns = 1, nswbands
1534 : !
1535 : ! If 2 or more configurations share identical properties at a given level k,
1536 : ! the properties (at level k) are computed once and copied to
1537 : ! all the configurations for efficiency.
1538 : !
1539 1180526130 : do isn = is0+1, is1
1540 639361524 : jj = icond(isn,km1,i)
1541 12787230480 : exptdn(:,k,jj,i) = exptdn(:,k,j,i)
1542 12787230480 : rdndif(:,k,jj,i) = rdndif(:,k,j,i)
1543 13169618886 : tdntot(:,k,jj,i) = tdntot(:,k,j,i)
1544 : end do
1545 :
1546 : !
1547 : ! end do l0 = 1, nuniqd(k,i)
1548 : !
1549 : end do
1550 : !
1551 : ! end do k = 1, pverp
1552 : !
1553 : end do
1554 : !
1555 : ! end do i = 1, Nday
1556 : !
1557 : end do
1558 : !CSD$ END PARALLEL
1559 : !
1560 : ! Solve for properties involving upward propagation of radiation.
1561 : ! The bulk properties are:
1562 : !
1563 : ! (1. rupdif Ref to dif rad for layers below
1564 : ! (2. rupdir Ref to dir rad for layers below
1565 : !
1566 : ! Specify surface boundary conditions (surface albedos)
1567 : !
1568 :
1569 : ! nswbands, 0:pverp, nconfgmax, pcols
1570 391783 : rupdir = 0._r8
1571 391783 : rupdif = 0._r8
1572 6272383 : do i = 1, Nday
1573 118003783 : do ns = 1, nswbands
1574 830740610 : rupdir(ns,pverp,1:nconfig(i),i) = albdir(i,ns)
1575 836621210 : rupdif(ns,pverp,1:nconfig(i),i) = albdif(i,ns)
1576 : end do
1577 : end do
1578 :
1579 6272383 : do i = 1, Nday
1580 165048583 : do k = pver, 0, -1
1581 966025267 : do l0 = 1, nuniqu(k,i)
1582 801368467 : is0 = istrtu(l0,k,i)
1583 801368467 : is1 = istrtu(l0+1,k,i)-1
1584 :
1585 801368467 : j = iconu(is0,k,i)
1586 :
1587 : !
1588 : ! If cloud in layer, use cloudy layer radiative properties (ccon == 1)
1589 : ! If clear layer, use clear-sky layer radiative properties (ccon /= 1)
1590 : !
1591 801368467 : if ( ccon(j,k,i) == 1 ) then
1592 1663980180 : Ttdif(:) = tdif(:,i,k)
1593 1663980180 : Trdif(:) = rdif(:,i,k)
1594 1663980180 : Trdir(:) = rdir(:,i,k)
1595 1663980180 : Ttdir(:) = tdir(:,i,k)
1596 1663980180 : Texplay(:) = explay(:,i,k)
1597 : else
1598 14363389160 : Ttdif(:) = tdifc(:,i,k)
1599 14363389160 : Trdif(:) = rdifc(:,i,k)
1600 14363389160 : Trdir(:) = rdirc(:,i,k)
1601 14363389160 : Ttdir(:) = tdirc(:,i,k)
1602 14363389160 : Texplay(:) = explayc(:,i,k)
1603 : end if
1604 :
1605 16027369340 : do ns = 1, nswbands
1606 15226000873 : xrupd = rupdif(ns,k+1,j,i)
1607 15226000873 : xrups = rupdir(ns,k+1,j,i)
1608 :
1609 : !
1610 : ! If cloud in layer, use cloudy layer radiative properties (ccon == 1)
1611 : ! If clear layer, use clear-sky layer radiative properties (ccon /= 1)
1612 : !
1613 15226000873 : yexpt = Texplay(ns)
1614 15226000873 : yrupd = Trdif(ns)
1615 15226000873 : ytupd = Ttdif(ns)
1616 :
1617 15226000873 : rdenom = 1._r8/( 1._r8 - yrupd*xrupd)
1618 15226000873 : tdnmexp = (Ttdir(ns)-yexpt)
1619 15226000873 : rdirexp = xrups*yexpt
1620 :
1621 15226000873 : zrupd = yrupd + xrupd*(ytupd**2)*rdenom
1622 15226000873 : zrups = Trdir(ns) + ytupd*(rdirexp + xrupd*tdnmexp)*rdenom
1623 :
1624 15226000873 : rupdif(ns,k,j,i) = zrupd
1625 16027369340 : rupdir(ns,k,j,i) = zrups
1626 : end do ! ns = 1, nswbands
1627 : !
1628 : ! If 2 or more configurations share identical properties at a given level k,
1629 : ! the properties (at level k) are computed once and copied to
1630 : ! all the configurations for efficiency.
1631 : !
1632 1180526130 : do isn = is0+1, is1
1633 220381463 : jj = iconu(isn,k,i)
1634 4407629260 : rupdif(:,k,jj,i) = rupdif(:,k,j,i)
1635 5208997727 : rupdir(:,k,jj,i) = rupdir(:,k,j,i)
1636 : end do
1637 :
1638 : !
1639 : ! end do l0 = 1, nuniqu(k,i)
1640 : !
1641 : end do
1642 : !
1643 : ! end do k = pver,0,-1
1644 : !
1645 : end do
1646 : !
1647 : ! end do i = 1, Nday
1648 : !
1649 : end do
1650 : !
1651 : !----------------------------------------------------------------------
1652 : !
1653 : ! STEP 3
1654 : !
1655 : ! Compute up and down fluxes for each interface k. This requires
1656 : ! adding up the contributions from all possible permutations
1657 : ! of streams in all max-overlap regions, weighted by the
1658 : ! product of the fractional areas of the streams in each region
1659 : ! (the random overlap assumption). The adding principle has been
1660 : ! used in step 2 to combine the bulk radiative properties
1661 : ! above and below the interface.
1662 : !
1663 :
1664 : !
1665 : ! Initialize the fluxes
1666 : !
1667 391783 : fluxup = 0.0_r8
1668 391783 : fluxdn = 0.0_r8
1669 :
1670 6272383 : do i = 1, Nday
1671 44114973 : do iconfig = 1, nconfig(i)
1672 37842590 : xwgt = v_wgtv(iconfig,i)
1673 :
1674 1103315710 : do k = 0, pverp
1675 21229692990 : do ns = 1, nswbands
1676 20132257880 : xexpt = exptdn(ns,k,iconfig,i)
1677 20132257880 : xtdnt = tdntot(ns,k,iconfig,i)
1678 20132257880 : xrdnd = rdndif(ns,k,iconfig,i)
1679 20132257880 : xrupd = rupdif(ns,k,iconfig,i)
1680 20132257880 : xrups = rupdir(ns,k,iconfig,i)
1681 : !
1682 : ! Flux computation
1683 : !
1684 20132257880 : rdenom = 1._r8/(1._r8 - xrdnd * xrupd)
1685 :
1686 : fluxup(ns,k,i) = fluxup(ns,k,i) + xwgt * &
1687 20132257880 : ((xexpt * xrups + (xtdnt - xexpt) * xrupd) * rdenom)
1688 : fluxdn(ns,k,i) = fluxdn(ns,k,i) + xwgt * &
1689 21191850400 : (xexpt + (xtdnt - xexpt + xexpt * xrups * xrdnd) * rdenom)
1690 : end do ! do ns = 1, nswbands
1691 : end do
1692 : !
1693 : ! End do iconfig = 1, nconfig(i)
1694 : !
1695 : end do
1696 : !
1697 : ! End do iconfig = 1, Nday
1698 : !
1699 : end do
1700 : !
1701 : ! Normalize by total area covered by cloud configurations included
1702 : ! in solution
1703 : !
1704 6272383 : do i = 1, Nday
1705 170929183 : do k = 0, pverp
1706 3299016600 : do ns = 1, nswbands
1707 3128479200 : fluxup(ns,k,i)=fluxup(ns,k,i) / totwgt(i)
1708 3293136000 : fluxdn(ns,k,i)=fluxdn(ns,k,i) / totwgt(i)
1709 : end do ! do i = 1, nday
1710 : end do ! do k = 0, pverp
1711 : end do ! do i = 1, nday
1712 :
1713 :
1714 : !
1715 : ! Initialize the direct-beam flux at surface
1716 : !
1717 118003783 : wexptdn(:,1:Nday) = 0.0_r8
1718 :
1719 7835660 : do ns = 1,nswbands
1720 7443877 : wgtint = nirwgt(ns)
1721 :
1722 :
1723 119175277 : do i=1,Nday
1724 838184487 : do iconfig = 1, nconfig(i)
1725 : !
1726 : ! Note: exptdn can be directly indexed by iconfig at k=pverp.
1727 : !
1728 830740610 : wexptdn(ns,i) = wexptdn(ns,i) + v_wgtv(iconfig,i) * exptdn(ns,pverp,iconfig,i)
1729 : end do
1730 : end do
1731 :
1732 119175277 : do i=1,Nday
1733 111731400 : wexptdn(ns,i) = wexptdn(ns,i) / totwgt(i)
1734 : !
1735 : ! Monochromatic computation completed; accumulate in totals
1736 : !
1737 111731400 : if ( do_spctrl_scaling ) then
1738 0 : solflx(i) = solin(i)*frcsol(ns)*psf(ns)*sfac(ns)
1739 : else
1740 111731400 : solflx(i) = solin(i)*frcsol(ns)*psf(ns)
1741 : endif
1742 111731400 : fsnt(i) = fsnt(i) + solflx(i)*(fluxdn(ns,1,i) - fluxup(ns,1,i))
1743 111731400 : fsntoa(i)= fsntoa(i) + solflx(i)*(fluxdn(ns,0,i) - fluxup(ns,0,i))
1744 111731400 : fsutoa(i)= fsutoa(i) + solflx(i)*(fluxup(ns,0,i))
1745 111731400 : fsns(i) = fsns(i) + solflx(i)*(fluxdn(ns,pverp,i)-fluxup(ns,pverp,i))
1746 111731400 : fsdtoa(i) = fsdtoa(i) + solflx(i)
1747 111731400 : fswup(i,0) = fswup(i,0) + solflx(i)*fluxup(ns,0,i)
1748 111731400 : fswdn(i,0) = fswdn(i,0) + solflx(i)*fluxdn(ns,0,i)
1749 : !
1750 : ! Down spectral fluxes need to be in mks; thus the .001 conversion factors
1751 : !
1752 111731400 : if (wavmid(ns) < 0.7_r8) then
1753 52925400 : sols(i) = sols(i) + wexptdn(ns,i)*solflx(i)*0.001_r8
1754 52925400 : solsd(i) = solsd(i)+(fluxdn(ns,pverp,i)-wexptdn(ns,i))*solflx(i)*0.001_r8
1755 : else
1756 58806000 : soll(i) = soll(i) + wexptdn(ns,i)*solflx(i)*0.001_r8
1757 58806000 : solld(i) = solld(i)+(fluxdn(ns,pverp,i)-wexptdn(ns,i))*solflx(i)*0.001_r8
1758 58806000 : fsnrtoaq(i) = fsnrtoaq(i) + solflx(i)*(fluxdn(ns,0,i) - fluxup(ns,0,i))
1759 : end if
1760 119175277 : fsnirtoa(i) = fsnirtoa(i) + wgtint*solflx(i)*(fluxdn(ns,0,i) - fluxup(ns,0,i))
1761 :
1762 : !
1763 : ! End do i=1,Nday
1764 : !
1765 : end do
1766 :
1767 :
1768 208428556 : do k=0,pver
1769 3225176356 : do i=1,Nday
1770 : !
1771 : ! Compute flux divergence in each layer using the interface up and down
1772 : ! fluxes:
1773 : !
1774 3016747800 : kp1 = k+1
1775 3016747800 : flxdiv = (fluxdn(ns,k,i) - fluxdn(ns,kp1,i)) + (fluxup(ns,kp1,i) - fluxup(ns,k,i))
1776 3016747800 : totfld(i,k) = totfld(i,k) + solflx(i)*flxdiv
1777 3016747800 : fswdn(i,kp1) = fswdn(i,kp1) + solflx(i)*fluxdn(ns,kp1,i)
1778 3016747800 : fswup(i,kp1) = fswup(i,kp1) + solflx(i)*fluxup(ns,kp1,i)
1779 3016747800 : fns(i,kp1) = fswdn(i,kp1) - fswup(i,kp1)
1780 3217732479 : if (single_column.and.scm_crm_mode) then
1781 0 : fus(i,kp1)=fswup(i,kp1)
1782 0 : fds(i,kp1)=fswdn(i,kp1)
1783 : endif
1784 : end do
1785 : end do
1786 : !
1787 : ! Perform clear-sky calculation
1788 : !
1789 :
1790 119175277 : exptdnc(1:Nday,0) = 1.0_r8
1791 119175277 : rdndifc(1:Nday,0) = 0.0_r8
1792 119175277 : tdntotc(1:Nday,0) = 1.0_r8
1793 119175277 : rupdirc(1:Nday,pverp) = albdir(1:Nday,ns)
1794 119175277 : rupdifc(1:Nday,pverp) = albdif(1:Nday,ns)
1795 :
1796 :
1797 208428556 : do k = 1, pverp
1798 3225176356 : do i=1,Nday
1799 3016747800 : km1 = k - 1
1800 3016747800 : xexpt = exptdnc(i,km1)
1801 3016747800 : xrdnd = rdndifc(i,km1)
1802 3016747800 : yrdnd = rdifc(ns,i,km1)
1803 3016747800 : ytdnd = tdifc(ns,i,km1)
1804 :
1805 3016747800 : exptdnc(i,k) = xexpt*explayc(ns,i,km1)
1806 :
1807 3016747800 : rdenom = 1._r8/(1._r8 - yrdnd*xrdnd)
1808 3016747800 : rdirexp = rdirc(ns,i,km1)*xexpt
1809 3016747800 : tdnmexp = tdntotc(i,km1) - xexpt
1810 :
1811 : tdntotc(i,k) = xexpt*tdirc(ns,i,km1) + ytdnd*(tdnmexp + xrdnd*rdirexp)* &
1812 3016747800 : rdenom
1813 3217732479 : rdndifc(i,k) = yrdnd + xrdnd*(ytdnd**2)*rdenom
1814 : !
1815 : ! End do i=1,Nday
1816 : !
1817 : end do
1818 : end do
1819 :
1820 208428556 : do k=pver,0,-1
1821 3225176356 : do i=1,Nday
1822 3016747800 : xrupd = rupdifc(i,k+1)
1823 3016747800 : yexpt = explayc(ns,i,k)
1824 3016747800 : yrupd = rdifc(ns,i,k)
1825 3016747800 : ytupd = tdifc(ns,i,k)
1826 :
1827 3016747800 : rdenom = 1._r8/( 1._r8 - yrupd*xrupd)
1828 :
1829 3016747800 : rupdirc(i,k) = rdirc(ns,i,k) + ytupd*(rupdirc(i,k+1)*yexpt + &
1830 3016747800 : xrupd*(tdirc(ns,i,k)-yexpt))*rdenom
1831 3217732479 : rupdifc(i,k) = yrupd + xrupd*ytupd**2*rdenom
1832 : !
1833 : ! End do i=1,Nday
1834 : !
1835 : end do
1836 : end do
1837 :
1838 215872433 : do k=0,pverp
1839 3344351633 : do i=1,Nday
1840 3128479200 : rdenom = 1._r8/(1._r8 - rdndifc(i,k)*rupdifc(i,k))
1841 : fluxup(ns,k,i) = (exptdnc(i,k)*rupdirc(i,k) + (tdntotc(i,k)-exptdnc(i,k))*rupdifc(i,k))* &
1842 3128479200 : rdenom
1843 : fluxdn(ns,k,i) = exptdnc(i,k) + &
1844 : (tdntotc(i,k) - exptdnc(i,k) + exptdnc(i,k)*rupdirc(i,k)*rdndifc(i,k))* &
1845 3336907756 : rdenom
1846 : !
1847 : ! End do i=1,Nday
1848 : !
1849 : end do
1850 : end do
1851 :
1852 208428556 : do k=0,pver
1853 3225176356 : do i=1,Nday
1854 : !
1855 : ! Compute flux divergence in each layer using the interface up and down
1856 : ! fluxes:
1857 : !
1858 3016747800 : kp1 = k+1
1859 3016747800 : flxdiv = (fluxdn(ns,k,i) - fluxdn(ns,kp1,i)) + (fluxup(ns,kp1,i) - fluxup(ns,k,i))
1860 3217732479 : totfldc(i,k) = totfldc(i,k) + solflx(i)*flxdiv
1861 : end do
1862 : end do
1863 :
1864 119175277 : do i=1,Nday
1865 111731400 : fsntc(i) = fsntc(i)+solflx(i)*(fluxdn(ns,1,i)-fluxup(ns,1,i))
1866 111731400 : fsntoac(i) = fsntoac(i)+solflx(i)*(fluxdn(ns,0,i)-fluxup(ns,0,i))
1867 111731400 : fsnsc(i) = fsnsc(i)+solflx(i)*(fluxdn(ns,pverp,i)-fluxup(ns,pverp,i))
1868 111731400 : fsdsc(i) = fsdsc(i)+solflx(i)*(fluxdn(ns,pverp,i))
1869 111731400 : fsnrtoac(i) = fsnrtoac(i)+wgtint*solflx(i)*(fluxdn(ns,0,i)-fluxup(ns,0,i))
1870 119175277 : if (single_column.and.scm_crm_mode) then
1871 0 : do k = 1,pverp
1872 0 : fusc(i,k)=fusc(i,k) + solflx(i) * fluxup(ns,k,i)
1873 0 : fdsc(i,k)=fdsc(i,k) + solflx(i) * fluxdn(ns,k,i)
1874 : enddo
1875 : endif
1876 : !
1877 : ! End do i=1,Nday
1878 : !
1879 : end do
1880 :
1881 208820339 : do k = 1,pverp
1882 3225176356 : do i=1,Nday
1883 3217732479 : fcns(i,k)=fcns(i,k) + solflx(i)*(fluxdn(ns,k,i)-fluxup(ns,k,i))
1884 : enddo
1885 : enddo
1886 : !
1887 : ! End of clear sky calculation
1888 : !
1889 :
1890 : !
1891 : ! End of spectral interval loop
1892 : !
1893 : end do
1894 :
1895 6272383 : do i=1,Nday
1896 :
1897 : !
1898 : ! Compute solar heating rate (J/kg/s)
1899 : !
1900 158776200 : do k=1,pver
1901 152895600 : qrs(i,k) = -1.E-4_r8*gravit*totfld(i,k)/(pint(i,k) - pint(i,k+1))
1902 158776200 : qrsc(i,k) = -1.E-4_r8*gravit*totfldc(i,k)/(pint(i,k) - pint(i,k+1))
1903 : end do
1904 : !
1905 : ! Set the downwelling flux at the surface
1906 : !
1907 6272383 : fsds(i) = fswdn(i,pverp)
1908 : !
1909 : ! End do i=1,Nday
1910 : !
1911 : end do
1912 : !
1913 : ! Rearrange output arrays.
1914 : !
1915 : ! intent(inout)
1916 : !
1917 391783 : call ExpDayNite(pmxrgn, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1918 391783 : call ExpDayNite(nmxrgn, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1919 : !
1920 : ! intent(out)
1921 : !
1922 391783 : call ExpDayNite(solin, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1923 391783 : call ExpDayNite(qrs, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
1924 391783 : call ExpDayNite(qrsc, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
1925 391783 : call ExpDayNite(fns, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1926 391783 : call ExpDayNite(fcns, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1927 391783 : call ExpDayNite(fsns, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1928 391783 : call ExpDayNite(fsnt, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1929 391783 : call ExpDayNite(fsntoa, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1930 391783 : call ExpDayNite(fsutoa, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1931 391783 : call ExpDayNite(fsds, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1932 391783 : call ExpDayNite(fsnsc, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1933 391783 : call ExpDayNite(fsdsc, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1934 391783 : call ExpDayNite(fsntc, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1935 391783 : call ExpDayNite(fsdtoa, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1936 391783 : call ExpDayNite(fsntoac, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1937 391783 : call ExpDayNite(sols, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1938 391783 : call ExpDayNite(soll, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1939 391783 : call ExpDayNite(solsd, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1940 391783 : call ExpDayNite(solld, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1941 391783 : call ExpDayNite(fsnirtoa, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1942 391783 : call ExpDayNite(fsnrtoac, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1943 391783 : call ExpDayNite(fsnrtoaq, Nday, IdxDay, Nnite, IdxNite, 1, pcols)
1944 391783 : call ExpDayNite(tauxcl_out, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
1945 391783 : call ExpDayNite(tauxci_out, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
1946 :
1947 : ! these outfld calls don't work for spmd only outfield in scm mode (nonspmd)
1948 391783 : if (single_column.and.scm_crm_mode) then
1949 : ! Following outputs added for CRM
1950 0 : call ExpDayNite(fus, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1951 0 : call ExpDayNite(fds, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1952 0 : call ExpDayNite(fusc, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1953 0 : call ExpDayNite(fdsc, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
1954 0 : call outfld('FUS ', fus*1.e-3_r8, pcols, lchnk)
1955 0 : call outfld('FDS ', fds*1.e-3_r8, pcols, lchnk)
1956 0 : call outfld('FUSC ', fusc*1.e-3_r8, pcols, lchnk)
1957 0 : call outfld('FDSC ', fdsc*1.e-3_r8, pcols, lchnk)
1958 : endif
1959 : ! write(iulog, '(a, x, i3)') 'radcswmx : exiting, chunk identifier', lchnk
1960 :
1961 : return
1962 749232 : end subroutine radcswmx
1963 :
1964 : !-------------------------------------------------------------------------------
1965 :
1966 7443877 : subroutine raddedmx(coszrs ,ndayc ,abh2o , &
1967 : abo3 ,abco2 ,abo2 ,uh2o ,uo3 , &
1968 : uco2 ,uo2 ,trayoslp,pflx ,ns , &
1969 : tauxcl ,wcl ,gcl ,fcl ,tauxci , &
1970 : wci ,gci ,fci ,aer_tau ,aer_tau_w, &
1971 : aer_tau_w_g, aer_tau_w_f ,rdir ,rdif ,tdir , &
1972 : tdif ,explay ,rdirc ,rdifc ,tdirc , &
1973 : tdifc ,explayc )
1974 : !-----------------------------------------------------------------------
1975 : !
1976 : ! Purpose:
1977 : ! Computes layer reflectivities and transmissivities, from the top down
1978 : ! to the surface using the delta-Eddington solutions for each layer
1979 : !
1980 : ! Method:
1981 : ! For more details , see Briegleb, Bruce P., 1992: Delta-Eddington
1982 : ! Approximation for Solar Radiation in the NCAR Community Climate Model,
1983 : ! Journal of Geophysical Research, Vol 97, D7, pp7603-7612).
1984 : !
1985 : ! Modified for maximum/random cloud overlap by Bill Collins and John
1986 : ! Truesdale
1987 : !
1988 : ! Author: Bill Collins
1989 : !
1990 : !-----------------------------------------------------------------------
1991 :
1992 : !
1993 : ! Minimum total transmission below which no layer computation are done:
1994 : !
1995 : real(r8) trmin ! Minimum total transmission allowed
1996 : real(r8) wray ! Rayleigh single scatter albedo
1997 : real(r8) gray ! Rayleigh asymmetry parameter
1998 : real(r8) fray ! Rayleigh forward scattered fraction
1999 :
2000 : parameter (trmin = 1.e-3_r8)
2001 : parameter (wray = 0.999999_r8)
2002 : parameter (gray = 0.0_r8)
2003 : parameter (fray = 0.1_r8)
2004 : !
2005 : !------------------------------Arguments--------------------------------
2006 : !
2007 : ! Input arguments
2008 : !
2009 : real(r8), intent(in) :: coszrs(pcols) ! Cosine zenith angle
2010 : real(r8), intent(in) :: trayoslp ! Tray/sslp
2011 : real(r8), intent(in) :: pflx(pcols,0:pverp) ! Interface pressure
2012 : real(r8), intent(in) :: abh2o ! Absorption coefficiant for h2o
2013 : real(r8), intent(in) :: abo3 ! Absorption coefficiant for o3
2014 : real(r8), intent(in) :: abco2 ! Absorption coefficiant for co2
2015 : real(r8), intent(in) :: abo2 ! Absorption coefficiant for o2
2016 : real(r8), intent(in) :: uh2o(pcols,0:pver) ! Layer absorber amount of h2o
2017 : real(r8), intent(in) :: uo3(pcols,0:pver) ! Layer absorber amount of o3
2018 : real(r8), intent(in) :: uco2(pcols,0:pver) ! Layer absorber amount of co2
2019 : real(r8), intent(in) :: uo2(pcols,0:pver) ! Layer absorber amount of o2
2020 : real(r8), intent(in) :: tauxcl(pcols,0:pver) ! Cloud extinction optical depth (liquid)
2021 : real(r8), intent(in) :: wcl(pcols,0:pver) ! Cloud single scattering albedo (liquid)
2022 : real(r8), intent(in) :: gcl(pcols,0:pver) ! Cloud asymmetry parameter (liquid)
2023 : real(r8), intent(in) :: fcl(pcols,0:pver) ! Cloud forward scattered fraction (liquid)
2024 : real(r8), intent(in) :: tauxci(pcols,0:pver) ! Cloud extinction optical depth (ice)
2025 : real(r8), intent(in) :: wci(pcols,0:pver) ! Cloud single scattering albedo (ice)
2026 : real(r8), intent(in) :: gci(pcols,0:pver) ! Cloud asymmetry parameter (ice)
2027 : real(r8), intent(in) :: fci(pcols,0:pver) ! Cloud forward scattered fraction (ice)
2028 : real(r8), intent(inout) :: aer_tau(pcols,0:pver) ! Aerosol extinction optical depth
2029 : real(r8), intent(inout) :: aer_tau_w(pcols,0:pver) ! Aerosol single scattering albedo * tau
2030 : real(r8), intent(inout) :: aer_tau_w_g(pcols,0:pver) ! Aerosol asymmetry parameter * w * t
2031 : real(r8), intent(inout) :: aer_tau_w_f(pcols,0:pver) ! Aerosol forward scattered fraction * w * tau
2032 :
2033 : integer, intent(in) :: ndayc ! Number of daylight columns
2034 : integer, intent(in) :: ns ! Index of spectral interval
2035 : !
2036 : ! Input/Output arguments
2037 : !
2038 : ! Following variables are defined for each layer; 0 refers to extra
2039 : ! layer above top of model:
2040 : !
2041 : real(r8), intent(inout) :: rdir(nswbands,pcols,0:pver) ! Layer reflectivity to direct rad
2042 : real(r8), intent(inout) :: rdif(nswbands,pcols,0:pver) ! Layer reflectivity to diffuse rad
2043 : real(r8), intent(inout) :: tdir(nswbands,pcols,0:pver) ! Layer transmission to direct rad
2044 : real(r8), intent(inout) :: tdif(nswbands,pcols,0:pver) ! Layer transmission to diffuse rad
2045 : real(r8), intent(inout) :: explay(nswbands,pcols,0:pver) ! Solar beam exp transm for layer
2046 : !
2047 : ! Corresponding quantities for clear-skies
2048 : !
2049 : real(r8), intent(inout) :: rdirc(nswbands,pcols,0:pver) ! Clear layer reflec. to direct rad
2050 : real(r8), intent(inout) :: rdifc(nswbands,pcols,0:pver) ! Clear layer reflec. to diffuse rad
2051 : real(r8), intent(inout) :: tdirc(nswbands,pcols,0:pver) ! Clear layer trans. to direct rad
2052 : real(r8), intent(inout) :: tdifc(nswbands,pcols,0:pver) ! Clear layer trans. to diffuse rad
2053 : real(r8), intent(inout) :: explayc(nswbands,pcols,0:pver)! Solar beam exp transm clear layer
2054 : !
2055 : !---------------------------Local variables-----------------------------
2056 : !
2057 : integer i ! Column indices
2058 : integer k ! Level index
2059 : integer nn ! Index of column loops (max=ndayc)
2060 :
2061 : real(r8) taugab(pcols) ! Layer total gas absorption optical depth
2062 : real(r8) tauray(pcols) ! Layer rayleigh optical depth
2063 : real(r8) taucsc ! Layer cloud scattering optical depth
2064 : real(r8) tautot ! Total layer optical depth
2065 : real(r8) wtot ! Total layer single scatter albedo
2066 : real(r8) gtot ! Total layer asymmetry parameter
2067 : real(r8) ftot ! Total layer forward scatter fraction
2068 : real(r8) wtau ! rayleigh layer scattering optical depth
2069 : real(r8) wt ! layer total single scattering albedo
2070 : real(r8) ts ! layer scaled extinction optical depth
2071 : real(r8) ws ! layer scaled single scattering albedo
2072 : real(r8) gs ! layer scaled asymmetry parameter
2073 : !
2074 : !---------------------------Statement functions-------------------------
2075 : !
2076 : ! Statement functions and other local variables
2077 : !
2078 : real(r8) alpha ! Term in direct reflect and transmissivity
2079 : real(r8) gamma ! Term in direct reflect and transmissivity
2080 : real(r8) el ! Term in alpha,gamma,n,u
2081 : real(r8) taus ! Scaled extinction optical depth
2082 : real(r8) omgs ! Scaled single particle scattering albedo
2083 : real(r8) asys ! Scaled asymmetry parameter
2084 : real(r8) u ! Term in diffuse reflect and
2085 : ! transmissivity
2086 : real(r8) n ! Term in diffuse reflect and
2087 : ! transmissivity
2088 : real(r8) lm ! Temporary for el
2089 : real(r8) ne ! Temporary for n
2090 : real(r8) w ! Dummy argument for statement function
2091 : real(r8) uu ! Dummy argument for statement function
2092 : real(r8) g ! Dummy argument for statement function
2093 : real(r8) e ! Dummy argument for statement function
2094 : real(r8) f ! Dummy argument for statement function
2095 : real(r8) t ! Dummy argument for statement function
2096 : real(r8) et ! Dummy argument for statement function
2097 : !
2098 : ! Intermediate terms for delta-eddington solution
2099 : !
2100 : real(r8) alp ! Temporary for alpha
2101 : real(r8) gam ! Temporary for gamma
2102 : real(r8) ue ! Temporary for u
2103 : real(r8) arg ! Exponential argument
2104 : real(r8) extins ! Extinction
2105 : real(r8) amg ! Alp - gam
2106 : real(r8) apg ! Alp + gam
2107 : !
2108 : ! ssa <=1 limit for aerosol
2109 : !
2110 : real(r8) :: w_limited(pcols,0:pver) ! Aerosol ssa (limited to < 0.999999)
2111 : real(r8) :: aer_g_limit(pcols,0:pver) ! Aerosol tau_w_g (limited ssa)
2112 : real(r8) :: aer_f_limit(pcols,0:pver) ! Aerosol tau_w_f (limited ssa)
2113 : !
2114 : alpha(w,uu,g,e) = .75_r8*w*uu*((1._r8 + g*(1._r8-w))/(1._r8 - e*e*uu*uu))
2115 : gamma(w,uu,g,e) = .50_r8*w*((3._r8*g*(1._r8-w)*uu*uu + 1._r8)/(1._r8-e*e*uu*uu))
2116 : el(w,g) = sqrt(3._r8*(1._r8-w)*(1._r8 - w*g))
2117 : taus(w,f,t) = (1._r8 - w*f)*t
2118 : omgs(w,f) = (1._r8 - f)*w/(1._r8 - w*f)
2119 : asys(g,f) = (g - f)/(1._r8 - f)
2120 : u(w,g,e) = 1.5_r8*(1._r8 - w*g)/e
2121 : n(uu,et) = ((uu+1._r8)*(uu+1._r8)/et ) - ((uu-1._r8)*(uu-1._r8)*et)
2122 : !
2123 : !-----------------------------------------------------------------------
2124 : !
2125 : ! Compute layer radiative properties
2126 : !
2127 : ! Compute radiative properties (reflectivity and transmissivity for
2128 : ! direct and diffuse radiation incident from above, under clear
2129 : ! and cloudy conditions) and transmission of direct radiation
2130 : ! (under clear and cloudy conditions) for each layer.
2131 : !
2132 208428556 : do k=0,pver
2133 3225176356 : do i=1,ndayc
2134 3016747800 : if(aer_tau(i,k) > 0._r8) then !where(aer_tau > 0._r8)
2135 0 : aer_g_limit(i,k) = aer_tau_w_g(i,k) / aer_tau_w(i,k)
2136 0 : aer_f_limit(i,k) = aer_tau_w_f(i,k) / aer_tau_w(i,k)
2137 0 : aer_tau_w(i,k) = aer_tau(i,k) * min(aer_tau_w(i,k)/aer_tau(i,k) , 0.999999_r8)
2138 : else
2139 3016747800 : aer_tau(i,k) = 0._r8
2140 3016747800 : aer_tau_w(i,k) = 0._r8
2141 3016747800 : aer_g_limit(i,k) = 0._r8
2142 3016747800 : aer_f_limit(i,k) = 0._r8
2143 : endif
2144 3016747800 : aer_tau_w_g(i,k) = aer_tau_w(i,k) * aer_g_limit(i,k)
2145 3217732479 : aer_tau_w_f(i,k) = aer_tau_w(i,k) * aer_f_limit(i,k)
2146 : enddo
2147 : enddo
2148 :
2149 208428556 : do k=0,pver
2150 3225176356 : do i=1,ndayc
2151 3016747800 : tauray(i) = trayoslp*(pflx(i,k+1)-pflx(i,k))
2152 3016747800 : taugab(i) = abh2o*uh2o(i,k) + abo3*uo3(i,k) + abco2*uco2(i,k) + abo2*uo2(i,k)
2153 3016747800 : tautot = tauxcl(i,k) + tauxci(i,k) + tauray(i) + taugab(i) + aer_tau(i,k)
2154 3016747800 : taucsc = tauxcl(i,k)*wcl(i,k) + tauxci(i,k)*wci(i,k) + aer_tau_w(i,k)
2155 3016747800 : wtau = wray*tauray(i)
2156 3016747800 : wt = wtau + taucsc
2157 3016747800 : wtot = wt/tautot
2158 : gtot = (wtau*gray + gcl(i,k)*wcl(i,k)*tauxcl(i,k) &
2159 3016747800 : + gci(i,k)*wci(i,k)*tauxci(i,k) + aer_tau_w_g(i,k))/wt
2160 : ftot = (wtau*fray + fcl(i,k)*wcl(i,k)*tauxcl(i,k) &
2161 3016747800 : + fci(i,k)*wci(i,k)*tauxci(i,k) + aer_tau_w_f(i,k))/wt
2162 3016747800 : ts = taus(wtot,ftot,tautot)
2163 3016747800 : ws = omgs(wtot,ftot)
2164 3016747800 : gs = asys(gtot,ftot)
2165 3016747800 : lm = el(ws,gs)
2166 3016747800 : alp = alpha(ws,coszrs(i),gs,lm)
2167 3016747800 : gam = gamma(ws,coszrs(i),gs,lm)
2168 3016747800 : ue = u(ws,gs,lm)
2169 : !
2170 : ! Limit argument of exponential to 25, in case lm very large:
2171 : !
2172 3016747800 : arg = min(lm*ts,25._r8)
2173 3016747800 : extins = exp(-arg)
2174 3016747800 : ne = n(ue,extins)
2175 3016747800 : rdif(ns,i,k) = (ue+1._r8)*(ue-1._r8)*(1._r8/extins - extins)/ne
2176 3016747800 : tdif(ns,i,k) = 4._r8*ue/ne
2177 : !
2178 : ! Limit argument of exponential to 25, in case coszrs is very small:
2179 : !
2180 3016747800 : arg = min(ts/coszrs(i),25._r8)
2181 3016747800 : explay(ns,i,k) = exp(-arg)
2182 3016747800 : apg = alp + gam
2183 3016747800 : amg = alp - gam
2184 3016747800 : rdir(ns,i,k) = amg*(tdif(ns,i,k)*explay(ns,i,k)-1._r8) + apg*rdif(ns,i,k)
2185 3016747800 : tdir(ns,i,k) = apg*tdif(ns,i,k) + (amg*rdif(ns,i,k)-(apg-1._r8))*explay(ns,i,k)
2186 : !
2187 : ! Under rare conditions, reflectivies and transmissivities can be
2188 : ! negative; zero out any negative values
2189 : !
2190 3016747800 : rdir(ns,i,k) = max(rdir(ns,i,k),0.0_r8)
2191 3016747800 : tdir(ns,i,k) = max(tdir(ns,i,k),0.0_r8)
2192 3016747800 : rdif(ns,i,k) = max(rdif(ns,i,k),0.0_r8)
2193 3016747800 : tdif(ns,i,k) = max(tdif(ns,i,k),0.0_r8)
2194 : !
2195 : ! Clear-sky calculation
2196 : !
2197 3217732479 : if (tauxcl(i,k) == 0.0_r8 .and. tauxci(i,k) == 0.0_r8) then
2198 :
2199 2281020946 : rdirc(ns,i,k) = rdir(ns,i,k)
2200 2281020946 : tdirc(ns,i,k) = tdir(ns,i,k)
2201 2281020946 : rdifc(ns,i,k) = rdif(ns,i,k)
2202 2281020946 : tdifc(ns,i,k) = tdif(ns,i,k)
2203 2281020946 : explayc(ns,i,k) = explay(ns,i,k)
2204 : else
2205 735726854 : tautot = tauray(i) + taugab(i) + aer_tau(i,k)
2206 735726854 : taucsc = aer_tau_w(i,k)
2207 : !
2208 : ! wtau already computed for all-sky
2209 : !
2210 735726854 : wt = wtau + taucsc
2211 735726854 : wtot = wt/tautot
2212 735726854 : gtot = (wtau*gray + aer_tau_w_g(i,k))/wt
2213 735726854 : ftot = (wtau*fray + aer_tau_w_f(i,k))/wt
2214 735726854 : ts = taus(wtot,ftot,tautot)
2215 735726854 : ws = omgs(wtot,ftot)
2216 735726854 : gs = asys(gtot,ftot)
2217 735726854 : lm = el(ws,gs)
2218 735726854 : alp = alpha(ws,coszrs(i),gs,lm)
2219 735726854 : gam = gamma(ws,coszrs(i),gs,lm)
2220 735726854 : ue = u(ws,gs,lm)
2221 : !
2222 : ! Limit argument of exponential to 25, in case lm very large:
2223 : !
2224 735726854 : arg = min(lm*ts,25._r8)
2225 735726854 : extins = exp(-arg)
2226 735726854 : ne = n(ue,extins)
2227 735726854 : rdifc(ns,i,k) = (ue+1._r8)*(ue-1._r8)*(1._r8/extins - extins)/ne
2228 735726854 : tdifc(ns,i,k) = 4._r8*ue/ne
2229 : !
2230 : ! Limit argument of exponential to 25, in case coszrs is very small:
2231 : !
2232 735726854 : arg = min(ts/coszrs(i),25._r8)
2233 735726854 : explayc(ns,i,k) = exp(-arg)
2234 735726854 : apg = alp + gam
2235 735726854 : amg = alp - gam
2236 : rdirc(ns,i,k) = amg*(tdifc(ns,i,k)*explayc(ns,i,k)-1._r8)+ &
2237 735726854 : apg*rdifc(ns,i,k)
2238 : tdirc(ns,i,k) = apg*tdifc(ns,i,k) + (amg*rdifc(ns,i,k) - (apg-1._r8))* &
2239 735726854 : explayc(ns,i,k)
2240 : !
2241 : ! Under rare conditions, reflectivies and transmissivities can be
2242 : ! negative; zero out any negative values
2243 : !
2244 735726854 : rdirc(ns,i,k) = max(rdirc(ns,i,k),0.0_r8)
2245 735726854 : tdirc(ns,i,k) = max(tdirc(ns,i,k),0.0_r8)
2246 735726854 : rdifc(ns,i,k) = max(rdifc(ns,i,k),0.0_r8)
2247 735726854 : tdifc(ns,i,k) = max(tdifc(ns,i,k),0.0_r8)
2248 : end if
2249 : end do
2250 : end do
2251 :
2252 749232 : end subroutine raddedmx
2253 :
2254 : !-------------------------------------------------------------------------------
2255 :
2256 1536 : subroutine radsw_init(gravx)
2257 : !-----------------------------------------------------------------------
2258 : !
2259 : ! Purpose:
2260 : ! Initialize various constants for radiation scheme; note that
2261 : ! the radiation scheme uses cgs units.
2262 : !
2263 : ! Author: W. Collins (H2O parameterization) and J. Kiehl
2264 : !
2265 : !-----------------------------------------------------------------------
2266 : !
2267 : ! Input arguments
2268 : !
2269 : real(r8), intent(in) :: gravx ! Acceleration of gravity (MKS)
2270 :
2271 : real(r8), parameter :: ref_tsi = 1367._r8 ! value supplied by Dan Marsh -- see solvar_woods.F90
2272 : !
2273 : !-----------------------------------------------------------------------
2274 : !
2275 : ! Set general radiation consts; convert to cgs units where appropriate:
2276 : !
2277 1536 : gravit = 100._r8*gravx
2278 1536 : rga = 1._r8/gravit
2279 1536 : sslp = 1.013250e6_r8
2280 :
2281 1536 : end subroutine radsw_init
2282 :
2283 : !-------------------------------------------------------------------------------
2284 :
2285 : end module radsw
|
