diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_cloud_optics.F90 b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_cloud_optics.F90
new file mode 100644
index 0000000000000000000000000000000000000000..0fa20eb824abb7c52908766d5829644651e31143
--- /dev/null
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_cloud_optics.F90
@@ -0,0 +1,491 @@
+! radiation_cloud_optics.F90 - Computing cloud optical properties
+!
+! (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
+
+module radiation_cloud_optics
+
+ implicit none
+ public
+
+contains
+
+ ! Provides elemental function "delta_eddington_scat_od"
+#include "radiation_delta_eddington.h"
+
+ !---------------------------------------------------------------------
+ ! Load cloud scattering data; this subroutine delegates to one
+ ! in radiation_cloud_optics_data.F90, but checks the size of
+ ! what is returned
+ subroutine setup_cloud_optics(config)
+
+ use parkind1, only : jprb
+ use yomhook, only : lhook, dr_hook
+
+ use radiation_io, only : nulerr, radiation_abort
+ use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
+ & IIceModelBaran2016, IIceModelBaran2017, &
+ & IIceModelYi, &
+ & ILiquidModelSOCRATES, ILiquidModelSlingo
+ use radiation_cloud_optics_data, only : cloud_optics_type
+ use radiation_ice_optics_fu, only : NIceOpticsCoeffsFuSW, &
+ & NIceOpticsCoeffsFuLW
+ use radiation_ice_optics_baran, only : NIceOpticsCoeffsBaran, &
+ & NIceOpticsCoeffsBaran2016
+ use radiation_ice_optics_baran2017, only : NIceOpticsCoeffsBaran2017, &
+ & NIceOpticsGeneralCoeffsBaran2017
+ use radiation_ice_optics_yi, only : NIceOpticsCoeffsYiSW, &
+ & NIceOpticsCoeffsYiLW
+ use radiation_liquid_optics_socrates, only : NLiqOpticsCoeffsSOCRATES
+ use radiation_liquid_optics_slingo, only : NLiqOpticsCoeffsSlingoSW, &
+ & NLiqOpticsCoeffsLindnerLiLW
+
+ type(config_type), intent(inout) :: config
+
+ real(jprb) :: hook_handle
+
+ if (lhook) call dr_hook('radiation_cloud_optics:setup_cloud_optics',0,hook_handle)
+
+ call config%cloud_optics%setup(trim(config%liq_optics_file_name), &
+ & trim(config%ice_optics_file_name), iverbose=config%iverbosesetup)
+
+ ! Check liquid coefficients
+ if (size(config%cloud_optics%liq_coeff_lw, 1) /= config%n_bands_lw) then
+ write(nulerr,'(a,i0,a,i0,a)') &
+ & '*** Error: number of longwave bands for droplets (', &
+ & size(config%cloud_optics%liq_coeff_lw, 1), &
+ & ') does not match number for gases (', config%n_bands_lw, ')'
+ call radiation_abort()
+ end if
+ if (size(config%cloud_optics%liq_coeff_sw, 1) /= config%n_bands_sw) then
+ write(nulerr,'(a,i0,a,i0,a)') &
+ & '*** Error: number of shortwave bands for droplets (', &
+ & size(config%cloud_optics%liq_coeff_sw, 1), &
+ & ') does not match number for gases (', config%n_bands_sw, ')'
+ call radiation_abort()
+ end if
+
+ if (config%i_liq_model == ILiquidModelSOCRATES) then
+ if (size(config%cloud_optics%liq_coeff_lw, 2) /= NLiqOpticsCoeffsSOCRATES) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of liquid cloud optical coefficients (', &
+ & size(config%cloud_optics%liq_coeff_lw, 2), &
+ & ') does not match number expected (', NLiqOpticsCoeffsSOCRATES,')'
+ call radiation_abort()
+ end if
+ else if (config%i_liq_model == ILiquidModelSlingo) then
+ if (size(config%cloud_optics%liq_coeff_sw, 2) /= NLiqOpticsCoeffsSlingoSW) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of shortwave liquid cloud optical coefficients (', &
+ & size(config%cloud_optics%liq_coeff_sw, 2), &
+ & ') does not match number expected (', NLiqOpticsCoeffsSlingoSW,')'
+ call radiation_abort()
+ end if
+ if (size(config%cloud_optics%liq_coeff_lw, 2) /= NLiqOpticsCoeffsLindnerLiLW) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of longwave liquid cloud optical coefficients (', &
+ & size(config%cloud_optics%liq_coeff_lw, 2), &
+ & ') does not match number expected (', NLiqOpticsCoeffsLindnerLiLw,')'
+ call radiation_abort()
+ end if
+ end if
+
+ ! Check ice coefficients
+ if (size(config%cloud_optics%ice_coeff_lw, 1) /= config%n_bands_lw) then
+ write(nulerr,'(a,i0,a,i0,a)') &
+ & '*** Error: number of longwave bands for ice particles (', &
+ & size(config%cloud_optics%ice_coeff_lw, 1), &
+ & ') does not match number for gases (', config%n_bands_lw, ')'
+ call radiation_abort()
+ end if
+ if (size(config%cloud_optics%ice_coeff_sw, 1) /= config%n_bands_sw) then
+ write(nulerr,'(a,i0,a,i0,a)') &
+ & '*** Error: number of shortwave bands for ice particles (', &
+ & size(config%cloud_optics%ice_coeff_sw, 1), &
+ & ') does not match number for gases (', config%n_bands_sw, ')'
+ call radiation_abort()
+ end if
+
+ if (config%i_ice_model == IIceModelFu) then
+ if (size(config%cloud_optics%ice_coeff_lw, 2) &
+ & /= NIceOpticsCoeffsFuLW) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of LW ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_lw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsFuLW,')'
+ call radiation_abort()
+ end if
+ if (size(config%cloud_optics%ice_coeff_sw, 2) &
+ & /= NIceOpticsCoeffsFuSW) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of SW ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_sw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsFuSW,')'
+ call radiation_abort()
+ end if
+ else if (config%i_ice_model == IIceModelBaran &
+ & .and. size(config%cloud_optics%ice_coeff_lw, 2) &
+ & /= NIceOpticsCoeffsBaran) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_lw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsBaran,')'
+ call radiation_abort()
+ else if (config%i_ice_model == IIceModelBaran2016 &
+ & .and. size(config%cloud_optics%ice_coeff_lw, 2) &
+ & /= NIceOpticsCoeffsBaran2016) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_lw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsBaran2016,')'
+ call radiation_abort()
+ else if (config%i_ice_model == IIceModelBaran2017) then
+ if (size(config%cloud_optics%ice_coeff_lw, 2) &
+ & /= NIceOpticsCoeffsBaran2017) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_lw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsBaran2017,')'
+ call radiation_abort()
+ else if (.not. allocated(config%cloud_optics%ice_coeff_gen)) then
+ write(nulerr,'(a)') &
+ & '*** Error: coeff_gen needed for Baran-2017 ice optics parameterization'
+ call radiation_abort()
+ else if (size(config%cloud_optics%ice_coeff_gen) &
+ & /= NIceOpticsGeneralCoeffsBaran2017) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of general ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_gen), &
+ & ') does not match number expected (', NIceOpticsGeneralCoeffsBaran2017,')'
+ call radiation_abort()
+ end if
+ else if (config%i_ice_model == IIceModelYi) then
+ if (size(config%cloud_optics%ice_coeff_lw, 2) &
+ & /= NIceOpticsCoeffsYiLW) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of LW ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_lw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsYiLW,')'
+ call radiation_abort()
+ end if
+ if (size(config%cloud_optics%ice_coeff_sw, 2) &
+ & /= NIceOpticsCoeffsYiSW) then
+ write(nulerr,'(a,i0,a,i0,a,i0,a)') &
+ & '*** Error: number of SW ice-particle optical coefficients (', &
+ & size(config%cloud_optics%ice_coeff_sw, 2), &
+ & ') does not match number expected (', NIceOpticsCoeffsYiSW,')'
+ call radiation_abort()
+ end if
+ end if
+
+ if (lhook) call dr_hook('radiation_cloud_optics:setup_cloud_optics',1,hook_handle)
+
+ end subroutine setup_cloud_optics
+
+
+ !---------------------------------------------------------------------
+ ! Compute cloud optical properties
+ subroutine cloud_optics(nlev,istartcol,iendcol, &
+ & config, thermodynamics, cloud, &
+ & od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
+ & od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
+
+ use parkind1, only : jprb
+ use yomhook, only : lhook, dr_hook
+
+ use radiation_io, only : nulout, nulerr, radiation_abort
+ use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
+ & IIceModelBaran2016, IIceModelBaran2017, &
+ & IIceModelYi, &
+ & ILiquidModelSOCRATES, ILiquidModelSlingo
+ use radiation_thermodynamics, only : thermodynamics_type
+ use radiation_cloud, only : cloud_type
+ use radiation_constants, only : AccelDueToGravity
+ use radiation_cloud_optics_data, only : cloud_optics_type
+ use radiation_ice_optics_fu, only : calc_ice_optics_fu_sw, &
+ & calc_ice_optics_fu_lw
+ use radiation_ice_optics_baran, only : calc_ice_optics_baran, &
+ & calc_ice_optics_baran2016
+ use radiation_ice_optics_baran2017, only : calc_ice_optics_baran2017
+ use radiation_ice_optics_yi, only : calc_ice_optics_yi_sw, &
+ & calc_ice_optics_yi_lw
+ use radiation_liquid_optics_socrates, only:calc_liq_optics_socrates
+ use radiation_liquid_optics_slingo, only:calc_liq_optics_slingo, &
+ & calc_liq_optics_lindner_li
+
+ integer, intent(in) :: nlev ! number of model levels
+ integer, intent(in) :: istartcol, iendcol ! range of columns to process
+ type(config_type), intent(in), target :: config
+ type(thermodynamics_type),intent(in) :: thermodynamics
+ type(cloud_type), intent(in) :: cloud
+
+ ! Layer optical depth, single scattering albedo and g factor of
+ ! clouds in each longwave band, where the latter two
+ ! variables are only defined if cloud longwave scattering is
+ ! enabled (otherwise both are treated as zero).
+ real(jprb), dimension(config%n_bands_lw,nlev,istartcol:iendcol), intent(out) :: &
+ & od_lw_cloud
+ real(jprb), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol), &
+ & intent(out) :: ssa_lw_cloud, g_lw_cloud
+
+ ! Layer optical depth, single scattering albedo and g factor of
+ ! clouds in each shortwave band
+ real(jprb), dimension(config%n_bands_sw,nlev,istartcol:iendcol), intent(out) :: &
+ & od_sw_cloud, ssa_sw_cloud, g_sw_cloud
+
+ ! Longwave and shortwave optical depth, scattering optical depth
+ ! and asymmetry factor, for liquid and ice in all bands but a
+ ! single cloud layer
+ real(jprb), dimension(config%n_bands_lw) :: &
+ & od_lw_liq, scat_od_lw_liq, g_lw_liq, &
+ & od_lw_ice, scat_od_lw_ice, g_lw_ice
+ real(jprb), dimension(config%n_bands_sw) :: &
+ & od_sw_liq, scat_od_sw_liq, g_sw_liq, &
+ & od_sw_ice, scat_od_sw_ice, g_sw_ice
+
+ ! In-cloud water path of cloud liquid or ice (i.e. liquid or ice
+ ! gridbox-mean water path divided by cloud fraction); kg m-2
+ real(jprb) :: lwp_in_cloud, iwp_in_cloud
+
+ ! Full-level temperature (K)
+ real(jprb) :: temperature
+
+ ! Factor to convert gridbox-mean mixing ratio to in-cloud water
+ ! path
+ real(jprb) :: factor
+
+ ! Pointer to the cloud optics coefficients for brevity of
+ ! access
+ type(cloud_optics_type), pointer :: ho
+
+ integer :: jcol, jlev
+
+ real(jprb) :: hook_handle
+
+ if (lhook) call dr_hook('radiation_cloud_optics:cloud_optics',0,hook_handle)
+
+ if (config%iverbose >= 2) then
+ write(nulout,'(a)') 'Computing cloud absorption/scattering properties'
+ end if
+
+ ho => config%cloud_optics
+
+ ! Array-wise assignment
+ od_lw_cloud = 0.0_jprb
+ od_sw_cloud = 0.0_jprb
+ ssa_sw_cloud = 0.0_jprb
+ g_sw_cloud = 0.0_jprb
+ if (config%do_lw_cloud_scattering) then
+ ssa_lw_cloud = 0.0_jprb
+ g_lw_cloud = 0.0_jprb
+ end if
+
+ do jlev = 1,nlev
+ do jcol = istartcol,iendcol
+ ! Only do anything if cloud is present (assume that
+ ! cloud%crop_cloud_fraction has already been called)
+ if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
+
+ ! Compute in-cloud liquid and ice water path
+ if (config%is_homogeneous) then
+ ! Homogeneous solvers assume cloud fills the box
+ ! horizontally, so we don't divide by cloud fraction
+ factor = ( thermodynamics%pressure_hl(jcol,jlev+1) &
+ & -thermodynamics%pressure_hl(jcol,jlev ) ) &
+ & / AccelDueToGravity
+ else
+ factor = ( thermodynamics%pressure_hl(jcol,jlev+1) &
+ & -thermodynamics%pressure_hl(jcol,jlev ) ) &
+ & / (AccelDueToGravity * cloud%fraction(jcol,jlev))
+ end if
+ lwp_in_cloud = factor * cloud%q_liq(jcol,jlev)
+ iwp_in_cloud = factor * cloud%q_ice(jcol,jlev)
+
+ ! Only compute liquid properties if liquid cloud is
+ ! present
+ if (lwp_in_cloud > 0.0_jprb) then
+ if (config%i_liq_model == ILiquidModelSOCRATES) then
+ ! Compute longwave properties
+ call calc_liq_optics_socrates(config%n_bands_lw, &
+ & config%cloud_optics%liq_coeff_lw, &
+ & lwp_in_cloud, cloud%re_liq(jcol,jlev), &
+ & od_lw_liq, scat_od_lw_liq, g_lw_liq)
+ ! Compute shortwave properties
+ call calc_liq_optics_socrates(config%n_bands_sw, &
+ & config%cloud_optics%liq_coeff_sw, &
+ & lwp_in_cloud, cloud%re_liq(jcol,jlev), &
+ & od_sw_liq, scat_od_sw_liq, g_sw_liq)
+ else if (config%i_liq_model == ILiquidModelSlingo) then
+ ! Compute longwave properties
+ call calc_liq_optics_lindner_li(config%n_bands_lw, &
+ & config%cloud_optics%liq_coeff_lw, &
+ & lwp_in_cloud, cloud%re_liq(jcol,jlev), &
+ & od_lw_liq, scat_od_lw_liq, g_lw_liq)
+ ! Compute shortwave properties
+ call calc_liq_optics_slingo(config%n_bands_sw, &
+ & config%cloud_optics%liq_coeff_sw, &
+ & lwp_in_cloud, cloud%re_liq(jcol,jlev), &
+ & od_sw_liq, scat_od_sw_liq, g_sw_liq)
+ else
+ write(nulerr,*) '*** Error: Unknown liquid model with code', &
+ & config%i_liq_model
+ call radiation_abort()
+ end if
+
+ if (.not. config%do_sw_delta_scaling_with_gases) then
+ ! Delta-Eddington scaling in the shortwave only
+ call delta_eddington_scat_od(od_sw_liq, scat_od_sw_liq, g_sw_liq)
+ end if
+ else
+ ! Liquid not present: set properties to zero
+ od_lw_liq = 0.0_jprb
+ scat_od_lw_liq = 0.0_jprb
+ g_lw_liq = 0.0_jprb
+
+ od_sw_liq = 0.0_jprb
+ scat_od_sw_liq = 0.0_jprb
+ g_sw_liq = 0.0_jprb
+ end if ! Liquid present
+
+ ! Only compute ice properties if ice cloud is present
+ if (iwp_in_cloud > 0.0_jprb) then
+ if (config%i_ice_model == IIceModelBaran) then
+ ! Compute longwave properties
+ call calc_ice_optics_baran(config%n_bands_lw, &
+ & config%cloud_optics%ice_coeff_lw, &
+ & iwp_in_cloud, cloud%q_ice(jcol,jlev), &
+ & od_lw_ice, scat_od_lw_ice, g_lw_ice)
+ ! Compute shortwave properties
+ call calc_ice_optics_baran(config%n_bands_sw, &
+ & config%cloud_optics%ice_coeff_sw, &
+ & iwp_in_cloud, cloud%q_ice(jcol,jlev), &
+ & od_sw_ice, scat_od_sw_ice, g_sw_ice)
+ else if (config%i_ice_model == IIceModelBaran2016) then
+ temperature = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev) &
+ & +thermodynamics%temperature_hl(jcol,jlev+1))
+ ! Compute longwave properties
+ call calc_ice_optics_baran2016(config%n_bands_lw, &
+ & config%cloud_optics%ice_coeff_lw, &
+ & iwp_in_cloud, cloud%q_ice(jcol,jlev), &
+ & temperature, &
+ & od_lw_ice, scat_od_lw_ice, g_lw_ice)
+ ! Compute shortwave properties
+ call calc_ice_optics_baran2016(config%n_bands_sw, &
+ & config%cloud_optics%ice_coeff_sw, &
+ & iwp_in_cloud, cloud%q_ice(jcol,jlev), &
+ & temperature, &
+ & od_sw_ice, scat_od_sw_ice, g_sw_ice)
+ else if (config%i_ice_model == IIceModelBaran2017) then
+ temperature = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev) &
+ & +thermodynamics%temperature_hl(jcol,jlev+1))
+ ! Compute longwave properties
+ call calc_ice_optics_baran2017(config%n_bands_lw, &
+ & config%cloud_optics%ice_coeff_gen, &
+ & config%cloud_optics%ice_coeff_lw, &
+ & iwp_in_cloud, cloud%q_ice(jcol,jlev), &
+ & temperature, &
+ & od_lw_ice, scat_od_lw_ice, g_lw_ice)
+ ! Compute shortwave properties
+ call calc_ice_optics_baran2017(config%n_bands_sw, &
+ & config%cloud_optics%ice_coeff_gen, &
+ & config%cloud_optics%ice_coeff_sw, &
+ & iwp_in_cloud, cloud%q_ice(jcol,jlev), &
+ & temperature, &
+ & od_sw_ice, scat_od_sw_ice, g_sw_ice)
+ else if (config%i_ice_model == IIceModelFu) then
+ ! Compute longwave properties
+ call calc_ice_optics_fu_lw(config%n_bands_lw, &
+ & config%cloud_optics%ice_coeff_lw, &
+ & iwp_in_cloud, cloud%re_ice(jcol,jlev), &
+ & od_lw_ice, scat_od_lw_ice, g_lw_ice)
+ if (config%do_fu_lw_ice_optics_bug) then
+ ! Reproduce bug in old IFS scheme
+ scat_od_lw_ice = od_lw_ice - scat_od_lw_ice
+ end if
+ ! Compute shortwave properties
+ call calc_ice_optics_fu_sw(config%n_bands_sw, &
+ & config%cloud_optics%ice_coeff_sw, &
+ & iwp_in_cloud, cloud%re_ice(jcol,jlev), &
+ & od_sw_ice, scat_od_sw_ice, g_sw_ice)
+ else if (config%i_ice_model == IIceModelYi) then
+ ! Compute longwave properties
+ call calc_ice_optics_yi_lw(config%n_bands_lw, &
+ & config%cloud_optics%ice_coeff_lw, &
+ & iwp_in_cloud, cloud%re_ice(jcol,jlev), &
+ & od_lw_ice, scat_od_lw_ice, g_lw_ice)
+ ! Compute shortwave properties
+ call calc_ice_optics_yi_sw(config%n_bands_sw, &
+ & config%cloud_optics%ice_coeff_sw, &
+ & iwp_in_cloud, cloud%re_ice(jcol,jlev), &
+ & od_sw_ice, scat_od_sw_ice, g_sw_ice)
+ else
+ write(nulerr,*) '*** Error: Unknown ice model with code', &
+ & config%i_ice_model
+ call radiation_abort()
+ end if
+
+ if (.not. config%do_sw_delta_scaling_with_gases) then
+ ! Delta-Eddington scaling in both longwave and shortwave
+ ! (assume that particles are larger than wavelength even
+ ! in longwave)
+ call delta_eddington_scat_od(od_sw_ice, scat_od_sw_ice, g_sw_ice)
+ end if
+
+ call delta_eddington_scat_od(od_lw_ice, scat_od_lw_ice, g_lw_ice)
+ else
+ ! Ice not present: set properties to zero
+ od_lw_ice = 0.0_jprb
+ scat_od_lw_ice = 0.0_jprb
+ g_lw_ice = 0.0_jprb
+
+ od_sw_ice = 0.0_jprb
+ scat_od_sw_ice = 0.0_jprb
+ g_sw_ice = 0.0_jprb
+ end if ! Ice present
+
+ ! Combine liquid and ice
+ if (config%do_lw_cloud_scattering) then
+ od_lw_cloud(:,jlev,jcol) = od_lw_liq + od_lw_ice
+ where (scat_od_lw_liq+scat_od_lw_ice > 0.0_jprb)
+ g_lw_cloud(:,jlev,jcol) = (g_lw_liq * scat_od_lw_liq &
+ & + g_lw_ice * scat_od_lw_ice) &
+ & / (scat_od_lw_liq+scat_od_lw_ice)
+ elsewhere
+ g_lw_cloud(:,jlev,jcol) = 0.0_jprb
+ end where
+ ssa_lw_cloud(:,jlev,jcol) = (scat_od_lw_liq + scat_od_lw_ice) &
+ & / (od_lw_liq + od_lw_ice)
+ else
+ ! If longwave scattering is to be neglected then the
+ ! best approximation is to set the optical depth equal
+ ! to the absorption optical depth
+ od_lw_cloud(:,jlev,jcol) = od_lw_liq - scat_od_lw_liq &
+ & + od_lw_ice - scat_od_lw_ice
+ end if
+ od_sw_cloud(:,jlev,jcol) = od_sw_liq + od_sw_ice
+ g_sw_cloud(:,jlev,jcol) = (g_sw_liq * scat_od_sw_liq &
+ & + g_sw_ice * scat_od_sw_ice) &
+ & / (scat_od_sw_liq + scat_od_sw_ice)
+ ssa_sw_cloud(:,jlev,jcol) &
+ & = (scat_od_sw_liq + scat_od_sw_ice) / (od_sw_liq + od_sw_ice)
+ end if ! Cloud present
+ end do ! Loop over column
+ end do ! Loop over level
+
+ if (lhook) call dr_hook('radiation_cloud_optics:cloud_optics',1,hook_handle)
+
+ end subroutine cloud_optics
+
+end module radiation_cloud_optics
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_config.F90 b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_config.F90
new file mode 100644
index 0000000000000000000000000000000000000000..2a77f3a3ee06265dd224dac37cb6015cff0e2565
--- /dev/null
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_config.F90
@@ -0,0 +1,1980 @@
+! 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
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_delta_eddington.h b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_delta_eddington.h
new file mode 100644
index 0000000000000000000000000000000000000000..04f43bafde36958c84e50a5662009d817975b3a2
--- /dev/null
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_delta_eddington.h
@@ -0,0 +1,93 @@
+! radiation_delta_eddington.h - Delta-Eddington scaling
+!
+! (C) Copyright 2015- 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
+!
+! This file is intended to be included inside a module to ensure that
+! these simple functions may be inlined
+
+!---------------------------------------------------------------------
+! Perform in-place delta-Eddington scaling of the phase function
+elemental subroutine delta_eddington(od, ssa, g)
+
+ use parkind1, only : jprb
+
+ ! Total optical depth, single scattering albedo and asymmetry
+ ! factor
+ real(jprb), intent(inout) :: od, ssa, g
+
+ ! Fraction of the phase function deemed to be in the forward lobe
+ ! and therefore treated as if it is not scattered at all
+ real(jprb) :: f
+
+ f = g*g
+ od = od * (1.0_jprb - ssa*f)
+ ssa = ssa * (1.0_jprb - f) / (1.0_jprb - ssa*f)
+ g = g / (1.0_jprb + g)
+
+end subroutine delta_eddington
+
+
+!---------------------------------------------------------------------
+! Perform in-place delta-Eddington scaling of the phase function, but
+! using extensive variables (i.e. the scattering optical depth,
+! scat_od, rather than the single-scattering albedo, and the
+! scattering-optical-depth-multiplied-by-asymmetry-factor, scat_od_g,
+! rather than the asymmetry factor.
+elemental subroutine delta_eddington_extensive(od, scat_od, scat_od_g)
+
+ use parkind1, only : jprb
+
+ ! Total optical depth, scattering optical depth and asymmetry factor
+ ! multiplied by the scattering optical depth
+ real(jprb), intent(inout) :: od, scat_od, scat_od_g
+
+ ! Fraction of the phase function deemed to be in the forward lobe
+ ! and therefore treated as if it is not scattered at all
+ real(jprb) :: f, g
+
+ if (scat_od > 0.0_jprb) then
+ g = scat_od_g / scat_od
+ else
+ g = 0.0
+ end if
+
+ f = g*g
+ od = od - scat_od * f
+ scat_od = scat_od * (1.0_jprb - f)
+ scat_od_g = scat_od * g / (1.0_jprb + g)
+
+end subroutine delta_eddington_extensive
+
+
+!---------------------------------------------------------------------
+! Perform in-place delta-Eddington scaling of the phase function,
+! using the scattering optical depth rather than the single scattering
+! albedo
+elemental subroutine delta_eddington_scat_od(od, scat_od, g)
+
+ use parkind1, only : jprb
+
+ ! Total optical depth, scattering optical depth and asymmetry factor
+ real(jprb), intent(inout) :: od, scat_od, g
+
+ ! Fraction of the phase function deemed to be in the forward lobe
+ ! and therefore treated as if it is not scattered at all
+ real(jprb) :: f
+
+ f = g*g
+ od = od - scat_od * f
+ scat_od = scat_od * (1.0_jprb - f)
+ g = g / (1.0_jprb + g)
+
+end subroutine delta_eddington_scat_od
+