Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

radiation_homogeneous_sw.F90 @ 5097

Last change on this file since 5097 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

Inclusion of the ecrad directory containing the sources of the ECRAD code
- interface routine : radiation_scheme.F90
Adaptation of compilation scripts :
- compilation under CPP key CPP_ECRAD
- compilation with option "-rad ecard" or "-ecard true"
- The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
Runing LMDZ simulations with ecrad, you need :
- logical key iflag_rrtm = 2 in physiq.def
- namelist_ecrad (DefLists?)
- the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...

File size: 15.8 KB

Line
1	! radiation_homogeneous_sw.F90 - Shortwave homogeneous-column (no cloud fraction) 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	! 2019-01-14 R. Hogan Save spectral flux profile if required
20
21	module radiation_homogeneous_sw
22
23	public
24
25	contains
26
27	! Provides elemental function "delta_eddington"
28	#include "radiation_delta_eddington.h"
29
30	!---------------------------------------------------------------------
31	! Shortwave homogeneous solver, in which clouds are assumed to fill
32	! the gridbox horizontally
33	subroutine solver_homogeneous_sw(nlev,istartcol,iendcol, &
34	& config, single_level, cloud, &
35	& od, ssa, g, od_cloud, ssa_cloud, g_cloud, &
36	& albedo_direct, albedo_diffuse, incoming_sw, &
37	& flux)
38
39	use parkind1, only : jprb
40	use yomhook, only : lhook, dr_hook
41
42	use radiation_config, only : config_type
43	use radiation_single_level, only : single_level_type
44	use radiation_cloud, only : cloud_type
45	use radiation_flux, only : flux_type, &
46	& indexed_sum_profile, add_indexed_sum_profile
47	use radiation_two_stream, only : calc_two_stream_gammas_sw, &
48	& calc_reflectance_transmittance_sw
49	use radiation_constants, only : Pi, GasConstantDryAir, &
50	& AccelDueToGravity
51	use radiation_adding_ica_sw, only : adding_ica_sw
52
53	implicit none
54
55	! Inputs
56	integer, intent(in) :: nlev ! number of model levels
57	integer, intent(in) :: istartcol, iendcol ! range of columns to process
58	type(config_type), intent(in) :: config
59	type(single_level_type), intent(in) :: single_level
60	type(cloud_type), intent(in) :: cloud
61
62	! Gas and aerosol optical depth, single-scattering albedo and
63	! asymmetry factor at each shortwave g-point
64	real(jprb), intent(in), dimension(config%n_g_sw, nlev, istartcol:iendcol) :: &
65	& od, ssa, g
66
67	! Cloud and precipitation optical depth, single-scattering albedo and
68	! asymmetry factor in each shortwave band
69	real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
70	& od_cloud, ssa_cloud, g_cloud
71
72	! Direct and diffuse surface albedos, and the incoming shortwave
73	! flux into a plane perpendicular to the incoming radiation at
74	! top-of-atmosphere in each of the shortwave g points
75	real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
76	& albedo_direct, albedo_diffuse, incoming_sw
77
78	! Output
79	type(flux_type), intent(inout):: flux
80
81	! Local variables
82
83	! Cosine of solar zenith angle
84	real(jprb) :: cos_sza
85
86	! Diffuse reflectance and transmittance for each layer
87	real(jprb), dimension(config%n_g_sw, nlev) :: reflectance, transmittance
88
89	! Fraction of direct beam scattered by a layer into the upwelling
90	! or downwelling diffuse streams
91	real(jprb), dimension(config%n_g_sw, nlev) :: ref_dir, trans_dir_diff
92
93	! Transmittance for the direct beam in clear and all skies
94	real(jprb), dimension(config%n_g_sw, nlev) :: trans_dir_dir
95
96	! Fluxes per g point
97	real(jprb), dimension(config%n_g_sw, nlev+1) :: flux_up, flux_dn_diffuse, flux_dn_direct
98
99	! Combined gas+aerosol+cloud optical depth, single scattering
100	! albedo and asymmetry factor
101	real(jprb), dimension(config%n_g_sw) :: od_total, ssa_total, g_total
102
103	! Two-stream coefficients
104	real(jprb), dimension(config%n_g_sw) :: gamma1, gamma2, gamma3
105
106	! Optical depth of cloud in g-point space
107	real(jprb), dimension(config%n_g_sw) :: od_cloud_g
108
109	! Is there any cloud in the profile?
110	logical :: is_cloudy_profile
111
112	! Number of g points
113	integer :: ng
114
115	! Loop indices for level and column
116	integer :: jlev, jcol
117
118	real(jprb) :: hook_handle
119
120	if (lhook) call dr_hook('radiation_homogeneous_sw:solver_homogeneous_sw',0,hook_handle)
121
122	ng = config%n_g_sw
123
124	! Loop through columns
125	do jcol = istartcol,iendcol
126	! Only perform calculation if sun above the horizon
127	if (single_level%cos_sza(jcol) > 0.0_jprb) then
128
129	cos_sza = single_level%cos_sza(jcol)
130
131	! Is there any cloud in the profile?
132	is_cloudy_profile = .false.
133	do jlev = 1,nlev
134	if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
135	is_cloudy_profile = .true.
136	exit
137	end if
138	end do
139
140	! If clear-sky fluxes need to be computed then we first
141	! compute the reflectance and transmittance of all layers,
142	! neglecting clouds. If clear-sky fluxes are not required then
143	! we only do the clear-sky layers since these will be needed
144	! when we come to do the total-sky fluxes.
145	if (.not. config%do_sw_delta_scaling_with_gases) then
146	! Delta-Eddington scaling has already been performed to the
147	! aerosol part of od, ssa and g
148	do jlev = 1,nlev
149	if (config%do_clear .or. cloud%fraction(jcol,jlev) &
150	& < config%cloud_fraction_threshold) then
151	call calc_two_stream_gammas_sw(ng, cos_sza, &
152	& ssa(:,jlev,jcol), g(:,jlev,jcol), &
153	& gamma1, gamma2, gamma3)
154	call calc_reflectance_transmittance_sw(ng, &
155	& cos_sza, &
156	& od(:,jlev,jcol), ssa(:,jlev,jcol), &
157	& gamma1, gamma2, gamma3, &
158	& reflectance(:,jlev), transmittance(:,jlev), &
159	& ref_dir(:,jlev), trans_dir_diff(:,jlev), &
160	& trans_dir_dir(:,jlev) )
161
162	end if
163	end do
164	else
165	! Apply delta-Eddington scaling to the aerosol-gas mixture
166	do jlev = 1,nlev
167	if (config%do_clear .or. cloud%fraction(jcol,jlev) &
168	& < config%cloud_fraction_threshold) then
169	od_total = od(:,jlev,jcol)
170	ssa_total = ssa(:,jlev,jcol)
171	g_total = g(:,jlev,jcol)
172	call delta_eddington(od_total, ssa_total, g_total)
173	call calc_two_stream_gammas_sw(ng, &
174	& cos_sza, ssa_total, g_total, &
175	& gamma1, gamma2, gamma3)
176	call calc_reflectance_transmittance_sw(ng, &
177	& cos_sza, od_total, ssa_total, &
178	& gamma1, gamma2, gamma3, &
179	& reflectance(:,jlev), transmittance(:,jlev), &
180	& ref_dir(:,jlev), trans_dir_diff(:,jlev), &
181	& trans_dir_dir(:,jlev) )
182	end if
183	end do
184	end if
185
186	if (config%do_clear) then
187	! Use adding method to compute fluxes
188	call adding_ica_sw(ng, nlev, incoming_sw(:,jcol), &
189	& albedo_diffuse(:,jcol), albedo_direct(:,jcol), &
190	& spread(cos_sza,1,ng), reflectance, transmittance, ref_dir, trans_dir_diff, &
191	& trans_dir_dir, flux_up, flux_dn_diffuse, flux_dn_direct)
192
193	! Sum over g-points to compute and save clear-sky broadband
194	! fluxes
195	flux%sw_up_clear(jcol,:) = sum(flux_up,1)
196	if (allocated(flux%sw_dn_direct_clear)) then
197	flux%sw_dn_direct_clear(jcol,:) &
198	& = sum(flux_dn_direct,1)
199	flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
200	& + flux%sw_dn_direct_clear(jcol,:)
201	else
202	flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
203	& + sum(flux_dn_direct,1)
204	end if
205	! Store spectral downwelling fluxes at surface
206	flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
207	flux%sw_dn_direct_surf_clear_g(:,jcol) = flux_dn_direct(:,nlev+1)
208
209	! Save the spectral fluxes if required
210	if (config%do_save_spectral_flux) then
211	call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_sw, &
212	& flux%sw_up_clear_band(:,jcol,:))
213	call indexed_sum_profile(flux_dn_direct, config%i_spec_from_reordered_g_sw, &
214	& flux%sw_dn_clear_band(:,jcol,:))
215	if (allocated(flux%sw_dn_direct_clear_band)) then
216	flux%sw_dn_direct_clear_band(:,jcol,:) &
217	& = flux%sw_dn_clear_band(:,jcol,:)
218	end if
219	call add_indexed_sum_profile(flux_dn_diffuse, &
220	& config%i_spec_from_reordered_g_sw, &
221	& flux%sw_dn_clear_band(:,jcol,:))
222	end if
223
224	end if ! Do clear-sky calculations
225
226	! Now the total-sky calculation. If this is a clear profile
227	! and clear-sky fluxes have been calculated then we can simply
228	! copy over the clear-sky fluxes, otherwise we need to compute
229	! fluxes now.
230	if (is_cloudy_profile .or. .not. config%do_clear) then
231	do jlev = 1,nlev
232	! Compute combined gas+aerosol+cloud optical properties;
233	! note that for clear layers, the reflectance and
234	! transmittance have already been calculated
235	if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
236	od_cloud_g = od_cloud(config%i_band_from_reordered_g_sw,jlev,jcol)
237	od_total = od(:,jlev,jcol) + od_cloud_g
238	ssa_total = 0.0_jprb
239	g_total = 0.0_jprb
240	where (od_total > 0.0_jprb)
241	ssa_total = (ssa(:,jlev,jcol)*od(:,jlev,jcol) &
242	& + ssa_cloud(config%i_band_from_reordered_g_sw,jlev,jcol) &
243	& * od_cloud_g) &
244	& / od_total
245	end where
246	where (ssa_total > 0.0_jprb .and. od_total > 0.0_jprb)
247	g_total = (g(:,jlev,jcol)ssa(:,jlev,jcol)od(:,jlev,jcol) &
248	& + g_cloud(config%i_band_from_reordered_g_sw,jlev,jcol) &
249	& * ssa_cloud(config%i_band_from_reordered_g_sw,jlev,jcol) &
250	& * od_cloud_g) &
251	& / (ssa_total*od_total)
252	end where
253
254	! Apply delta-Eddington scaling to the cloud-aerosol-gas
255	! mixture
256	if (config%do_sw_delta_scaling_with_gases) then
257	call delta_eddington(od_total, ssa_total, g_total)
258	end if
259
260	! Compute cloudy-sky reflectance, transmittance etc at
261	! each model level
262	call calc_two_stream_gammas_sw(ng, &
263	& cos_sza, ssa_total, g_total, &
264	& gamma1, gamma2, gamma3)
265	call calc_reflectance_transmittance_sw(ng, &
266	& cos_sza, od_total, ssa_total, &
267	& gamma1, gamma2, gamma3, &
268	& reflectance(:,jlev), transmittance(:,jlev), &
269	& ref_dir(:,jlev), trans_dir_diff(:,jlev), &
270	& trans_dir_dir(:,jlev) )
271
272	end if
273	end do
274
275	! Use adding method to compute fluxes for an overcast sky
276	call adding_ica_sw(ng, nlev, incoming_sw(:,jcol), &
277	& albedo_diffuse(:,jcol), albedo_direct(:,jcol), &
278	& spread(cos_sza,1,ng), reflectance, transmittance, ref_dir, trans_dir_diff, &
279	& trans_dir_dir, flux_up, flux_dn_diffuse, flux_dn_direct)
280
281	! Store overcast broadband fluxes
282	flux%sw_up(jcol,:) = sum(flux_up,1)
283	if (allocated(flux%sw_dn_direct)) then
284	flux%sw_dn_direct(jcol,:) = sum(flux_dn_direct,1)
285	flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
286	& + flux%sw_dn_direct(jcol,:)
287	else
288	flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
289	& + sum(flux_dn_direct,1)
290	end if
291
292	! Likewise for surface spectral fluxes
293	flux%sw_dn_diffuse_surf_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
294	flux%sw_dn_direct_surf_g(:,jcol) = flux_dn_direct(:,nlev+1)
295
296	! Save the spectral fluxes if required
297	if (config%do_save_spectral_flux) then
298	call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_sw, &
299	& flux%sw_up_band(:,jcol,:))
300	call indexed_sum_profile(flux_dn_direct, config%i_spec_from_reordered_g_sw, &
301	& flux%sw_dn_band(:,jcol,:))
302	if (allocated(flux%sw_dn_direct_band)) then
303	flux%sw_dn_direct_band(:,jcol,:) &
304	& = flux%sw_dn_band(:,jcol,:)
305	end if
306	call add_indexed_sum_profile(flux_dn_diffuse, &
307	& config%i_spec_from_reordered_g_sw, &
308	& flux%sw_dn_band(:,jcol,:))
309	end if
310
311	else
312	! No cloud in profile and clear-sky fluxes already
313	! calculated: copy them over
314	flux%sw_up(jcol,:) = flux%sw_up_clear(jcol,:)
315	flux%sw_dn(jcol,:) = flux%sw_dn_clear(jcol,:)
316	if (allocated(flux%sw_dn_direct)) then
317	flux%sw_dn_direct(jcol,:) = flux%sw_dn_direct_clear(jcol,:)
318	end if
319	flux%sw_dn_diffuse_surf_g(:,jcol) = flux%sw_dn_diffuse_surf_clear_g(:,jcol)
320	flux%sw_dn_direct_surf_g(:,jcol) = flux%sw_dn_direct_surf_clear_g(:,jcol)
321
322	if (config%do_save_spectral_flux) then
323	flux%sw_up_band(:,jcol,:) = flux%sw_up_clear_band(:,jcol,:)
324	flux%sw_dn_band(:,jcol,:) = flux%sw_dn_clear_band(:,jcol,:)
325	if (allocated(flux%sw_dn_direct_band)) then
326	flux%sw_dn_direct_band(:,jcol,:) = flux%sw_dn_direct_clear_band(:,jcol,:)
327	end if
328	end if
329
330	end if ! Cloud is present in profile
331
332	else
333	! Set fluxes to zero if sun is below the horizon
334	flux%sw_up(jcol,:) = 0.0_jprb
335	flux%sw_dn(jcol,:) = 0.0_jprb
336	if (allocated(flux%sw_dn_direct)) then
337	flux%sw_dn_direct(jcol,:) = 0.0_jprb
338	end if
339	flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
340	flux%sw_dn_direct_surf_g(:,jcol) = 0.0_jprb
341
342	if (config%do_clear) then
343	flux%sw_up_clear(jcol,:) = 0.0_jprb
344	flux%sw_dn_clear(jcol,:) = 0.0_jprb
345	if (allocated(flux%sw_dn_direct_clear)) then
346	flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
347	end if
348	flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
349	flux%sw_dn_direct_surf_clear_g(:,jcol) = 0.0_jprb
350	end if
351
352	if (config%do_save_spectral_flux) then
353	flux%sw_dn_band(:,jcol,:) = 0.0_jprb
354	flux%sw_up_band(:,jcol,:) = 0.0_jprb
355	if (allocated(flux%sw_dn_direct_band)) then
356	flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb
357	end if
358	if (config%do_clear) then
359	flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb
360	flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb
361	if (allocated(flux%sw_dn_direct_clear_band)) then
362	flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb
363	end if
364	end if
365	end if
366
367	end if ! sun above horizon
368	end do
369
370	if (lhook) call dr_hook('radiation_homogeneous_sw:solver_homogeneous_sw',1,hook_handle)
371
372	end subroutine solver_homogeneous_sw
373
374	end module radiation_homogeneous_sw

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

Download in other formats:

Original Format