Context Navigation

radiation_tripleclouds_sw.F90 @ 4773

Last change on this file since 4773 was 4773, checked in by idelkadi, 5 months ago

Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)

Implementation of the double call of Ecrad in LMDZ

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