Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

radiation_tripleclouds_sw.F90

Last change on this file was 4946, checked in by idelkadi, 7 months ago
The addition of explicit loops with the "omp simd reduction" directive to solve the slowness problem linked to the "sum" command (svn4938 revesion) led to non-reproducibility in MPI and mixed MPI-OMP modes in the case of Tripleclouds and Mcica solvers. We return to the svn4848 versions of the radiation_tripleclouds_w.F90 and radiation_mcica_w.F90 routines.
File size: 27.2 KB

Line
1	! radiation_tripleclouds_sw.F90 - Shortwave "Tripleclouds" solver
2	!
3	! (C) Copyright 2016- ECMWF.
4	!
5	! This software is licensed under the terms of the Apache Licence Version 2.0
6	! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7	!
8	! In applying this licence, ECMWF does not waive the privileges and immunities
9	! granted to it by virtue of its status as an intergovernmental organisation
10	! nor does it submit to any jurisdiction.
11	!
12	! Author: Robin Hogan
13	! Email: r.j.hogan@ecmwf.int
14	!
15	! Modifications
16	! 2017-04-11 R. Hogan Receive albedos at g-points
17	! 2017-04-22 R. Hogan Store surface fluxes at all g-points
18	! 2017-10-23 R. Hogan Renamed single-character variables
19	! 2018-10-08 R. Hogan Call calc_region_properties
20	! 2019-01-02 R. Hogan Fixed problem of do_save_spectral_flux .and. .not. do_sw_direct
21	! 2020-09-18 R. Hogan Replaced some array expressions with loops for speed
22	! 2021-10-01 P. Ukkonen Performance optimizations: batched computations
23
24	module radiation_tripleclouds_sw
25
26	public
27
28	contains
29	! Provides elemental function "delta_eddington"
30	#include "radiation_delta_eddington.h"
31
32	! Small routine for scaling cloud optical depth in the cloudy
33	! regions
34	#include "radiation_optical_depth_scaling.h"
35
36	!---------------------------------------------------------------------
37	! This module contains just one subroutine, the shortwave
38	! "Tripleclouds" solver in which cloud inhomogeneity is treated by
39	! dividing each model level into three regions, one clear and two
40	! cloudy (with differing optical depth). This approach was described
41	! by Shonk and Hogan (2008).
42	subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
43	& config, single_level, cloud, &
44	& od, ssa, g, od_cloud, ssa_cloud, g_cloud, &
45	& albedo_direct, albedo_diffuse, incoming_sw, &
46	& flux)
47
48	use parkind1, only : jprb
49	use yomhook, only : lhook, dr_hook, jphook
50
51	! use radiation_io, only : nulout
52	use radiation_config, only : config_type, IPdfShapeGamma
53	use radiation_single_level, only : single_level_type
54	use radiation_cloud, only : cloud_type
55	use radiation_regions, only : calc_region_properties
56	use radiation_overlap, only : calc_overlap_matrices
57	use radiation_flux, only : flux_type, &
58	& indexed_sum, add_indexed_sum
59	use radiation_matrix, only : singlemat_x_vec
60	use radiation_two_stream, only : calc_ref_trans_sw
61
62	implicit none
63
64	! Number of regions
65	integer, parameter :: nregions = 3
66
67	! Inputs
68	integer, intent(in) :: nlev ! number of model levels
69	integer, intent(in) :: istartcol, iendcol ! range of columns to process
70	type(config_type), intent(in) :: config
71	type(single_level_type), intent(in) :: single_level
72	type(cloud_type), intent(in) :: cloud
73
74	! Gas and aerosol optical depth, single-scattering albedo and
75	! asymmetry factor at each shortwave g-point
76	! real(jprb), intent(in), dimension(istartcol:iendcol,nlev,config%n_g_sw) :: &
77	real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: &
78	& od, ssa, g
79
80	! Cloud and precipitation optical depth, single-scattering albedo and
81	! asymmetry factor in each shortwave band
82	real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
83	& od_cloud, ssa_cloud, g_cloud
84
85	! Optical depth, single scattering albedo and asymmetry factor in
86	! each g-point of each cloudy region including gas, aerosol and
87	! clouds
88	real(jprb), dimension(config%n_g_sw,2:nregions) &
89	& :: od_total, ssa_total, g_total
90
91	! Direct and diffuse surface albedos, and the incoming shortwave
92	! flux into a plane perpendicular to the incoming radiation at
93	! top-of-atmosphere in each of the shortwave g points
94	real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
95	& albedo_direct, albedo_diffuse, incoming_sw
96
97	! Output
98	type(flux_type), intent(inout):: flux
99
100	! Local variables
101
102	! The area fractions of each region
103	real(jprb) :: region_fracs(1:nregions,nlev,istartcol:iendcol)
104
105	! The scaling used for the optical depth in the cloudy regions
106	real(jprb) :: od_scaling(2:nregions,nlev,istartcol:iendcol)
107
108	! Directional overlap matrices defined at all layer interfaces
109	! including top-of-atmosphere and the surface
110	real(jprb), dimension(nregions,nregions,nlev+1, &
111	& istartcol:iendcol) :: u_matrix, v_matrix
112
113	! Diffuse reflection and transmission matrices in the cloudy
114	! regions of each layer
115	real(jprb), dimension(config%n_g_sw, 2:nregions, nlev) &
116	& :: reflectance, transmittance
117
118	! Terms translating the direct flux entering the layer from above
119	! to the reflected radiation exiting upwards (ref_dir) and the
120	! scattered radiation exiting downwards (trans_dir_diff), along with the
121	! direct unscattered transmission matrix (trans_dir_dir).
122	real(jprb), dimension(config%n_g_sw, 2:nregions, nlev) &
123	& :: ref_dir, trans_dir_diff, trans_dir_dir
124
125	! As above but for the clear regions; clear and cloudy layers are
126	! separated out so that calc_ref_trans_sw can be called on the
127	! entire clear-sky atmosphere at once
128	real(jprb), dimension(config%n_g_sw, nlev) &
129	& :: reflectance_clear, transmittance_clear, &
130	& ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear
131
132	! Total albedo of the atmosphere/surface just above a layer
133	! interface with respect to downwelling diffuse and direct
134	! (respectively) radiation at that interface, where level index =
135	! 1 corresponds to the top-of-atmosphere
136	real(jprb), dimension(config%n_g_sw, nregions, nlev+1) &
137	& :: total_albedo, total_albedo_direct
138
139	! ...equivalent values for clear-skies
140	real(jprb), dimension(config%n_g_sw, nlev+1) &
141	& :: total_albedo_clear, total_albedo_clear_direct
142
143	! Total albedo of the atmosphere just below a layer interface
144	real(jprb), dimension(config%n_g_sw, nregions) &
145	& :: total_albedo_below, total_albedo_below_direct
146
147	! Direct downwelling flux below and above an interface between
148	! layers into a plane perpendicular to the direction of the sun
149	real(jprb), dimension(config%n_g_sw, nregions) :: direct_dn
150	! Diffuse equivalents
151	real(jprb), dimension(config%n_g_sw, nregions) :: flux_dn, flux_up
152
153	! ...clear-sky equivalent (no distinction between "above/below")
154	real(jprb), dimension(config%n_g_sw) &
155	& :: direct_dn_clear, flux_dn_clear, flux_up_clear
156
157	! Clear-sky equivalent, but actually its reciprocal to replace
158	! some divisions by multiplications
159	real(jprb), dimension(config%n_g_sw, nregions) :: inv_denom
160
161	! Identify clear-sky layers, with pseudo layers for outer space
162	! and below the ground, both treated as single-region clear skies
163	logical :: is_clear_sky_layer(0:nlev+1)
164
165	! Scattering optical depth of gas+aerosol and of cloud
166	real(jprb) :: scat_od, scat_od_cloud
167
168	real(jprb) :: mu0
169
170	integer :: jcol, jlev, jg, jreg, iband, jreg2, ng
171
172	real(jphook) :: hook_handle
173
174	if (lhook) call dr_hook('radiation_tripleclouds_sw:solver_tripleclouds_sw',0,hook_handle)
175
176	! --------------------------------------------------------
177	! Section 1: Prepare general variables and arrays
178	! --------------------------------------------------------
179	! Copy array dimensions to local variables for convenience
180	ng = config%n_g_sw
181
182	! Compute the wavelength-independent region fractions and
183	! optical-depth scalings
184	call calc_region_properties(nlev,nregions,istartcol,iendcol, &
185	& config%i_cloud_pdf_shape == IPdfShapeGamma, &
186	& cloud%fraction, cloud%fractional_std, region_fracs, &
187	& od_scaling, config%cloud_fraction_threshold)
188
189	! Compute wavelength-independent overlap matrices u_matrix and
190	! v_matrix
191	call calc_overlap_matrices(nlev,nregions,istartcol,iendcol, &
192	& region_fracs, cloud%overlap_param, &
193	& u_matrix, v_matrix, &
194	& decorrelation_scaling=config%cloud_inhom_decorr_scaling, &
195	& cloud_fraction_threshold=config%cloud_fraction_threshold, &
196	& use_beta_overlap=config%use_beta_overlap, &
197	& cloud_cover=flux%cloud_cover_sw)
198
199	! Main loop over columns
200	do jcol = istartcol, iendcol
201	! --------------------------------------------------------
202	! Section 2: Prepare column-specific variables and arrays
203	! --------------------------------------------------------
204
205	! Copy local cosine of the solar zenith angle
206	mu0 = single_level%cos_sza(jcol)
207
208	! Skip profile if sun is too low in the sky
209	if (mu0 < 1.0e-10_jprb) then
210	flux%sw_dn(jcol,:) = 0.0_jprb
211	flux%sw_up(jcol,:) = 0.0_jprb
212	if (allocated(flux%sw_dn_direct)) then
213	flux%sw_dn_direct(jcol,:) = 0.0_jprb
214	end if
215	if (config%do_clear) then
216	flux%sw_dn_clear(jcol,:) = 0.0_jprb
217	flux%sw_up_clear(jcol,:) = 0.0_jprb
218	if (allocated(flux%sw_dn_direct_clear)) then
219	flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
220	end if
221	end if
222
223	if (config%do_save_spectral_flux) then
224	flux%sw_dn_band(:,jcol,:) = 0.0_jprb
225	flux%sw_up_band(:,jcol,:) = 0.0_jprb
226	if (allocated(flux%sw_dn_direct_band)) then
227	flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb
228	end if
229	if (config%do_clear) then
230	flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb
231	flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb
232	if (allocated(flux%sw_dn_direct_clear_band)) then
233	flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb
234	end if
235	end if
236	end if
237
238	flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
239	flux%sw_dn_direct_surf_g(:,jcol) = 0.0_jprb
240	if (config%do_clear) then
241	flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
242	flux%sw_dn_direct_surf_clear_g(:,jcol) = 0.0_jprb
243	end if
244
245	cycle
246	end if ! sun is below the horizon
247
248	! At this point mu0 >= 1.0e-10
249
250	! Define which layers contain cloud; assume that
251	! cloud%crop_cloud_fraction has already been called
252	is_clear_sky_layer = .true.
253	do jlev = 1,nlev
254	if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
255	is_clear_sky_layer(jlev) = .false.
256	end if
257	end do
258
259	! --------------------------------------------------------
260	! Section 3: Loop over layers to compute reflectance and transmittance
261	! --------------------------------------------------------
262	! In this section the reflectance, transmittance and sources
263	! are computed for each layer
264
265	! Clear-sky quantities for all layers at once
266	call calc_ref_trans_sw(ng*nlev, &
267	& mu0, od(:,:,jcol), ssa(:,:,jcol), g(:,:,jcol), &
268	& reflectance_clear, transmittance_clear, &
269	& ref_dir_clear, trans_dir_diff_clear, &
270	& trans_dir_dir_clear)
271
272	! Cloudy layers
273	do jlev = 1,nlev ! Start at top-of-atmosphere
274	if (.not. is_clear_sky_layer(jlev)) then
275	do jreg = 2,nregions
276	do jg = 1,ng
277	! Mapping from g-point to band
278	iband = config%i_band_from_reordered_g_sw(jg)
279	scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol)
280	scat_od_cloud = od_cloud(iband,jlev,jcol) &
281	& * ssa_cloud(iband,jlev,jcol) * od_scaling(jreg,jlev,jcol)
282	! Add scaled cloud optical depth to clear-sky value
283	od_total(jg,jreg) = od(jg,jlev,jcol) &
284	& + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
285	! Compute single-scattering albedo and asymmetry
286	! factor of gas-cloud combination
287	ssa_total(jg,jreg) = (scat_od+scat_od_cloud) &
288	& / od_total(jg,jreg)
289	g_total(jg,jreg) = (scat_od*g(jg,jlev,jcol) &
290	& + scat_od_cloud * g_cloud(iband,jlev,jcol)) &
291	& / (scat_od + scat_od_cloud)
292	end do
293	end do
294
295	if (config%do_sw_delta_scaling_with_gases) then
296	! Apply delta-Eddington scaling to the aerosol-gas(-cloud)
297	! mixture
298	call delta_eddington(od_total, ssa_total, g_total)
299	end if
300	! Both cloudy regions at once
301	call calc_ref_trans_sw(ng*(nregions-1), &
302	& mu0, od_total, ssa_total, g_total, &
303	& reflectance(:,:,jlev), transmittance(:,:,jlev), &
304	& ref_dir(:,:,jlev), trans_dir_diff(:,:,jlev), &
305	& trans_dir_dir(:,:,jlev))
306	end if
307	end do
308
309	! --------------------------------------------------------
310	! Section 4: Compute total albedos
311	! --------------------------------------------------------
312
313	total_albedo(:,:,:) = 0.0_jprb
314	total_albedo_direct(:,:,:) = 0.0_jprb
315
316	! Copy surface albedos in clear-sky region
317	do jg = 1,ng
318	total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol)
319	total_albedo_direct(jg,1,nlev+1) &
320	& = mu0 * albedo_direct(jg,jcol)
321	end do
322
323	! If there is cloud in the lowest layer then we need the albedos
324	! underneath
325	if (.not. is_clear_sky_layer(nlev)) then
326	do jreg = 2,nregions
327	total_albedo(:,jreg,nlev+1) = total_albedo(:,1,nlev+1)
328	total_albedo_direct(:,jreg,nlev+1) = total_albedo_direct(:,1,nlev+1)
329	end do
330	end if
331
332	if (config%do_clear) then
333	total_albedo_clear(:,nlev+1) = total_albedo(:,1,nlev+1)
334	total_albedo_clear_direct(:,nlev+1) = total_albedo_direct(:,1,nlev+1)
335	end if
336
337	! Work up from the surface computing the total albedo of the
338	! atmosphere below that point using the adding method
339	do jlev = nlev,1,-1
340
341	total_albedo_below = 0.0_jprb
342	total_albedo_below_direct = 0.0_jprb
343
344	if (config%do_clear) then
345	! For clear-skies there is no need to consider "above" and
346	! "below" quantities since with no cloud overlap to worry
347	! about, these are the same
348	do jg = 1,ng
349	inv_denom(jg,1) = 1.0_jprb &
350	& / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance_clear(jg,jlev))
351	total_albedo_clear(jg,jlev) = reflectance_clear(jg,jlev) &
352	& + transmittance_clear(jg,jlev) * transmittance_clear(jg,jlev) &
353	& * total_albedo_clear(jg,jlev+1) * inv_denom(jg,1)
354
355	total_albedo_clear_direct(jg,jlev) = ref_dir_clear(jg,jlev) &
356	& + (trans_dir_dir_clear(jg,jlev)*total_albedo_clear_direct(jg,jlev+1) &
357	& +trans_dir_diff_clear(jg,jlev)*total_albedo_clear(jg,jlev+1)) &
358	& * transmittance_clear(jg,jlev) * inv_denom(jg,1)
359	end do
360	end if
361
362	! All-sky fluxes: first the clear region
363	do jg = 1,ng
364	inv_denom(jg,1) = 1.0_jprb &
365	& / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance_clear(jg,jlev))
366	total_albedo_below(jg,1) = reflectance_clear(jg,jlev) &
367	& + transmittance_clear(jg,jlev) * transmittance_clear(jg,jlev) &
368	& * total_albedo(jg,1,jlev+1) * inv_denom(jg,1)
369	total_albedo_below_direct(jg,1) = ref_dir_clear(jg,jlev) &
370	& + (trans_dir_dir_clear(jg,jlev)*total_albedo_direct(jg,1,jlev+1) &
371	& +trans_dir_diff_clear(jg,jlev)*total_albedo(jg,1,jlev+1)) &
372	& * transmittance_clear(jg,jlev) * inv_denom(jg,1)
373	end do
374
375	! Then the cloudy regions if any cloud is present in this layer
376	if (.not. is_clear_sky_layer(jlev)) then
377	do jreg = 2,nregions
378	do jg = 1,ng
379	inv_denom(jg,jreg) = 1.0_jprb / (1.0_jprb &
380	& - total_albedo(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev))
381	total_albedo_below(jg,jreg) = reflectance(jg,jreg,jlev) &
382	& + transmittance(jg,jreg,jlev) * transmittance(jg,jreg,jlev) &
383	& * total_albedo(jg,jreg,jlev+1) * inv_denom(jg,jreg)
384	total_albedo_below_direct(jg,jreg) = ref_dir(jg,jreg,jlev) &
385	& + (trans_dir_dir(jg,jreg,jlev)*total_albedo_direct(jg,jreg,jlev+1) &
386	& +trans_dir_diff(jg,jreg,jlev)*total_albedo(jg,jreg,jlev+1)) &
387	& * transmittance(jg,jreg,jlev) * inv_denom(jg,jreg)
388	end do
389	end do
390	end if
391
392	! Account for cloud overlap when converting albedo below a
393	! layer interface to the equivalent values just above
394	if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then
395	total_albedo(:,:,jlev) = total_albedo_below(:,:)
396	total_albedo_direct(:,:,jlev) = total_albedo_below_direct(:,:)
397	else
398	! Use overlap matrix and exclude "anomalous" horizontal
399	! transport described by Shonk & Hogan (2008). Therefore,
400	! the operation we perform is essentially diag(total_albedo)
401	! = matmul(transpose(v_matrix)), diag(total_albedo_below)).
402	do jreg = 1,nregions
403	do jreg2 = 1,nregions
404	total_albedo(:,jreg,jlev) &
405	& = total_albedo(:,jreg,jlev) &
406	& + total_albedo_below(:,jreg2) &
407	& * v_matrix(jreg2,jreg,jlev,jcol)
408	total_albedo_direct(:,jreg,jlev) &
409	& = total_albedo_direct(:,jreg,jlev) &
410	& + total_albedo_below_direct(:,jreg2) &
411	& * v_matrix(jreg2,jreg,jlev,jcol)
412	end do
413	end do
414
415	end if
416
417	end do ! Reverse loop over levels
418
419	! --------------------------------------------------------
420	! Section 5: Compute fluxes
421	! --------------------------------------------------------
422
423	! Top-of-atmosphere fluxes into the regions of the top-most
424	! layer, zero since we assume no diffuse downwelling
425	flux_dn = 0.0_jprb
426	! Direct downwelling flux (into a plane perpendicular to the
427	! sun) entering the top of each region in the top-most layer
428	do jreg = 1,nregions
429	direct_dn(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
430	end do
431	flux_up = direct_dn*total_albedo_direct(:,:,1)
432
433	if (config%do_clear) then
434	flux_dn_clear = 0.0_jprb
435	direct_dn_clear(:) = incoming_sw(:,jcol)
436	flux_up_clear = direct_dn_clear*total_albedo_clear_direct(:,1)
437	end if
438
439	! Store TOA spectral fluxes
440	flux%sw_up_toa_g(:,jcol) = sum(flux_up,2)
441	flux%sw_dn_toa_g(:,jcol) = incoming_sw(:,jcol)*mu0
442	if (config%do_clear) then
443	flux%sw_up_toa_clear_g(:,jcol) = flux_up_clear
444	end if
445
446	! Store the TOA broadband fluxes
447	flux%sw_up(jcol,1) = sum(sum(flux_up,1))
448	flux%sw_dn(jcol,1) = mu0 * sum(sum(direct_dn,1))
449	if (allocated(flux%sw_dn_direct)) then
450	flux%sw_dn_direct(jcol,1) = flux%sw_dn(jcol,1)
451	end if
452	if (config%do_clear) then
453	flux%sw_up_clear(jcol,1) = sum(flux_up_clear)
454	flux%sw_dn_clear(jcol,1) = mu0 * sum(direct_dn_clear)
455	if (allocated(flux%sw_dn_direct_clear)) then
456	flux%sw_dn_direct_clear(jcol,1) = flux%sw_dn_clear(jcol,1)
457	end if
458	end if
459
460	! Save the spectral fluxes if required; some redundancy here as
461	! the TOA downwelling flux is the same in clear and cloudy skies
462	if (config%do_save_spectral_flux) then
463	call indexed_sum(sum(flux_up,2), &
464	& config%i_spec_from_reordered_g_sw, &
465	& flux%sw_up_band(:,jcol,1))
466	call indexed_sum(sum(direct_dn,2), &
467	& config%i_spec_from_reordered_g_sw, &
468	& flux%sw_dn_band(:,jcol,1))
469	flux%sw_dn_band(:,jcol,1) = &
470	& mu0 * flux%sw_dn_band(:,jcol,1)
471	if (allocated(flux%sw_dn_direct_band)) then
472	flux%sw_dn_direct_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1)
473	end if
474	call add_indexed_sum(sum(flux_dn,2), &
475	& config%i_spec_from_reordered_g_sw, &
476	& flux%sw_dn_band(:,jcol,1))
477	if (config%do_clear) then
478	call indexed_sum(flux_up_clear, &
479	& config%i_spec_from_reordered_g_sw, &
480	& flux%sw_up_clear_band(:,jcol,1))
481	call indexed_sum(direct_dn_clear, &
482	& config%i_spec_from_reordered_g_sw, &
483	& flux%sw_dn_clear_band(:,jcol,1))
484	flux%sw_dn_clear_band(:,jcol,1) &
485	& = mu0 * flux%sw_dn_clear_band(:,jcol,1)
486	if (allocated(flux%sw_dn_direct_clear_band)) then
487	flux%sw_dn_direct_clear_band(:,jcol,1) &
488	& = flux%sw_dn_clear_band(:,jcol,1)
489	end if
490	call add_indexed_sum(flux_dn_clear, &
491	& config%i_spec_from_reordered_g_sw, &
492	& flux%sw_dn_clear_band(:,jcol,1))
493	end if
494	end if
495
496	! Final loop back down through the atmosphere to compute fluxes
497	do jlev = 1,nlev
498	if (config%do_clear) then
499	do jg = 1,ng
500	flux_dn_clear(jg) = (transmittance_clear(jg,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) &
501	& * (trans_dir_dir_clear(jg,jlev)total_albedo_clear_direct(jg,jlev+1)reflectance_clear(jg,jlev) &
502	& + trans_dir_diff_clear(jg,jlev) )) &
503	& / (1.0_jprb - reflectance_clear(jg,jlev)*total_albedo_clear(jg,jlev+1))
504	direct_dn_clear(jg) = trans_dir_dir_clear(jg,jlev)*direct_dn_clear(jg)
505	flux_up_clear(jg) = direct_dn_clear(jg)*total_albedo_clear_direct(jg,jlev+1) &
506	& + flux_dn_clear(jg)*total_albedo_clear(jg,jlev+1)
507	end do
508	end if
509
510	! All-sky fluxes: first the clear region...
511	do jg = 1,ng
512	flux_dn(jg,1) = (transmittance_clear(jg,jlev)*flux_dn(jg,1) + direct_dn(jg,1) &
513	& * (trans_dir_dir_clear(jg,jlev)total_albedo_direct(jg,1,jlev+1)reflectance_clear(jg,jlev) &
514	& + trans_dir_diff_clear(jg,jlev) )) &
515	& / (1.0_jprb - reflectance_clear(jg,jlev)*total_albedo(jg,1,jlev+1))
516	direct_dn(jg,1) = trans_dir_dir_clear(jg,jlev)*direct_dn(jg,1)
517	flux_up(jg,1) = direct_dn(jg,1)*total_albedo_direct(jg,1,jlev+1) &
518	& + flux_dn(jg,1)*total_albedo(jg,1,jlev+1)
519	end do
520
521	! ...then the cloudy regions if there are any
522	if (is_clear_sky_layer(jlev)) then
523	flux_dn(:,2:nregions) = 0.0_jprb
524	flux_up(:,2:nregions) = 0.0_jprb
525	direct_dn(:,2:nregions)= 0.0_jprb
526	else
527	do jreg = 2,nregions
528	do jg = 1,ng
529	flux_dn(jg,jreg) = (transmittance(jg,jreg,jlev)*flux_dn(jg,jreg) + direct_dn(jg,jreg) &
530	& * (trans_dir_dir(jg,jreg,jlev)total_albedo_direct(jg,jreg,jlev+1)reflectance(jg,jreg,jlev) &
531	& + trans_dir_diff(jg,jreg,jlev) )) &
532	& / (1.0_jprb - reflectance(jg,jreg,jlev)*total_albedo(jg,jreg,jlev+1))
533	direct_dn(jg,jreg) = trans_dir_dir(jg,jreg,jlev)*direct_dn(jg,jreg)
534	flux_up(jg,jreg) = direct_dn(jg,jreg)*total_albedo_direct(jg,jreg,jlev+1) &
535	& + flux_dn(jg,jreg)*total_albedo(jg,jreg,jlev+1)
536	end do
537	end do
538	end if
539
540	if (.not. (is_clear_sky_layer(jlev) &
541	& .and. is_clear_sky_layer(jlev+1))) then
542	! Account for overlap rules in translating fluxes just above
543	! a layer interface to the values just below
544	flux_dn = singlemat_x_vec(ng,ng,nregions, &
545	& v_matrix(:,:,jlev+1,jcol), flux_dn)
546	direct_dn = singlemat_x_vec(ng,ng,nregions, &
547	& v_matrix(:,:,jlev+1,jcol), direct_dn)
548	end if ! Otherwise the fluxes in each region are the same so
549	! nothing to do
550
551	! Store the broadband fluxes
552	flux%sw_up(jcol,jlev+1) = sum(sum(flux_up,1))
553	if (allocated(flux%sw_dn_direct)) then
554	flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1))
555	flux%sw_dn(jcol,jlev+1) &
556	& = flux%sw_dn_direct(jcol,jlev+1) + sum(sum(flux_dn,1))
557	else
558	flux%sw_dn(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1)) + sum(sum(flux_dn,1))
559	end if
560	if (config%do_clear) then
561	flux%sw_up_clear(jcol,jlev+1) = sum(flux_up_clear)
562	if (allocated(flux%sw_dn_direct_clear)) then
563	flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear)
564	flux%sw_dn_clear(jcol,jlev+1) &
565	& = flux%sw_dn_direct_clear(jcol,jlev+1) + sum(flux_dn_clear)
566	else
567	flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear) &
568	& + sum(flux_dn_clear)
569	end if
570	end if
571
572	! Save the spectral fluxes if required
573	if (config%do_save_spectral_flux) then
574	call indexed_sum(sum(flux_up,2), &
575	& config%i_spec_from_reordered_g_sw, &
576	& flux%sw_up_band(:,jcol,jlev+1))
577	call indexed_sum(sum(direct_dn,2), &
578	& config%i_spec_from_reordered_g_sw, &
579	& flux%sw_dn_band(:,jcol,jlev+1))
580	flux%sw_dn_band(:,jcol,jlev+1) = &
581	& mu0 * flux%sw_dn_band(:,jcol,jlev+1)
582	if (allocated(flux%sw_dn_direct_band)) then
583	flux%sw_dn_direct_band(:,jcol,jlev+1) &
584	& = flux%sw_dn_band(:,jcol,jlev+1)
585	end if
586	call add_indexed_sum(sum(flux_dn,2), &
587	& config%i_spec_from_reordered_g_sw, &
588	& flux%sw_dn_band(:,jcol,jlev+1))
589	if (config%do_clear) then
590	call indexed_sum(flux_up_clear, &
591	& config%i_spec_from_reordered_g_sw, &
592	& flux%sw_up_clear_band(:,jcol,jlev+1))
593	call indexed_sum(direct_dn_clear, &
594	& config%i_spec_from_reordered_g_sw, &
595	& flux%sw_dn_clear_band(:,jcol,jlev+1))
596	flux%sw_dn_clear_band(:,jcol,jlev+1) = &
597	& mu0 * flux%sw_dn_clear_band(:,jcol,jlev+1)
598	if (allocated(flux%sw_dn_direct_clear_band)) then
599	flux%sw_dn_direct_clear_band(:,jcol,jlev+1) &
600	& = flux%sw_dn_clear_band(:,jcol,jlev+1)
601	end if
602	call add_indexed_sum(flux_dn_clear, &
603	& config%i_spec_from_reordered_g_sw, &
604	& flux%sw_dn_clear_band(:,jcol,jlev+1))
605	end if
606	end if
607
608	end do ! Final loop over levels
609
610	! Store surface spectral fluxes, if required (after the end of
611	! the final loop over levels, the current values of these arrays
612	! will be the surface values)
613	flux%sw_dn_diffuse_surf_g(:,jcol) = sum(flux_dn,2)
614	flux%sw_dn_direct_surf_g(:,jcol) = mu0 * sum(direct_dn,2)
615	if (config%do_clear) then
616	flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_clear
617	flux%sw_dn_direct_surf_clear_g(:,jcol) = mu0 * direct_dn_clear
618	end if
619
620	end do ! Loop over columns
621
622	if (lhook) call dr_hook('radiation_tripleclouds_sw:solver_tripleclouds_sw',1,hook_handle)
623
624	end subroutine solver_tripleclouds_sw
625
626	end module radiation_tripleclouds_sw

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format