Context Navigation

radiation_config.F90 @ 5213

Last change on this file since 5213 was 4489, checked in by lguez, 20 months ago
Merge LMDZ_ECRad branch back into trunk!
File size: 86.6 KB

Line
1	! radiation_config.F90 - Derived type to configure the radiation scheme
2	!
3	! (C) Copyright 2014- 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-07-22 R. Hogan Added Yi et al. ice optics model
17	! 2017-10-23 R. Hogan Renamed single-character variables
18	! 2018-03-15 R. Hogan Added logicals controlling surface spectral treatment
19	! 2018-08-29 R. Hogan Added monochromatic single-scattering albedo / asymmetry factor
20	! 2018-09-03 R. Hogan Added min_cloud_effective_size
21	! 2018-09-04 R. Hogan Added encroachment_scaling
22	! 2018-09-13 R. Hogan Added IEncroachmentFractal
23	! 2019-01-02 R. Hogan Added Cloudless solvers
24	! 2019-01-14 R. Hogan Added out_of_bounds_[1,2,3]d for checker routines
25	! 2019-01-18 R. Hogan Added albedo weighting
26	! 2019-02-03 R. Hogan Added ability to fix out-of-physical-bounds inputs
27	! 2019-02-10 R. Hogan Renamed "encroachment" to "entrapment"
28	! 2020-05-18 R. Hogan Moved out_of_bounds_* to radiation_check.F90
29	! 2021-07-04 R. Hogan Numerous changes for ecCKD and general cloud/aerosol optics
30	!
31	! Note: The aim is for ecRad in the IFS to be as similar as possible
32	! to the offline version, so if you make any changes to this or any
33	! files in this directory, please inform Robin Hogan.
34	!
35
36	module radiation_config
37
38	use parkind1, only : jprb
39
40	use radiation_cloud_optics_data, only : cloud_optics_type
41	use radiation_general_cloud_optics_data, only : general_cloud_optics_type
42	use radiation_aerosol_optics_data, only : aerosol_optics_type
43	use radiation_pdf_sampler, only : pdf_sampler_type
44	use radiation_cloud_cover, only : OverlapName, &
45	& IOverlapMaximumRandom, IOverlapExponentialRandom, IOverlapExponential
46	use radiation_ecckd, only : ckd_model_type
47
48	implicit none
49	public
50
51	! Configuration codes: use C-style enumerators to avoid having to
52	! remember the numbers
53
54	! Solvers: can be specified for longwave and shortwave
55	! independently, except for "Homogeneous", which must be the same
56	! for both
57	enum, bind(c)
58	enumerator ISolverCloudless, ISolverHomogeneous, ISolverMcICA, &
59	& ISolverSpartacus, ISolverTripleclouds
60	end enum
61	character(len=*), parameter :: SolverName(0:4) = (/ 'Cloudless ', &
62	& 'Homogeneous ', &
63	& 'McICA ', &
64	& 'SPARTACUS ', &
65	& 'Tripleclouds' /)
66
67	! SPARTACUS shortwave solver can treat the reflection of radiation
68	! back up into different regions in various ways
69	enum, bind(c)
70	enumerator &
71	& IEntrapmentZero, & ! No entrapment, as Tripleclouds
72	& IEntrapmentEdgeOnly, & ! Only radiation passed through cloud edge is horizontally homogenized
73	& IEntrapmentExplicit, & ! Estimate horiz migration dist, account for fractal clouds
74	& IEntrapmentExplicitNonFractal, & ! As above but ignore fractal nature of clouds
75	& IEntrapmentMaximum ! Complete horizontal homogenization within regions (old SPARTACUS assumption)
76	end enum
77
78	! Names available in the radiation namelist for variable
79	! sw_entrapment_name
80	character(len=*), parameter :: EntrapmentName(0:4) = [ 'Zero ', &
81	& 'Edge-only ', &
82	& 'Explicit ', &
83	& 'Non-fractal', &
84	& 'Maximum ' ]
85	! For backwards compatibility, the radiation namelist also supports
86	! the equivalent variable sw_encroachment_name with the following
87	! names
88	character(len=*), parameter :: EncroachmentName(0:4) = [ 'Zero ', &
89	& 'Minimum ', &
90	& 'Fractal ', &
91	& 'Computed', &
92	& 'Maximum ' ]
93
94	! Two-stream models
95	! This is not configurable at run-time
96
97	! Gas models
98	enum, bind(c)
99	enumerator IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD
100	end enum
101	character(len=*), parameter :: GasModelName(0:2) = (/ 'Monochromatic', &
102	& 'RRTMG-IFS ', &
103	& 'ECCKD '/)
104
105	! Hydrometeor scattering models
106	enum, bind(c)
107	enumerator ILiquidModelMonochromatic, &
108	& ILiquidModelSOCRATES, ILiquidModelSlingo
109	end enum
110	character(len=*), parameter :: LiquidModelName(0:2) = (/ 'Monochromatic', &
111	& 'SOCRATES ', &
112	& 'Slingo ' /)
113
114	enum, bind(c)
115	enumerator IIceModelMonochromatic, IIceModelFu, &
116	& IIceModelBaran, IIceModelBaran2016, IIceModelBaran2017, &
117	& IIceModelYi
118	end enum
119	character(len=*), parameter :: IceModelName(0:5) = (/ 'Monochromatic', &
120	& 'Fu-IFS ', &
121	& 'Baran ', &
122	& 'Baran2016 ', &
123	& 'Baran2017 ', &
124	& 'Yi ' /)
125
126	! Cloud PDF distribution shapes
127	enum, bind(c)
128	enumerator IPdfShapeLognormal, IPdfShapeGamma
129	end enum
130	character(len=*), parameter :: PdfShapeName(0:1) = (/ 'Lognormal', &
131	& 'Gamma ' /)
132
133	! Maximum number of different aerosol types that can be provided
134	integer, parameter :: NMaxAerosolTypes = 256
135
136	! Maximum number of different cloud types that can be provided
137	integer, parameter :: NMaxCloudTypes = 12
138
139	! Maximum number of shortwave albedo and longwave emissivity
140	! intervals
141	integer, parameter :: NMaxAlbedoIntervals = 256
142
143	! Length of string buffer for printing config information
144	integer, parameter :: NPrintStringLen = 60
145
146	!---------------------------------------------------------------------
147	! Derived type containing all the configuration information needed
148	! to run the radiation scheme. The intention is that this is fixed
149	! for a given model run. The parameters are to list first those
150	! quantities that can be set directly by the user, for example using a
151	! namelist, and second those quantities that are computed afterwards
152	! from the user-supplied numbers, especially the details of the gas
153	! optics model.
154	type config_type
155	! USER-CONFIGURABLE PARAMETERS
156
157	! Override default solar spectrum
158	logical :: use_spectral_solar_scaling = .false.
159
160	! Directory in which gas, cloud and aerosol data files are to be
161	! found
162	character(len=511) :: directory_name = '.'
163
164	! If this is true then support arbitrary hydrometeor types (not
165	! just ice and liquid) and arbitrary spectral discretization (not
166	! just RRTMG). It is required that this is true if the ecCKD gas
167	! optics model is selected. General cloud optics has only been
168	! available from ecRad version 1.5.
169	logical :: use_general_cloud_optics = .true.
170
171	! If this is true then support aerosol properties at an arbitrary
172	! spectral discretization (not just RRTMG). It is required that
173	! this is true if the ecCKD gas optics model is selected.
174	logical :: use_general_aerosol_optics = .true.
175
176	! Cloud is deemed to be present in a layer if cloud fraction
177	! exceeds this value
178	real(jprb) :: cloud_fraction_threshold = 1.0e-6_jprb
179	! ...and total cloud water mixing ratio exceeds this value
180	real(jprb) :: cloud_mixing_ratio_threshold = 1.0e-9_jprb
181
182	! Overlap scheme
183	integer :: i_overlap_scheme = IOverlapExponentialRandom
184
185	! Use the Shonk et al. (2010) "beta" overlap parameter, rather
186	! than the "alpha" overlap parameter of Hogan and Illingworth
187	! (2000)?
188	logical :: use_beta_overlap = .false.
189
190	! Use a more vectorizable McICA cloud generator, at the expense of
191	! more random numbers being generated? This is the default on NEC
192	! SX.
193	#ifdef __SX__
194	logical :: use_vectorizable_generator = .true.
195	#else
196	logical :: use_vectorizable_generator = .false.
197	#endif
198
199	! Shape of sub-grid cloud water PDF
200	integer :: i_cloud_pdf_shape = IPdfShapeGamma
201
202	! The ratio of the overlap decorrelation length for cloud
203	! inhomogeneities to the overlap decorrelation length for cloud
204	! boundaries. Observations suggest this has a value of 0.5
205	! (e.g. from the decorrelation lengths of Hogan and Illingworth
206	! 2003 and Hogan and Illingworth 2000).
207	real(jprb) :: cloud_inhom_decorr_scaling = 0.5_jprb
208
209	! Factor controlling how much of the cloud edge length interfaces
210	! directly between the clear-sky region (region a) and the
211	! optically thick cloudy region (region c). If Lxy is the length
212	! of the interfaces between regions x and y, and Lab and Lbc have
213	! been computed already, then
214	! Lac=clear_to_thick_fraction*min(Lab,Lbc).
215	real(jprb) :: clear_to_thick_fraction = 0.0_jprb
216
217	! Factor allowing lateral transport when the sun is close to
218	! overhead; consider atand(overhead_sun_factor) to be the number
219	! of degrees that the sun angle is perturbed from zenith for the
220	! purposes of computing lateral transport. A value of up to 0.1
221	! seems to be necessary to account for the fact that some forward
222	! scattered radiation is treated as unscattered by delta-Eddington
223	! scaling; therefore it ought to have the chance to escape.
224	real(jprb) :: overhead_sun_factor = 0.0_jprb
225
226	! Minimum gas optical depth in a single layer at any wavelength,
227	! for stability
228	real(jprb) :: min_gas_od_lw = 1.0e-15_jprb
229	real(jprb) :: min_gas_od_sw = 0.0_jprb
230
231	! Maximum gas optical depth in a layer before that g-point will
232	! not be considered for 3D treatment: a limit is required to avoid
233	! expensive computation of matrix exponentials on matrices with
234	! large elements
235	real(jprb) :: max_gas_od_3d = 8.0_jprb
236
237	! Maximum total optical depth of a cloudy region for stability:
238	! optical depth will be capped at this value in the SPARTACUS
239	! solvers
240	real(jprb) :: max_cloud_od = 16.0_jprb
241
242	! How much longwave scattering is included?
243	logical :: do_lw_cloud_scattering = .true.
244	logical :: do_lw_aerosol_scattering = .true.
245
246	! Number of regions used to describe clouds and clear skies. A
247	! value of 2 means one clear and one cloudy region, so clouds are
248	! horizontally homogeneous, while a value of 3 means two cloudy
249	! regions with different optical depth, thereby representing
250	! inhomogeneity via the Shonk & Hogan (2008) "Tripleclouds"
251	! method.
252	integer :: nregions = 3
253
254	! Code specifying the solver to be used: use the enumerations
255	! defined above
256	integer :: i_solver_sw = ISolverMcICA
257	integer :: i_solver_lw = ISolverMcICA
258
259	! Do shortwave delta-Eddington scaling on the cloud-aerosol-gas
260	! mixture (as in the original IFS scheme), rather than the more
261	! correct approach of separately scaling the cloud and aerosol
262	! scattering properties before merging with gases. Note that
263	! .true. is not compatible with the SPARTACUS solver.
264	logical :: do_sw_delta_scaling_with_gases = .false.
265
266	! Codes describing the gas model
267	integer :: i_gas_model = IGasModelIFSRRTMG
268
269	! Optics if i_gas_model==IGasModelMonochromatic.
270	! The wavelength to use for the Planck function in metres. If this
271	! is positive then the output longwave fluxes will be in units of
272	! W m-2 um-1. If this is zero or negative (the default) then
273	! sigma*T^4 will be used and the output longwave fluxes will be in
274	! W m-2.
275	real(jprb) :: mono_lw_wavelength = -1.0_jprb
276	! Total zenith optical depth of the atmosphere in the longwave and
277	! shortwave, distributed vertically according to the pressure.
278	! Default is zero.
279	real(jprb) :: mono_lw_total_od = 0.0_jprb
280	real(jprb) :: mono_sw_total_od = 0.0_jprb
281	! Single-scattering albedo and asymmetry factor: values typical
282	! for liquid clouds with effective radius of 10 microns, at (SW)
283	! 0.55 micron wavelength and (LW) 10.7 microns wavelength
284	real(jprb) :: mono_sw_single_scattering_albedo = 0.999999_jprb
285	real(jprb) :: mono_sw_asymmetry_factor = 0.86_jprb
286	real(jprb) :: mono_lw_single_scattering_albedo = 0.538_jprb
287	real(jprb) :: mono_lw_asymmetry_factor = 0.925_jprb
288
289	! Codes describing particle scattering models
290	integer :: i_liq_model = ILiquidModelSOCRATES
291	integer :: i_ice_model = IIceModelBaran
292
293	! The mapping from albedo/emissivity intervals to SW/LW bands can
294	! either be done by finding the interval containing the central
295	! wavenumber of the band (nearest neighbour), or by a weighting
296	! according to the spectral overlap of each interval with each
297	! band
298	logical :: do_nearest_spectral_sw_albedo = .false.
299	logical :: do_nearest_spectral_lw_emiss = .false.
300
301	! User-defined monotonically increasing wavelength bounds (m)
302	! between input surface albedo/emissivity intervals. Implicitly
303	! the first interval starts at zero and the last ends at
304	! infinity. These must be set with define_sw_albedo_intervals and
305	! define_lw_emiss_intervals.
306	real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1) = -1.0_jprb
307	real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) = -1.0_jprb
308
309	! The index to the surface albedo/emissivity intervals for each of
310	! the wavelength bounds specified in sw_albedo_wavelength_bound
311	! and lw_emiss_wavelength_bound
312	integer :: i_sw_albedo_index(NMaxAlbedoIntervals) = 0
313	integer :: i_lw_emiss_index (NMaxAlbedoIntervals) = 0
314
315	! Do we compute longwave and/or shortwave radiation?
316	logical :: do_lw = .true.
317	logical :: do_sw = .true.
318
319	! Do we compute clear-sky fluxes and/or solar direct fluxes?
320	logical :: do_clear = .true.
321	logical :: do_sw_direct = .true.
322
323	! Do we include 3D effects?
324	logical :: do_3d_effects = .true.
325
326	character(len=511) :: cloud_type_name(NMaxCloudTypes) = ["","","","","","","","","","","",""]
327	! &
328	! & = ["mie_droplet ", &
329	! & "baum-general-habit-mixture_ice"]
330
331	! Spectral averaging method to use with generalized cloud optics;
332	! see Edwards & Slingo (1996) for definition. Experimentation
333	! with ecRad suggests that "thick" averaging is more accurate for
334	! both liquid and ice clouds.
335	logical :: use_thick_cloud_spectral_averaging(NMaxCloudTypes) &
336	& = [.true.,.true.,.true.,.true.,.true.,.true., &
337	& .true.,.true.,.true.,.true.,.true.,.true.]
338
339	! To what extent do we include "entrapment" effects in the
340	! SPARTACUS solver? This essentially means that in a situation
341	! like this
342	!
343	! 000111
344	! 222222
345	!
346	! Radiation downwelling from region 1 may be reflected back into
347	! region 0 due to some degree of homogenization of the radiation
348	! in region 2. Hogan and Shonk (2013) referred to this as
349	! "anomalous horizontal transport" for a 1D model, although for 3D
350	! calculations it is desirable to include at least some of it. The
351	! options are described by the IEntrapment* parameters above.
352	integer :: i_3d_sw_entrapment = IEntrapmentExplicit
353
354	! In the longwave, the equivalent process it either "on" (like
355	! maximum entrapment) or "off" (like zero entrapment):
356	logical :: do_3d_lw_multilayer_effects = .false.
357
358	! Do we account for the effective emissivity of the side of
359	! clouds?
360	logical :: do_lw_side_emissivity = .true.
361
362	! The 3D transfer rate "X" is such that if transport out of a
363	! region was the only process occurring then by the base of a
364	! layer only exp(-X) of the original flux would remain in that
365	! region. The transfer rate computed geometrically can be very
366	! high for the clear-sky regions in layers with high cloud
367	! fraction. For stability reasons it is necessary to provide a
368	! maximum possible 3D transfer rate.
369	real(jprb) :: max_3d_transfer_rate = 10.0_jprb
370
371	! It has also sometimes been found necessary to set a minimum
372	! cloud effective size for stability (metres)
373	real(jprb) :: min_cloud_effective_size = 100.0_jprb
374
375	! Given a horizontal migration distance, there is still
376	! uncertainty about how much entrapment occurs associated with how
377	! one assumes cloud boundaries line up in adjacent layers. This
378	! factor can be varied between 0.0 (the boundaries line up to the
379	! greatest extent possible given the overlap parameter) and 1.0
380	! (the boundaries line up to the minimum extent possible). In the
381	! Hogan et al. entrapment paper it is referred to as the overhang
382	! factor zeta, and a value of 0 matches the Monte Carlo
383	! calculations best.
384	real(jprb) :: overhang_factor = 0.0_jprb
385
386	! By default, the Meador & Weaver (1980) expressions are used
387	! instead of the matrix exponential whenever 3D effects can be
388	! neglected (e.g. cloud-free layers or clouds with infinitely
389	! large effective cloud size), but setting the following to true
390	! uses the matrix exponential everywhere, enabling the two
391	! methods to be compared. Note that Meador & Weaver will still be
392	! used for very optically thick g points where the matrix
393	! exponential can produce incorrect results.
394	logical :: use_expm_everywhere = .false.
395
396	! Aerosol descriptors: aerosol_type_mapping must be of length
397	! n_aerosol_types, and contains 0 if that type is to be ignored,
398	! positive numbers to map on to the indices of hydrophobic
399	! aerosols in the aerosol optics configuration file, and negative
400	! numbers to map on to (the negative of) the indices of
401	! hydrophilic aerosols in the configuration file.
402	logical :: use_aerosols = .false.
403	integer :: n_aerosol_types = 0
404	integer :: i_aerosol_type_map(NMaxAerosolTypes)
405
406	! Save the gas and cloud optical properties for each g point in
407	! "radiative_properties.nc"?
408	logical :: do_save_radiative_properties = .false.
409
410	! Save the flux profiles in each band?
411	logical :: do_save_spectral_flux = .false.
412
413	! Save the surface downwelling shortwave fluxes in each band?
414	logical :: do_surface_sw_spectral_flux = .true.
415
416	! Compute the longwave derivatives needed to apply the approximate
417	! radiation updates of Hogan and Bozzo (2015)
418	logical :: do_lw_derivatives = .false.
419
420	! Save the flux profiles in each g-point (overrides
421	! do_save_spectral_flux if TRUE)?
422	logical :: do_save_gpoint_flux = .false.
423
424	! In the IFS environment, setting up RRTM has already been done
425	! so not needed to do it again
426	logical :: do_setup_ifsrrtm = .true.
427
428	! In the IFS environment the old scheme has a bug in the Fu
429	! longwave ice optics whereby the single scattering albedo is one
430	! minus what it should be. Unfortunately fixing it makes
431	! forecasts worse. Setting the following to true reproduces the
432	! bug.
433	logical :: do_fu_lw_ice_optics_bug = .false.
434
435	! Control verbosity: 0=none (no output to standard output; write
436	! to standard error only if an error occurs), 1=warning, 2=info,
437	! 3=progress, 4=detailed, 5=debug. Separate settings for the
438	! setup of the scheme and the execution of it.
439	integer :: iverbosesetup = 2
440	integer :: iverbose = 1
441
442	! Are we doing radiative transfer in complex surface canopies
443	! (streets/vegetation), in which case tailored downward fluxes are
444	! needed at the top of the canopy?
445	logical :: do_canopy_fluxes_sw = .false.
446	logical :: do_canopy_fluxes_lw = .false.
447	! If so, do we use the full spectrum as in the atmosphere, or just
448	! the reduced spectrum in which the shortwave albedo and longwave
449	! emissivity are provided?
450	logical :: use_canopy_full_spectrum_sw = .false.
451	logical :: use_canopy_full_spectrum_lw = .false.
452	! Do we treat gas radiative transfer in streets/vegetation?
453	logical :: do_canopy_gases_sw = .false.
454	logical :: do_canopy_gases_lw = .false.
455
456	! Optics file names for overriding the ones generated from the
457	! other options. If these remain empty then the generated names
458	! will be used (see the "consolidate_config" routine below). If
459	! the user assigns one of these and it starts with a '/' character
460	! then that will be used instead. If the user assigns one and it
461	! doesn't start with a '/' character then it will be prepended by
462	! the contents of directory_name.
463	character(len=511) :: ice_optics_override_file_name = ''
464	character(len=511) :: liq_optics_override_file_name = ''
465	character(len=511) :: aerosol_optics_override_file_name = ''
466	character(len=511) :: gas_optics_sw_override_file_name = ''
467	character(len=511) :: gas_optics_lw_override_file_name = ''
468
469	! Optionally override the look-up table file for the cloud-water
470	! PDF used by the McICA solver
471	character(len=511) :: cloud_pdf_override_file_name = ''
472
473	! Do we compute cloud, aerosol and surface optical properties per
474	! g point? Not available with RRTMG gas optics model.
475	logical :: do_cloud_aerosol_per_sw_g_point = .true.
476	logical :: do_cloud_aerosol_per_lw_g_point = .true.
477
478	! Do we weight the mapping from surface emissivity/albedo to
479	! g-point/band weighting by a reference Planck function (more
480	! accurate) or weighting each wavenumber equally (less accurate
481	! but consistent with IFS Cycle 48r1 and earlier)?
482	logical :: do_weighted_surface_mapping = .true.
483
484	! COMPUTED PARAMETERS
485
486	! Users of this library should not edit these parameters directly;
487	! they are set by the "consolidate" routine
488
489	! Has "consolidate" been called?
490	logical :: is_consolidated = .false.
491
492	! Fraction of each g point in each wavenumber interval,
493	! dimensioned (n_wav_frac_[l\|s]w, n_g_[l\|s]w)
494	real(jprb), allocatable, dimension(:,:) :: g_frac_sw, g_frac_lw
495
496	! If the nearest surface albedo/emissivity interval is to be used
497	! for each SW/LW band then the following arrays will be allocated
498	! to the length of the number of bands and contain the index to
499	! the relevant interval
500	integer, allocatable, dimension(:) :: i_albedo_from_band_sw
501	integer, allocatable, dimension(:) :: i_emiss_from_band_lw
502
503	! ...alternatively, this matrix dimensioned
504	! (n_albedo_intervals,n_bands_sw) providing the weights needed for
505	! computing the albedo in each ecRad band from the albedo in each
506	! native albedo band - see radiation_single_level.F90
507	real(jprb), allocatable, dimension(:,:) :: sw_albedo_weights
508	! ...and similarly in the longwave, dimensioned
509	! (n_emiss_intervals,n_bands_lw)
510	real(jprb), allocatable, dimension(:,:) :: lw_emiss_weights
511
512	! Arrays of length the number of g-points that convert from
513	! g-point to the band index
514	integer, allocatable, dimension(:) :: i_band_from_g_lw
515	integer, allocatable, dimension(:) :: i_band_from_g_sw
516
517	! We allow for the possibility for g-points to be ordered in terms
518	! of likely absorption (weakest to strongest) across the shortwave
519	! or longwave spectrum, in order that in SPARTACUS we select only
520	! the first n g-points that will not have too large an absorption,
521	! and therefore matrix exponentials that are both finite and not
522	! too expensive to compute. The following two arrays map the
523	! reordered g-points to the original ones.
524	integer, allocatable, dimension(:) :: i_g_from_reordered_g_lw
525	integer, allocatable, dimension(:) :: i_g_from_reordered_g_sw
526
527	! The following map the reordered g-points to the bands
528	integer, allocatable, dimension(:) :: i_band_from_reordered_g_lw
529	integer, allocatable, dimension(:) :: i_band_from_reordered_g_sw
530
531	! The following map the reordered g-points to the spectral
532	! information being saved: if do_save_gpoint_flux==TRUE then this
533	! will map on to the original g points, but if only
534	! do_save_spectral_flux==TRUE then this will map on to the bands
535	integer, pointer, dimension(:) :: i_spec_from_reordered_g_lw
536	integer, pointer, dimension(:) :: i_spec_from_reordered_g_sw
537
538	! Number of spectral intervals used for the canopy radiative
539	! transfer calculation; they are either equal to
540	! n_albedo_intervals/n_emiss_intervals or n_g_sw/n_g_lw
541	integer :: n_canopy_bands_sw = 1
542	integer :: n_canopy_bands_lw = 1
543
544	! Data structures containing gas optics description in the case of
545	! ecCKD
546	type(ckd_model_type) :: gas_optics_sw, gas_optics_lw
547
548	! Data structure containing cloud scattering data
549	type(cloud_optics_type) :: cloud_optics
550
551	! Number of general cloud types, default liquid and ice
552	integer :: n_cloud_types = 2
553
554	! List of data structures (one per cloud type) containing cloud
555	! scattering data
556	type(general_cloud_optics_type), allocatable :: cloud_optics_sw(:)
557	type(general_cloud_optics_type), allocatable :: cloud_optics_lw(:)
558
559	! Data structure containing aerosol scattering data
560	type(aerosol_optics_type) :: aerosol_optics
561
562	! Object for sampling from a gamma or lognormal distribution
563	type(pdf_sampler_type) :: pdf_sampler
564
565	! Optics file names
566	character(len=511) :: ice_optics_file_name, &
567	& liq_optics_file_name, &
568	& aerosol_optics_file_name, &
569	& gas_optics_sw_file_name, &
570	& gas_optics_lw_file_name
571
572	! McICA PDF look-up table file name
573	character(len=511) :: cloud_pdf_file_name
574
575	! Number of gpoints and bands in the shortwave and longwave - set
576	! to zero as will be set properly later
577	integer :: n_g_sw = 0, n_g_lw = 0
578	integer :: n_bands_sw = 0, n_bands_lw = 0
579
580	! Number of spectral points to save (equal either to the number of
581	! g points or the number of bands
582	integer :: n_spec_sw = 0, n_spec_lw = 0
583
584	! Number of wavenumber intervals used to describe the mapping from
585	! g-points to wavenumber space
586	integer :: n_wav_frac_sw = 0, n_wav_frac_lw = 0
587
588	! Dimensions to store variables that are only needed if longwave
589	! scattering is included. "n_g_lw_if_scattering" is equal to
590	! "n_g_lw" if aerosols are allowed to scatter in the longwave,
591	! and zero otherwise. "n_bands_lw_if_scattering" is equal to
592	! "n_bands_lw" if clouds are allowed to scatter in the longwave,
593	! and zero otherwise.
594	integer :: n_g_lw_if_scattering = 0, n_bands_lw_if_scattering = 0
595
596	! Treat clouds as horizontally homogeneous within the gribox
597	logical :: is_homogeneous = .false.
598
599	! If the solvers are both "Cloudless" then we don't need to do any
600	! cloud processing
601	logical :: do_clouds = .true.
602
603	contains
604	procedure :: read => read_config_from_namelist
605	procedure :: consolidate => consolidate_config
606	procedure :: set => set_config
607	procedure :: print => print_config
608	procedure :: get_sw_weights
609	procedure :: define_sw_albedo_intervals
610	procedure :: define_lw_emiss_intervals
611	procedure :: set_aerosol_wavelength_mono
612	procedure :: consolidate_sw_albedo_intervals
613	procedure :: consolidate_lw_emiss_intervals
614
615	end type config_type
616
617	! procedure, private :: print_logical, print_real, print_int
618
619	contains
620
621
622	!---------------------------------------------------------------------
623	! This subroutine reads configuration data from a namelist file, and
624	! anything that is not in the namelists will be set to default
625	! values. If optional output argument "is_success" is present, then
626	! on error (e.g. missing file) it will be set to .false.; if this
627	! argument is missing then on error the program will be aborted. You
628	! may either specify the file_name or the unit of an open file to
629	! read, but not both.
630	subroutine read_config_from_namelist(this, file_name, unit, is_success)
631
632	use yomhook, only : lhook, dr_hook
633	use radiation_io, only : nulout, nulerr, nulrad, radiation_abort
634
635	class(config_type), intent(inout) :: this
636	character(*), intent(in), optional :: file_name
637	integer, intent(in), optional :: unit
638	logical, intent(out), optional :: is_success
639
640	integer :: iosopen, iosread ! Status after calling open and read
641
642	! The following variables are read from the namelists and map
643	! directly onto members of the config_type derived type
644
645	! To be read from the radiation_config namelist
646	logical :: do_sw, do_lw, do_clear, do_sw_direct
647	logical :: do_3d_effects, use_expm_everywhere, use_aerosols
648	logical :: use_general_cloud_optics, use_general_aerosol_optics
649	logical :: do_lw_side_emissivity
650	logical :: do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug
651	logical :: do_lw_aerosol_scattering, do_lw_cloud_scattering
652	logical :: do_save_radiative_properties, do_save_spectral_flux
653	logical :: do_save_gpoint_flux, do_surface_sw_spectral_flux
654	logical :: use_beta_overlap, do_lw_derivatives, use_vectorizable_generator
655	logical :: do_sw_delta_scaling_with_gases
656	logical :: do_canopy_fluxes_sw, do_canopy_fluxes_lw
657	logical :: use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw
658	logical :: do_canopy_gases_sw, do_canopy_gases_lw
659	logical :: do_cloud_aerosol_per_sw_g_point, do_cloud_aerosol_per_lw_g_point
660	logical :: do_weighted_surface_mapping
661	integer :: n_regions, iverbose, iverbosesetup, n_aerosol_types
662	real(jprb):: mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od
663	real(jprb):: mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo
664	real(jprb):: mono_lw_asymmetry_factor, mono_sw_asymmetry_factor
665	real(jprb):: cloud_inhom_decorr_scaling, cloud_fraction_threshold
666	real(jprb):: clear_to_thick_fraction, max_gas_od_3d, max_cloud_od
667	real(jprb):: cloud_mixing_ratio_threshold, overhead_sun_factor
668	real(jprb):: max_3d_transfer_rate, min_cloud_effective_size
669	real(jprb):: overhang_factor, encroachment_scaling
670	character(511) :: directory_name, aerosol_optics_override_file_name
671	character(511) :: liq_optics_override_file_name, ice_optics_override_file_name
672	character(511) :: cloud_pdf_override_file_name
673	character(511) :: gas_optics_sw_override_file_name, gas_optics_lw_override_file_name
674	character(63) :: liquid_model_name, ice_model_name, gas_model_name
675	character(63) :: sw_solver_name, lw_solver_name, overlap_scheme_name
676	character(63) :: sw_entrapment_name, sw_encroachment_name, cloud_pdf_shape_name
677	character(len=511) :: cloud_type_name(NMaxCloudTypes) = ["","","","","","","","","","","",""]
678	logical :: use_thick_cloud_spectral_averaging(NMaxCloudTypes) &
679	& = [.false.,.false.,.false.,.false.,.false.,.false., &
680	& .false.,.false.,.false.,.false.,.false.,.false.]
681	integer :: i_aerosol_type_map(NMaxAerosolTypes) ! More than 256 is an error
682
683	logical :: do_nearest_spectral_sw_albedo
684	logical :: do_nearest_spectral_lw_emiss
685	real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1)
686	real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1)
687	integer :: i_sw_albedo_index(NMaxAlbedoIntervals)
688	integer :: i_lw_emiss_index (NMaxAlbedoIntervals)
689
690	integer :: iunit ! Unit number of namelist file
691
692	namelist /radiation/ do_sw, do_lw, do_sw_direct, &
693	& do_3d_effects, do_lw_side_emissivity, do_clear, &
694	& do_save_radiative_properties, sw_entrapment_name, sw_encroachment_name, &
695	& do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug, &
696	& do_save_spectral_flux, do_save_gpoint_flux, &
697	& do_surface_sw_spectral_flux, do_lw_derivatives, &
698	& do_lw_aerosol_scattering, do_lw_cloud_scattering, &
699	& n_regions, directory_name, gas_model_name, &
700	& ice_optics_override_file_name, liq_optics_override_file_name, &
701	& aerosol_optics_override_file_name, cloud_pdf_override_file_name, &
702	& gas_optics_sw_override_file_name, gas_optics_lw_override_file_name, &
703	& liquid_model_name, ice_model_name, max_3d_transfer_rate, &
704	& min_cloud_effective_size, overhang_factor, encroachment_scaling, &
705	& use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw, &
706	& do_canopy_fluxes_sw, do_canopy_fluxes_lw, &
707	& do_canopy_gases_sw, do_canopy_gases_lw, &
708	& use_general_cloud_optics, use_general_aerosol_optics, &
709	& do_sw_delta_scaling_with_gases, overlap_scheme_name, &
710	& sw_solver_name, lw_solver_name, use_beta_overlap, use_vectorizable_generator, &
711	& use_expm_everywhere, iverbose, iverbosesetup, &
712	& cloud_inhom_decorr_scaling, cloud_fraction_threshold, &
713	& clear_to_thick_fraction, max_gas_od_3d, max_cloud_od, &
714	& cloud_mixing_ratio_threshold, overhead_sun_factor, &
715	& n_aerosol_types, i_aerosol_type_map, use_aerosols, &
716	& mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od, &
717	& mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo, &
718	& mono_lw_asymmetry_factor, mono_sw_asymmetry_factor, &
719	& cloud_pdf_shape_name, cloud_type_name, use_thick_cloud_spectral_averaging, &
720	& do_nearest_spectral_sw_albedo, do_nearest_spectral_lw_emiss, &
721	& sw_albedo_wavelength_bound, lw_emiss_wavelength_bound, &
722	& i_sw_albedo_index, i_lw_emiss_index, &
723	& do_cloud_aerosol_per_lw_g_point, &
724	& do_cloud_aerosol_per_sw_g_point, do_weighted_surface_mapping
725
726	real(jprb) :: hook_handle
727
728	if (lhook) call dr_hook('radiation_config:read',0,hook_handle)
729
730	! Copy default values from the original structure
731	do_sw = this%do_sw
732	do_lw = this%do_lw
733	do_sw_direct = this%do_sw_direct
734	do_3d_effects = this%do_3d_effects
735	do_3d_lw_multilayer_effects = this%do_3d_lw_multilayer_effects
736	do_lw_side_emissivity = this%do_lw_side_emissivity
737	do_clear = this%do_clear
738	do_lw_aerosol_scattering = this%do_lw_aerosol_scattering
739	do_lw_cloud_scattering = this%do_lw_cloud_scattering
740	do_sw_delta_scaling_with_gases = this%do_sw_delta_scaling_with_gases
741	do_fu_lw_ice_optics_bug = this%do_fu_lw_ice_optics_bug
742	do_canopy_fluxes_sw = this%do_canopy_fluxes_sw
743	do_canopy_fluxes_lw = this%do_canopy_fluxes_lw
744	use_canopy_full_spectrum_sw = this%use_canopy_full_spectrum_sw
745	use_canopy_full_spectrum_lw = this%use_canopy_full_spectrum_lw
746	do_canopy_gases_sw = this%do_canopy_gases_sw
747	do_canopy_gases_lw = this%do_canopy_gases_lw
748	n_regions = this%nregions
749	directory_name = this%directory_name
750	cloud_pdf_override_file_name = this%cloud_pdf_override_file_name
751	liq_optics_override_file_name = this%liq_optics_override_file_name
752	ice_optics_override_file_name = this%ice_optics_override_file_name
753	aerosol_optics_override_file_name = this%aerosol_optics_override_file_name
754	gas_optics_sw_override_file_name = this%gas_optics_sw_override_file_name
755	gas_optics_lw_override_file_name = this%gas_optics_lw_override_file_name
756	use_expm_everywhere = this%use_expm_everywhere
757	use_aerosols = this%use_aerosols
758	do_save_radiative_properties = this%do_save_radiative_properties
759	do_save_spectral_flux = this%do_save_spectral_flux
760	do_save_gpoint_flux = this%do_save_gpoint_flux
761	do_lw_derivatives = this%do_lw_derivatives
762	do_surface_sw_spectral_flux = this%do_surface_sw_spectral_flux
763	iverbose = this%iverbose
764	iverbosesetup = this%iverbosesetup
765	use_general_cloud_optics = this%use_general_cloud_optics
766	use_general_aerosol_optics = this%use_general_aerosol_optics
767	cloud_fraction_threshold = this%cloud_fraction_threshold
768	cloud_mixing_ratio_threshold = this%cloud_mixing_ratio_threshold
769	use_beta_overlap = this%use_beta_overlap
770	use_vectorizable_generator = this%use_vectorizable_generator
771	cloud_inhom_decorr_scaling = this%cloud_inhom_decorr_scaling
772	clear_to_thick_fraction = this%clear_to_thick_fraction
773	overhead_sun_factor = this%overhead_sun_factor
774	max_gas_od_3d = this%max_gas_od_3d
775	max_cloud_od = this%max_cloud_od
776	max_3d_transfer_rate = this%max_3d_transfer_rate
777	min_cloud_effective_size = this%min_cloud_effective_size
778	cloud_type_name = this%cloud_type_name
779	use_thick_cloud_spectral_averaging = this%use_thick_cloud_spectral_averaging
780
781	overhang_factor = this%overhang_factor
782	encroachment_scaling = -1.0_jprb
783	gas_model_name = '' !DefaultGasModelName
784	liquid_model_name = '' !DefaultLiquidModelName
785	ice_model_name = '' !DefaultIceModelName
786	sw_solver_name = '' !DefaultSwSolverName
787	lw_solver_name = '' !DefaultLwSolverName
788	sw_entrapment_name = ''
789	sw_encroachment_name = ''
790	overlap_scheme_name = ''
791	cloud_pdf_shape_name = ''
792	n_aerosol_types = this%n_aerosol_types
793	mono_lw_wavelength = this%mono_lw_wavelength
794	mono_lw_total_od = this%mono_lw_total_od
795	mono_sw_total_od = this%mono_sw_total_od
796	mono_lw_single_scattering_albedo = this%mono_lw_single_scattering_albedo
797	mono_sw_single_scattering_albedo = this%mono_sw_single_scattering_albedo
798	mono_lw_asymmetry_factor = this%mono_lw_asymmetry_factor
799	mono_sw_asymmetry_factor = this%mono_sw_asymmetry_factor
800	i_aerosol_type_map = this%i_aerosol_type_map
801	do_nearest_spectral_sw_albedo = this%do_nearest_spectral_sw_albedo
802	do_nearest_spectral_lw_emiss = this%do_nearest_spectral_lw_emiss
803	sw_albedo_wavelength_bound = this%sw_albedo_wavelength_bound
804	lw_emiss_wavelength_bound = this%lw_emiss_wavelength_bound
805	i_sw_albedo_index = this%i_sw_albedo_index
806	i_lw_emiss_index = this%i_lw_emiss_index
807	do_cloud_aerosol_per_lw_g_point = this%do_cloud_aerosol_per_lw_g_point
808	do_cloud_aerosol_per_sw_g_point = this%do_cloud_aerosol_per_sw_g_point
809	do_weighted_surface_mapping = this%do_weighted_surface_mapping
810
811	if (present(file_name) .and. present(unit)) then
812	write(nulerr,'(a)') '*** Error: cannot specify both file_name and unit in call to config_type%read'
813	call radiation_abort('Radiation configuration error')
814	else if (.not. present(file_name) .and. .not. present(unit)) then
815	write(nulerr,'(a)') '*** Error: neither file_name nor unit specified in call to config_type%read'
816	call radiation_abort('Radiation configuration error')
817	end if
818
819	if (present(file_name)) then
820	! Open the namelist file
821	iunit = nulrad
822	open(unit=iunit, iostat=iosopen, file=trim(file_name))
823	else
824	! Assume that iunit represents and open file
825	iosopen = 0
826	iunit = unit
827	end if
828
829	if (iosopen /= 0) then
830	! An error occurred opening the file
831	if (present(is_success)) then
832	is_success = .false.
833	! We now continue the subroutine so that the default values
834	! are placed in the config structure
835	else
836	write(nulerr,'(a,a,a)') '*** Error: namelist file "', &
837	& trim(file_name), '" not found'
838	call radiation_abort('Radiation configuration error')
839	end if
840	else
841
842	! This version exits correctly, but provides less information
843	! about how the namelist was incorrect
844	read(unit=iunit, iostat=iosread, nml=radiation)
845
846	! Depending on compiler this version provides more information
847	! about the error in the namelist
848	!read(unit=iunit, nml=radiation)
849
850	if (iosread /= 0) then
851	! An error occurred reading the file
852	if (present(is_success)) then
853	is_success = .false.
854	! We now continue the subroutine so that the default values
855	! are placed in the config structure
856	else if (present(file_name)) then
857	write(nulerr,'(a,a,a)') '*** Error reading namelist "radiation" from file "', &
858	& trim(file_name), '"'
859	close(unit=iunit)
860	call radiation_abort('Radiation configuration error')
861	else
862	write(nulerr,'(a,i0)') '*** Error reading namelist "radiation" from unit ', &
863	& iunit
864	call radiation_abort('Radiation configuration error')
865	end if
866	end if
867
868	if (present(file_name)) then
869	close(unit=iunit)
870	end if
871	end if
872
873	! Copy namelist data into configuration object
874
875	! Start with verbosity levels, which should be within limits
876	if (iverbose < 0) then
877	iverbose = 0
878	end if
879	this%iverbose = iverbose
880
881	if (iverbosesetup < 0) then
882	iverbosesetup = 0
883	end if
884	this%iverbosesetup = iverbosesetup
885
886	this%do_lw = do_lw
887	this%do_sw = do_sw
888	this%do_clear = do_clear
889	this%do_sw_direct = do_sw_direct
890	this%do_3d_effects = do_3d_effects
891	this%do_3d_lw_multilayer_effects = do_3d_lw_multilayer_effects
892	this%do_lw_side_emissivity = do_lw_side_emissivity
893	this%use_expm_everywhere = use_expm_everywhere
894	this%use_aerosols = use_aerosols
895	this%do_lw_cloud_scattering = do_lw_cloud_scattering
896	this%do_lw_aerosol_scattering = do_lw_aerosol_scattering
897	this%nregions = n_regions
898	this%do_surface_sw_spectral_flux = do_surface_sw_spectral_flux
899	this%do_sw_delta_scaling_with_gases = do_sw_delta_scaling_with_gases
900	this%do_fu_lw_ice_optics_bug = do_fu_lw_ice_optics_bug
901	this%do_canopy_fluxes_sw = do_canopy_fluxes_sw
902	this%do_canopy_fluxes_lw = do_canopy_fluxes_lw
903	this%use_canopy_full_spectrum_sw = use_canopy_full_spectrum_sw
904	this%use_canopy_full_spectrum_lw = use_canopy_full_spectrum_lw
905	this%do_canopy_gases_sw = do_canopy_gases_sw
906	this%do_canopy_gases_lw = do_canopy_gases_lw
907	this%mono_lw_wavelength = mono_lw_wavelength
908	this%mono_lw_total_od = mono_lw_total_od
909	this%mono_sw_total_od = mono_sw_total_od
910	this%mono_lw_single_scattering_albedo = mono_lw_single_scattering_albedo
911	this%mono_sw_single_scattering_albedo = mono_sw_single_scattering_albedo
912	this%mono_lw_asymmetry_factor = mono_lw_asymmetry_factor
913	this%mono_sw_asymmetry_factor = mono_sw_asymmetry_factor
914	this%use_beta_overlap = use_beta_overlap
915	this%use_vectorizable_generator = use_vectorizable_generator
916	this%cloud_inhom_decorr_scaling = cloud_inhom_decorr_scaling
917	this%clear_to_thick_fraction = clear_to_thick_fraction
918	this%overhead_sun_factor = overhead_sun_factor
919	this%max_gas_od_3d = max_gas_od_3d
920	this%max_cloud_od = max_cloud_od
921	this%max_3d_transfer_rate = max_3d_transfer_rate
922	this%min_cloud_effective_size = max(1.0e-6_jprb, min_cloud_effective_size)
923	this%cloud_type_name = cloud_type_name
924	this%use_thick_cloud_spectral_averaging = use_thick_cloud_spectral_averaging
925	if (encroachment_scaling >= 0.0_jprb) then
926	this%overhang_factor = encroachment_scaling
927	if (iverbose >= 1) then
928	write(nulout, '(a)') 'Warning: radiation namelist parameter "encroachment_scaling" is deprecated: use "overhang_factor"'
929	end if
930	else
931	this%overhang_factor = overhang_factor
932	end if
933	this%directory_name = directory_name
934	this%cloud_pdf_override_file_name = cloud_pdf_override_file_name
935	this%liq_optics_override_file_name = liq_optics_override_file_name
936	this%ice_optics_override_file_name = ice_optics_override_file_name
937	this%aerosol_optics_override_file_name = aerosol_optics_override_file_name
938	this%gas_optics_sw_override_file_name = gas_optics_sw_override_file_name
939	this%gas_optics_lw_override_file_name = gas_optics_lw_override_file_name
940	this%use_general_cloud_optics = use_general_cloud_optics
941	this%use_general_aerosol_optics = use_general_aerosol_optics
942	this%cloud_fraction_threshold = cloud_fraction_threshold
943	this%cloud_mixing_ratio_threshold = cloud_mixing_ratio_threshold
944	this%n_aerosol_types = n_aerosol_types
945	this%do_save_radiative_properties = do_save_radiative_properties
946	this%do_lw_derivatives = do_lw_derivatives
947	this%do_save_spectral_flux = do_save_spectral_flux
948	this%do_save_gpoint_flux = do_save_gpoint_flux
949	this%do_nearest_spectral_sw_albedo = do_nearest_spectral_sw_albedo
950	this%do_nearest_spectral_lw_emiss = do_nearest_spectral_lw_emiss
951	this%sw_albedo_wavelength_bound = sw_albedo_wavelength_bound
952	this%lw_emiss_wavelength_bound = lw_emiss_wavelength_bound
953	this%i_sw_albedo_index = i_sw_albedo_index
954	this%i_lw_emiss_index = i_lw_emiss_index
955	this%do_cloud_aerosol_per_lw_g_point = do_cloud_aerosol_per_lw_g_point
956	this%do_cloud_aerosol_per_sw_g_point = do_cloud_aerosol_per_sw_g_point
957	this%do_weighted_surface_mapping = do_weighted_surface_mapping
958
959	if (do_save_gpoint_flux) then
960	! Saving the fluxes every g-point overrides saving as averaged
961	! in a band, but this%do_save_spectral_flux needs to be TRUE as
962	! it is tested inside the solver routines to decide whether to
963	! save anything
964	this%do_save_spectral_flux = .true.
965	end if
966
967	! Determine liquid optics model
968	call get_enum_code(liquid_model_name, LiquidModelName, &
969	& 'liquid_model_name', this%i_liq_model)
970
971	! Determine ice optics model
972	call get_enum_code(ice_model_name, IceModelName, &
973	& 'ice_model_name', this%i_ice_model)
974
975	! Determine gas optics model
976	call get_enum_code(gas_model_name, GasModelName, &
977	& 'gas_model_name', this%i_gas_model)
978
979	! Determine solvers
980	call get_enum_code(sw_solver_name, SolverName, &
981	& 'sw_solver_name', this%i_solver_sw)
982	call get_enum_code(lw_solver_name, SolverName, &
983	& 'lw_solver_name', this%i_solver_lw)
984
985	if (len_trim(sw_encroachment_name) > 1) then
986	call get_enum_code(sw_encroachment_name, EncroachmentName, &
987	& 'sw_encroachment_name', this%i_3d_sw_entrapment)
988	write(nulout, '(a)') 'Warning: radiation namelist string "sw_encroachment_name" is deprecated: use "sw_entrapment_name"'
989	else
990	call get_enum_code(sw_entrapment_name, EntrapmentName, &
991	& 'sw_entrapment_name', this%i_3d_sw_entrapment)
992	end if
993
994	! Determine overlap scheme
995	call get_enum_code(overlap_scheme_name, OverlapName, &
996	& 'overlap_scheme_name', this%i_overlap_scheme)
997
998	! Determine cloud PDF shape
999	call get_enum_code(cloud_pdf_shape_name, PdfShapeName, &
1000	& 'cloud_pdf_shape_name', this%i_cloud_pdf_shape)
1001
1002	this%i_aerosol_type_map = 0
1003	if (this%use_aerosols) then
1004	this%i_aerosol_type_map(1:n_aerosol_types) &
1005	& = i_aerosol_type_map(1:n_aerosol_types)
1006	end if
1007
1008	! Will clouds be used at all?
1009	if ((this%do_sw .and. this%i_solver_sw /= ISolverCloudless) &
1010	& .or. (this%do_lw .and. this%i_solver_lw /= ISolverCloudless)) then
1011	this%do_clouds = .true.
1012	else
1013	this%do_clouds = .false.
1014	end if
1015
1016	if (this%i_gas_model == IGasModelIFSRRTMG &
1017	& .and. (this%use_general_cloud_optics &
1018	& .or. this%use_general_aerosol_optics)) then
1019	if (this%do_sw .and. this%do_cloud_aerosol_per_sw_g_point) then
1020	write(nulout,'(a)') 'Warning: RRTMG SW only supports cloud/aerosol/surface optical properties per band, not per g-point'
1021	this%do_cloud_aerosol_per_sw_g_point = .false.
1022	end if
1023	if (this%do_lw .and. this%do_cloud_aerosol_per_lw_g_point) then
1024	write(nulout,'(a)') 'Warning: RRTMG LW only supports cloud/aerosol/surface optical properties per band, not per g-point'
1025	this%do_cloud_aerosol_per_lw_g_point = .false.
1026	end if
1027	end if
1028
1029
1030	! Normal subroutine exit
1031	if (present(is_success)) then
1032	is_success = .true.
1033	end if
1034
1035	if (lhook) call dr_hook('radiation_config:read',1,hook_handle)
1036
1037	end subroutine read_config_from_namelist
1038
1039
1040	!---------------------------------------------------------------------
1041	! This routine is called by radiation_interface:setup_radiation and
1042	! it converts the user specified options into some more specific
1043	! data such as data file names
1044	subroutine consolidate_config(this)
1045
1046	use parkind1, only : jprd
1047	use yomhook, only : lhook, dr_hook
1048	use radiation_io, only : nulout, nulerr, radiation_abort
1049
1050	class(config_type), intent(inout) :: this
1051
1052	real(jprb) :: hook_handle
1053
1054	if (lhook) call dr_hook('radiation_config:consolidate',0,hook_handle)
1055
1056	! Check consistency of models
1057	if (this%do_canopy_fluxes_sw .and. .not. this%do_surface_sw_spectral_flux) then
1058	if (this%iverbosesetup >= 1) then
1059	write(nulout,'(a)') 'Warning: turning on do_surface_sw_spectral_flux as required by do_canopy_fluxes_sw'
1060	end if
1061	this%do_surface_sw_spectral_flux = .true.
1062	end if
1063
1064	! Will clouds be used at all?
1065	if ((this%do_sw .and. this%i_solver_sw /= ISolverCloudless) &
1066	& .or. (this%do_lw .and. this%i_solver_lw /= ISolverCloudless)) then
1067	this%do_clouds = .true.
1068	else
1069	this%do_clouds = .false.
1070	end if
1071
1072	! SPARTACUS only works with Exp-Ran overlap scheme
1073	if (( this%i_solver_sw == ISolverSPARTACUS &
1074	& .or. this%i_solver_lw == ISolverSPARTACUS &
1075	& .or. this%i_solver_sw == ISolverTripleclouds &
1076	& .or. this%i_solver_lw == ISolverTripleclouds) &
1077	& .and. this%i_overlap_scheme /= IOverlapExponentialRandom) then
1078	write(nulerr,'(a)') '*** Error: SPARTACUS/Tripleclouds solvers can only do Exponential-Random overlap'
1079	call radiation_abort('Radiation configuration error')
1080	end if
1081
1082	if (jprb < jprd .and. this%iverbosesetup >= 1 &
1083	& .and. (this%i_solver_sw == ISolverSPARTACUS &
1084	& .or. this%i_solver_lw == ISolverSPARTACUS)) then
1085	write(nulout,'(a)') 'Warning: the SPARTACUS solver may be unstable in single precision'
1086	end if
1087
1088	! If ecCKD gas optics model is being used set relevant file names
1089	if (this%i_gas_model == IGasModelECCKD) then
1090
1091	! This gas optics model requires the general cloud and
1092	! aerosol optics settings
1093	if (.not. this%use_general_cloud_optics) then
1094	write(nulerr,'(a)') '*** Error: ecCKD gas optics model requires general cloud optics'
1095	call radiation_abort('Radiation configuration error')
1096	end if
1097	if (.not. this%use_general_aerosol_optics) then
1098	write(nulerr,'(a)') '*** Error: ecCKD gas optics model requires general aerosol optics'
1099	call radiation_abort('Radiation configuration error')
1100	end if
1101
1102	if (len_trim(this%gas_optics_sw_override_file_name) > 0) then
1103	if (this%gas_optics_sw_override_file_name(1:1) == '/') then
1104	this%gas_optics_sw_file_name = trim(this%gas_optics_sw_override_file_name)
1105	else
1106	this%gas_optics_sw_file_name = trim(this%directory_name) &
1107	& // '/' // trim(this%gas_optics_sw_override_file_name)
1108	end if
1109	else
1110	! In the IFS, the gas optics files should be specified in
1111	! ifs/module/radiation_setup.F90, not here
1112	this%gas_optics_sw_file_name = trim(this%directory_name) &
1113	& // "/ecckd-1.0_sw_climate_rgb-32b_ckd-definition.nc"
1114	end if
1115
1116	if (len_trim(this%gas_optics_lw_override_file_name) > 0) then
1117	if (this%gas_optics_lw_override_file_name(1:1) == '/') then
1118	this%gas_optics_lw_file_name = trim(this%gas_optics_lw_override_file_name)
1119	else
1120	this%gas_optics_lw_file_name = trim(this%directory_name) &
1121	& // '/' // trim(this%gas_optics_lw_override_file_name)
1122	end if
1123	else
1124	! In the IFS, the gas optics files should be specified in
1125	! ifs/module/radiation_setup.F90, not here
1126	this%gas_optics_lw_file_name = trim(this%directory_name) &
1127	& // "/ecckd-1.0_lw_climate_fsck-32b_ckd-definition.nc"
1128	end if
1129
1130	end if
1131
1132	! Set aerosol optics file name
1133	if (len_trim(this%aerosol_optics_override_file_name) > 0) then
1134	if (this%aerosol_optics_override_file_name(1:1) == '/') then
1135	this%aerosol_optics_file_name = trim(this%aerosol_optics_override_file_name)
1136	else
1137	this%aerosol_optics_file_name = trim(this%directory_name) &
1138	& // '/' // trim(this%aerosol_optics_override_file_name)
1139	end if
1140	else
1141	! In the IFS, the aerosol optics file should be specified in
1142	! ifs/module/radiation_setup.F90, not here
1143	if (this%use_general_aerosol_optics) then
1144	this%aerosol_optics_file_name &
1145	& = trim(this%directory_name) // "/aerosol_ifs_48R1.nc"
1146	else
1147	this%aerosol_optics_file_name &
1148	& = trim(this%directory_name) // "/aerosol_ifs_rrtm_46R1_with_NI_AM.nc"
1149	end if
1150	end if
1151
1152	! Set liquid optics file name
1153	if (len_trim(this%liq_optics_override_file_name) > 0) then
1154	if (this%liq_optics_override_file_name(1:1) == '/') then
1155	this%liq_optics_file_name = trim(this%liq_optics_override_file_name)
1156	else
1157	this%liq_optics_file_name = trim(this%directory_name) &
1158	& // '/' // trim(this%liq_optics_override_file_name)
1159	end if
1160	else if (this%i_liq_model == ILiquidModelSOCRATES) then
1161	this%liq_optics_file_name &
1162	& = trim(this%directory_name) // "/socrates_droplet_scattering_rrtm.nc"
1163	else if (this%i_liq_model == ILiquidModelSlingo) then
1164	this%liq_optics_file_name &
1165	& = trim(this%directory_name) // "/slingo_droplet_scattering_rrtm.nc"
1166	end if
1167
1168	! Set ice optics file name
1169	if (len_trim(this%ice_optics_override_file_name) > 0) then
1170	if (this%ice_optics_override_file_name(1:1) == '/') then
1171	this%ice_optics_file_name = trim(this%ice_optics_override_file_name)
1172	else
1173	this%ice_optics_file_name = trim(this%directory_name) &
1174	& // '/' // trim(this%ice_optics_override_file_name)
1175	end if
1176	else if (this%i_ice_model == IIceModelFu) then
1177	this%ice_optics_file_name &
1178	& = trim(this%directory_name) // "/fu_ice_scattering_rrtm.nc"
1179	else if (this%i_ice_model == IIceModelBaran) then
1180	this%ice_optics_file_name &
1181	& = trim(this%directory_name) // "/baran_ice_scattering_rrtm.nc"
1182	else if (this%i_ice_model == IIceModelBaran2016) then
1183	this%ice_optics_file_name &
1184	& = trim(this%directory_name) // "/baran2016_ice_scattering_rrtm.nc"
1185	else if (this%i_ice_model == IIceModelBaran2017) then
1186	this%ice_optics_file_name &
1187	& = trim(this%directory_name) // "/baran2017_ice_scattering_rrtm.nc"
1188	else if (this%i_ice_model == IIceModelYi) then
1189	this%ice_optics_file_name &
1190	& = trim(this%directory_name) // "/yi_ice_scattering_rrtm.nc"
1191	end if
1192
1193	! Set cloud-water PDF look-up table file name
1194	if (len_trim(this%cloud_pdf_override_file_name) > 0) then
1195	if (this%cloud_pdf_override_file_name(1:1) == '/') then
1196	this%cloud_pdf_file_name = trim(this%cloud_pdf_override_file_name)
1197	else
1198	this%cloud_pdf_file_name = trim(this%directory_name) &
1199	& // '/' // trim(this%cloud_pdf_override_file_name)
1200	end if
1201	elseif (this%i_cloud_pdf_shape == IPdfShapeLognormal) then
1202	this%cloud_pdf_file_name = trim(this%directory_name) // "/mcica_lognormal.nc"
1203	else
1204	this%cloud_pdf_file_name = trim(this%directory_name) // "/mcica_gamma.nc"
1205	end if
1206
1207	! Aerosol data
1208	if (this%n_aerosol_types < 0 &
1209	& .or. this%n_aerosol_types > NMaxAerosolTypes) then
1210	write(nulerr,'(a,i0)') '*** Error: number of aerosol types must be between 0 and ', &
1211	& NMaxAerosolTypes
1212	call radiation_abort('Radiation configuration error')
1213	end if
1214
1215	if (this%use_aerosols .and. this%n_aerosol_types == 0) then
1216	if (this%iverbosesetup >= 2) then
1217	write(nulout, '(a)') 'Aerosols on but n_aerosol_types=0: optical properties to be computed outside ecRad'
1218	end if
1219	end if
1220
1221	! In the monochromatic case we need to override the liquid, ice
1222	! and aerosol models to ensure compatibility
1223	if (this%i_gas_model == IGasModelMonochromatic) then
1224	this%i_liq_model = ILiquidModelMonochromatic
1225	this%i_ice_model = IIceModelMonochromatic
1226	this%use_aerosols = .false.
1227	end if
1228
1229	! McICA solver currently can't store full profiles of spectral fluxes
1230	if (this%i_solver_sw == ISolverMcICA) then
1231	this%do_save_spectral_flux = .false.
1232	end if
1233
1234	if (this%i_solver_sw == ISolverSPARTACUS .and. this%do_sw_delta_scaling_with_gases) then
1235	write(nulerr,'(a)') '*** Error: SW delta-Eddington scaling with gases not possible with SPARTACUS solver'
1236	call radiation_abort('Radiation configuration error')
1237	end if
1238
1239	if ((this%do_lw .and. this%do_sw) .and. &
1240	& ( ( this%i_solver_sw == ISolverHomogeneous &
1241	& .and. this%i_solver_lw /= ISolverHomogeneous) &
1242	& .or. ( this%i_solver_sw /= ISolverHomogeneous &
1243	& .and. this%i_solver_lw == ISolverHomogeneous) &
1244	& ) ) then
1245	write(nulerr,'(a)') '*** Error: if one solver is "Homogeneous" then the other must be'
1246	call radiation_abort('Radiation configuration error')
1247	end if
1248
1249	! Set is_homogeneous if the active solvers are homogeneous, since
1250	! this affects how "in-cloud" water contents are computed
1251	if ( (this%do_sw .and. this%i_solver_sw == ISolverHomogeneous) &
1252	& .or. (this%do_lw .and. this%i_solver_lw == ISolverHomogeneous)) then
1253	this%is_homogeneous = .true.
1254	end if
1255
1256	this%is_consolidated = .true.
1257
1258	if (lhook) call dr_hook('radiation_config:consolidate',1,hook_handle)
1259
1260	end subroutine consolidate_config
1261
1262
1263	!---------------------------------------------------------------------
1264	! This subroutine sets members of the configuration object via
1265	! optional arguments, and any member not specified is left
1266	! untouched. Therefore, this should be called after taking data from
1267	! the namelist.
1268	subroutine set_config(config, directory_name, &
1269	& do_lw, do_sw, &
1270	& do_lw_aerosol_scattering, do_lw_cloud_scattering, &
1271	& do_sw_direct)
1272
1273	class(config_type), intent(inout):: config
1274	character(len=*), intent(in), optional :: directory_name
1275	logical, intent(in), optional :: do_lw, do_sw
1276	logical, intent(in), optional :: do_lw_aerosol_scattering
1277	logical, intent(in), optional :: do_lw_cloud_scattering
1278	logical, intent(in), optional :: do_sw_direct
1279
1280	if (present(do_lw)) then
1281	config%do_lw = do_lw
1282	end if
1283
1284	if(present(do_sw)) then
1285	config%do_sw = do_sw
1286	end if
1287
1288	if (present(do_sw_direct)) then
1289	config%do_sw_direct = do_sw_direct
1290	end if
1291
1292	if (present(directory_name)) then
1293	config%directory_name = trim(directory_name)
1294	end if
1295
1296	if (present(do_lw_aerosol_scattering)) then
1297	config%do_lw_aerosol_scattering = .true.
1298	end if
1299
1300	if (present(do_lw_cloud_scattering)) then
1301	config%do_lw_cloud_scattering = .true.
1302	end if
1303
1304	end subroutine set_config
1305
1306
1307	!---------------------------------------------------------------------
1308	! Print configuration information to standard output
1309	subroutine print_config(this, iverbose)
1310
1311	use radiation_io, only : nulout
1312
1313	class(config_type), intent(in) :: this
1314
1315	integer, optional, intent(in) :: iverbose
1316	integer :: i_local_verbose
1317
1318	if (present(iverbose)) then
1319	i_local_verbose = iverbose
1320	else
1321	i_local_verbose = this%iverbose
1322	end if
1323
1324	if (i_local_verbose >= 2) then
1325	!---------------------------------------------------------------------
1326	write(nulout, '(a)') 'General settings:'
1327	write(nulout, '(a,a,a)') ' Data files expected in "', &
1328	& trim(this%directory_name), '"'
1329	call print_logical(' Clear-sky calculations are', 'do_clear', this%do_clear)
1330	call print_logical(' Saving intermediate radiative properties', &
1331	& 'do_save_radiative_properties', this%do_save_radiative_properties)
1332	call print_logical(' Saving spectral flux profiles', &
1333	& 'do_save_spectral_flux', this%do_save_spectral_flux)
1334	call print_enum(' Gas model is', GasModelName, 'i_gas_model', &
1335	& this%i_gas_model)
1336	call print_logical(' Aerosols are', 'use_aerosols', this%use_aerosols)
1337	if (this%use_aerosols) then
1338	call print_logical(' General aerosol optics', &
1339	& 'use_general_aerosol_optics', this%use_general_aerosol_optics)
1340	end if
1341	if (this%do_clouds) then
1342	write(nulout,'(a)') ' Clouds are ON'
1343	else
1344	write(nulout,'(a)') ' Clouds are OFF'
1345	end if
1346	if (this%do_sw) then
1347	call print_logical(' Do cloud/aerosol/surface SW properties per g-point', &
1348	& 'do_cloud_aerosol_per_sw_g_point', this%do_cloud_aerosol_per_sw_g_point)
1349	end if
1350	if (this%do_lw) then
1351	call print_logical(' Do cloud/aerosol/surface LW properties per g-point', &
1352	& 'do_cloud_aerosol_per_lw_g_point', this%do_cloud_aerosol_per_lw_g_point)
1353	end if
1354
1355	!---------------------------------------------------------------------
1356	write(nulout, '(a)') 'Surface settings:'
1357	if (this%do_sw) then
1358	call print_logical(' Saving surface shortwave spectral fluxes', &
1359	& 'do_surface_sw_spectral_flux', this%do_surface_sw_spectral_flux)
1360	call print_logical(' Saving surface shortwave fluxes in abledo bands', &
1361	& 'do_canopy_fluxes_sw', this%do_canopy_fluxes_sw)
1362	end if
1363	if (this%do_lw) then
1364	call print_logical(' Saving surface longwave fluxes in emissivity bands', &
1365	& 'do_canopy_fluxes_lw', this%do_canopy_fluxes_lw)
1366	call print_logical(' Longwave derivative calculation is', &
1367	& 'do_lw_derivatives',this%do_lw_derivatives)
1368	end if
1369	if (this%do_sw) then
1370	call print_logical(' Nearest-neighbour spectral albedo mapping', &
1371	& 'do_nearest_spectral_sw_albedo', this%do_nearest_spectral_sw_albedo)
1372	end if
1373	if (this%do_lw) then
1374	call print_logical(' Nearest-neighbour spectral emissivity mapping', &
1375	& 'do_nearest_spectral_lw_emiss', this%do_nearest_spectral_lw_emiss)
1376	end if
1377	call print_logical(' Planck-weighted surface albedo/emiss mapping', &
1378	& 'do_weighted_surface_mapping', this%do_weighted_surface_mapping)
1379
1380	!---------------------------------------------------------------------
1381	if (this%do_clouds) then
1382	write(nulout, '(a)') 'Cloud settings:'
1383	call print_real(' Cloud fraction threshold', &
1384	& 'cloud_fraction_threshold', this%cloud_fraction_threshold)
1385	call print_real(' Cloud mixing-ratio threshold', &
1386	& 'cloud_mixing_ratio_threshold', this%cloud_mixing_ratio_threshold)
1387	call print_logical(' General cloud optics', &
1388	& 'use_general_cloud_optics', this%use_general_cloud_optics)
1389	if (.not. this%use_general_cloud_optics) then
1390	call print_enum(' Liquid optics scheme is', LiquidModelName, &
1391	& 'i_liq_model',this%i_liq_model)
1392	call print_enum(' Ice optics scheme is', IceModelName, &
1393	& 'i_ice_model',this%i_ice_model)
1394	if (this%i_ice_model == IIceModelFu) then
1395	call print_logical(' Longwave ice optics bug in Fu scheme is', &
1396	& 'do_fu_lw_ice_optics_bug',this%do_fu_lw_ice_optics_bug)
1397	end if
1398	end if
1399	call print_enum(' Cloud overlap scheme is', OverlapName, &
1400	& 'i_overlap_scheme',this%i_overlap_scheme)
1401	call print_logical(' Use "beta" overlap parameter is', &
1402	& 'use_beta_overlap', this%use_beta_overlap)
1403	call print_enum(' Cloud PDF shape is', PdfShapeName, &
1404	& 'i_cloud_pdf_shape',this%i_cloud_pdf_shape)
1405	call print_real(' Cloud inhom decorrelation scaling', &
1406	& 'cloud_inhom_decorr_scaling', this%cloud_inhom_decorr_scaling)
1407	end if
1408
1409	!---------------------------------------------------------------------
1410	write(nulout, '(a)') 'Solver settings:'
1411	if (this%do_sw) then
1412	call print_enum(' Shortwave solver is', SolverName, &
1413	& 'i_solver_sw', this%i_solver_sw)
1414
1415	if (this%i_gas_model == IGasModelMonochromatic) then
1416	call print_real(' Shortwave atmospheric optical depth', &
1417	& 'mono_sw_total_od', this%mono_sw_total_od)
1418	call print_real(' Shortwave particulate single-scattering albedo', &
1419	& 'mono_sw_single_scattering_albedo', &
1420	& this%mono_sw_single_scattering_albedo)
1421	call print_real(' Shortwave particulate asymmetry factor', &
1422	& 'mono_sw_asymmetry_factor', &
1423	& this%mono_sw_asymmetry_factor)
1424	end if
1425	call print_logical(' Shortwave delta scaling after merge with gases', &
1426	& 'do_sw_delta_scaling_with_gases', &
1427	& this%do_sw_delta_scaling_with_gases)
1428	else
1429	call print_logical(' Shortwave calculations are','do_sw',this%do_sw)
1430	end if
1431
1432	if (this%do_lw) then
1433	call print_enum(' Longwave solver is', SolverName, 'i_solver_lw', &
1434	& this%i_solver_lw)
1435
1436	if (this%i_gas_model == IGasModelMonochromatic) then
1437	if (this%mono_lw_wavelength > 0.0_jprb) then
1438	call print_real(' Longwave effective wavelength (m)', &
1439	& 'mono_lw_wavelength', this%mono_lw_wavelength)
1440	else
1441	write(nulout,'(a)') ' Longwave fluxes are broadband (mono_lw_wavelength<=0)'
1442	end if
1443	call print_real(' Longwave atmospheric optical depth', &
1444	& 'mono_lw_total_od', this%mono_lw_total_od)
1445	call print_real(' Longwave particulate single-scattering albedo', &
1446	& 'mono_lw_single_scattering_albedo', &
1447	& this%mono_lw_single_scattering_albedo)
1448	call print_real(' Longwave particulate asymmetry factor', &
1449	& 'mono_lw_asymmetry_factor', &
1450	& this%mono_lw_asymmetry_factor)
1451	end if
1452	call print_logical(' Longwave cloud scattering is', &
1453	& 'do_lw_cloud_scattering',this%do_lw_cloud_scattering)
1454	call print_logical(' Longwave aerosol scattering is', &
1455	& 'do_lw_aerosol_scattering',this%do_lw_aerosol_scattering)
1456	else
1457	call print_logical(' Longwave calculations are','do_lw', this%do_lw)
1458	end if
1459
1460	if (this%i_solver_sw == ISolverSpartacus &
1461	& .or. this%i_solver_lw == ISolverSpartacus) then
1462	write(nulout, '(a)') ' SPARTACUS options:'
1463	call print_integer(' Number of regions', 'n_regions', this%nregions)
1464	call print_real(' Max cloud optical depth per layer', &
1465	& 'max_cloud_od', this%max_cloud_od)
1466	call print_enum(' Shortwave entrapment is', EntrapmentName, &
1467	& 'i_3d_sw_entrapment', this%i_3d_sw_entrapment)
1468	call print_logical(' Multilayer longwave horizontal transport is', &
1469	'do_3d_lw_multilayer_effects', this%do_3d_lw_multilayer_effects)
1470	call print_logical(' Use matrix exponential everywhere is', &
1471	& 'use_expm_everywhere', this%use_expm_everywhere)
1472	call print_logical(' 3D effects are', 'do_3d_effects', &
1473	& this%do_3d_effects)
1474
1475	if (this%do_3d_effects) then
1476	call print_logical(' Longwave side emissivity parameterization is', &
1477	& 'do_lw_side_emissivity', this%do_lw_side_emissivity)
1478	call print_real(' Clear-to-thick edge fraction is', &
1479	& 'clear_to_thick_fraction', this%clear_to_thick_fraction)
1480	call print_real(' Overhead sun factor is', &
1481	& 'overhead_sun_factor', this%overhead_sun_factor)
1482	call print_real(' Max gas optical depth for 3D effects', &
1483	& 'max_gas_od_3d', this%max_gas_od_3d)
1484	call print_real(' Max 3D transfer rate', &
1485	& 'max_3d_transfer_rate', this%max_3d_transfer_rate)
1486	call print_real(' Min cloud effective size (m)', &
1487	& 'min_cloud_effective_size', this%min_cloud_effective_size)
1488	call print_real(' Overhang factor', &
1489	& 'overhang_factor', this%overhang_factor)
1490	end if
1491
1492	else if (this%i_solver_sw == ISolverMcICA &
1493	& .or. this%i_solver_lw == ISolverMcICA) then
1494	call print_logical(' Use vectorizable McICA cloud generator', &
1495	& 'use_vectorizable_generator', this%use_vectorizable_generator)
1496	end if
1497
1498	end if
1499
1500	end subroutine print_config
1501
1502
1503
1504	!---------------------------------------------------------------------
1505	! In order to estimate UV and photosynthetically active radiation,
1506	! we need weighted sum of fluxes considering wavelength range
1507	! required. This routine returns information for how to correctly
1508	! weight output spectral fluxes for a range of input wavelengths.
1509	! Note that this is approximate; internally it may be assumed that
1510	! the energy is uniformly distributed in wavenumber space, for
1511	! example. If the character string "weighting_name" is present, and
1512	! iverbose>=2, then information on the weighting will be provided on
1513	! nulout.
1514	subroutine get_sw_weights(this, wavelength1, wavelength2, &
1515	& nweights, iband, weight, weighting_name)
1516
1517	use parkind1, only : jprb
1518	use radiation_io, only : nulout, nulerr, radiation_abort
1519	use radiation_spectral_definition, only : SolarReferenceTemperature
1520
1521	class(config_type), intent(in) :: this
1522	! Range of wavelengths to get weights for (m)
1523	real(jprb), intent(in) :: wavelength1, wavelength2
1524	! Output number of weights needed
1525	integer, intent(out) :: nweights
1526	! Only write to the first nweights of these arrays: they contain
1527	! the indices to the non-zero bands, and the weight in each of
1528	! those bands
1529	integer, intent(out) :: iband(:)
1530	real(jprb), intent(out) :: weight(:)
1531	character(len=*), optional, intent(in) :: weighting_name
1532
1533	real(jprb), allocatable :: mapping(:,:)
1534
1535	! Internally we deal with wavenumber
1536	real(jprb) :: wavenumber1, wavenumber2 ! cm-1
1537
1538	real(jprb) :: wavenumber1_band, wavenumber2_band ! cm-1
1539
1540	integer :: jband ! Loop index for spectral band
1541
1542	if (this%n_bands_sw <= 0) then
1543	write(nulerr,'(a)') '*** Error: get_sw_weights called before number of shortwave bands set'
1544	call radiation_abort('Radiation configuration error')
1545	end if
1546
1547	! Convert wavelength range (m) to wavenumber (cm-1)
1548	wavenumber1 = 0.01_jprb / wavelength2
1549	wavenumber2 = 0.01_jprb / wavelength1
1550
1551	call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(SolarReferenceTemperature, &
1552	& [wavelength1, wavelength2], [1, 2, 3], mapping, &
1553	& use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point), use_fluxes=.true.)
1554
1555	! "mapping" now contains a 3*nband matrix, where mapping(2,:)
1556	! contains the weights of interest. We now find the non-zero weights
1557	nweights = 0
1558	do jband = 1,size(mapping,2)
1559	if (mapping(2,jband) > 0.0_jprb) then
1560	nweights = nweights+1
1561	iband(nweights) = jband;
1562	weight(nweights) = mapping(2,jband)
1563	end if
1564	end do
1565
1566	if (nweights == 0) then
1567	write(nulerr,'(a,e8.4,a,e8.4,a)') '*** Error: wavelength range ', &
1568	& wavelength1, ' to ', wavelength2, ' m is outside shortwave band'
1569	call radiation_abort('Radiation configuration error')
1570	else if (this%iverbosesetup >= 2 .and. present(weighting_name)) then
1571	write(nulout,'(a,a,a,f6.0,a,f6.0,a)') 'Spectral weights for ', &
1572	& weighting_name, ' (', wavenumber1, ' to ', &
1573	& wavenumber2, ' cm-1):'
1574	if (this%do_cloud_aerosol_per_sw_g_point) then
1575	do jband = 1, nweights
1576	write(nulout, '(a,i0,a,f8.4)') ' Shortwave g point ', iband(jband), ': ', weight(jband)
1577	end do
1578	else
1579	do jband = 1, nweights
1580	wavenumber1_band = this%gas_optics_sw%spectral_def%wavenumber1_band(iband(jband))
1581	wavenumber2_band = this%gas_optics_sw%spectral_def%wavenumber2_band(iband(jband))
1582	write(nulout, '(a,i0,a,f6.0,a,f6.0,a,f8.4)') ' Shortwave band ', &
1583	& iband(jband), ' (', wavenumber1_band, ' to ', &
1584	& wavenumber2_band, ' cm-1): ', weight(jband)
1585	end do
1586	end if
1587	end if
1588
1589	end subroutine get_sw_weights
1590
1591
1592	!---------------------------------------------------------------------
1593	! The input shortwave surface albedo coming in is likely to be in
1594	! different spectral intervals to the gas model in the radiation
1595	! scheme. We assume that the input albedo is defined within
1596	! "ninterval" spectral intervals covering the wavelength range 0 to
1597	! infinity, but allow for the possibility that two intervals may be
1598	! indexed back to the same albedo band.
1599	subroutine define_sw_albedo_intervals(this, ninterval, wavelength_bound, &
1600	& i_intervals, do_nearest)
1601
1602	use radiation_io, only : nulerr, radiation_abort
1603	use radiation_spectral_definition, only : SolarReferenceTemperature
1604
1605	class(config_type), intent(inout) :: this
1606	! Number of spectral intervals in which albedo is defined
1607	integer, intent(in) :: ninterval
1608	! Monotonically increasing wavelength bounds between intervals,
1609	! not including the outer bounds (which are assumed to be zero and
1610	! infinity)
1611	real(jprb), intent(in) :: wavelength_bound(ninterval-1)
1612	! The albedo indices corresponding to each interval
1613	integer, intent(in) :: i_intervals(ninterval)
1614	logical, optional, intent(in) :: do_nearest
1615
1616	if (ninterval > NMaxAlbedoIntervals) then
1617	write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, &
1618	& ' albedo intervals exceeds maximum of ', NMaxAlbedoIntervals
1619	call radiation_abort('Radiation configuration error')
1620	end if
1621
1622	if (present(do_nearest)) then
1623	this%do_nearest_spectral_sw_albedo = do_nearest
1624	else
1625	this%do_nearest_spectral_sw_albedo = .false.
1626	end if
1627	this%sw_albedo_wavelength_bound(1:ninterval-1) = wavelength_bound(1:ninterval-1)
1628	this%sw_albedo_wavelength_bound(ninterval:) = -1.0_jprb
1629	this%i_sw_albedo_index(1:ninterval) = i_intervals(1:ninterval)
1630	this%i_sw_albedo_index(ninterval+1:) = 0
1631
1632	! If this routine is called before setup_radiation then the
1633	! spectral intervals are not yet known
1634	! consolidate_sw_albedo_intervals is called later. Otherwise it
1635	! is called immediately and overwrites any existing mapping.
1636	if (this%is_consolidated) then
1637	call this%consolidate_sw_albedo_intervals
1638	end if
1639
1640	end subroutine define_sw_albedo_intervals
1641
1642
1643	!---------------------------------------------------------------------
1644	! As define_sw_albedo_intervals but for longwave emissivity
1645	subroutine define_lw_emiss_intervals(this, ninterval, wavelength_bound, &
1646	& i_intervals, do_nearest)
1647
1648	use radiation_io, only : nulerr, radiation_abort
1649	use radiation_spectral_definition, only : TerrestrialReferenceTemperature
1650
1651	class(config_type), intent(inout) :: this
1652	! Number of spectral intervals in which emissivity is defined
1653	integer, intent(in) :: ninterval
1654	! Monotonically increasing wavelength bounds between intervals,
1655	! not including the outer bounds (which are assumed to be zero and
1656	! infinity)
1657	real(jprb), intent(in) :: wavelength_bound(ninterval-1)
1658	! The emissivity indices corresponding to each interval
1659	integer, intent(in) :: i_intervals(ninterval)
1660	logical, optional, intent(in) :: do_nearest
1661
1662	if (ninterval > NMaxAlbedoIntervals) then
1663	write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, &
1664	& ' emissivity intervals exceeds maximum of ', NMaxAlbedoIntervals
1665	call radiation_abort('Radiation configuration error')
1666	end if
1667
1668	if (present(do_nearest)) then
1669	this%do_nearest_spectral_lw_emiss = do_nearest
1670	else
1671	this%do_nearest_spectral_lw_emiss = .false.
1672	end if
1673	this%lw_emiss_wavelength_bound(1:ninterval-1) = wavelength_bound(1:ninterval-1)
1674	this%lw_emiss_wavelength_bound(ninterval:) = -1.0_jprb
1675	this%i_lw_emiss_index(1:ninterval) = i_intervals(1:ninterval)
1676	this%i_lw_emiss_index(ninterval+1:) = 0
1677
1678	if (this%is_consolidated) then
1679	call this%consolidate_lw_emiss_intervals
1680	end if
1681
1682	end subroutine define_lw_emiss_intervals
1683
1684
1685	!---------------------------------------------------------------------
1686	! Set the wavelengths (m) at which monochromatic aerosol properties
1687	! are required. This routine must be called before consolidation of
1688	! settings.
1689	subroutine set_aerosol_wavelength_mono(this, wavelength_mono)
1690
1691	use radiation_io, only : nulerr, radiation_abort
1692
1693	class(config_type), intent(inout) :: this
1694	real(jprb), intent(in) :: wavelength_mono(:)
1695
1696	if (this%is_consolidated) then
1697	write(nulerr,'(a)') '*** Errror: set_aerosol_wavelength_mono must be called before setup_radiation'
1698	call radiation_abort('Radiation configuration error')
1699	end if
1700
1701	if (allocated(this%aerosol_optics%wavelength_mono)) then
1702	deallocate(this%aerosol_optics%wavelength_mono)
1703	end if
1704	allocate(this%aerosol_optics%wavelength_mono(size(wavelength_mono)))
1705	this%aerosol_optics%wavelength_mono = wavelength_mono
1706
1707	end subroutine set_aerosol_wavelength_mono
1708
1709
1710	!---------------------------------------------------------------------
1711	! Consolidate the surface shortwave albedo intervals with the
1712	! band/g-point intervals
1713	subroutine consolidate_sw_albedo_intervals(this)
1714
1715	use radiation_io, only : nulout
1716	use radiation_spectral_definition, only : SolarReferenceTemperature
1717
1718	class(config_type), intent(inout) :: this
1719
1720	integer :: ninterval, jint, jband
1721
1722	! Count the number of albedo/emissivity intervals
1723	ninterval = 0
1724	do jint = 1,NMaxAlbedoIntervals
1725	if (this%i_sw_albedo_index(jint) > 0) then
1726	ninterval = jint
1727	else
1728	exit
1729	end if
1730	end do
1731
1732	if (ninterval < 1) then
1733	! The user has not specified shortwave albedo bands - assume
1734	! only one
1735	ninterval = 1
1736	this%i_sw_albedo_index(1) = 1
1737	this%i_sw_albedo_index(2:) = 0
1738	if (this%use_canopy_full_spectrum_sw) then
1739	this%n_canopy_bands_sw = this%n_g_sw
1740	else
1741	this%n_canopy_bands_sw = 1
1742	end if
1743	else
1744	if (this%use_canopy_full_spectrum_sw) then
1745	this%n_canopy_bands_sw = this%n_g_sw
1746	else
1747	this%n_canopy_bands_sw = maxval(this%i_sw_albedo_index(1:ninterval))
1748	end if
1749	end if
1750
1751	if (this%do_weighted_surface_mapping) then
1752	call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(SolarReferenceTemperature, &
1753	& this%sw_albedo_wavelength_bound(1:ninterval-1), this%i_sw_albedo_index(1:ninterval), &
1754	& this%sw_albedo_weights, use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point))
1755	else
1756	! Weight each wavenumber equally as in IFS Cycles 48 and earlier
1757	call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(-1.0_jprb, &
1758	& this%sw_albedo_wavelength_bound(1:ninterval-1), this%i_sw_albedo_index(1:ninterval), &
1759	& this%sw_albedo_weights, use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point))
1760	end if
1761
1762	! Legacy method uses input band with largest weight
1763	if (this%do_nearest_spectral_sw_albedo) then
1764	allocate(this%i_albedo_from_band_sw(this%n_bands_sw))
1765	this%i_albedo_from_band_sw = maxloc(this%sw_albedo_weights, dim=1)
1766	end if
1767
1768	if (this%iverbosesetup >= 2) then
1769	write(nulout, '(a)') 'Surface shortwave albedo'
1770	if (.not. this%do_nearest_spectral_sw_albedo) then
1771	call this%gas_optics_sw%spectral_def%print_mapping_from_bands(this%sw_albedo_weights, &
1772	& use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point))
1773	else if (ninterval <= 1) then
1774	write(nulout, '(a)') 'All shortwave bands will use the same albedo'
1775	else
1776	write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_albedo_from_band_sw), &
1777	& ' shortwave intervals to albedo intervals:'
1778	do jband = 1,size(this%i_albedo_from_band_sw)
1779	write(nulout,'(a,i0)',advance='no') ' ', this%i_albedo_from_band_sw(jband)
1780	end do
1781	write(nulout, '()')
1782	end if
1783	end if
1784
1785	end subroutine consolidate_sw_albedo_intervals
1786
1787
1788	!---------------------------------------------------------------------
1789	! Consolidate the surface longwave emissivity intervals with the
1790	! band/g-point intervals
1791	subroutine consolidate_lw_emiss_intervals(this)
1792
1793	use radiation_io, only : nulout
1794	use radiation_spectral_definition, only : TerrestrialReferenceTemperature
1795
1796	class(config_type), intent(inout) :: this
1797
1798	integer :: ninterval, jint, jband
1799
1800	! Count the number of albedo/emissivity intervals
1801	ninterval = 0
1802	do jint = 1,NMaxAlbedoIntervals
1803	if (this%i_lw_emiss_index(jint) > 0) then
1804	ninterval = jint
1805	else
1806	exit
1807	end if
1808	end do
1809
1810	if (ninterval < 1) then
1811	! The user has not specified longwave emissivity bands - assume
1812	! only one
1813	ninterval = 1
1814	this%i_lw_emiss_index(1) = 1
1815	this%i_lw_emiss_index(2:) = 0
1816	if (this%use_canopy_full_spectrum_sw) then
1817	this%n_canopy_bands_lw = this%n_g_lw
1818	else
1819	this%n_canopy_bands_lw = 1
1820	end if
1821	else
1822	if (this%use_canopy_full_spectrum_lw) then
1823	this%n_canopy_bands_lw = this%n_g_lw
1824	else
1825	this%n_canopy_bands_lw = maxval(this%i_lw_emiss_index(1:ninterval))
1826	end if
1827	end if
1828
1829	if (this%do_weighted_surface_mapping) then
1830	call this%gas_optics_lw%spectral_def%calc_mapping_from_bands(TerrestrialReferenceTemperature, &
1831	& this%lw_emiss_wavelength_bound(1:ninterval-1), this%i_lw_emiss_index(1:ninterval), &
1832	& this%lw_emiss_weights, use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point))
1833	else
1834	! Weight each wavenumber equally as in IFS Cycles 48 and earlier
1835	call this%gas_optics_lw%spectral_def%calc_mapping_from_bands(-1.0_jprb, &
1836	& this%lw_emiss_wavelength_bound(1:ninterval-1), this%i_lw_emiss_index(1:ninterval), &
1837	& this%lw_emiss_weights, use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point))
1838	end if
1839
1840	! Legacy method uses input band with largest weight
1841	if (this%do_nearest_spectral_lw_emiss) then
1842	allocate(this%i_emiss_from_band_lw(this%n_bands_lw))
1843	this%i_emiss_from_band_lw = maxloc(this%lw_emiss_weights, dim=1)
1844	end if
1845
1846	if (this%iverbosesetup >= 2) then
1847	write(nulout, '(a)') 'Surface longwave emissivity'
1848	if (.not. this%do_nearest_spectral_lw_emiss) then
1849	call this%gas_optics_lw%spectral_def%print_mapping_from_bands(this%lw_emiss_weights, &
1850	& use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point))
1851	else if (ninterval <= 1) then
1852	write(nulout, '(a)') 'All longwave bands will use the same emissivty'
1853	else
1854	write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_emiss_from_band_lw), &
1855	& ' longwave intervals to emissivity intervals:'
1856	do jband = 1,size(this%i_emiss_from_band_lw)
1857	write(nulout,'(a,i0)',advance='no') ' ', this%i_emiss_from_band_lw(jband)
1858	end do
1859	write(nulout, '()')
1860	end if
1861	end if
1862
1863	end subroutine consolidate_lw_emiss_intervals
1864
1865
1866	!---------------------------------------------------------------------
1867	! Return the 0-based index for str in enum_str, or abort if it is
1868	! not found
1869	subroutine get_enum_code(str, enum_str, var_name, icode)
1870
1871	use radiation_io, only : nulerr, radiation_abort
1872
1873	character(len=*), intent(in) :: str
1874	character(len=*), intent(in) :: enum_str(0:)
1875	character(len=*), intent(in) :: var_name
1876	integer, intent(out) :: icode
1877
1878	integer :: jc
1879	logical :: is_not_found
1880
1881	! If string is empty then we don't modify icode but assume it has
1882	! a sensible default value
1883	if (len_trim(str) > 1) then
1884	is_not_found = .true.
1885
1886	do jc = 0,size(enum_str)-1
1887	if (trim(str) == trim(enum_str(jc))) then
1888	icode = jc
1889	is_not_found = .false.
1890	exit
1891	end if
1892	end do
1893	if (is_not_found) then
1894	write(nulerr,'(a,a,a,a,a)',advance='no') '*** Error: ', trim(var_name), &
1895	& ' must be one of: "', enum_str(0), '"'
1896	do jc = 1,size(enum_str)-1
1897	write(nulerr,'(a,a,a)',advance='no') ', "', trim(enum_str(jc)), '"'
1898	end do
1899	write(nulerr,'(a)') ''
1900	call radiation_abort('Radiation configuration error')
1901	end if
1902	end if
1903
1904	end subroutine get_enum_code
1905
1906
1907	!---------------------------------------------------------------------
1908	! Print one line of information: logical
1909	subroutine print_logical(message_str, name, val)
1910	use radiation_io, only : nulout
1911	character(len=*), intent(in) :: message_str
1912	character(len=*), intent(in) :: name
1913	logical, intent(in) :: val
1914	character(4) :: on_or_off
1915	character(NPrintStringLen) :: str
1916	if (val) then
1917	on_or_off = ' ON '
1918	else
1919	on_or_off = ' OFF'
1920	end if
1921	write(str, '(a,a4)') message_str, on_or_off
1922	write(nulout,'(a,a,a,a,l1,a)') str, ' (', name, '=', val,')'
1923	end subroutine print_logical
1924
1925
1926	!---------------------------------------------------------------------
1927	! Print one line of information: integer
1928	subroutine print_integer(message_str, name, val)
1929	use radiation_io, only : nulout
1930	character(len=*), intent(in) :: message_str
1931	character(len=*), intent(in) :: name
1932	integer, intent(in) :: val
1933	character(NPrintStringLen) :: str
1934	write(str, '(a,a,i0)') message_str, ' = ', val
1935	write(nulout,'(a,a,a,a)') str, ' (', name, ')'
1936	end subroutine print_integer
1937
1938
1939	!---------------------------------------------------------------------
1940	! Print one line of information: real
1941	subroutine print_real(message_str, name, val)
1942	use parkind1, only : jprb
1943	use radiation_io, only : nulout
1944	character(len=*), intent(in) :: message_str
1945	character(len=*), intent(in) :: name
1946	real(jprb), intent(in) :: val
1947	character(NPrintStringLen) :: str
1948	write(str, '(a,a,g8.3)') message_str, ' = ', val
1949	write(nulout,'(a,a,a,a)') str, ' (', name, ')'
1950	end subroutine print_real
1951
1952
1953	!---------------------------------------------------------------------
1954	! Print one line of information: enum
1955	subroutine print_enum(message_str, enum_str, name, val)
1956	use radiation_io, only : nulout
1957	character(len=*), intent(in) :: message_str
1958	character(len=*), intent(in) :: enum_str(0:)
1959	character(len=*), intent(in) :: name
1960	integer, intent(in) :: val
1961	character(NPrintStringLen) :: str
1962	write(str, '(a,a,a,a)') message_str, ' "', trim(enum_str(val)), '"'
1963	write(nulout,'(a,a,a,a,i0,a)') str, ' (', name, '=', val,')'
1964	end subroutine print_enum
1965
1966	end module radiation_config

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

Download in other formats:

Original Format