source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_config.F90 @ 4629

Last change on this file since 4629 was 4489, checked in by lguez, 19 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
36module 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
619contains
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
1966end module radiation_config
Note: See TracBrowser for help on using the repository browser.