Context Navigation

radiation_tripleclouds_sw.F90 @ 4936

Last change on this file since 4936 was 4853, checked in by idelkadi, 8 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
interface routines between lmdz and ecrad grouped in a new "lmdz" directory
double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
clean-up in the read.F routine: delete unitules arguments
modification of default gas model in namelist (default: ECCKD)

File size: 28.6 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(config%n_g_sw,nlev,istartcol:iendcol) &
77	& :: od, ssa, g
78
79	! Cloud and precipitation optical depth, single-scattering albedo and
80	! asymmetry factor in each shortwave band
81	real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) &
82	& :: od_cloud, ssa_cloud, g_cloud
83
84	! Optical depth, single scattering albedo and asymmetry factor in
85	! each g-point of each cloudy region including gas, aerosol and
86	! clouds
87	real(jprb), dimension(config%n_g_sw,2:nregions) &
88	& :: od_total, ssa_total, g_total
89
90	! Direct and diffuse surface albedos, and the incoming shortwave
91	! flux into a plane perpendicular to the incoming radiation at
92	! top-of-atmosphere in each of the shortwave g points
93	real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) &
94	& :: albedo_direct, albedo_diffuse, incoming_sw
95
96	! Output
97	type(flux_type), intent(inout):: flux
98
99	! Local variables
100
101	! The area fractions of each region
102	real(jprb) :: region_fracs(1:nregions,nlev,istartcol:iendcol)
103
104	! The scaling used for the optical depth in the cloudy regions
105	real(jprb) :: od_scaling(2:nregions,nlev,istartcol:iendcol)
106
107	! Directional overlap matrices defined at all layer interfaces
108	! including top-of-atmosphere and the surface
109	real(jprb), dimension(nregions,nregions,nlev+1, &
110	& istartcol:iendcol) :: u_matrix, v_matrix
111
112	! Diffuse reflection and transmission matrices in the cloudy
113	! regions of each layer
114	real(jprb), dimension(config%n_g_sw, 2:nregions, nlev) &
115	& :: reflectance, transmittance
116
117	! Terms translating the direct flux entering the layer from above
118	! to the reflected radiation exiting upwards (ref_dir) and the
119	! scattered radiation exiting downwards (trans_dir_diff), along with the
120	! direct unscattered transmission matrix (trans_dir_dir).
121	real(jprb), dimension(config%n_g_sw, 2:nregions, nlev) &
122	& :: ref_dir, trans_dir_diff, trans_dir_dir
123
124	! As above but for the clear regions; clear and cloudy layers are
125	! separated out so that calc_ref_trans_sw can be called on the
126	! entire clear-sky atmosphere at once
127	real(jprb), dimension(config%n_g_sw, nlev) &
128	& :: reflectance_clear, transmittance_clear, &
129	& ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear
130
131	! Total albedo of the atmosphere/surface just above a layer
132	! interface with respect to downwelling diffuse and direct
133	! (respectively) radiation at that interface, where level index =
134	! 1 corresponds to the top-of-atmosphere
135	real(jprb), dimension(config%n_g_sw, nregions, nlev+1) &
136	& :: total_albedo, total_albedo_direct
137
138	! ...equivalent values for clear-skies
139	real(jprb), dimension(config%n_g_sw, nlev+1) &
140	& :: total_albedo_clear, total_albedo_clear_direct
141
142	! Total albedo of the atmosphere just below a layer interface
143	real(jprb), dimension(config%n_g_sw, nregions) &
144	& :: total_albedo_below, total_albedo_below_direct
145
146	! Direct downwelling flux below and above an interface between
147	! layers into a plane perpendicular to the direction of the sun
148	real(jprb), dimension(config%n_g_sw, nregions) :: direct_dn
149	! Diffuse equivalents
150	real(jprb), dimension(config%n_g_sw, nregions) :: flux_dn, flux_up
151
152	! ...clear-sky equivalent (no distinction between "above/below")
153	real(jprb), dimension(config%n_g_sw) &
154	& :: direct_dn_clear, flux_dn_clear, flux_up_clear
155
156	! Clear-sky equivalent, but actually its reciprocal to replace
157	! some divisions by multiplications
158	real(jprb), dimension(config%n_g_sw, nregions) :: inv_denom
159
160	! Identify clear-sky layers, with pseudo layers for outer space
161	! and below the ground, both treated as single-region clear skies
162	logical :: is_clear_sky_layer(0:nlev+1)
163
164	! Scattering optical depth of gas+aerosol and of cloud
165	real(jprb) :: scat_od, scat_od_cloud
166
167	! Temporaries to speed up summations
168	real(jprb) :: sum_dn_diff, sum_dn_dir, sum_up
169
170	! Local cosine of solar zenith angle
171	real(jprb) :: mu0
172
173	integer :: jcol, jlev, jg, jreg, iband, jreg2, ng
174
175	real(jphook) :: hook_handle
176
177	if (lhook) call dr_hook('radiation_tripleclouds_sw:solver_tripleclouds_sw',0,hook_handle)
178
179	! --------------------------------------------------------
180	! Section 1: Prepare general variables and arrays
181	! --------------------------------------------------------
182	! Copy array dimensions to local variables for convenience
183	ng = config%n_g_sw
184
185	! Compute the wavelength-independent region fractions and
186	! optical-depth scalings
187	call calc_region_properties(nlev,nregions,istartcol,iendcol, &
188	& config%i_cloud_pdf_shape == IPdfShapeGamma, &
189	& cloud%fraction, cloud%fractional_std, region_fracs, &
190	& od_scaling, config%cloud_fraction_threshold)
191
192	! Compute wavelength-independent overlap matrices u_matrix and
193	! v_matrix
194	call calc_overlap_matrices(nlev,nregions,istartcol,iendcol, &
195	& region_fracs, cloud%overlap_param, &
196	& u_matrix, v_matrix, &
197	& decorrelation_scaling=config%cloud_inhom_decorr_scaling, &
198	& cloud_fraction_threshold=config%cloud_fraction_threshold, &
199	& use_beta_overlap=config%use_beta_overlap, &
200	& cloud_cover=flux%cloud_cover_sw)
201
202	! Main loop over columns
203	do jcol = istartcol, iendcol
204	! --------------------------------------------------------
205	! Section 2: Prepare column-specific variables and arrays
206	! --------------------------------------------------------
207
208	! Copy local cosine of the solar zenith angle
209	mu0 = single_level%cos_sza(jcol)
210
211	! Skip profile if sun is too low in the sky
212	if (mu0 < 1.0e-10_jprb) then
213	flux%sw_dn(jcol,:) = 0.0_jprb
214	flux%sw_up(jcol,:) = 0.0_jprb
215	if (allocated(flux%sw_dn_direct)) then
216	flux%sw_dn_direct(jcol,:) = 0.0_jprb
217	end if
218	if (config%do_clear) then
219	flux%sw_dn_clear(jcol,:) = 0.0_jprb
220	flux%sw_up_clear(jcol,:) = 0.0_jprb
221	if (allocated(flux%sw_dn_direct_clear)) then
222	flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
223	end if
224	end if
225
226	if (config%do_save_spectral_flux) then
227	flux%sw_dn_band(:,jcol,:) = 0.0_jprb
228	flux%sw_up_band(:,jcol,:) = 0.0_jprb
229	if (allocated(flux%sw_dn_direct_band)) then
230	flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb
231	end if
232	if (config%do_clear) then
233	flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb
234	flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb
235	if (allocated(flux%sw_dn_direct_clear_band)) then
236	flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb
237	end if
238	end if
239	end if
240
241	flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
242	flux%sw_dn_direct_surf_g(:,jcol) = 0.0_jprb
243	if (config%do_clear) then
244	flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
245	flux%sw_dn_direct_surf_clear_g(:,jcol) = 0.0_jprb
246	end if
247
248	cycle
249	end if ! sun is below the horizon
250
251	! At this point mu0 >= 1.0e-10
252
253	! Define which layers contain cloud; assume that
254	! cloud%crop_cloud_fraction has already been called
255	is_clear_sky_layer = .true.
256	do jlev = 1,nlev
257	if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
258	is_clear_sky_layer(jlev) = .false.
259	end if
260	end do
261
262	! --------------------------------------------------------
263	! Section 3: Loop over layers to compute reflectance and transmittance
264	! --------------------------------------------------------
265	! In this section the reflectance, transmittance and sources
266	! are computed for each layer
267
268	! Clear-sky quantities for all layers at once
269	call calc_ref_trans_sw(ng*nlev, &
270	& mu0, od(:,:,jcol), ssa(:,:,jcol), g(:,:,jcol), &
271	& reflectance_clear, transmittance_clear, &
272	& ref_dir_clear, trans_dir_diff_clear, &
273	& trans_dir_dir_clear)
274
275	! Cloudy layers
276	do jlev = 1,nlev ! Start at top-of-atmosphere
277	if (.not. is_clear_sky_layer(jlev)) then
278	do jreg = 2,nregions
279	do jg = 1,ng
280	! Mapping from g-point to band
281	iband = config%i_band_from_reordered_g_sw(jg)
282	scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol)
283	scat_od_cloud = od_cloud(iband,jlev,jcol) &
284	& * ssa_cloud(iband,jlev,jcol) * od_scaling(jreg,jlev,jcol)
285	! Add scaled cloud optical depth to clear-sky value
286	od_total(jg,jreg) = od(jg,jlev,jcol) &
287	& + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
288	! Compute single-scattering albedo and asymmetry
289	! factor of gas-cloud combination
290	ssa_total(jg,jreg) = (scat_od+scat_od_cloud) &
291	& / od_total(jg,jreg)
292	g_total(jg,jreg) = (scat_od*g(jg,jlev,jcol) &
293	& + scat_od_cloud * g_cloud(iband,jlev,jcol)) &
294	& / (scat_od + scat_od_cloud)
295	end do
296	end do
297
298	if (config%do_sw_delta_scaling_with_gases) then
299	! Apply delta-Eddington scaling to the aerosol-gas(-cloud)
300	! mixture
301	call delta_eddington(od_total, ssa_total, g_total)
302	end if
303	! Both cloudy regions at once
304	call calc_ref_trans_sw(ng*(nregions-1), &
305	& mu0, od_total, ssa_total, g_total, &
306	& reflectance(:,:,jlev), transmittance(:,:,jlev), &
307	& ref_dir(:,:,jlev), trans_dir_diff(:,:,jlev), &
308	& trans_dir_dir(:,:,jlev))
309	end if
310	end do
311
312	! --------------------------------------------------------
313	! Section 4: Compute total albedos
314	! --------------------------------------------------------
315
316	total_albedo(:,:,:) = 0.0_jprb
317	total_albedo_direct(:,:,:) = 0.0_jprb
318
319	! Copy surface albedos in clear-sky region
320	do jg = 1,ng
321	total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol)
322	total_albedo_direct(jg,1,nlev+1) &
323	& = mu0 * albedo_direct(jg,jcol)
324	end do
325
326	! If there is cloud in the lowest layer then we need the albedos
327	! underneath
328	if (.not. is_clear_sky_layer(nlev)) then
329	do jreg = 2,nregions
330	total_albedo(:,jreg,nlev+1) = total_albedo(:,1,nlev+1)
331	total_albedo_direct(:,jreg,nlev+1) = total_albedo_direct(:,1,nlev+1)
332	end do
333	end if
334
335	if (config%do_clear) then
336	total_albedo_clear(:,nlev+1) = total_albedo(:,1,nlev+1)
337	total_albedo_clear_direct(:,nlev+1) = total_albedo_direct(:,1,nlev+1)
338	end if
339
340	! Work up from the surface computing the total albedo of the
341	! atmosphere below that point using the adding method
342	do jlev = nlev,1,-1
343
344	total_albedo_below = 0.0_jprb
345	total_albedo_below_direct = 0.0_jprb
346
347	if (config%do_clear) then
348	! For clear-skies there is no need to consider "above" and
349	! "below" quantities since with no cloud overlap to worry
350	! about, these are the same
351	do jg = 1,ng
352	inv_denom(jg,1) = 1.0_jprb &
353	& / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance_clear(jg,jlev))
354	total_albedo_clear(jg,jlev) = reflectance_clear(jg,jlev) &
355	& + transmittance_clear(jg,jlev) * transmittance_clear(jg,jlev) &
356	& * total_albedo_clear(jg,jlev+1) * inv_denom(jg,1)
357
358	total_albedo_clear_direct(jg,jlev) = ref_dir_clear(jg,jlev) &
359	& + (trans_dir_dir_clear(jg,jlev)*total_albedo_clear_direct(jg,jlev+1) &
360	& +trans_dir_diff_clear(jg,jlev)*total_albedo_clear(jg,jlev+1)) &
361	& * transmittance_clear(jg,jlev) * inv_denom(jg,1)
362	end do
363	end if
364
365	! All-sky fluxes: first the clear region
366	do jg = 1,ng
367	inv_denom(jg,1) = 1.0_jprb &
368	& / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance_clear(jg,jlev))
369	total_albedo_below(jg,1) = reflectance_clear(jg,jlev) &
370	& + transmittance_clear(jg,jlev) * transmittance_clear(jg,jlev) &
371	& * total_albedo(jg,1,jlev+1) * inv_denom(jg,1)
372	total_albedo_below_direct(jg,1) = ref_dir_clear(jg,jlev) &
373	& + (trans_dir_dir_clear(jg,jlev)*total_albedo_direct(jg,1,jlev+1) &
374	& +trans_dir_diff_clear(jg,jlev)*total_albedo(jg,1,jlev+1)) &
375	& * transmittance_clear(jg,jlev) * inv_denom(jg,1)
376	end do
377
378	! Then the cloudy regions if any cloud is present in this layer
379	if (.not. is_clear_sky_layer(jlev)) then
380	do jreg = 2,nregions
381	do jg = 1,ng
382	inv_denom(jg,jreg) = 1.0_jprb / (1.0_jprb &
383	& - total_albedo(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev))
384	total_albedo_below(jg,jreg) = reflectance(jg,jreg,jlev) &
385	& + transmittance(jg,jreg,jlev) * transmittance(jg,jreg,jlev) &
386	& * total_albedo(jg,jreg,jlev+1) * inv_denom(jg,jreg)
387	total_albedo_below_direct(jg,jreg) = ref_dir(jg,jreg,jlev) &
388	& + (trans_dir_dir(jg,jreg,jlev)*total_albedo_direct(jg,jreg,jlev+1) &
389	& +trans_dir_diff(jg,jreg,jlev)*total_albedo(jg,jreg,jlev+1)) &
390	& * transmittance(jg,jreg,jlev) * inv_denom(jg,jreg)
391	end do
392	end do
393	end if
394
395	! Account for cloud overlap when converting albedo below a
396	! layer interface to the equivalent values just above
397	if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then
398	total_albedo(:,:,jlev) = total_albedo_below(:,:)
399	total_albedo_direct(:,:,jlev) = total_albedo_below_direct(:,:)
400	else
401	! Use overlap matrix and exclude "anomalous" horizontal
402	! transport described by Shonk & Hogan (2008). Therefore,
403	! the operation we perform is essentially diag(total_albedo)
404	! = matmul(transpose(v_matrix)), diag(total_albedo_below)).
405	do jreg = 1,nregions
406	do jreg2 = 1,nregions
407	total_albedo(:,jreg,jlev) &
408	& = total_albedo(:,jreg,jlev) &
409	& + total_albedo_below(:,jreg2) &
410	& * v_matrix(jreg2,jreg,jlev,jcol)
411	total_albedo_direct(:,jreg,jlev) &
412	& = total_albedo_direct(:,jreg,jlev) &
413	& + total_albedo_below_direct(:,jreg2) &
414	& * v_matrix(jreg2,jreg,jlev,jcol)
415	end do
416	end do
417
418	end if
419
420	end do ! Reverse loop over levels
421
422	! --------------------------------------------------------
423	! Section 5: Compute fluxes
424	! --------------------------------------------------------
425
426	! Top-of-atmosphere fluxes into the regions of the top-most
427	! layer, zero since we assume no diffuse downwelling
428	flux_dn = 0.0_jprb
429	! Direct downwelling flux (into a plane perpendicular to the
430	! sun) entering the top of each region in the top-most layer
431	do jreg = 1,nregions
432	direct_dn(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
433	end do
434	flux_up = direct_dn*total_albedo_direct(:,:,1)
435
436	if (config%do_clear) then
437	flux_dn_clear = 0.0_jprb
438	direct_dn_clear(:) = incoming_sw(:,jcol)
439	flux_up_clear = direct_dn_clear*total_albedo_clear_direct(:,1)
440	end if
441
442	! Store TOA spectral fluxes
443	flux%sw_up_toa_g(:,jcol) = sum(flux_up,2)
444	flux%sw_dn_toa_g(:,jcol) = incoming_sw(:,jcol)*mu0
445	if (config%do_clear) then
446	flux%sw_up_toa_clear_g(:,jcol) = flux_up_clear
447	end if
448
449	! Store the TOA broadband fluxes, noting that there is no
450	! diffuse downwelling at TOA. The intrinsic "sum" command has
451	! been found to be very slow; better performance is found on
452	! x86-64 architecture with explicit loops and the "omp simd
453	! reduction" directive.
454	sum_up = 0.0_jprb
455	sum_dn_dir = 0.0_jprb
456	do jreg = 1,nregions
457	!$omp simd reduction(+:sum_up, sum_dn_dir)
458	do jg = 1,ng
459	sum_up = sum_up + flux_up(jg,jreg)
460	sum_dn_dir = sum_dn_dir + direct_dn(jg,jreg)
461	end do
462	end do
463	flux%sw_up(jcol,1) = sum_up
464	flux%sw_dn(jcol,1) = mu0 * sum_dn_dir
465	if (allocated(flux%sw_dn_direct)) then
466	flux%sw_dn_direct(jcol,1) = flux%sw_dn(jcol,1)
467	end if
468	if (config%do_clear) then
469	sum_up = 0.0_jprb
470	sum_dn_dir = 0.0_jprb
471	!$omp simd reduction(+:sum_up, sum_dn_dir)
472	do jg = 1,ng
473	sum_up = sum_up + flux_up_clear(jg)
474	sum_dn_dir = sum_dn_dir + direct_dn_clear(jg)
475	end do
476	flux%sw_up_clear(jcol,1) = sum_up
477	flux%sw_dn_clear(jcol,1) = mu0 * sum_dn_dir
478	if (allocated(flux%sw_dn_direct_clear)) then
479	flux%sw_dn_direct_clear(jcol,1) = flux%sw_dn_clear(jcol,1)
480	end if
481	end if
482
483	! Save the spectral fluxes if required; some redundancy here as
484	! the TOA downwelling flux is the same in clear and cloudy skies
485	if (config%do_save_spectral_flux) then
486	call indexed_sum(sum(flux_up,2), &
487	& config%i_spec_from_reordered_g_sw, &
488	& flux%sw_up_band(:,jcol,1))
489	call indexed_sum(sum(direct_dn,2), &
490	& config%i_spec_from_reordered_g_sw, &
491	& flux%sw_dn_band(:,jcol,1))
492	flux%sw_dn_band(:,jcol,1) = mu0 * flux%sw_dn_band(:,jcol,1)
493	if (allocated(flux%sw_dn_direct_band)) then
494	flux%sw_dn_direct_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1)
495	end if
496	call add_indexed_sum(sum(flux_dn,2), &
497	& config%i_spec_from_reordered_g_sw, &
498	& flux%sw_dn_band(:,jcol,1))
499	if (config%do_clear) then
500	call indexed_sum(flux_up_clear, &
501	& config%i_spec_from_reordered_g_sw, &
502	& flux%sw_up_clear_band(:,jcol,1))
503	call indexed_sum(direct_dn_clear, &
504	& config%i_spec_from_reordered_g_sw, &
505	& flux%sw_dn_clear_band(:,jcol,1))
506	flux%sw_dn_clear_band(:,jcol,1) &
507	& = mu0 * flux%sw_dn_clear_band(:,jcol,1)
508	if (allocated(flux%sw_dn_direct_clear_band)) then
509	flux%sw_dn_direct_clear_band(:,jcol,1) &
510	& = flux%sw_dn_clear_band(:,jcol,1)
511	end if
512	call add_indexed_sum(flux_dn_clear, &
513	& config%i_spec_from_reordered_g_sw, &
514	& flux%sw_dn_clear_band(:,jcol,1))
515	end if
516	end if
517
518	! Final loop back down through the atmosphere to compute fluxes
519	do jlev = 1,nlev
520	if (config%do_clear) then
521	do jg = 1,ng
522	flux_dn_clear(jg) = (transmittance_clear(jg,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) &
523	& * (trans_dir_dir_clear(jg,jlev)total_albedo_clear_direct(jg,jlev+1)reflectance_clear(jg,jlev) &
524	& + trans_dir_diff_clear(jg,jlev) )) &
525	& / (1.0_jprb - reflectance_clear(jg,jlev)*total_albedo_clear(jg,jlev+1))
526	direct_dn_clear(jg) = trans_dir_dir_clear(jg,jlev)*direct_dn_clear(jg)
527	flux_up_clear(jg) = direct_dn_clear(jg)*total_albedo_clear_direct(jg,jlev+1) &
528	& + flux_dn_clear(jg)*total_albedo_clear(jg,jlev+1)
529	end do
530	end if
531
532	! All-sky fluxes: first the clear region...
533	do jg = 1,ng
534	flux_dn(jg,1) = (transmittance_clear(jg,jlev)*flux_dn(jg,1) + direct_dn(jg,1) &
535	& * (trans_dir_dir_clear(jg,jlev)total_albedo_direct(jg,1,jlev+1)reflectance_clear(jg,jlev) &
536	& + trans_dir_diff_clear(jg,jlev) )) &
537	& / (1.0_jprb - reflectance_clear(jg,jlev)*total_albedo(jg,1,jlev+1))
538	direct_dn(jg,1) = trans_dir_dir_clear(jg,jlev)*direct_dn(jg,1)
539	flux_up(jg,1) = direct_dn(jg,1)*total_albedo_direct(jg,1,jlev+1) &
540	& + flux_dn(jg,1)*total_albedo(jg,1,jlev+1)
541	end do
542
543	! ...then the cloudy regions if there are any
544	if (is_clear_sky_layer(jlev)) then
545	flux_dn(:,2:nregions) = 0.0_jprb
546	flux_up(:,2:nregions) = 0.0_jprb
547	direct_dn(:,2:nregions)= 0.0_jprb
548	else
549	do jreg = 2,nregions
550	do jg = 1,ng
551	flux_dn(jg,jreg) = (transmittance(jg,jreg,jlev)*flux_dn(jg,jreg) + direct_dn(jg,jreg) &
552	& * (trans_dir_dir(jg,jreg,jlev)total_albedo_direct(jg,jreg,jlev+1)reflectance(jg,jreg,jlev) &
553	& + trans_dir_diff(jg,jreg,jlev) )) &
554	& / (1.0_jprb - reflectance(jg,jreg,jlev)*total_albedo(jg,jreg,jlev+1))
555	direct_dn(jg,jreg) = trans_dir_dir(jg,jreg,jlev)*direct_dn(jg,jreg)
556	flux_up(jg,jreg) = direct_dn(jg,jreg)*total_albedo_direct(jg,jreg,jlev+1) &
557	& + flux_dn(jg,jreg)*total_albedo(jg,jreg,jlev+1)
558	end do
559	end do
560	end if
561
562	if (.not. (is_clear_sky_layer(jlev) &
563	& .and. is_clear_sky_layer(jlev+1))) then
564	! Account for overlap rules in translating fluxes just above
565	! a layer interface to the values just below
566	flux_dn = singlemat_x_vec(ng,ng,nregions, &
567	& v_matrix(:,:,jlev+1,jcol), flux_dn)
568	direct_dn = singlemat_x_vec(ng,ng,nregions, &
569	& v_matrix(:,:,jlev+1,jcol), direct_dn)
570	end if ! Otherwise the fluxes in each region are the same so
571	! nothing to do
572
573	! Store the broadband fluxes. The intrinsic "sum" command has
574	! been found to be very slow; better performance is found on
575	! x86-64 architecture with explicit loops and the "omp simd
576	! reduction" directive.
577	sum_up = 0.0_jprb
578	sum_dn_dir = 0.0_jprb
579	sum_dn_diff = 0.0_jprb
580	do jreg = 1,nregions
581	!$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
582	do jg = 1,ng
583	sum_up = sum_up + flux_up(jg,jreg)
584	sum_dn_diff = sum_dn_diff + flux_dn(jg,jreg)
585	sum_dn_dir = sum_dn_dir + direct_dn(jg,jreg)
586	end do
587	end do
588	flux%sw_up(jcol,jlev+1) = sum_up
589	flux%sw_dn(jcol,jlev+1) = mu0 * sum_dn_dir + sum_dn_diff
590	if (allocated(flux%sw_dn_direct)) then
591	flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum_dn_dir
592	end if
593	if (config%do_clear) then
594	sum_up = 0.0_jprb
595	sum_dn_dir = 0.0_jprb
596	sum_dn_diff = 0.0_jprb
597	!$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
598	do jg = 1,ng
599	sum_up = sum_up + flux_up_clear(jg)
600	sum_dn_diff = sum_dn_diff + flux_dn_clear(jg)
601	sum_dn_dir = sum_dn_dir + direct_dn_clear(jg)
602	end do
603	flux%sw_up_clear(jcol,jlev+1) = sum_up
604	flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum_dn_dir + sum_dn_diff
605	if (allocated(flux%sw_dn_direct_clear)) then
606	flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum_dn_dir
607	end if
608	end if
609
610	! Save the spectral fluxes if required
611	if (config%do_save_spectral_flux) then
612	call indexed_sum(sum(flux_up,2), &
613	& config%i_spec_from_reordered_g_sw, &
614	& flux%sw_up_band(:,jcol,jlev+1))
615	call indexed_sum(sum(direct_dn,2), &
616	& config%i_spec_from_reordered_g_sw, &
617	& flux%sw_dn_band(:,jcol,jlev+1))
618	flux%sw_dn_band(:,jcol,jlev+1) = &
619	& mu0 * flux%sw_dn_band(:,jcol,jlev+1)
620	if (allocated(flux%sw_dn_direct_band)) then
621	flux%sw_dn_direct_band(:,jcol,jlev+1) &
622	& = flux%sw_dn_band(:,jcol,jlev+1)
623	end if
624	call add_indexed_sum(sum(flux_dn,2), &
625	& config%i_spec_from_reordered_g_sw, &
626	& flux%sw_dn_band(:,jcol,jlev+1))
627	if (config%do_clear) then
628	call indexed_sum(flux_up_clear, &
629	& config%i_spec_from_reordered_g_sw, &
630	& flux%sw_up_clear_band(:,jcol,jlev+1))
631	call indexed_sum(direct_dn_clear, &
632	& config%i_spec_from_reordered_g_sw, &
633	& flux%sw_dn_clear_band(:,jcol,jlev+1))
634	flux%sw_dn_clear_band(:,jcol,jlev+1) = &
635	& mu0 * flux%sw_dn_clear_band(:,jcol,jlev+1)
636	if (allocated(flux%sw_dn_direct_clear_band)) then
637	flux%sw_dn_direct_clear_band(:,jcol,jlev+1) &
638	& = flux%sw_dn_clear_band(:,jcol,jlev+1)
639	end if
640	call add_indexed_sum(flux_dn_clear, &
641	& config%i_spec_from_reordered_g_sw, &
642	& flux%sw_dn_clear_band(:,jcol,jlev+1))
643	end if
644	end if
645	end do ! Final loop over levels
646
647	! Store surface spectral fluxes, if required (after the end of
648	! the final loop over levels, the current values of these arrays
649	! will be the surface values)
650	flux%sw_dn_diffuse_surf_g(:,jcol) = sum(flux_dn,2)
651	flux%sw_dn_direct_surf_g(:,jcol) = mu0 * sum(direct_dn,2)
652	if (config%do_clear) then
653	flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_clear
654	flux%sw_dn_direct_surf_clear_g(:,jcol) = mu0 * direct_dn_clear
655	end if
656
657	end do ! Loop over columns
658
659	if (lhook) call dr_hook('radiation_tripleclouds_sw:solver_tripleclouds_sw',1,hook_handle)
660
661	end subroutine solver_tripleclouds_sw
662
663	end module radiation_tripleclouds_sw

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

Download in other formats:

Original Format