diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/data/YBUR_ice_scattering_rrtm.nc b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YBUR_ice_scattering_rrtm.nc
new file mode 100644
index 0000000000000000000000000000000000000000..cfcaf7ff2ecada6346d89231a840851ffa998bc8
Binary files /dev/null and b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YBUR_ice_scattering_rrtm.nc differ
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/data/YCOL_ice_scattering_rrtm.nc b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YCOL_ice_scattering_rrtm.nc
new file mode 100644
index 0000000000000000000000000000000000000000..b51cdc499d6352c8d79a690fce41bef49eb95706
Binary files /dev/null and b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YCOL_ice_scattering_rrtm.nc differ
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/data/YDRO_ice_scattering_rrtm.nc b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YDRO_ice_scattering_rrtm.nc
new file mode 100644
index 0000000000000000000000000000000000000000..f6412d1704477f6d1fe5cfda9d8068ec59c7f394
Binary files /dev/null and b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YDRO_ice_scattering_rrtm.nc differ
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/data/YHBU_ice_scattering_rrtm.nc b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YHBU_ice_scattering_rrtm.nc
new file mode 100644
index 0000000000000000000000000000000000000000..127f166ddd6073d17afeb31ec962628965571cd1
Binary files /dev/null and b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YHBU_ice_scattering_rrtm.nc differ
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/data/YHCO_ice_scattering_rrtm.nc b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YHCO_ice_scattering_rrtm.nc
new file mode 100644
index 0000000000000000000000000000000000000000..3b3ad6fef7c98f08bc9d2811d298c2dc0bc1911b
Binary files /dev/null and b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YHCO_ice_scattering_rrtm.nc differ
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/data/YPLA_ice_scattering_rrtm.nc b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YPLA_ice_scattering_rrtm.nc
new file mode 100644
index 0000000000000000000000000000000000000000..92baaf1e9b331dfc351a01b1841f2bd813bd07d3
Binary files /dev/null and b/src/LIB/RAD/ecrad-1.4.0_mnh/data/YPLA_ice_scattering_rrtm.nc differ
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/ifs/radiation_setup.F90 b/src/LIB/RAD/ecrad-1.4.0_mnh/ifs/radiation_setup.F90
index 0677d3f94c66fd2ab89ac6efee68a4e993669f38..69d2d44fc5601fa4672b0192a98d00c428d52483 100644
--- a/src/LIB/RAD/ecrad-1.4.0_mnh/ifs/radiation_setup.F90
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/ifs/radiation_setup.F90
@@ -40,7 +40,7 @@ MODULE RADIATION_SETUP
USE radiation_config, ONLY : config_type, &
& ISolverMcICA, ISolverSpartacus, &
& ILiquidModelSlingo, ILiquidModelSOCRATES, &
- & IIceModelFu, IIceModelBaran, &
+ & IIceModelFu, IIceModelBaran, IIceModelShapes, &
& IOverlapExponential
USE MODD_PARAM_ECRAD_n , ONLY : rad_config
@@ -161,6 +161,8 @@ CONTAINS
rad_config%i_ice_model = IIceModelFu
ELSEIF (NICEOPT == 4) THEN
rad_config%i_ice_model = IIceModelBaran
+ ELSEIF (NICEOPT == 7) THEN
+ rad_config%i_ice_model = IIceModelShapes
ELSE
WRITE(NULERR,'(a,i0)') 'Unavailable ice optics model in modular radiation scheme: NICEOPT=', &
& NICEOPT
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
index 0fa20eb824abb7c52908766d5829644651e31143..f636fd604bac0b890f3dbaa2916c6a88b56b7cd4 100644
--- 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
@@ -18,6 +18,7 @@
module radiation_cloud_optics
implicit none
+
public
contains
@@ -37,7 +38,7 @@ contains
use radiation_io, only : nulerr, radiation_abort
use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
& IIceModelBaran2016, IIceModelBaran2017, &
- & IIceModelYi, &
+ & IIceModelYi,IIceModelShapes, &
& ILiquidModelSOCRATES, ILiquidModelSlingo
use radiation_cloud_optics_data, only : cloud_optics_type
use radiation_ice_optics_fu, only : NIceOpticsCoeffsFuSW, &
@@ -48,6 +49,8 @@ contains
& NIceOpticsGeneralCoeffsBaran2017
use radiation_ice_optics_yi, only : NIceOpticsCoeffsYiSW, &
& NIceOpticsCoeffsYiLW
+ use radiation_ice_optics_shapes, only : NIceOpticsCoeffsShapesSW, &
+ & NIceOpticsCoeffsShapesLW
use radiation_liquid_optics_socrates, only : NLiqOpticsCoeffsSOCRATES
use radiation_liquid_optics_slingo, only : NLiqOpticsCoeffsSlingoSW, &
& NLiqOpticsCoeffsLindnerLiLW
@@ -188,6 +191,23 @@ contains
& ') does not match number expected (', NIceOpticsCoeffsYiSW,')'
call radiation_abort()
end if
+ else if (config%i_ice_model == IIceModelShapes) then
+ if (size(config%cloud_optics%ice_coeff_lw, 2) &
+ & /= NIceOpticsCoeffsShapesLW) 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 (', NIceOpticsCoeffsShapesLW,')'
+ call radiation_abort()
+ end if
+ if (size(config%cloud_optics%ice_coeff_sw, 2) &
+ & /= NIceOpticsCoeffsShapesSW) 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 (', NIceOpticsCoeffsShapesSW,')'
+ call radiation_abort()
+ end if
end if
if (lhook) call dr_hook('radiation_cloud_optics:setup_cloud_optics',1,hook_handle)
@@ -208,7 +228,7 @@ contains
use radiation_io, only : nulout, nulerr, radiation_abort
use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
& IIceModelBaran2016, IIceModelBaran2017, &
- & IIceModelYi, &
+ & IIceModelYi,IIceModelShapes, &
& ILiquidModelSOCRATES, ILiquidModelSlingo
use radiation_thermodynamics, only : thermodynamics_type
use radiation_cloud, only : cloud_type
@@ -221,6 +241,8 @@ contains
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_ice_optics_shapes, only : calc_ice_optics_shapes_sw, &
+ & calc_ice_optics_shapes_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
@@ -344,10 +366,12 @@ contains
call radiation_abort()
end if
+ ! Delta-Eddington scaling in the shortwave only
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
+ !call delta_eddington_scat_od(od_lw_liq, scat_od_lw_liq, g_lw_liq)
+
else
! Liquid not present: set properties to zero
od_lw_liq = 0.0_jprb
@@ -427,6 +451,17 @@ contains
& 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 if (config%i_ice_model == IIceModelShapes) then
+ ! Compute longwave properties
+ call calc_ice_optics_shapes_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_shapes_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)
@@ -436,14 +471,14 @@ contains
call radiation_abort()
end if
+ ! Delta-Eddington scaling in both longwave and shortwave
+ ! (assume that particles are larger than wavelength even
+ ! in longwave)
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
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
index 2a77f3a3ee06265dd224dac37cb6015cff0e2565..0f641bb4012cc8cf73955dd206e029fbc9632970 100644
--- 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
@@ -32,7 +32,7 @@
!
module radiation_config
-
+ USE MODD_PARAM_LIMA, ONLY : CPRISTINE_ICE_LIMA
use parkind1, only : jprb
use radiation_cloud_optics_data, only : cloud_optics_type
@@ -69,7 +69,6 @@ module radiation_config
& 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
@@ -110,14 +109,15 @@ module radiation_config
enum, bind(c)
enumerator IIceModelMonochromatic, IIceModelFu, &
& IIceModelBaran, IIceModelBaran2016, IIceModelBaran2017, &
- & IIceModelYi
+ & IIceModelYi, IIceModelShapes
end enum
- character(len=*), parameter :: IceModelName(0:5) = (/ 'Monochromatic', &
+ character(len=*), parameter :: IceModelName(0:6) = (/ 'Monochromatic', &
& 'Fu-IFS ', &
& 'Baran ', &
& 'Baran2016 ', &
& 'Baran2017 ', &
- & 'Yi ' /)
+ & 'Yi ', &
+ & 'Shapes ' /)
! Cloud PDF distribution shapes
enum, bind(c)
@@ -238,7 +238,7 @@ module radiation_config
! 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
+ ! integer :: i_cloud_model
! Optics if i_gas_model==IGasModelMonochromatic.
! The wavelength to use for the Planck function in metres. If this
@@ -274,10 +274,10 @@ module radiation_config
! 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.
+ ! 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
@@ -419,8 +419,8 @@ module radiation_config
! 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) :: 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
@@ -428,8 +428,7 @@ module radiation_config
character(len=511) :: cloud_pdf_override_file_name = ''
! Has "consolidate" been called?
- logical :: is_consolidated = .false.
-
+ logical :: is_consolidated = .false.
! COMPUTED PARAMETERS
! Users of this library should not edit these parameters directly;
! they are set by the "consolidate" routine
@@ -441,6 +440,7 @@ module radiation_config
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
@@ -502,7 +502,7 @@ module radiation_config
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
@@ -601,7 +601,7 @@ contains
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.
+ 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)
@@ -1006,6 +1006,9 @@ contains
else if (this%i_ice_model == IIceModelYi) then
this%ice_optics_file_name &
& = trim(this%directory_name) // "/yi_ice_scattering_rrtm.nc"
+ else if (this%i_ice_model == IIceModelShapes) then
+ this%ice_optics_file_name &
+ & = trim(this%directory_name) // "/" // CPRISTINE_ICE_LIMA // "_ice_scattering_rrtm.nc"
end if
! Set cloud-water PDF look-up table file name
@@ -1326,7 +1329,7 @@ contains
if (this%n_bands_sw <= 0) then
write(nulerr,'(a)') '*** Error: get_sw_weights called before number of shortwave bands set'
- call radiation_abort()
+ call radiation_abort()
end if
! Convert wavelength range (m) to wavenumber (cm-1)
@@ -1976,5 +1979,5 @@ contains
end function out_of_bounds_3d
-
+
end module radiation_config
diff --git a/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_ice_optics_shapes.F90 b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_ice_optics_shapes.F90
new file mode 100644
index 0000000000000000000000000000000000000000..011f2666e762e97b6c081ef98f0b1710620277db
--- /dev/null
+++ b/src/LIB/RAD/ecrad-1.4.0_mnh/radiation/radiation_ice_optics_shapes.F90
@@ -0,0 +1,142 @@
+! radiation_ice_optics_shapes.F90 -
+!
+! (C) Copyright 2017- 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.
+!
+! Authors: Marie Taufour
+! Email: marie.taufour@aero.obs-mip.fr
+!
+! The reference for this ice optics parameterization is ....
+
+module radiation_ice_optics_shapes
+
+ implicit none
+ public
+
+ ! The number of ice coefficients depends on the parameterization
+ integer, parameter :: NIceOpticsCoeffsShapesSW = 69
+ integer, parameter :: NIceOpticsCoeffsShapesLW = 69
+
+ integer, parameter :: NSingleCoeffs = 23
+
+contains
+
+ !---------------------------------------------------------------------
+ ! Compute shortwave ice-particle scattering properties using Yi et
+ ! al. (2013) parameterization
+ subroutine calc_ice_optics_shapes_sw(nb, coeff, ice_wp, &
+ & re, od, scat_od, g)
+
+ use parkind1, only : jprb, jpim
+ !use yomhook, only : lhook, dr_hook
+
+ ! Number of bands
+ integer, intent(in) :: nb
+ ! Coefficients read from a data file
+ real(jprb), intent(in) :: coeff(:,:)
+ ! Ice water path (kg m-2)
+ real(jprb), intent(in) :: ice_wp
+ ! Effective radius (m)
+ real(jprb), intent(in) :: re
+ ! Total optical depth, scattering optical depth and asymmetry factor
+ real(jprb), intent(out) :: od(nb), scat_od(nb), g(nb)
+
+ ! Yi's effective diameter (microns)
+ real(jprb) :: de_um
+ ! Ice water path in g m-2
+ real (jprb) :: iwp_gm_2
+ ! LUT temp variables
+ real(jprb) :: wts_1, wts_2
+ integer(jpim) :: lu_idx
+ real(kind=jprb), parameter :: lu_scale = 0.2_jprb
+ real(kind=jprb), parameter :: lu_offset = 1.0_jprb
+ !real(jprb) :: hook_handle
+
+ !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_yi_sw',0,hook_handle)
+
+ ! Convert to effective diameter using the relationship in the IFS
+ !de_um = re * (1.0e6_jprb / 0.64952_jprb)
+ de_um = re * 2.0e6_jprb
+
+ ! limit de_um to validity of LUT
+ de_um = max(de_um,10.0_jprb)
+ de_um = min(de_um,119.99_jprb) !avoid greater than or equal to 120 um
+
+ iwp_gm_2 = ice_wp * 1000.0_jprb
+
+ lu_idx = floor(de_um * lu_scale - lu_offset)
+ wts_2 = (de_um * lu_scale - lu_offset) - lu_idx
+ wts_1 = 1.0_jprb - wts_2
+ od = 0.001_jprb * iwp_gm_2 * &
+ & ( wts_1 * coeff(1:nb,lu_idx) + wts_2 * coeff(1:nb,lu_idx+1) )
+ scat_od = od * &
+ & ( wts_1 * coeff(1:nb,lu_idx+NSingleCoeffs) + wts_2 * coeff(1:nb,lu_idx+NSingleCoeffs+1) )
+ g = wts_1 * coeff(1:nb,lu_idx+2*NSingleCoeffs) + wts_2 * coeff(1:nb,lu_idx+2*NSingleCoeffs+1)
+
+ !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_yi_sw',1,hook_handle)
+
+ end subroutine calc_ice_optics_shapes_sw
+
+
+ !---------------------------------------------------------------------
+ ! Compute longwave ice-particle scattering properties using .....
+ subroutine calc_ice_optics_shapes_lw(nb, coeff, ice_wp, &
+ & re, od, scat_od, g)
+
+ use parkind1, only : jprb, jpim
+ !use yomhook, only : lhook, dr_hook
+
+ ! Number of bands
+ integer, intent(in) :: nb
+ ! Coefficients read from a data file
+ real(jprb), intent(in) :: coeff(:,:)
+ ! Ice water path (kg m-2)
+ real(jprb), intent(in) :: ice_wp
+ ! Effective radius (m)
+ real(jprb), intent(in) :: re
+ ! Total optical depth, scattering optical depth and asymmetry factor
+ real(jprb), intent(out) :: od(nb), scat_od(nb), g(nb)
+
+ ! Yi's effective diameter (microns)
+ real(jprb) :: de_um
+ ! Ice water path in g m-2
+ real (jprb) :: iwp_gm_2
+ ! LUT temp variables
+ real(jprb) :: wts_1, wts_2
+ integer(jpim) :: lu_idx
+ real(kind=jprb), parameter :: lu_scale = 0.2_jprb
+ real(kind=jprb), parameter :: lu_offset = 1.0_jprb
+ !real(jprb) :: hook_handle
+
+ !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_yi_sw',0,hook_handle)
+
+ ! Convert to effective diameter using the relationship in the IFS
+ !de_um = re * (1.0e6_jprb / 0.64952_jprb)
+ de_um = re * 2.0e6_jprb
+
+ ! limit de_um to validity of LUT
+ de_um = max(de_um,10.0_jprb)
+ de_um = min(de_um,119.99_jprb) !avoid greater than or equal to 120 um
+
+ iwp_gm_2 = ice_wp * 1000.0_jprb
+
+ lu_idx = floor(de_um * lu_scale - lu_offset)
+ wts_2 = (de_um * lu_scale - lu_offset) - lu_idx
+ wts_1 = 1.0_jprb - wts_2
+ od = 0.001_jprb * iwp_gm_2 * &
+ & ( wts_1 * coeff(1:nb,lu_idx) + wts_2 * coeff(1:nb,lu_idx+1) )
+ scat_od = od * &
+ & ( wts_1 * coeff(1:nb,lu_idx+NSingleCoeffs) + wts_2 * coeff(1:nb,lu_idx+NSingleCoeffs+1) )
+ g = wts_1 * coeff(1:nb,lu_idx+2*NSingleCoeffs) + wts_2 * coeff(1:nb,lu_idx+2*NSingleCoeffs+1)
+
+ !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_yi_lw',1,hook_handle)
+
+ end subroutine calc_ice_optics_shapes_lw
+
+end module radiation_ice_optics_shapes