! radiation_config.F90 - Derived type to configure the radiation scheme ! ! (C) Copyright 2014- ECMWF. ! ! This software is licensed under the terms of the Apache Licence Version 2.0 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. ! ! In applying this licence, ECMWF does not waive the privileges and immunities ! granted to it by virtue of its status as an intergovernmental organisation ! nor does it submit to any jurisdiction. ! ! Author: Robin Hogan ! Email: r.j.hogan@ecmwf.int ! ! Modifications ! 2017-07-22 R. Hogan Added Yi et al. ice optics model ! 2017-10-23 R. Hogan Renamed single-character variables ! 2018-03-15 R. Hogan Added logicals controlling surface spectral treatment ! 2018-08-29 R. Hogan Added monochromatic single-scattering albedo / asymmetry factor ! 2018-09-03 R. Hogan Added min_cloud_effective_size ! 2018-09-04 R. Hogan Added encroachment_scaling ! 2018-09-13 R. Hogan Added IEncroachmentFractal ! 2019-01-02 R. Hogan Added Cloudless solvers ! 2019-01-14 R. Hogan Added out_of_bounds_[1,2,3]d for checker routines ! 2019-01-18 R. Hogan Added albedo weighting ! 2019-02-03 R. Hogan Added ability to fix out-of-physical-bounds inputs ! 2019-02-10 R. Hogan Renamed "encroachment" to "entrapment" ! ! Note: The aim is for ecRad in the IFS to be as similar as possible ! to the offline version, so if you make any changes to this or any ! files in this directory, please inform Robin Hogan. ! module radiation_config use parkind1, only : jprb use radiation_cloud_optics_data, only : cloud_optics_type use radiation_aerosol_optics_data, only : aerosol_optics_type use radiation_pdf_sampler, only : pdf_sampler_type use radiation_cloud_cover, only : OverlapName, & & IOverlapMaximumRandom, IOverlapExponentialRandom, IOverlapExponential implicit none public ! Configuration codes: use C-style enumerators to avoid having to ! remember the numbers ! Solvers: can be specified for longwave and shortwave ! independently, except for "Homogeneous", which must be the same ! for both enum, bind(c) enumerator ISolverCloudless, ISolverHomogeneous, ISolverMcICA, & & ISolverSpartacus, ISolverTripleclouds end enum character(len=*), parameter :: SolverName(0:4) = (/ 'Cloudless ', & & 'Homogeneous ', & & 'McICA ', & & 'SPARTACUS ', & & 'Tripleclouds' /) ! SPARTACUS shortwave solver can treat the reflection of radiation ! back up into different regions in various ways enum, bind(c) enumerator & & IEntrapmentZero, & ! No entrapment, as Tripleclouds & IEntrapmentEdgeOnly, & ! Only radiation passed through cloud edge is horizontally homogenized & IEntrapmentExplicit, & ! Estimate horiz migration dist, account for fractal clouds & IEntrapmentExplicitNonFractal, & ! As above but ignore fractal nature of clouds & IEntrapmentMaximum ! Complete horizontal homogenization within regions (old SPARTACUS assumption) end enum ! Names available in the radiation namelist for variable ! sw_entrapment_name character(len=*), parameter :: EntrapmentName(0:4) = [ 'Zero ', & & 'Edge-only ', & & 'Explicit ', & & 'Non-fractal', & & 'Maximum ' ] ! For backwards compatibility, the radiation namelist also supports ! the equivalent variable sw_encroachment_name with the following ! names character(len=*), parameter :: EncroachmentName(0:4) = [ 'Zero ', & & 'Minimum ', & & 'Fractal ', & & 'Computed', & & 'Maximum ' ] ! Two-stream models ! This is not configurable at run-time ! Gas models enum, bind(c) enumerator IGasModelMonochromatic, IGasModelIFSRRTMG end enum character(len=*), parameter :: GasModelName(0:1) = (/ 'Monochromatic', & & 'RRTMG-IFS ' /) ! Hydrometeor scattering models enum, bind(c) enumerator ILiquidModelMonochromatic, & & ILiquidModelSOCRATES, ILiquidModelSlingo end enum character(len=*), parameter :: LiquidModelName(0:2) = (/ 'Monochromatic', & & 'SOCRATES ', & & 'Slingo ' /) enum, bind(c) enumerator IIceModelMonochromatic, IIceModelFu, & & IIceModelBaran, IIceModelBaran2016, IIceModelBaran2017, & & IIceModelYi end enum character(len=*), parameter :: IceModelName(0:5) = (/ 'Monochromatic', & & 'Fu-IFS ', & & 'Baran ', & & 'Baran2016 ', & & 'Baran2017 ', & & 'Yi ' /) ! Cloud PDF distribution shapes enum, bind(c) enumerator IPdfShapeLognormal, IPdfShapeGamma end enum character(len=*), parameter :: PdfShapeName(0:1) = (/ 'Lognormal', & & 'Gamma ' /) ! Maximum number of different aerosol types that can be provided integer, parameter :: NMaxAerosolTypes = 256 ! Maximum number of shortwave albedo and longwave emissivity ! intervals integer, parameter :: NMaxAlbedoIntervals = 256 ! Length of string buffer for printing config information integer, parameter :: NPrintStringLen = 60 !--------------------------------------------------------------------- ! Derived type containing all the configuration information needed ! to run the radiation scheme. The intention is that this is fixed ! for a given model run. The parameters are to list first those ! quantities that can be set directly by the user, for example using a ! namelist, and second those quantities that are computed afterwards ! from the user-supplied numbers, especially the details of the gas ! optics model. type config_type ! USER-CONFIGURABLE PARAMETERS ! Override default solar spectrum logical :: use_spectral_solar_scaling = .false. ! Directory in which gas, cloud and aerosol data files are to be ! found character(len=511) :: directory_name = '.' ! Cloud is deemed to be present in a layer if cloud fraction ! exceeds this value real(jprb) :: cloud_fraction_threshold = 1.0e-6_jprb ! ...and total cloud water mixing ratio exceeds this value real(jprb) :: cloud_mixing_ratio_threshold = 1.0e-9_jprb ! Overlap scheme integer :: i_overlap_scheme = IOverlapExponentialRandom ! Use the Shonk et al. (2010) "beta" overlap parameter, rather ! than the "alpha" overlap parameter of Hogan and Illingworth ! (2000)? logical :: use_beta_overlap = .false. ! Shape of sub-grid cloud water PDF integer :: i_cloud_pdf_shape = IPdfShapeGamma ! The ratio of the overlap decorrelation length for cloud ! inhomogeneities to the overlap decorrelation length for cloud ! boundaries. Observations suggest this has a value of 0.5 ! (e.g. from the decorrelation lengths of Hogan and Illingworth ! 2003 and Hogan and Illingworth 2000). real(jprb) :: cloud_inhom_decorr_scaling = 0.5_jprb ! Factor controlling how much of the cloud edge length interfaces ! directly between the clear-sky region (region a) and the ! optically thick cloudy region (region c). If Lxy is the length ! of the interfaces between regions x and y, and Lab and Lbc have ! been computed already, then ! Lac=clear_to_thick_fraction*min(Lab,Lbc). real(jprb) :: clear_to_thick_fraction = 0.0_jprb ! Factor allowing lateral transport when the sun is close to ! overhead; consider atand(overhead_sun_factor) to be the number ! of degrees that the sun angle is perturbed from zenith for the ! purposes of computing lateral transport. A value of up to 0.1 ! seems to be necessary to account for the fact that some forward ! scattered radiation is treated as unscattered by delta-Eddington ! scaling; therefore it ought to have the chance to escape. real(jprb) :: overhead_sun_factor = 0.0_jprb ! Minimum gas optical depth in a single layer at any wavelength, ! for stability real(jprb) :: min_gas_od_lw = 1.0e-15_jprb real(jprb) :: min_gas_od_sw = 0.0_jprb ! Maximum gas optical depth in a layer before that g-point will ! not be considered for 3D treatment: a limit is required to avoid ! expensive computation of matrix exponentials on matrices with ! large elements real(jprb) :: max_gas_od_3d = 8.0_jprb ! Maximum total optical depth of a cloudy region for stability: ! optical depth will be capped at this value in the SPARTACUS ! solvers real(jprb) :: max_cloud_od = 16.0_jprb ! How much longwave scattering is included? logical :: do_lw_cloud_scattering = .true. logical :: do_lw_aerosol_scattering = .true. ! Number of regions used to describe clouds and clear skies. A ! value of 2 means one clear and one cloudy region, so clouds are ! horizontally homogeneous, while a value of 3 means two cloudy ! regions with different optical depth, thereby representing ! inhomogeneity via the Shonk & Hogan (2008) "Tripleclouds" ! method. integer :: nregions = 3 ! Code specifying the solver to be used: use the enumerations ! defined above integer :: i_solver_sw = ISolverMcICA integer :: i_solver_lw = ISolverMcICA ! Do shortwave delta-Eddington scaling on the cloud-aerosol-gas ! mixture (as in the original IFS scheme), rather than the more ! correct approach of separately scaling the cloud and aerosol ! scattering properties before merging with gases. Note that ! .true. is not compatible with the SPARTACUS solver. logical :: do_sw_delta_scaling_with_gases = .false. ! Codes describing the gas and cloud scattering models to use, the ! latter of which is currently not used integer :: i_gas_model = IGasModelIFSRRTMG ! integer :: i_cloud_model ! Optics if i_gas_model==IGasModelMonochromatic. ! The wavelength to use for the Planck function in metres. If this ! is positive then the output longwave fluxes will be in units of ! W m-2 um-1. If this is zero or negative (the default) then ! sigma*T^4 will be used and the output longwave fluxes will be in ! W m-2. real(jprb) :: mono_lw_wavelength = -1.0_jprb ! Total zenith optical depth of the atmosphere in the longwave and ! shortwave, distributed vertically according to the pressure. ! Default is zero. real(jprb) :: mono_lw_total_od = 0.0_jprb real(jprb) :: mono_sw_total_od = 0.0_jprb ! Single-scattering albedo and asymmetry factor: values typical ! for liquid clouds with effective radius of 10 microns, at (SW) ! 0.55 micron wavelength and (LW) 10.7 microns wavelength real(jprb) :: mono_sw_single_scattering_albedo = 0.999999_jprb real(jprb) :: mono_sw_asymmetry_factor = 0.86_jprb real(jprb) :: mono_lw_single_scattering_albedo = 0.538_jprb real(jprb) :: mono_lw_asymmetry_factor = 0.925_jprb ! Codes describing particle scattering models integer :: i_liq_model = ILiquidModelSOCRATES integer :: i_ice_model = IIceModelBaran ! The mapping from albedo/emissivity intervals to SW/LW bands can ! either be done by finding the interval containing the central ! wavenumber of the band (nearest neighbour), or by a weighting ! according to the spectral overlap of each interval with each ! band logical :: do_nearest_spectral_sw_albedo = .true. logical :: do_nearest_spectral_lw_emiss = .true. ! User-defined monotonically increasing wavelength bounds (m) ! between input surface albedo/emissivity intervals. Implicitly ! the first interval starts at zero and the last ends at infinity. real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1) = -1.0_jprb real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) = -1.0_jprb ! The index to the surface albedo/emissivity intervals for each of ! the wavelength bounds specified in sw_albedo_wavelength_bound ! and lw_emiss_wavelength_bound integer :: i_sw_albedo_index(NMaxAlbedoIntervals) = 0 integer :: i_lw_emiss_index (NMaxAlbedoIntervals) = 0 ! Do we compute longwave and/or shortwave radiation? logical :: do_lw = .true. logical :: do_sw = .true. ! Do we compute clear-sky fluxes and/or solar direct fluxes? logical :: do_clear = .true. logical :: do_sw_direct = .true. ! Do we include 3D effects? logical :: do_3d_effects = .true. ! To what extent do we include "entrapment" effects in the ! SPARTACUS solver? This essentially means that in a situation ! like this ! ! 000111 ! 222222 ! ! Radiation downwelling from region 1 may be reflected back into ! region 0 due to some degree of homogenization of the radiation ! in region 2. Hogan and Shonk (2013) referred to this as ! "anomalous horizontal transport" for a 1D model, although for 3D ! calculations it is desirable to include at least some of it. The ! options are described by the IEntrapment* parameters above. integer :: i_3d_sw_entrapment = IEntrapmentExplicit ! In the longwave, the equivalent process it either "on" (like ! maximum entrapment) or "off" (like zero entrapment): logical :: do_3d_lw_multilayer_effects = .false. ! Do we account for the effective emissivity of the side of ! clouds? logical :: do_lw_side_emissivity = .true. ! The 3D transfer rate "X" is such that if transport out of a ! region was the only process occurring then by the base of a ! layer only exp(-X) of the original flux would remain in that ! region. The transfer rate computed geometrically can be very ! high for the clear-sky regions in layers with high cloud ! fraction. For stability reasons it is necessary to provide a ! maximum possible 3D transfer rate. real(jprb) :: max_3d_transfer_rate = 10.0_jprb ! It has also sometimes been found necessary to set a minimum ! cloud effective size for stability (metres) real(jprb) :: min_cloud_effective_size = 100.0_jprb ! Given a horizontal migration distance, there is still ! uncertainty about how much entrapment occurs associated with how ! one assumes cloud boundaries line up in adjacent layers. This ! factor can be varied between 0.0 (the boundaries line up to the ! greatest extent possible given the overlap parameter) and 1.0 ! (the boundaries line up to the minimum extent possible). In the ! Hogan et al. entrapment paper it is referred to as the overhang ! factor zeta, and a value of 0 matches the Monte Carlo ! calculations best. real(jprb) :: overhang_factor = 0.0_jprb ! By default, the Meador & Weaver (1980) expressions are used ! instead of the matrix exponential whenever 3D effects can be ! neglected (e.g. cloud-free layers or clouds with infinitely ! large effective cloud size), but setting the following to true ! uses the matrix exponential everywhere, enabling the two ! methods to be compared. Note that Meador & Weaver will still be ! used for very optically thick g points where the matrix ! exponential can produce incorrect results. logical :: use_expm_everywhere = .false. ! Aerosol descriptors: aerosol_type_mapping must be of length ! n_aerosol_types, and contains 0 if that type is to be ignored, ! positive numbers to map on to the indices of hydrophobic ! aerosols in the aerosol optics configuration file, and negative ! numbers to map on to (the negative of) the indices of ! hydrophilic aerosols in the configuration file. logical :: use_aerosols = .false. integer :: n_aerosol_types = 0 integer :: i_aerosol_type_map(NMaxAerosolTypes) ! Save the gas and cloud optical properties for each g point in ! "radiative_properties.nc"? logical :: do_save_radiative_properties = .false. ! Save the flux profiles in each band? logical :: do_save_spectral_flux = .false. ! Save the surface downwelling shortwave fluxes in each band? logical :: do_surface_sw_spectral_flux = .true. ! Compute the longwave derivatives needed to apply the approximate ! radiation updates of Hogan and Bozzo (2015) logical :: do_lw_derivatives = .false. ! Save the flux profiles in each g-point (overrides ! do_save_spectral_flux if TRUE)? logical :: do_save_gpoint_flux = .false. ! In the IFS environment, setting up RRTM has already been done ! so not needed to do it again logical :: do_setup_ifsrrtm = .true. ! In the IFS environment the old scheme has a bug in the Fu ! longwave ice optics whereby the single scattering albedo is one ! minus what it should be. Unfortunately fixing it makes ! forecasts worse. Setting the following to true reproduces the ! bug. logical :: do_fu_lw_ice_optics_bug = .false. ! Control verbosity: 0=none (no output to standard output; write ! to standard error only if an error occurs), 1=warning, 2=info, ! 3=progress, 4=detailed, 5=debug. Separate settings for the ! setup of the scheme and the execution of it. integer :: iverbosesetup = 2 integer :: iverbose = 1 ! Are we doing radiative transfer in complex surface canopies ! (streets/vegetation), in which case tailored downward fluxes are ! needed at the top of the canopy? logical :: do_canopy_fluxes_sw = .false. logical :: do_canopy_fluxes_lw = .false. ! If so, do we use the full spectrum as in the atmosphere, or just ! the reduced spectrum in which the shortwave albedo and longwave ! emissivity are provided? logical :: use_canopy_full_spectrum_sw = .false. logical :: use_canopy_full_spectrum_lw = .false. ! Do we treat gas radiative transfer in streets/vegetation? logical :: do_canopy_gases_sw = .false. logical :: do_canopy_gases_lw = .false. ! Optics file names for overriding the ones generated from the ! other options. If these remain empty then the generated names ! will be used (see the "consolidate_config" routine below). If ! the user assigns one of these and it starts with a '/' character ! then that will be used instead. If the user assigns one and it ! doesn't start with a '/' character then it will be prepended by ! the contents of directory_name. character(len=511) :: ice_optics_override_file_name = '' character(len=511) :: liq_optics_override_file_name = '' character(len=511) :: aerosol_optics_override_file_name = '' ! Optionally override the look-up table file for the cloud-water ! PDF used by the McICA solver character(len=511) :: cloud_pdf_override_file_name = '' ! Has "consolidate" been called? logical :: is_consolidated = .false. ! COMPUTED PARAMETERS ! Users of this library should not edit these parameters directly; ! they are set by the "consolidate" routine ! Wavenumber range for each band, in cm-1, which will be allocated ! to be of length n_bands_sw or n_bands_lw real(jprb), allocatable, dimension(:) :: wavenumber1_sw real(jprb), allocatable, dimension(:) :: wavenumber2_sw real(jprb), allocatable, dimension(:) :: wavenumber1_lw real(jprb), allocatable, dimension(:) :: wavenumber2_lw ! If the nearest surface albedo/emissivity interval is to be used ! for each SW/LW band then the following arrays will be allocated ! to the length of the number of bands and contain the index to ! the relevant interval integer, allocatable, dimension(:) :: i_albedo_from_band_sw integer, allocatable, dimension(:) :: i_emiss_from_band_lw ! ...alternatively, this matrix dimensioned ! (n_albedo_intervals,n_bands_sw) providing the weights needed for ! computing the albedo in each ecRad band from the albedo in each ! native albedo band - see radiation_single_level.F90 real(jprb), allocatable, dimension(:,:) :: sw_albedo_weights ! ...and similarly in the longwave, dimensioned ! (n_emiss_intervals,n_bands_lw) real(jprb), allocatable, dimension(:,:) :: lw_emiss_weights ! Arrays of length the number of g-points that convert from ! g-point to the band index integer, allocatable, dimension(:) :: i_band_from_g_lw integer, allocatable, dimension(:) :: i_band_from_g_sw ! We allow for the possibility for g-points to be ordered in terms ! of likely absorption (weakest to strongest) across the shortwave ! or longwave spectrum, in order that in SPARTACUS we select only ! the first n g-points that will not have too large an absorption, ! and therefore matrix exponentials that are both finite and not ! too expensive to compute. The following two arrays map the ! reordered g-points to the original ones. integer, allocatable, dimension(:) :: i_g_from_reordered_g_lw integer, allocatable, dimension(:) :: i_g_from_reordered_g_sw ! The following map the reordered g-points to the bands integer, allocatable, dimension(:) :: i_band_from_reordered_g_lw integer, allocatable, dimension(:) :: i_band_from_reordered_g_sw ! The following map the reordered g-points to the spectral ! information being saved: if do_save_gpoint_flux==TRUE then this ! will map on to the original g points, but if only ! do_save_spectral_flux==TRUE then this will map on to the bands integer, pointer, dimension(:) :: i_spec_from_reordered_g_lw integer, pointer, dimension(:) :: i_spec_from_reordered_g_sw ! Number of spectral intervals used for the canopy radiative ! transfer calculation; they are either equal to ! n_albedo_intervals/n_emiss_intervals or n_g_sw/n_g_lw integer :: n_canopy_bands_sw = 1 integer :: n_canopy_bands_lw = 1 ! Data structure containing cloud scattering data type(cloud_optics_type) :: cloud_optics ! Data structure containing aerosol scattering data type(aerosol_optics_type) :: aerosol_optics ! Object for sampling from a gamma or lognormal distribution type(pdf_sampler_type) :: pdf_sampler ! Optics file names character(len=511) :: ice_optics_file_name, & & liq_optics_file_name, & & aerosol_optics_file_name ! McICA PDF look-up table file name character(len=511) :: cloud_pdf_file_name ! Number of gpoints and bands in the shortwave and longwave - set ! to zero as will be set properly later integer :: n_g_sw = 0, n_g_lw = 0 integer :: n_bands_sw = 0, n_bands_lw = 0 ! Number of spectral points to save (equal either to the number of ! g points or the number of bands integer :: n_spec_sw = 0, n_spec_lw = 0 ! Dimensions to store variables that are only needed if longwave ! scattering is included. "n_g_lw_if_scattering" is equal to ! "n_g_lw" if aerosols are allowed to scatter in the longwave, ! and zero otherwise. "n_bands_lw_if_scattering" is equal to ! "n_bands_lw" if clouds are allowed to scatter in the longwave, ! and zero otherwise. integer :: n_g_lw_if_scattering = 0, n_bands_lw_if_scattering = 0 ! Treat clouds as horizontally homogeneous within the gribox logical :: is_homogeneous = .false. ! If the solvers are both "Cloudless" then we don't need to do any ! cloud processing logical :: do_clouds = .true. contains procedure :: read => read_config_from_namelist procedure :: consolidate => consolidate_config procedure :: set => set_config procedure :: print => print_config procedure :: get_sw_weights procedure :: define_sw_albedo_intervals procedure :: define_lw_emiss_intervals procedure :: consolidate_intervals end type config_type ! procedure, private :: print_logical, print_real, print_int contains !--------------------------------------------------------------------- ! This subroutine reads configuration data from a namelist file, and ! anything that is not in the namelists will be set to default ! values. If optional output argument "is_success" is present, then ! on error (e.g. missing file) it will be set to .false.; if this ! argument is missing then on error the program will be aborted. You ! may either specify the file_name or the unit of an open file to ! read, but not both. subroutine read_config_from_namelist(this, file_name, unit, is_success) use yomhook, only : lhook, dr_hook use radiation_io, only : nulout, nulerr, nulrad, radiation_abort class(config_type), intent(inout) :: this character(*), intent(in), optional :: file_name integer, intent(in), optional :: unit logical, intent(out), optional :: is_success integer :: iosopen, iosread ! Status after calling open and read ! The following variables are read from the namelists and map ! directly onto members of the config_type derived type ! To be read from the radiation_config namelist logical :: do_sw, do_lw, do_clear, do_sw_direct logical :: do_3d_effects, use_expm_everywhere, use_aerosols logical :: do_lw_side_emissivity logical :: do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug logical :: do_lw_aerosol_scattering, do_lw_cloud_scattering logical :: do_save_radiative_properties, do_save_spectral_flux logical :: do_save_gpoint_flux, do_surface_sw_spectral_flux logical :: use_beta_overlap, do_lw_derivatives logical :: do_sw_delta_scaling_with_gases logical :: do_canopy_fluxes_sw, do_canopy_fluxes_lw logical :: use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw logical :: do_canopy_gases_sw, do_canopy_gases_lw integer :: n_regions, iverbose, iverbosesetup, n_aerosol_types real(jprb):: mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od real(jprb):: mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo real(jprb):: mono_lw_asymmetry_factor, mono_sw_asymmetry_factor real(jprb):: cloud_inhom_decorr_scaling, cloud_fraction_threshold real(jprb):: clear_to_thick_fraction, max_gas_od_3d, max_cloud_od real(jprb):: cloud_mixing_ratio_threshold, overhead_sun_factor real(jprb):: max_3d_transfer_rate, min_cloud_effective_size real(jprb):: overhang_factor, encroachment_scaling character(511) :: directory_name, aerosol_optics_override_file_name character(511) :: liq_optics_override_file_name, ice_optics_override_file_name character(511) :: cloud_pdf_override_file_name character(63) :: liquid_model_name, ice_model_name, gas_model_name character(63) :: sw_solver_name, lw_solver_name, overlap_scheme_name character(63) :: sw_entrapment_name, sw_encroachment_name, cloud_pdf_shape_name integer :: i_aerosol_type_map(NMaxAerosolTypes) ! More than 256 is an error logical :: do_nearest_spectral_sw_albedo = .true. logical :: do_nearest_spectral_lw_emiss = .true. real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1) real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) integer :: i_sw_albedo_index(NMaxAlbedoIntervals) integer :: i_lw_emiss_index (NMaxAlbedoIntervals) integer :: iunit ! Unit number of namelist file namelist /radiation/ do_sw, do_lw, do_sw_direct, & & do_3d_effects, do_lw_side_emissivity, do_clear, & & do_save_radiative_properties, sw_entrapment_name, sw_encroachment_name, & & do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug, & & do_save_spectral_flux, do_save_gpoint_flux, & & do_surface_sw_spectral_flux, do_lw_derivatives, & & do_lw_aerosol_scattering, do_lw_cloud_scattering, & & n_regions, directory_name, gas_model_name, & & ice_optics_override_file_name, liq_optics_override_file_name, & & aerosol_optics_override_file_name, cloud_pdf_override_file_name, & & liquid_model_name, ice_model_name, max_3d_transfer_rate, & & min_cloud_effective_size, overhang_factor, encroachment_scaling, & & use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw, & & do_canopy_fluxes_sw, do_canopy_fluxes_lw, & & do_canopy_gases_sw, do_canopy_gases_lw, & & do_sw_delta_scaling_with_gases, overlap_scheme_name, & & sw_solver_name, lw_solver_name, use_beta_overlap, & & use_expm_everywhere, iverbose, iverbosesetup, & & cloud_inhom_decorr_scaling, cloud_fraction_threshold, & & clear_to_thick_fraction, max_gas_od_3d, max_cloud_od, & & cloud_mixing_ratio_threshold, overhead_sun_factor, & & n_aerosol_types, i_aerosol_type_map, use_aerosols, & & mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od, & & mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo, & & mono_lw_asymmetry_factor, mono_sw_asymmetry_factor, & & cloud_pdf_shape_name, & & do_nearest_spectral_sw_albedo, do_nearest_spectral_lw_emiss, & & sw_albedo_wavelength_bound, lw_emiss_wavelength_bound, & & i_sw_albedo_index, i_lw_emiss_index real(jprb) :: hook_handle if (lhook) call dr_hook('radiation_config:read',0,hook_handle) ! Copy default values from the original structure do_sw = this%do_sw do_lw = this%do_lw do_sw_direct = this%do_sw_direct do_3d_effects = this%do_3d_effects do_3d_lw_multilayer_effects = this%do_3d_lw_multilayer_effects do_lw_side_emissivity = this%do_lw_side_emissivity do_clear = this%do_clear do_lw_aerosol_scattering = this%do_lw_aerosol_scattering do_lw_cloud_scattering = this%do_lw_cloud_scattering do_sw_delta_scaling_with_gases = this%do_sw_delta_scaling_with_gases do_fu_lw_ice_optics_bug = this%do_fu_lw_ice_optics_bug do_canopy_fluxes_sw = this%do_canopy_fluxes_sw do_canopy_fluxes_lw = this%do_canopy_fluxes_lw use_canopy_full_spectrum_sw = this%use_canopy_full_spectrum_sw use_canopy_full_spectrum_lw = this%use_canopy_full_spectrum_lw do_canopy_gases_sw = this%do_canopy_gases_sw do_canopy_gases_lw = this%do_canopy_gases_lw n_regions = this%nregions directory_name = this%directory_name cloud_pdf_override_file_name = this%cloud_pdf_override_file_name liq_optics_override_file_name = this%liq_optics_override_file_name ice_optics_override_file_name = this%ice_optics_override_file_name aerosol_optics_override_file_name = this%aerosol_optics_override_file_name use_expm_everywhere = this%use_expm_everywhere use_aerosols = this%use_aerosols do_save_radiative_properties = this%do_save_radiative_properties do_save_spectral_flux = this%do_save_spectral_flux do_save_gpoint_flux = this%do_save_gpoint_flux do_lw_derivatives = this%do_lw_derivatives do_surface_sw_spectral_flux = this%do_surface_sw_spectral_flux iverbose = this%iverbose iverbosesetup = this%iverbosesetup cloud_fraction_threshold = this%cloud_fraction_threshold cloud_mixing_ratio_threshold = this%cloud_mixing_ratio_threshold use_beta_overlap = this%use_beta_overlap cloud_inhom_decorr_scaling = this%cloud_inhom_decorr_scaling clear_to_thick_fraction = this%clear_to_thick_fraction overhead_sun_factor = this%overhead_sun_factor max_gas_od_3d = this%max_gas_od_3d max_cloud_od = this%max_cloud_od max_3d_transfer_rate = this%max_3d_transfer_rate min_cloud_effective_size = this%min_cloud_effective_size overhang_factor = this%overhang_factor encroachment_scaling = -1.0_jprb gas_model_name = '' !DefaultGasModelName liquid_model_name = '' !DefaultLiquidModelName ice_model_name = '' !DefaultIceModelName sw_solver_name = '' !DefaultSwSolverName lw_solver_name = '' !DefaultLwSolverName sw_entrapment_name = '' sw_encroachment_name = '' overlap_scheme_name = '' cloud_pdf_shape_name = '' n_aerosol_types = this%n_aerosol_types mono_lw_wavelength = this%mono_lw_wavelength mono_lw_total_od = this%mono_lw_total_od mono_sw_total_od = this%mono_sw_total_od mono_lw_single_scattering_albedo = this%mono_lw_single_scattering_albedo mono_sw_single_scattering_albedo = this%mono_sw_single_scattering_albedo mono_lw_asymmetry_factor = this%mono_lw_asymmetry_factor mono_sw_asymmetry_factor = this%mono_sw_asymmetry_factor i_aerosol_type_map = this%i_aerosol_type_map do_nearest_spectral_sw_albedo = this%do_nearest_spectral_sw_albedo do_nearest_spectral_lw_emiss = this%do_nearest_spectral_lw_emiss sw_albedo_wavelength_bound = this%sw_albedo_wavelength_bound lw_emiss_wavelength_bound = this%lw_emiss_wavelength_bound i_sw_albedo_index = this%i_sw_albedo_index i_lw_emiss_index = this%i_lw_emiss_index if (present(file_name) .and. present(unit)) then write(nulerr,'(a)') '*** Error: cannot specify both file_name and unit in call to config_type%read' call radiation_abort('Radiation configuration error') else if (.not. present(file_name) .and. .not. present(unit)) then write(nulerr,'(a)') '*** Error: neither file_name nor unit specified in call to config_type%read' call radiation_abort('Radiation configuration error') end if if (present(file_name)) then ! Open the namelist file iunit = nulrad open(unit=iunit, iostat=iosopen, file=trim(file_name)) else ! Assume that iunit represents and open file iosopen = 0 iunit = unit end if if (iosopen /= 0) then ! An error occurred opening the file if (present(is_success)) then is_success = .false. ! We now continue the subroutine so that the default values ! are placed in the config structure else write(nulerr,'(a,a,a)') '*** Error: namelist file "', & & trim(file_name), '" not found' call radiation_abort('Radiation configuration error') end if else read(unit=iunit, iostat=iosread, nml=radiation) if (iosread /= 0) then ! An error occurred reading the file if (present(is_success)) then is_success = .false. ! We now continue the subroutine so that the default values ! are placed in the config structure else if (present(file_name)) then write(nulerr,'(a,a,a)') '*** Error reading namelist "radiation" from file "', & & trim(file_name), '"' close(unit=iunit) call radiation_abort('Radiation configuration error') else write(nulerr,'(a,i0)') '*** Error reading namelist "radiation" from unit ', & & iunit call radiation_abort('Radiation configuration error') end if end if if (present(file_name)) then close(unit=iunit) end if end if ! Copy namelist data into configuration object ! Start with verbosity levels, which should be within limits if (iverbose < 0) then iverbose = 0 end if this%iverbose = iverbose if (iverbosesetup < 0) then iverbosesetup = 0 end if this%iverbosesetup = iverbosesetup this%do_lw = do_lw this%do_sw = do_sw this%do_clear = do_clear this%do_sw_direct = do_sw_direct this%do_3d_effects = do_3d_effects this%do_3d_lw_multilayer_effects = do_3d_lw_multilayer_effects this%do_lw_side_emissivity = do_lw_side_emissivity this%use_expm_everywhere = use_expm_everywhere this%use_aerosols = use_aerosols this%do_lw_cloud_scattering = do_lw_cloud_scattering this%do_lw_aerosol_scattering = do_lw_aerosol_scattering this%nregions = n_regions this%do_surface_sw_spectral_flux = do_surface_sw_spectral_flux this%do_sw_delta_scaling_with_gases = do_sw_delta_scaling_with_gases this%do_fu_lw_ice_optics_bug = do_fu_lw_ice_optics_bug this%do_canopy_fluxes_sw = do_canopy_fluxes_sw this%do_canopy_fluxes_lw = do_canopy_fluxes_lw this%use_canopy_full_spectrum_sw = use_canopy_full_spectrum_sw this%use_canopy_full_spectrum_lw = use_canopy_full_spectrum_lw this%do_canopy_gases_sw = do_canopy_gases_sw this%do_canopy_gases_lw = do_canopy_gases_lw this%mono_lw_wavelength = mono_lw_wavelength this%mono_lw_total_od = mono_lw_total_od this%mono_sw_total_od = mono_sw_total_od this%mono_lw_single_scattering_albedo = mono_lw_single_scattering_albedo this%mono_sw_single_scattering_albedo = mono_sw_single_scattering_albedo this%mono_lw_asymmetry_factor = mono_lw_asymmetry_factor this%mono_sw_asymmetry_factor = mono_sw_asymmetry_factor this%use_beta_overlap = use_beta_overlap this%cloud_inhom_decorr_scaling = cloud_inhom_decorr_scaling this%clear_to_thick_fraction = clear_to_thick_fraction this%overhead_sun_factor = overhead_sun_factor this%max_gas_od_3d = max_gas_od_3d this%max_cloud_od = max_cloud_od this%max_3d_transfer_rate = max_3d_transfer_rate this%min_cloud_effective_size = max(1.0e-6_jprb, min_cloud_effective_size) if (encroachment_scaling >= 0.0_jprb) then this%overhang_factor = encroachment_scaling if (iverbose >= 1) then write(nulout, '(a)') 'Warning: radiation namelist parameter "encroachment_scaling" is deprecated: use "overhang_factor"' end if else this%overhang_factor = overhang_factor end if this%directory_name = directory_name this%cloud_pdf_override_file_name = cloud_pdf_override_file_name this%liq_optics_override_file_name = liq_optics_override_file_name this%ice_optics_override_file_name = ice_optics_override_file_name this%aerosol_optics_override_file_name = aerosol_optics_override_file_name this%cloud_fraction_threshold = cloud_fraction_threshold this%cloud_mixing_ratio_threshold = cloud_mixing_ratio_threshold this%n_aerosol_types = n_aerosol_types this%do_save_radiative_properties = do_save_radiative_properties this%do_lw_derivatives = do_lw_derivatives this%do_save_spectral_flux = do_save_spectral_flux this%do_save_gpoint_flux = do_save_gpoint_flux this%do_nearest_spectral_sw_albedo = do_nearest_spectral_sw_albedo this%do_nearest_spectral_lw_emiss = do_nearest_spectral_lw_emiss this%sw_albedo_wavelength_bound = sw_albedo_wavelength_bound this%lw_emiss_wavelength_bound = lw_emiss_wavelength_bound this%i_sw_albedo_index = i_sw_albedo_index this%i_lw_emiss_index = i_lw_emiss_index if (do_save_gpoint_flux) then ! Saving the fluxes every g-point overrides saving as averaged ! in a band, but this%do_save_spectral_flux needs to be TRUE as ! it is tested inside the solver routines to decide whether to ! save anything this%do_save_spectral_flux = .true. end if ! Determine liquid optics model call get_enum_code(liquid_model_name, LiquidModelName, & & 'liquid_model_name', this%i_liq_model) ! Determine ice optics model call get_enum_code(ice_model_name, IceModelName, & & 'ice_model_name', this%i_ice_model) ! Determine gas optics model call get_enum_code(gas_model_name, GasModelName, & & 'gas_model_name', this%i_gas_model) ! Determine solvers call get_enum_code(sw_solver_name, SolverName, & & 'sw_solver_name', this%i_solver_sw) call get_enum_code(lw_solver_name, SolverName, & & 'lw_solver_name', this%i_solver_lw) if (len_trim(sw_encroachment_name) > 1) then call get_enum_code(sw_encroachment_name, EncroachmentName, & & 'sw_encroachment_name', this%i_3d_sw_entrapment) write(nulout, '(a)') 'Warning: radiation namelist string "sw_encroachment_name" is deprecated: use "sw_entrapment_name"' else call get_enum_code(sw_entrapment_name, EntrapmentName, & & 'sw_entrapment_name', this%i_3d_sw_entrapment) end if ! Determine overlap scheme call get_enum_code(overlap_scheme_name, OverlapName, & & 'overlap_scheme_name', this%i_overlap_scheme) ! Determine cloud PDF shape call get_enum_code(cloud_pdf_shape_name, PdfShapeName, & & 'cloud_pdf_shape_name', this%i_cloud_pdf_shape) this%i_aerosol_type_map = 0 if (this%use_aerosols) then this%i_aerosol_type_map(1:n_aerosol_types) & & = i_aerosol_type_map(1:n_aerosol_types) end if ! Will clouds be used at all? if ((this%do_sw .and. this%i_solver_sw /= ISolverCloudless) & & .or. (this%do_lw .and. this%i_solver_lw /= ISolverCloudless)) then this%do_clouds = .true. else this%do_clouds = .false. end if ! Normal subroutine exit if (present(is_success)) then is_success = .true. end if if (lhook) call dr_hook('radiation_config:read',1,hook_handle) end subroutine read_config_from_namelist !--------------------------------------------------------------------- ! This routine is called by radiation_interface:setup_radiation and ! it converts the user specified options into some more specific ! data such as data file names subroutine consolidate_config(this) use yomhook, only : lhook, dr_hook use radiation_io, only : nulout, nulerr, radiation_abort class(config_type), intent(inout) :: this real(jprb) :: hook_handle if (lhook) call dr_hook('radiation_config:consolidate',0,hook_handle) ! Check consistency of models if (this%do_canopy_fluxes_sw .and. .not. this%do_surface_sw_spectral_flux) then if (this%iverbosesetup >= 1) then write(nulout,'(a)') 'Warning: turning on do_surface_sw_spectral_flux as required by do_canopy_fluxes_sw' end if this%do_surface_sw_spectral_flux = .true. end if ! Will clouds be used at all? if ((this%do_sw .and. this%i_solver_sw /= ISolverCloudless) & & .or. (this%do_lw .and. this%i_solver_lw /= ISolverCloudless)) then this%do_clouds = .true. else this%do_clouds = .false. end if ! SPARTACUS only works with Exp-Ran overlap scheme if (( this%i_solver_sw == ISolverSPARTACUS & & .or. this%i_solver_lw == ISolverSPARTACUS & & .or. this%i_solver_sw == ISolverTripleclouds & & .or. this%i_solver_lw == ISolverTripleclouds) & & .and. this%i_overlap_scheme /= IOverlapExponentialRandom) then write(nulerr,'(a)') '*** Error: SPARTACUS/Tripleclouds solvers can only do Exponential-Random overlap' call radiation_abort('Radiation configuration error') end if ! Set aerosol optics file name if (len_trim(this%aerosol_optics_override_file_name) > 0) then if (this%aerosol_optics_override_file_name(1:1) == '/') then this%aerosol_optics_file_name = trim(this%aerosol_optics_override_file_name) else this%aerosol_optics_file_name = trim(this%directory_name) & & // '/' // trim(this%aerosol_optics_override_file_name) end if else ! In the IFS, the aerosol optics file should be specified in ! ifs/module/radiation_setup.F90, not here this%aerosol_optics_file_name & & = trim(this%directory_name) // "/aerosol_ifs_rrtm_45R2.nc" end if ! Set liquid optics file name if (len_trim(this%liq_optics_override_file_name) > 0) then if (this%liq_optics_override_file_name(1:1) == '/') then this%liq_optics_file_name = trim(this%liq_optics_override_file_name) else this%liq_optics_file_name = trim(this%directory_name) & & // '/' // trim(this%liq_optics_override_file_name) end if else if (this%i_liq_model == ILiquidModelSOCRATES) then this%liq_optics_file_name & & = trim(this%directory_name) // "/socrates_droplet_scattering_rrtm.nc" else if (this%i_liq_model == ILiquidModelSlingo) then this%liq_optics_file_name & & = trim(this%directory_name) // "/slingo_droplet_scattering_rrtm.nc" end if ! Set ice optics file name if (len_trim(this%ice_optics_override_file_name) > 0) then if (this%ice_optics_override_file_name(1:1) == '/') then this%ice_optics_file_name = trim(this%ice_optics_override_file_name) else this%ice_optics_file_name = trim(this%directory_name) & & // '/' // trim(this%ice_optics_override_file_name) end if else if (this%i_ice_model == IIceModelFu) then this%ice_optics_file_name & & = trim(this%directory_name) // "/fu_ice_scattering_rrtm.nc" else if (this%i_ice_model == IIceModelBaran) then this%ice_optics_file_name & & = trim(this%directory_name) // "/baran_ice_scattering_rrtm.nc" else if (this%i_ice_model == IIceModelBaran2016) then this%ice_optics_file_name & & = trim(this%directory_name) // "/baran2016_ice_scattering_rrtm.nc" else if (this%i_ice_model == IIceModelBaran2017) then this%ice_optics_file_name & & = trim(this%directory_name) // "/baran2017_ice_scattering_rrtm.nc" else if (this%i_ice_model == IIceModelYi) then this%ice_optics_file_name & & = trim(this%directory_name) // "/yi_ice_scattering_rrtm.nc" end if ! Set cloud-water PDF look-up table file name if (len_trim(this%cloud_pdf_override_file_name) > 0) then if (this%cloud_pdf_override_file_name(1:1) == '/') then this%cloud_pdf_file_name = trim(this%cloud_pdf_override_file_name) else this%cloud_pdf_file_name = trim(this%directory_name) & & // '/' // trim(this%cloud_pdf_override_file_name) end if elseif (this%i_cloud_pdf_shape == IPdfShapeLognormal) then this%cloud_pdf_file_name = trim(this%directory_name) // "/mcica_lognormal.nc" else this%cloud_pdf_file_name = trim(this%directory_name) // "/mcica_gamma.nc" end if ! Aerosol data if (this%n_aerosol_types < 0 & & .or. this%n_aerosol_types > NMaxAerosolTypes) then write(nulerr,'(a,i0)') '*** Error: number of aerosol types must be between 0 and ', & & NMaxAerosolTypes call radiation_abort('Radiation configuration error') end if if (this%use_aerosols .and. this%n_aerosol_types == 0) then if (this%iverbosesetup >= 2) then write(nulout, '(a)') 'Aerosols on but n_aerosol_types=0: optical properties to be computed outside ecRad' end if end if ! In the monochromatic case we need to override the liquid, ice ! and aerosol models to ensure compatibility if (this%i_gas_model == IGasModelMonochromatic) then this%i_liq_model = ILiquidModelMonochromatic this%i_ice_model = IIceModelMonochromatic this%use_aerosols = .false. end if ! McICA solver currently can't store full profiles of spectral fluxes if (this%i_solver_sw == ISolverMcICA) then this%do_save_spectral_flux = .false. end if if (this%i_solver_sw == ISolverSPARTACUS .and. this%do_sw_delta_scaling_with_gases) then write(nulerr,'(a)') '*** Error: SW delta-Eddington scaling with gases not possible with SPARTACUS solver' call radiation_abort('Radiation configuration error') end if if ((this%do_lw .and. this%do_sw) .and. & & ( ( this%i_solver_sw == ISolverHomogeneous & & .and. this%i_solver_lw /= ISolverHomogeneous) & & .or. ( this%i_solver_sw /= ISolverHomogeneous & & .and. this%i_solver_lw == ISolverHomogeneous) & & ) ) then write(nulerr,'(a)') '*** Error: if one solver is "Homogeneous" then the other must be' call radiation_abort('Radiation configuration error') end if ! Set is_homogeneous if the active solvers are homogeneous, since ! this affects how "in-cloud" water contents are computed if ( (this%do_sw .and. this%i_solver_sw == ISolverHomogeneous) & & .or. (this%do_lw .and. this%i_solver_lw == ISolverHomogeneous)) then this%is_homogeneous = .true. end if this%is_consolidated = .true. if (lhook) call dr_hook('radiation_config:consolidate',1,hook_handle) end subroutine consolidate_config !--------------------------------------------------------------------- ! This subroutine sets members of the configuration object via ! optional arguments, and any member not specified is left ! untouched. Therefore, this should be called after taking data from ! the namelist. subroutine set_config(config, directory_name, & & do_lw, do_sw, & & do_lw_aerosol_scattering, do_lw_cloud_scattering, & & do_sw_direct) class(config_type), intent(inout):: config character(len=*), intent(in), optional :: directory_name logical, intent(in), optional :: do_lw, do_sw logical, intent(in), optional :: do_lw_aerosol_scattering logical, intent(in), optional :: do_lw_cloud_scattering logical, intent(in), optional :: do_sw_direct if (present(do_lw)) then config%do_lw = do_lw end if if(present(do_sw)) then config%do_sw = do_sw end if if (present(do_sw_direct)) then config%do_sw_direct = do_sw_direct end if if (present(directory_name)) then config%directory_name = trim(directory_name) end if if (present(do_lw_aerosol_scattering)) then config%do_lw_aerosol_scattering = .true. end if if (present(do_lw_cloud_scattering)) then config%do_lw_cloud_scattering = .true. end if end subroutine set_config !--------------------------------------------------------------------- ! Print configuration information to standard output subroutine print_config(this, iverbose) use radiation_io, only : nulout class(config_type), intent(in) :: this integer, optional, intent(in) :: iverbose integer :: i_local_verbose if (present(iverbose)) then i_local_verbose = iverbose else i_local_verbose = this%iverbose end if if (i_local_verbose >= 2) then !--------------------------------------------------------------------- write(nulout, '(a)') 'General settings:' write(nulout, '(a,a,a)') ' Data files expected in "', & & trim(this%directory_name), '"' call print_logical(' Clear-sky calculations are', 'do_clear', this%do_clear) call print_logical(' Saving intermediate radiative properties', & & 'do_save_radiative_properties', this%do_save_radiative_properties) call print_logical(' Saving spectral flux profiles', & & 'do_save_spectral_flux', this%do_save_spectral_flux) call print_enum(' Gas model is', GasModelName, 'i_gas_model', & & this%i_gas_model) call print_logical(' Aerosols are', 'use_aerosols', this%use_aerosols) call print_logical(' Clouds are', 'do_clouds', this%do_clouds) !--------------------------------------------------------------------- write(nulout, '(a)') 'Surface settings:' if (this%do_sw) then call print_logical(' Saving surface shortwave spectral fluxes', & & 'do_surface_sw_spectral_flux', this%do_surface_sw_spectral_flux) call print_logical(' Saving surface shortwave fluxes in abledo bands', & & 'do_canopy_fluxes_sw', this%do_canopy_fluxes_sw) end if if (this%do_lw) then call print_logical(' Saving surface longwave fluxes in emissivity bands', & & 'do_canopy_fluxes_lw', this%do_canopy_fluxes_lw) call print_logical(' Longwave derivative calculation is', & & 'do_lw_derivatives',this%do_lw_derivatives) end if if (this%do_sw) then call print_logical(' Nearest-neighbour spectral albedo mapping', & & 'do_nearest_spectral_sw_albedo', this%do_nearest_spectral_sw_albedo) end if if (this%do_lw) then call print_logical(' Nearest-neighbour spectral emissivity mapping', & & 'do_nearest_spectral_lw_emiss', this%do_nearest_spectral_lw_emiss) end if !--------------------------------------------------------------------- if (this%do_clouds) then write(nulout, '(a)') 'Cloud settings:' call print_real(' Cloud fraction threshold', & & 'cloud_fraction_threshold', this%cloud_fraction_threshold) call print_real(' Cloud mixing-ratio threshold', & & 'cloud_mixing_ratio_threshold', this%cloud_mixing_ratio_threshold) call print_enum(' Liquid optics scheme is', LiquidModelName, & & 'i_liq_model',this%i_liq_model) call print_enum(' Ice optics scheme is', IceModelName, & & 'i_ice_model',this%i_ice_model) if (this%i_ice_model == IIceModelFu) then call print_logical(' Longwave ice optics bug in Fu scheme is', & & 'do_fu_lw_ice_optics_bug',this%do_fu_lw_ice_optics_bug) end if call print_enum(' Cloud overlap scheme is', OverlapName, & & 'i_overlap_scheme',this%i_overlap_scheme) call print_logical(' Use "beta" overlap parameter is', & & 'use_beta_overlap', this%use_beta_overlap) call print_enum(' Cloud PDF shape is', PdfShapeName, & & 'i_cloud_pdf_shape',this%i_cloud_pdf_shape) call print_real(' Cloud inhom decorrelation scaling', & & 'cloud_inhom_decorr_scaling', this%cloud_inhom_decorr_scaling) end if !--------------------------------------------------------------------- write(nulout, '(a)') 'Solver settings:' if (this%do_sw) then call print_enum(' Shortwave solver is', SolverName, & & 'i_solver_sw', this%i_solver_sw) if (this%i_gas_model == IGasModelMonochromatic) then call print_real(' Shortwave atmospheric optical depth', & & 'mono_sw_total_od', this%mono_sw_total_od) call print_real(' Shortwave particulate single-scattering albedo', & & 'mono_sw_single_scattering_albedo', & & this%mono_sw_single_scattering_albedo) call print_real(' Shortwave particulate asymmetry factor', & & 'mono_sw_asymmetry_factor', & & this%mono_sw_asymmetry_factor) end if call print_logical(' Shortwave delta scaling after merge with gases', & & 'do_sw_delta_scaling_with_gases', & & this%do_sw_delta_scaling_with_gases) else call print_logical(' Shortwave calculations are','do_sw',this%do_sw) end if if (this%do_lw) then call print_enum(' Longwave solver is', SolverName, 'i_solver_lw', & & this%i_solver_lw) if (this%i_gas_model == IGasModelMonochromatic) then if (this%mono_lw_wavelength > 0.0_jprb) then call print_real(' Longwave effective wavelength (m)', & & 'mono_lw_wavelength', this%mono_lw_wavelength) else write(nulout,'(a)') ' Longwave fluxes are broadband (mono_lw_wavelength<=0)' end if call print_real(' Longwave atmospheric optical depth', & & 'mono_lw_total_od', this%mono_lw_total_od) call print_real(' Longwave particulate single-scattering albedo', & & 'mono_lw_single_scattering_albedo', & & this%mono_lw_single_scattering_albedo) call print_real(' Longwave particulate asymmetry factor', & & 'mono_lw_asymmetry_factor', & & this%mono_lw_asymmetry_factor) end if call print_logical(' Longwave cloud scattering is', & & 'do_lw_cloud_scattering',this%do_lw_cloud_scattering) call print_logical(' Longwave aerosol scattering is', & & 'do_lw_aerosol_scattering',this%do_lw_aerosol_scattering) else call print_logical(' Longwave calculations are','do_lw', this%do_lw) end if if (this%i_solver_sw == ISolverSpartacus & & .or. this%i_solver_lw == ISolverSpartacus) then write(nulout, '(a)') ' SPARTACUS options:' call print_integer(' Number of regions', 'n_regions', this%nregions) call print_real(' Max cloud optical depth per layer', & & 'max_cloud_od', this%max_cloud_od) call print_enum(' Shortwave entrapment is', EntrapmentName, & & 'i_3d_sw_entrapment', this%i_3d_sw_entrapment) call print_logical(' Multilayer longwave horizontal transport is', & 'do_3d_lw_multilayer_effects', this%do_3d_lw_multilayer_effects) call print_logical(' Use matrix exponential everywhere is', & & 'use_expm_everywhere', this%use_expm_everywhere) call print_logical(' 3D effects are', 'do_3d_effects', & & this%do_3d_effects) if (this%do_3d_effects) then call print_logical(' Longwave side emissivity parameterization is', & & 'do_lw_side_emissivity', this%do_lw_side_emissivity) call print_real(' Clear-to-thick edge fraction is', & & 'clear_to_thick_fraction', this%clear_to_thick_fraction) call print_real(' Overhead sun factor is', & & 'overhead_sun_factor', this%overhead_sun_factor) call print_real(' Max gas optical depth for 3D effects', & & 'max_gas_od_3d', this%max_gas_od_3d) call print_real(' Max 3D transfer rate', & & 'max_3d_transfer_rate', this%max_3d_transfer_rate) call print_real(' Min cloud effective size (m)', & & 'min_cloud_effective_size', this%min_cloud_effective_size) call print_real(' Overhang factor', & & 'overhang_factor', this%overhang_factor) end if end if end if end subroutine print_config !--------------------------------------------------------------------- ! In order to estimate UV and photosynthetically active radiation, ! we need weighted sum of fluxes considering wavelength range ! required. This routine returns information for how to correctly ! weight output spectral fluxes for a range of input wavelengths. ! Note that this is approximate; internally it may be assumed that ! the energy is uniformly distributed in wavenumber space, for ! example. If the character string "weighting_name" is present, and ! iverbose>=2, then information on the weighting will be provided on ! nulout. subroutine get_sw_weights(this, wavelength1, wavelength2, & & nweights, iband, weight, weighting_name) use parkind1, only : jprb use radiation_io, only : nulout, nulerr, radiation_abort class(config_type), intent(in) :: this ! Range of wavelengths to get weights for (m) real(jprb), intent(in) :: wavelength1, wavelength2 ! Output number of weights needed integer, intent(out) :: nweights ! Only write to the first nweights of these arrays: they contain ! the indices to the non-zero bands, and the weight in each of ! those bands integer, intent(out) :: iband(:) real(jprb), intent(out) :: weight(:) character(len=*), optional, intent(in) :: weighting_name ! Internally we deal with wavenumber real(jprb) :: wavenumber1, wavenumber2 ! cm-1 integer :: jband ! Loop index for spectral band if (this%n_bands_sw <= 0) then write(nulerr,'(a)') '*** Error: get_sw_weights called before number of shortwave bands set' call radiation_abort() end if ! Convert wavelength range (m) to wavenumber (cm-1) wavenumber1 = 0.01_jprb / wavelength2 wavenumber2 = 0.01_jprb / wavelength1 nweights = 0 do jband = 1,this%n_bands_sw if (wavenumber1 < this%wavenumber2_sw(jband) & & .and. wavenumber2 > this%wavenumber1_sw(jband)) then nweights = nweights+1 iband(nweights) = jband weight(nweights) = (min(wavenumber2,this%wavenumber2_sw(jband)) & & - max(wavenumber1,this%wavenumber1_sw(jband))) & & / (this%wavenumber2_sw(jband) - this%wavenumber1_sw(jband)) end if end do if (nweights == 0) then write(nulerr,'(a,e8.4,a,e8.4,a)') '*** Error: wavelength range ', & & wavelength1, ' to ', wavelength2, ' m is outside shortwave band' call radiation_abort() else if (this%iverbosesetup >= 2 .and. present(weighting_name)) then write(nulout,'(a,a,a,f6.0,a,f6.0,a)') 'Spectral weights for ', & & weighting_name, ' (', wavenumber1, ' to ', & & wavenumber2, ' cm-1):' do jband = 1, nweights write(nulout, '(a,i0,a,f6.0,a,f6.0,a,f8.4)') ' Shortwave band ', & & iband(jband), ' (', this%wavenumber1_sw(iband(jband)), ' to ', & & this%wavenumber2_sw(iband(jband)), ' cm-1): ', weight(jband) end do end if end subroutine get_sw_weights !--------------------------------------------------------------------- ! The input shortwave surface albedo coming in is likely to be in ! different spectral intervals to the gas model in the radiation ! scheme. We assume that the input albedo is defined within ! "ninterval" spectral intervals covering the wavelength range 0 to ! infinity, but allow for the possibility that two intervals may be ! indexed back to the same albedo band. subroutine define_sw_albedo_intervals(this, ninterval, wavelength_bound, & & i_intervals, do_nearest) use radiation_io, only : nulerr, radiation_abort class(config_type), intent(inout) :: this ! Number of spectral intervals in which albedo is defined integer, intent(in) :: ninterval ! Monotonically increasing wavelength bounds between intervals, ! not including the outer bounds (which are assumed to be zero and ! infinity) real(jprb), intent(in) :: wavelength_bound(ninterval-1) ! The albedo indices corresponding to each interval integer, intent(in) :: i_intervals(ninterval) logical, optional, intent(in) :: do_nearest if (ninterval > NMaxAlbedoIntervals) then write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, & & ' albedo intervals exceeds maximum of ', NMaxAlbedoIntervals call radiation_abort(); end if if (present(do_nearest)) then this%do_nearest_spectral_sw_albedo = do_nearest else this%do_nearest_spectral_sw_albedo = .false. end if this%sw_albedo_wavelength_bound(1:ninterval-1) = wavelength_bound(1:ninterval-1) this%sw_albedo_wavelength_bound(ninterval:) = -1.0_jprb this%i_sw_albedo_index(1:ninterval) = i_intervals(1:ninterval) this%i_sw_albedo_index(ninterval+1:) = 0 if (this%is_consolidated) then call this%consolidate_intervals(.true., & & this%do_nearest_spectral_sw_albedo, & & this%sw_albedo_wavelength_bound, this%i_sw_albedo_index, & & this%wavenumber1_sw, this%wavenumber2_sw, & & this%i_albedo_from_band_sw, this%sw_albedo_weights) end if end subroutine define_sw_albedo_intervals !--------------------------------------------------------------------- ! As define_sw_albedo_intervals but for longwave emissivity subroutine define_lw_emiss_intervals(this, ninterval, wavelength_bound, & & i_intervals, do_nearest) use radiation_io, only : nulerr, radiation_abort class(config_type), intent(inout) :: this ! Number of spectral intervals in which emissivity is defined integer, intent(in) :: ninterval ! Monotonically increasing wavelength bounds between intervals, ! not including the outer bounds (which are assumed to be zero and ! infinity) real(jprb), intent(in) :: wavelength_bound(ninterval-1) ! The emissivity indices corresponding to each interval integer, intent(in) :: i_intervals(ninterval) logical, optional, intent(in) :: do_nearest if (ninterval > NMaxAlbedoIntervals) then write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, & & ' emissivity intervals exceeds maximum of ', NMaxAlbedoIntervals call radiation_abort(); end if if (present(do_nearest)) then this%do_nearest_spectral_lw_emiss = do_nearest else this%do_nearest_spectral_lw_emiss = .false. end if this%lw_emiss_wavelength_bound(1:ninterval-1) = wavelength_bound(1:ninterval-1) this%lw_emiss_wavelength_bound(ninterval:) = -1.0_jprb this%i_lw_emiss_index(1:ninterval) = i_intervals(1:ninterval) this%i_lw_emiss_index(ninterval+1:) = 0 if (this%is_consolidated) then call this%consolidate_intervals(.false., & & this%do_nearest_spectral_lw_emiss, & & this%lw_emiss_wavelength_bound, this%i_lw_emiss_index, & & this%wavenumber1_lw, this%wavenumber2_lw, & & this%i_emiss_from_band_lw, this%lw_emiss_weights) end if end subroutine define_lw_emiss_intervals !--------------------------------------------------------------------- ! This routine consolidates either the input shortwave albedo ! intervals with the shortwave bands, or the input longwave ! emissivity intervals with the longwave bands, depending on the ! arguments provided. subroutine consolidate_intervals(this, is_sw, do_nearest, & & wavelength_bound, i_intervals, wavenumber1, wavenumber2, & & i_mapping, weights) use radiation_io, only : nulout, nulerr, radiation_abort class(config_type), intent(inout) :: this ! Is this the shortwave? Otherwise longwave logical, intent(in) :: is_sw ! Do we find the nearest albedo interval to the centre of each ! band, or properly weight the contributions? This can be modified ! if there is only one albedo intervals. logical, intent(inout) :: do_nearest ! Monotonically increasing wavelength bounds between intervals, ! not including the outer bounds (which are assumed to be zero and ! infinity) real(jprb), intent(in) :: wavelength_bound(NMaxAlbedoIntervals-1) ! The albedo band indices corresponding to each interval integer, intent(in) :: i_intervals(NMaxAlbedoIntervals) ! Start and end wavenumber bounds for the ecRad bands (cm-1) real(jprb), intent(in) :: wavenumber1(:), wavenumber2(:) ! if do_nearest is TRUE then the result is expressed in i_mapping, ! which will be allocated to have the same length as wavenumber1, ! and contain the index of the albedo interval corresponding to ! that band integer, allocatable, intent(inout) :: i_mapping(:) ! ...otherwise the result is expressed in "weights", of ! size(n_intervals, n_bands) containing how much of each interval ! contributes to each band. real(jprb), allocatable, intent(inout) :: weights(:,:) ! Number and loop index of ecRad bands integer :: nband, jband ! Number and index of albedo/emissivity intervals integer :: ninterval, iinterval ! Sometimes an albedo or emissivity value will be used in more ! than one interval, so nvalue indicates how many values will ! actually be provided integer :: nvalue ! Wavenumber bounds of the albedo/emissivity interval real(jprb) :: wavenumber1_albedo, wavenumber2_albedo ! Reciprocal of the wavenumber range of the ecRad band real(jprb) :: recip_dwavenumber ! cm ! Midpoint/bound of wavenumber band real(jprb) :: wavenumber_mid, wavenumber_bound ! cm-1 nband = size(wavenumber1) ! Count the number of albedo/emissivity intervals ninterval = 0 do iinterval = 1,NMaxAlbedoIntervals if (i_intervals(iinterval) > 0) then ninterval = iinterval else exit end if end do if (ninterval < 2) then ! Zero or one albedo/emissivity intervals found, so we index all ! bands to one interval if (allocated(i_mapping)) then deallocate(i_mapping) end if allocate(i_mapping(nband)) i_mapping(:) = 1 do_nearest = .true. ninterval = 1 nvalue = 1 else ! Check wavelength is monotonically increasing do jband = 2,ninterval-1 if (wavelength_bound(jband) <= wavelength_bound(jband-1)) then if (is_sw) then write(nulerr, '(a,a)') '*** Error: wavelength bounds for shortwave albedo intervals ', & & 'must be monotonically increasing' else write(nulerr, '(a,a)') '*** Error: wavelength bounds for longwave emissivity intervals ', & & 'must be monotonically increasing' end if call radiation_abort() end if end do ! What is the maximum index, indicating the number of ! albedo/emissivity values to expect? nvalue = maxval(i_intervals(1:ninterval)) if (do_nearest) then ! Simpler nearest-neighbour mapping from band to ! albedo/emissivity interval if (allocated(i_mapping)) then deallocate(i_mapping) end if allocate(i_mapping(nband)) ! Loop over bands do jband = 1,nband ! Compute mid-point of band in wavenumber space (cm-1) wavenumber_mid = 0.5_jprb * (wavenumber1(jband) & & + wavenumber2(jband)) iinterval = 1 ! Convert wavelength (m) into wavenumber (cm-1) at the lower ! bound of the albedo interval wavenumber_bound = 0.01_jprb / wavelength_bound(iinterval) ! Find the albedo interval that has the largest overlap with ! the band; this approach assumes that the albedo intervals ! are larger than the spectral bands do while (wavenumber_bound >= wavenumber_mid & & .and. iinterval < ninterval) iinterval = iinterval + 1 if (iinterval < ninterval) then wavenumber_bound = 0.01_jprb / wavelength_bound(iinterval) else ! For the last interval there is no lower bound wavenumber_bound = 0.0_jprb end if end do ! Save the index of the band corresponding to the albedo ! interval and move onto the next band i_mapping(jband) = i_intervals(iinterval) end do else ! More accurate weighting if (allocated(weights)) then deallocate(weights) end if allocate(weights(nvalue,nband)) weights(:,:) = 0.0_jprb ! Loop over bands do jband = 1,nband recip_dwavenumber = 1.0_jprb / (wavenumber2(jband) & & - wavenumber1(jband)) ! Find the first overlapping albedo band iinterval = 1 ! Convert wavelength (m) into wavenumber (cm-1) at the lower ! bound (in wavenumber space) of the albedo/emissivty interval wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval) do while (wavenumber1_albedo >= wavenumber2(jband) & & .and. iinterval < ninterval) iinterval = iinterval + 1 wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval) end do wavenumber2_albedo = wavenumber2(jband) ! Add all overlapping bands do while (wavenumber2_albedo > wavenumber1(jband) & & .and. iinterval <= ninterval) weights(i_intervals(iinterval),jband) & & = weights(i_intervals(iinterval),jband) & & + recip_dwavenumber & & * (min(wavenumber2_albedo,wavenumber2(jband)) & & - max(wavenumber1_albedo,wavenumber1(jband))) wavenumber2_albedo = wavenumber1_albedo iinterval = iinterval + 1 if (iinterval < ninterval) then wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval) else wavenumber1_albedo = 0.0_jprb end if end do end do end if end if ! Define how many bands to use for reporting surface downwelling ! fluxes for canopy radiation scheme if (is_sw) then if (this%use_canopy_full_spectrum_sw) then this%n_canopy_bands_sw = this%n_g_sw else this%n_canopy_bands_sw = nvalue end if else if (this%use_canopy_full_spectrum_lw) then this%n_canopy_bands_lw = this%n_g_lw else this%n_canopy_bands_lw = nvalue end if end if if (this%iverbosesetup >= 2) then if (.not. do_nearest) then if (is_sw) then write(nulout, '(a,i0,a,i0,a)') 'Weighting of ', nvalue, ' albedo values in ', & & nband, ' shortwave bands (wavenumber ranges in cm-1):' else write(nulout, '(a,i0,a,i0,a)') 'Weighting of ', nvalue, ' emissivity values in ', & & nband, ' longwave bands (wavenumber ranges in cm-1):' end if do jband = 1,nband write(nulout,'(i6,a,i6,a)',advance='no') nint(wavenumber1(jband)), ' to', & & nint(wavenumber2(jband)), ':' do iinterval = 1,nvalue write(nulout,'(f5.2)',advance='no') weights(iinterval,jband) end do write(nulout, '()') end do else if (ninterval <= 1) then if (is_sw) then write(nulout, '(a)') 'All shortwave bands will use the same albedo' else write(nulout, '(a)') 'All longwave bands will use the same emissivty' end if else if (is_sw) then write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', nband, & & ' shortwave bands to albedo intervals:' else write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', nband, & & ' longwave bands to emissivity intervals:' end if do jband = 1,nband write(nulout,'(a,i0)',advance='no') ' ', i_mapping(jband) end do write(nulout, '()') end if end if end subroutine consolidate_intervals !--------------------------------------------------------------------- ! Return the 0-based index for str in enum_str, or abort if it is ! not found subroutine get_enum_code(str, enum_str, var_name, icode) use radiation_io, only : nulerr, radiation_abort character(len=*), intent(in) :: str character(len=*), intent(in) :: enum_str(0:) character(len=*), intent(in) :: var_name integer, intent(out) :: icode integer :: jc logical :: is_not_found ! If string is empty then we don't modify icode but assume it has ! a sensible default value if (len_trim(str) > 1) then is_not_found = .true. do jc = 0,size(enum_str)-1 if (trim(str) == trim(enum_str(jc))) then icode = jc is_not_found = .false. exit end if end do if (is_not_found) then write(nulerr,'(a,a,a,a,a)',advance='no') '*** Error: ', trim(var_name), & & ' must be one of: "', enum_str(0), '"' do jc = 1,size(enum_str)-1 write(nulerr,'(a,a,a)',advance='no') ', "', trim(enum_str(jc)), '"' end do write(nulerr,'(a)') '' call radiation_abort('Radiation configuration error') end if end if end subroutine get_enum_code !--------------------------------------------------------------------- ! Print one line of information: logical subroutine print_logical(message_str, name, val) use radiation_io, only : nulout character(len=*), intent(in) :: message_str character(len=*), intent(in) :: name logical, intent(in) :: val character(4) :: on_or_off character(NPrintStringLen) :: str if (val) then on_or_off = ' ON ' else on_or_off = ' OFF' end if write(str, '(a,a4)') message_str, on_or_off write(nulout,'(a,a,a,a,l1,a)') str, ' (', name, '=', val,')' end subroutine print_logical !--------------------------------------------------------------------- ! Print one line of information: integer subroutine print_integer(message_str, name, val) use radiation_io, only : nulout character(len=*), intent(in) :: message_str character(len=*), intent(in) :: name integer, intent(in) :: val character(NPrintStringLen) :: str write(str, '(a,a,i0)') message_str, ' = ', val write(nulout,'(a,a,a,a)') str, ' (', name, ')' end subroutine print_integer !--------------------------------------------------------------------- ! Print one line of information: real subroutine print_real(message_str, name, val) use parkind1, only : jprb use radiation_io, only : nulout character(len=*), intent(in) :: message_str character(len=*), intent(in) :: name real(jprb), intent(in) :: val character(NPrintStringLen) :: str write(str, '(a,a,g8.3)') message_str, ' = ', val write(nulout,'(a,a,a,a)') str, ' (', name, ')' end subroutine print_real !--------------------------------------------------------------------- ! Print one line of information: enum subroutine print_enum(message_str, enum_str, name, val) use radiation_io, only : nulout character(len=*), intent(in) :: message_str character(len=*), intent(in) :: enum_str(0:) character(len=*), intent(in) :: name integer, intent(in) :: val character(NPrintStringLen) :: str write(str, '(a,a,a,a)') message_str, ' "', trim(enum_str(val)), '"' write(nulout,'(a,a,a,a,i0,a)') str, ' (', name, '=', val,')' end subroutine print_enum !--------------------------------------------------------------------- ! Return .true. if 1D allocatable array "var" is out of physical ! range specified by boundmin and boundmax, and issue a warning. ! "do_fix" determines whether erroneous values are fixed to lie ! within the physical range. To check only a subset of the array, ! specify i1 and i2 for the range. function out_of_bounds_1d(var, var_name, boundmin, boundmax, do_fix, i1, i2) result (is_bad) use radiation_io, only : nulout real(jprb), allocatable, intent(inout) :: var(:) character(len=*), intent(in) :: var_name real(jprb), intent(in) :: boundmin, boundmax logical, intent(in) :: do_fix integer, optional, intent(in) :: i1, i2 logical :: is_bad real(jprb) :: varmin, varmax is_bad = .false. if (allocated(var)) then if (present(i1) .and. present(i2)) then varmin = minval(var(i1:i2)) varmax = maxval(var(i1:i2)) else varmin = minval(var) varmax = maxval(var) end if if (varmin < boundmin .or. varmax > boundmax) then write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') & & '*** Warning: ', var_name, ' range', varmin, ' to', varmax, & & ' is out of physical range', boundmin, 'to', boundmax is_bad = .true. if (do_fix) then if (present(i1) .and. present(i2)) then var(i1:i2) = max(boundmin, min(boundmax, var(i1:i2))) else var = max(boundmin, min(boundmax, var)) end if write(nulout,'(a)') ': corrected' else write(nulout,'(1x)') end if end if end if end function out_of_bounds_1d !--------------------------------------------------------------------- ! Return .true. if 2D allocatable array "var" is out of physical ! range specified by boundmin and boundmax, and issue a warning. To ! check only a subset of the array, specify i1 and i2 for the range ! of the first dimension and j1 and j2 for the range of the second. function out_of_bounds_2d(var, var_name, boundmin, boundmax, do_fix, & & i1, i2, j1, j2) result (is_bad) use radiation_io, only : nulout real(jprb), allocatable, intent(inout) :: var(:,:) character(len=*), intent(in) :: var_name real(jprb), intent(in) :: boundmin, boundmax logical, intent(in) :: do_fix integer, optional, intent(in) :: i1, i2, j1, j2 ! Local copies of indices integer :: ii1, ii2, jj1, jj2 logical :: is_bad real(jprb) :: varmin, varmax is_bad = .false. if (allocated(var)) then if (present(i1) .and. present(i2)) then ii1 = i1 ii2 = i2 else ii1 = lbound(var,1) ii2 = ubound(var,1) end if if (present(j1) .and. present(j2)) then jj1 = j1 jj2 = j2 else jj1 = lbound(var,2) jj2 = ubound(var,2) end if varmin = minval(var(ii1:ii2,jj1:jj2)) varmax = maxval(var(ii1:ii2,jj1:jj2)) if (varmin < boundmin .or. varmax > boundmax) then write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') & & '*** Warning: ', var_name, ' range', varmin, ' to', varmax,& & ' is out of physical range', boundmin, 'to', boundmax is_bad = .true. if (do_fix) then var(ii1:ii2,jj1:jj2) = max(boundmin, min(boundmax, var(ii1:ii2,jj1:jj2))) write(nulout,'(a)') ': corrected' else write(nulout,'(1x)') end if end if end if end function out_of_bounds_2d !--------------------------------------------------------------------- ! Return .true. if 3D allocatable array "var" is out of physical ! range specified by boundmin and boundmax, and issue a warning. To ! check only a subset of the array, specify i1 and i2 for the range ! of the first dimension, j1 and j2 for the second and k1 and k2 for ! the third. function out_of_bounds_3d(var, var_name, boundmin, boundmax, do_fix, & & i1, i2, j1, j2, k1, k2) result (is_bad) use radiation_io, only : nulout real(jprb), allocatable, intent(inout) :: var(:,:,:) character(len=*), intent(in) :: var_name real(jprb), intent(in) :: boundmin, boundmax logical, intent(in) :: do_fix integer, optional, intent(in) :: i1, i2, j1, j2, k1, k2 ! Local copies of indices integer :: ii1, ii2, jj1, jj2, kk1, kk2 logical :: is_bad real(jprb) :: varmin, varmax is_bad = .false. if (allocated(var)) then if (present(i1) .and. present(i2)) then ii1 = i1 ii2 = i2 else ii1 = lbound(var,1) ii2 = ubound(var,1) end if if (present(j1) .and. present(j2)) then jj1 = j1 jj2 = j2 else jj1 = lbound(var,2) jj2 = ubound(var,2) end if if (present(k1) .and. present(k2)) then kk1 = k1 kk2 = k2 else kk1 = lbound(var,3) kk2 = ubound(var,3) end if varmin = minval(var(ii1:ii2,jj1:jj2,kk1:kk2)) varmax = maxval(var(ii1:ii2,jj1:jj2,kk1:kk2)) if (varmin < boundmin .or. varmax > boundmax) then write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') & & '*** Warning: ', var_name, ' range', varmin, ' to', varmax,& & ' is out of physical range', boundmin, 'to', boundmax is_bad = .true. if (do_fix) then var(ii1:ii2,jj1:jj2,kk1:kk2) = max(boundmin, min(boundmax, & & var(ii1:ii2,jj1:jj2,kk1:kk2))) write(nulout,'(a)') ': corrected' else write(nulout,'(1x)') end if end if end if end function out_of_bounds_3d end module radiation_config