Ignore:
Timestamp:
Mar 31, 2023, 8:42:57 PM (20 months ago)
Author:
lguez
Message:

Merge LMDZ_ECRad branch back into trunk!

Location:
LMDZ6/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation_config.F90

    r4115 r4489  
    2626!   2019-02-03  R. Hogan  Added ability to fix out-of-physical-bounds inputs
    2727!   2019-02-10  R. Hogan  Renamed "encroachment" to "entrapment"
     28!   2020-05-18  R. Hogan  Moved out_of_bounds_* to radiation_check.F90
     29!   2021-07-04  R. Hogan  Numerous changes for ecCKD and general cloud/aerosol optics
    2830!
    2931! Note: The aim is for ecRad in the IFS to be as similar as possible
    … …  
    3739
    3840  use radiation_cloud_optics_data,   only : cloud_optics_type
     41  use radiation_general_cloud_optics_data,   only : general_cloud_optics_type
    3942  use radiation_aerosol_optics_data, only : aerosol_optics_type
    4043  use radiation_pdf_sampler,         only : pdf_sampler_type
    4144  use radiation_cloud_cover,         only : OverlapName, &
    4245       & IOverlapMaximumRandom, IOverlapExponentialRandom, IOverlapExponential
     46  use radiation_ecckd,               only : ckd_model_type
    4347
    4448  implicit none
    … …  
    7074       & IEntrapmentExplicitNonFractal, & ! As above but ignore fractal nature of clouds
    7175       & IEntrapmentMaximum ! Complete horizontal homogenization within regions (old SPARTACUS assumption)
    72 
    7376  end enum
    7477 
    … …  
    9497  ! Gas models
    9598  enum, bind(c)
    96      enumerator IGasModelMonochromatic, IGasModelIFSRRTMG
     99     enumerator IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD
    97100  end enum
    98   character(len=*), parameter :: GasModelName(0:1) = (/ 'Monochromatic', &
    99        &                                                'RRTMG-IFS    ' /)
     101  character(len=*), parameter :: GasModelName(0:2) = (/ 'Monochromatic', &
     102       &                                                'RRTMG-IFS    ', &
     103       &                                                'ECCKD        '/)
    100104
    101105  ! Hydrometeor scattering models
    … …  
    130134  integer, parameter :: NMaxAerosolTypes = 256
    131135
     136  ! Maximum number of different cloud types that can be provided
     137  integer, parameter :: NMaxCloudTypes = 12
     138
    132139  ! Maximum number of shortwave albedo and longwave emissivity
    133140  ! intervals
    … …  
    155162    character(len=511) :: directory_name = '.'
    156163
     164    ! If this is true then support arbitrary hydrometeor types (not
     165    ! just ice and liquid) and arbitrary spectral discretization (not
     166    ! just RRTMG). It is required that this is true if the ecCKD gas
     167    ! optics model is selected. General cloud optics has only been
     168    ! available from ecRad version 1.5.
     169    logical :: use_general_cloud_optics = .true.
     170
     171    ! If this is true then support aerosol properties at an arbitrary
     172    ! spectral discretization (not just RRTMG). It is required that
     173    ! this is true if the ecCKD gas optics model is selected.
     174    logical :: use_general_aerosol_optics = .true.
     175
    157176    ! Cloud is deemed to be present in a layer if cloud fraction
    158177    ! exceeds this value
    … …  
    168187    ! (2000)?
    169188    logical :: use_beta_overlap = .false.
     189
     190    ! Use a more vectorizable McICA cloud generator, at the expense of
     191    ! more random numbers being generated?  This is the default on NEC
     192    ! SX.
     193#ifdef __SX__
     194    logical :: use_vectorizable_generator = .true.
     195#else
     196    logical :: use_vectorizable_generator = .false.
     197#endif
    170198
    171199    ! Shape of sub-grid cloud water PDF
    … …  
    236264    logical :: do_sw_delta_scaling_with_gases = .false.
    237265
    238     ! Codes describing the gas and cloud scattering models to use, the
    239     ! latter of which is currently not used
     266    ! Codes describing the gas model
    240267    integer :: i_gas_model = IGasModelIFSRRTMG
    241     !     integer :: i_cloud_model
    242268
    243269    ! Optics if i_gas_model==IGasModelMonochromatic.
    … …  
    270296    ! according to the spectral overlap of each interval with each
    271297    ! band
    272     logical :: do_nearest_spectral_sw_albedo = .true.
    273     logical :: do_nearest_spectral_lw_emiss  = .true.
     298    logical :: do_nearest_spectral_sw_albedo = .false.
     299    logical :: do_nearest_spectral_lw_emiss  = .false.
    274300
    275301    ! User-defined monotonically increasing wavelength bounds (m)
    276302    ! between input surface albedo/emissivity intervals. Implicitly
    277     ! the first interval starts at zero and the last ends at infinity.
     303    ! the first interval starts at zero and the last ends at
     304    ! infinity. These must be set with define_sw_albedo_intervals and
     305    ! define_lw_emiss_intervals.
    278306    real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1) = -1.0_jprb
    279     real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1)  = -1.0_jprb
     307    real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) = -1.0_jprb
    280308
    281309    ! The index to the surface albedo/emissivity intervals for each of
    … …  
    296324    logical :: do_3d_effects = .true.
    297325   
     326    character(len=511) :: cloud_type_name(NMaxCloudTypes) = ["","","","","","","","","","","",""]
     327! &
     328!         &   = ["mie_droplet                   ", &
     329!         &      "baum-general-habit-mixture_ice"]
     330
     331    ! Spectral averaging method to use with generalized cloud optics;
     332    ! see Edwards & Slingo (1996) for definition.  Experimentation
     333    ! with ecRad suggests that "thick" averaging is more accurate for
     334    ! both liquid and ice clouds.
     335    logical :: use_thick_cloud_spectral_averaging(NMaxCloudTypes) &
     336         &  = [.true.,.true.,.true.,.true.,.true.,.true., &
     337         &     .true.,.true.,.true.,.true.,.true.,.true.]
     338
    298339    ! To what extent do we include "entrapment" effects in the
    299340    ! SPARTACUS solver? This essentially means that in a situation
    … …  
    420461    ! doesn't start with a '/' character then it will be prepended by
    421462    ! the contents of directory_name.
    422     character(len=511) :: ice_optics_override_file_name = ''
    423     character(len=511) :: liq_optics_override_file_name = ''
     463    character(len=511) :: ice_optics_override_file_name     = ''
     464    character(len=511) :: liq_optics_override_file_name     = ''
    424465    character(len=511) :: aerosol_optics_override_file_name = ''
     466    character(len=511) :: gas_optics_sw_override_file_name  = ''
     467    character(len=511) :: gas_optics_lw_override_file_name  = ''
    425468
    426469    ! Optionally override the look-up table file for the cloud-water
    … …  
    428471    character(len=511) :: cloud_pdf_override_file_name = ''
    429472
     473    ! Do we compute cloud, aerosol and surface optical properties per
     474    ! g point?  Not available with RRTMG gas optics model.
     475    logical :: do_cloud_aerosol_per_sw_g_point = .true.
     476    logical :: do_cloud_aerosol_per_lw_g_point = .true.
     477
     478    ! Do we weight the mapping from surface emissivity/albedo to
     479    ! g-point/band weighting by a reference Planck function (more
     480    ! accurate) or weighting each wavenumber equally (less accurate
     481    ! but consistent with IFS Cycle 48r1 and earlier)?
     482    logical :: do_weighted_surface_mapping = .true.
     483
     484    ! COMPUTED PARAMETERS
     485
     486    ! Users of this library should not edit these parameters directly;
     487    ! they are set by the "consolidate" routine
     488
    430489    ! Has "consolidate" been called? 
    431490    logical :: is_consolidated = .false.
    432491
    433     ! COMPUTED PARAMETERS
    434     ! Users of this library should not edit these parameters directly;
    435     ! they are set by the "consolidate" routine
    436 
    437     ! Wavenumber range for each band, in cm-1, which will be allocated
    438     ! to be of length n_bands_sw or n_bands_lw
    439     real(jprb), allocatable, dimension(:) :: wavenumber1_sw
    440     real(jprb), allocatable, dimension(:) :: wavenumber2_sw
    441     real(jprb), allocatable, dimension(:) :: wavenumber1_lw
    442     real(jprb), allocatable, dimension(:) :: wavenumber2_lw
     492    ! Fraction of each g point in each wavenumber interval,
     493    ! dimensioned (n_wav_frac_[l|s]w, n_g_[l|s]w)
     494    real(jprb), allocatable, dimension(:,:) :: g_frac_sw, g_frac_lw
    443495
    444496    ! If the nearest surface albedo/emissivity interval is to be used
    … …  
    490542    integer :: n_canopy_bands_lw = 1
    491543
     544    ! Data structures containing gas optics description in the case of
     545    ! ecCKD
     546    type(ckd_model_type)         :: gas_optics_sw, gas_optics_lw
     547
    492548    ! Data structure containing cloud scattering data
    493549    type(cloud_optics_type)      :: cloud_optics
     550
     551    ! Number of general cloud types, default liquid and ice
     552    integer :: n_cloud_types = 2
     553
     554    ! List of data structures (one per cloud type) containing cloud
     555    ! scattering data
     556    type(general_cloud_optics_type), allocatable :: cloud_optics_sw(:)
     557    type(general_cloud_optics_type), allocatable :: cloud_optics_lw(:)
    494558
    495559    ! Data structure containing aerosol scattering data
    … …  
    502566    character(len=511) :: ice_optics_file_name, &
    503567         &                liq_optics_file_name, &
    504          &                aerosol_optics_file_name
     568         &                aerosol_optics_file_name, &
     569         &                gas_optics_sw_file_name, &
     570         &                gas_optics_lw_file_name
    505571   
    506572    ! McICA PDF look-up table file name
    … …  
    515581    ! g points or the number of bands
    516582    integer :: n_spec_sw = 0, n_spec_lw = 0
     583
     584    ! Number of wavenumber intervals used to describe the mapping from
     585    ! g-points to wavenumber space
     586    integer :: n_wav_frac_sw = 0, n_wav_frac_lw = 0
    517587
    518588    ! Dimensions to store variables that are only needed if longwave
    … …  
    539609     procedure :: define_sw_albedo_intervals
    540610     procedure :: define_lw_emiss_intervals
    541      procedure :: consolidate_intervals
     611     procedure :: set_aerosol_wavelength_mono
     612     procedure :: consolidate_sw_albedo_intervals
     613     procedure :: consolidate_lw_emiss_intervals
    542614
    543615  end type config_type
    … …  
    574646    logical :: do_sw, do_lw, do_clear, do_sw_direct
    575647    logical :: do_3d_effects, use_expm_everywhere, use_aerosols
     648    logical :: use_general_cloud_optics, use_general_aerosol_optics
    576649    logical :: do_lw_side_emissivity
    577650    logical :: do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug
    … …  
    579652    logical :: do_save_radiative_properties, do_save_spectral_flux
    580653    logical :: do_save_gpoint_flux, do_surface_sw_spectral_flux
    581     logical :: use_beta_overlap, do_lw_derivatives
     654    logical :: use_beta_overlap, do_lw_derivatives, use_vectorizable_generator
    582655    logical :: do_sw_delta_scaling_with_gases
    583656    logical :: do_canopy_fluxes_sw, do_canopy_fluxes_lw
    584657    logical :: use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw
    585658    logical :: do_canopy_gases_sw, do_canopy_gases_lw
     659    logical :: do_cloud_aerosol_per_sw_g_point, do_cloud_aerosol_per_lw_g_point
     660    logical :: do_weighted_surface_mapping   
    586661    integer :: n_regions, iverbose, iverbosesetup, n_aerosol_types
    587662    real(jprb):: mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od
    … …  
    596671    character(511) :: liq_optics_override_file_name, ice_optics_override_file_name
    597672    character(511) :: cloud_pdf_override_file_name
     673    character(511) :: gas_optics_sw_override_file_name, gas_optics_lw_override_file_name
    598674    character(63)  :: liquid_model_name, ice_model_name, gas_model_name
    599675    character(63)  :: sw_solver_name, lw_solver_name, overlap_scheme_name
    600676    character(63)  :: sw_entrapment_name, sw_encroachment_name, cloud_pdf_shape_name
     677    character(len=511) :: cloud_type_name(NMaxCloudTypes) = ["","","","","","","","","","","",""]
     678    logical :: use_thick_cloud_spectral_averaging(NMaxCloudTypes) &
     679         &  = [.false.,.false.,.false.,.false.,.false.,.false., &
     680         &     .false.,.false.,.false.,.false.,.false.,.false.]
    601681    integer :: i_aerosol_type_map(NMaxAerosolTypes) ! More than 256 is an error
    602682
    603     logical :: do_nearest_spectral_sw_albedo = .true.
    604     logical :: do_nearest_spectral_lw_emiss  = .true.
     683    logical :: do_nearest_spectral_sw_albedo
     684    logical :: do_nearest_spectral_lw_emiss
    605685    real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1)
    606686    real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1)
    … …  
    609689
    610690    integer :: iunit ! Unit number of namelist file
    611 
    612     logical :: lldeb_conf = .false.
    613691
    614692    namelist /radiation/ do_sw, do_lw, do_sw_direct, &
    … …  
    622700         &  ice_optics_override_file_name, liq_optics_override_file_name, &
    623701         &  aerosol_optics_override_file_name, cloud_pdf_override_file_name, &
     702         &  gas_optics_sw_override_file_name, gas_optics_lw_override_file_name, &
    624703         &  liquid_model_name, ice_model_name, max_3d_transfer_rate, &
    625704         &  min_cloud_effective_size, overhang_factor, encroachment_scaling, &
    … …  
    627706         &  do_canopy_fluxes_sw, do_canopy_fluxes_lw, &
    628707         &  do_canopy_gases_sw, do_canopy_gases_lw, &
     708         &  use_general_cloud_optics, use_general_aerosol_optics, &
    629709         &  do_sw_delta_scaling_with_gases, overlap_scheme_name, &
    630          &  sw_solver_name, lw_solver_name, use_beta_overlap, &
     710         &  sw_solver_name, lw_solver_name, use_beta_overlap, use_vectorizable_generator, &
    631711         &  use_expm_everywhere, iverbose, iverbosesetup, &
    632712         &  cloud_inhom_decorr_scaling, cloud_fraction_threshold, &
    … …  
    637717         &  mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo, &
    638718         &  mono_lw_asymmetry_factor, mono_sw_asymmetry_factor, &
    639          &  cloud_pdf_shape_name, &
     719         &  cloud_pdf_shape_name, cloud_type_name, use_thick_cloud_spectral_averaging, &
    640720         &  do_nearest_spectral_sw_albedo, do_nearest_spectral_lw_emiss, &
    641721         &  sw_albedo_wavelength_bound, lw_emiss_wavelength_bound, &
    642          &  i_sw_albedo_index, i_lw_emiss_index
    643 
     722         &  i_sw_albedo_index, i_lw_emiss_index, &
     723         &  do_cloud_aerosol_per_lw_g_point, &
     724         &  do_cloud_aerosol_per_sw_g_point, do_weighted_surface_mapping
     725         
    644726    real(jprb) :: hook_handle
    645727
    … …  
    670752    ice_optics_override_file_name = this%ice_optics_override_file_name
    671753    aerosol_optics_override_file_name = this%aerosol_optics_override_file_name
     754    gas_optics_sw_override_file_name = this%gas_optics_sw_override_file_name
     755    gas_optics_lw_override_file_name = this%gas_optics_lw_override_file_name
    672756    use_expm_everywhere = this%use_expm_everywhere
    673757    use_aerosols = this%use_aerosols
    … …  
    679763    iverbose = this%iverbose
    680764    iverbosesetup = this%iverbosesetup
     765    use_general_cloud_optics = this%use_general_cloud_optics
     766    use_general_aerosol_optics = this%use_general_aerosol_optics
    681767    cloud_fraction_threshold = this%cloud_fraction_threshold
    682768    cloud_mixing_ratio_threshold = this%cloud_mixing_ratio_threshold
    683769    use_beta_overlap = this%use_beta_overlap
     770    use_vectorizable_generator = this%use_vectorizable_generator
    684771    cloud_inhom_decorr_scaling = this%cloud_inhom_decorr_scaling
    685772    clear_to_thick_fraction = this%clear_to_thick_fraction
    … …  
    689776    max_3d_transfer_rate = this%max_3d_transfer_rate
    690777    min_cloud_effective_size = this%min_cloud_effective_size
     778    cloud_type_name = this%cloud_type_name
     779    use_thick_cloud_spectral_averaging = this%use_thick_cloud_spectral_averaging
     780
    691781    overhang_factor = this%overhang_factor
    692782    encroachment_scaling = -1.0_jprb
    … …  
    715805    i_sw_albedo_index             = this%i_sw_albedo_index
    716806    i_lw_emiss_index              = this%i_lw_emiss_index
     807    do_cloud_aerosol_per_lw_g_point = this%do_cloud_aerosol_per_lw_g_point
     808    do_cloud_aerosol_per_sw_g_point = this%do_cloud_aerosol_per_sw_g_point
     809    do_weighted_surface_mapping   = this%do_weighted_surface_mapping
    717810
    718811    if (present(file_name) .and. present(unit)) then
    … …  
    746839      end if
    747840    else
     841
     842      ! This version exits correctly, but provides less information
     843      ! about how the namelist was incorrect
    748844      read(unit=iunit, iostat=iosread, nml=radiation)
     845
     846      ! Depending on compiler this version provides more information
     847      ! about the error in the namelist
     848      !read(unit=iunit, nml=radiation)
     849
    749850      if (iosread /= 0) then
    750851        ! An error occurred reading the file
    … …  
    812913    this%mono_sw_asymmetry_factor = mono_sw_asymmetry_factor
    813914    this%use_beta_overlap = use_beta_overlap
     915    this%use_vectorizable_generator = use_vectorizable_generator
    814916    this%cloud_inhom_decorr_scaling = cloud_inhom_decorr_scaling
    815917    this%clear_to_thick_fraction = clear_to_thick_fraction
    … …  
    819921    this%max_3d_transfer_rate = max_3d_transfer_rate
    820922    this%min_cloud_effective_size = max(1.0e-6_jprb, min_cloud_effective_size)
     923    this%cloud_type_name = cloud_type_name
     924    this%use_thick_cloud_spectral_averaging = use_thick_cloud_spectral_averaging
    821925    if (encroachment_scaling >= 0.0_jprb) then
    822926      this%overhang_factor = encroachment_scaling
    … …  
    832936    this%ice_optics_override_file_name = ice_optics_override_file_name
    833937    this%aerosol_optics_override_file_name = aerosol_optics_override_file_name
     938    this%gas_optics_sw_override_file_name = gas_optics_sw_override_file_name
     939    this%gas_optics_lw_override_file_name = gas_optics_lw_override_file_name
     940    this%use_general_cloud_optics      = use_general_cloud_optics
     941    this%use_general_aerosol_optics    = use_general_aerosol_optics
    834942    this%cloud_fraction_threshold = cloud_fraction_threshold
    835943    this%cloud_mixing_ratio_threshold = cloud_mixing_ratio_threshold
    … …  
    845953    this%i_sw_albedo_index             = i_sw_albedo_index
    846954    this%i_lw_emiss_index              = i_lw_emiss_index
    847 
    848 ! AI mars 2022
    849 if (lldeb_conf) then
    850 print*,'**************PARAMETRES DE CONFIGURATION OFFLINE*******************'
    851 print*,'config%iverbosesetup   = ', iverbosesetup
    852 print*,'config%do_lw   = ', do_lw
    853 print*,'config%do_sw   = ', do_sw
    854 print*,'config%do_clear   = ', do_clear
    855 print*,'config%do_sw_direct   = ', do_sw_direct
    856 print*,'config%do_3d_effects   = ', do_3d_effects
    857 print*,'config%do_3d_lw_multilayer_effects   = ', do_3d_lw_multilayer_effects
    858 print*,'config%do_lw_side_emissivity   = ', do_lw_side_emissivity
    859 print*,'config%use_expm_everywhere   = ', use_expm_everywhere
    860 print*,'config%use_aerosols   = ', use_aerosols
    861 print*,'config%do_lw_cloud_scattering   = ', do_lw_cloud_scattering
    862 print*,'config%do_lw_aerosol_scattering   = ', do_lw_aerosol_scattering
    863 print*,'config%nregions   = ', n_regions
    864 print*,'config%do_surface_sw_spectral_flux   = ', do_surface_sw_spectral_flux
    865 print*,'config%do_sw_delta_scaling_with_gases   = ', &
    866 do_sw_delta_scaling_with_gases
    867 print*,'config%do_fu_lw_ice_optics_bug   = ', do_fu_lw_ice_optics_bug
    868 print*,'config%do_canopy_fluxes_sw   = ', do_canopy_fluxes_sw
    869 print*,'config%do_canopy_fluxes_lw   = ', do_canopy_fluxes_lw
    870 print*,'config%use_canopy_full_spectrum_sw   = ', use_canopy_full_spectrum_sw
    871 print*,'config%use_canopy_full_spectrum_lw   = ', use_canopy_full_spectrum_lw
    872 print*,'config%do_canopy_gases_sw   = ', do_canopy_gases_sw
    873 print*,'config%do_canopy_gases_lw   = ', do_canopy_gases_lw
    874 print*,'config%mono_lw_wavelength   = ', mono_lw_wavelength
    875 print*,'config%mono_lw_total_od   = ', mono_lw_total_od
    876 print*,'config%mono_sw_total_od   = ', mono_sw_total_od
    877 print*,'config%mono_lw_single_scattering_albedo   = ', &
    878 mono_lw_single_scattering_albedo
    879 print*,'config%mono_sw_single_scattering_albedo   = ', &
    880 mono_sw_single_scattering_albedo
    881 print*,'config%mono_lw_asymmetry_factor   = ', mono_lw_asymmetry_factor
    882 print*,'config%mono_sw_asymmetry_factor   = ', mono_sw_asymmetry_factor
    883 print*,'config%use_beta_overlap   = ', use_beta_overlap
    884 print*,'config%cloud_inhom_decorr_scaling   = ', cloud_inhom_decorr_scaling
    885 print*,'config%clear_to_thick_fraction   = ', clear_to_thick_fraction
    886 print*,'config%overhead_sun_factor   = ', overhead_sun_factor
    887 print*,'config%max_gas_od_3d   = ', max_gas_od_3d
    888 print*,'config%max_cloud_od   = ', max_cloud_od
    889 print*,'config%max_3d_transfer_rate   = ', max_3d_transfer_rate
    890 print*,'config%min_cloud_effective_size   = ', &
    891 max(1.0e-6_jprb,min_cloud_effective_size)
    892 print*,'config%overhang_factor   = ', encroachment_scaling
    893 
    894 print*,'config%directory_name  = ',directory_name
    895 print*,'config%cloud_pdf_override_file_name  = ',cloud_pdf_override_file_name
    896 print*,'config%liq_optics_override_file_name  = ',liq_optics_override_file_name
    897 print*,'config%ice_optics_override_file_name  = ',ice_optics_override_file_name
    898 print*,'config%aerosol_optics_override_file_name  = ', &
    899 aerosol_optics_override_file_name
    900 print*,'config%cloud_fraction_threshold  = ',cloud_fraction_threshold
    901 print*,'config%cloud_mixing_ratio_threshold  = ',cloud_mixing_ratio_threshold
    902 print*,'config%n_aerosol_types  = ',n_aerosol_types
    903 print*,'config%do_save_radiative_properties  = ',do_save_radiative_properties
    904 print*,'config%do_lw_derivatives  = ',do_lw_derivatives
    905 print*,'config%do_save_spectral_flux  = ',do_save_spectral_flux
    906 print*,'config%do_save_gpoint_flux  = ',do_save_gpoint_flux
    907 print*,'config%do_nearest_spectral_sw_albedo  = ',do_nearest_spectral_sw_albedo
    908 print*,'config%do_nearest_spectral_lw_emiss   = ',do_nearest_spectral_lw_emiss
    909 print*,'config%sw_albedo_wavelength_bound     = ',sw_albedo_wavelength_bound
    910 print*,'config%lw_emiss_wavelength_bound      = ',lw_emiss_wavelength_bound
    911 print*,'config%i_sw_albedo_index              = ',i_sw_albedo_index
    912 print*,'config%i_lw_emiss_index               = ',i_lw_emiss_index
    913 print*,'************************************************************************'
    914 endif
     955    this%do_cloud_aerosol_per_lw_g_point = do_cloud_aerosol_per_lw_g_point
     956    this%do_cloud_aerosol_per_sw_g_point = do_cloud_aerosol_per_sw_g_point
     957    this%do_weighted_surface_mapping   = do_weighted_surface_mapping
     958
    915959    if (do_save_gpoint_flux) then
    916960      ! Saving the fluxes every g-point overrides saving as averaged
    … …  
    919963      ! save anything
    920964      this%do_save_spectral_flux = .true.
    921       print*,'config%do_save_spectral_flux = .true.'
    922965    end if
    923966
    … …  
    925968    call get_enum_code(liquid_model_name, LiquidModelName, &
    926969         &            'liquid_model_name', this%i_liq_model)
    927     print*,'config%i_liq_model =', this%i_liq_model
    928970
    929971    ! Determine ice optics model
    930972    call get_enum_code(ice_model_name, IceModelName, &
    931973         &            'ice_model_name', this%i_ice_model)
    932     print*,'config%i_ice_model =', this%i_ice_model
     974
    933975    ! Determine gas optics model
    934976    call get_enum_code(gas_model_name, GasModelName, &
    935977         &            'gas_model_name', this%i_gas_model)
    936     print*,'config%%i_gas_model = ', this%i_gas_model
    937978
    938979    ! Determine solvers
    939980    call get_enum_code(sw_solver_name, SolverName, &
    940981         &            'sw_solver_name', this%i_solver_sw)
    941     print*,'config%i_solver_sw = ', this%i_solver_sw
    942982    call get_enum_code(lw_solver_name, SolverName, &
    943983         &            'lw_solver_name', this%i_solver_lw)
    944     print*,'config%i_solver_lw = ', this%i_solver_lw
     984
    945985    if (len_trim(sw_encroachment_name) > 1) then
    946986      call get_enum_code(sw_encroachment_name, EncroachmentName, &
    … …  
    950990      call get_enum_code(sw_entrapment_name, EntrapmentName, &
    951991           &             'sw_entrapment_name', this%i_3d_sw_entrapment)
    952       print*,'config%i_3d_sw_entrapment = ', this%i_3d_sw_entrapment
    953992    end if
    954993
    … …  
    956995    call get_enum_code(overlap_scheme_name, OverlapName, &
    957996         &             'overlap_scheme_name', this%i_overlap_scheme)
    958     print*,'config%i_overlap_scheme = ', this%i_overlap_scheme
     997   
    959998    ! Determine cloud PDF shape
    960999    call get_enum_code(cloud_pdf_shape_name, PdfShapeName, &
    9611000         &             'cloud_pdf_shape_name', this%i_cloud_pdf_shape)
    962     print*,'config%i_cloud_pdf_shape = ', this%i_cloud_pdf_shape
     1001
    9631002    this%i_aerosol_type_map = 0
    9641003    if (this%use_aerosols) then
    9651004      this%i_aerosol_type_map(1:n_aerosol_types) &
    9661005           &  = i_aerosol_type_map(1:n_aerosol_types)
    967       print*,'config%i_aerosol_type_map = ', this%i_aerosol_type_map
    9681006    end if
    9691007
    … …  
    9751013      this%do_clouds = .false.
    9761014    end if
    977     print*,'config%do_clouds = ', this%do_clouds
     1015
     1016    if (this%i_gas_model == IGasModelIFSRRTMG &
     1017         & .and. (this%use_general_cloud_optics &
     1018         &        .or. this%use_general_aerosol_optics)) then
     1019      if (this%do_sw .and. this%do_cloud_aerosol_per_sw_g_point) then
     1020        write(nulout,'(a)') 'Warning: RRTMG SW only supports cloud/aerosol/surface optical properties per band, not per g-point'
     1021        this%do_cloud_aerosol_per_sw_g_point = .false.
     1022      end if
     1023      if (this%do_lw .and. this%do_cloud_aerosol_per_lw_g_point) then
     1024        write(nulout,'(a)') 'Warning: RRTMG LW only supports cloud/aerosol/surface optical properties per band, not per g-point'
     1025        this%do_cloud_aerosol_per_lw_g_point = .false.
     1026      end if
     1027    end if
     1028
    9781029
    9791030    ! Normal subroutine exit
    … …  
    9931044  subroutine consolidate_config(this)
    9941045
     1046    use parkind1,     only : jprd
    9951047    use yomhook,      only : lhook, dr_hook
    9961048    use radiation_io, only : nulout, nulerr, radiation_abort
    … …  
    10261078      write(nulerr,'(a)') '*** Error: SPARTACUS/Tripleclouds solvers can only do Exponential-Random overlap'
    10271079      call radiation_abort('Radiation configuration error')
     1080    end if
     1081
     1082    if (jprb < jprd .and. this%iverbosesetup >= 1 &
     1083         &  .and. (this%i_solver_sw == ISolverSPARTACUS &
     1084         &    .or. this%i_solver_lw == ISolverSPARTACUS)) then
     1085      write(nulout,'(a)') 'Warning: the SPARTACUS solver may be unstable in single precision'
     1086    end if
     1087
     1088    ! If ecCKD gas optics model is being used set relevant file names
     1089    if (this%i_gas_model == IGasModelECCKD) then
     1090
     1091      ! This gas optics model requires the general cloud and
     1092      ! aerosol optics settings
     1093      if (.not. this%use_general_cloud_optics) then
     1094        write(nulerr,'(a)') '*** Error: ecCKD gas optics model requires general cloud optics'
     1095        call radiation_abort('Radiation configuration error')
     1096      end if
     1097      if (.not. this%use_general_aerosol_optics) then
     1098        write(nulerr,'(a)') '*** Error: ecCKD gas optics model requires general aerosol optics'
     1099        call radiation_abort('Radiation configuration error')
     1100      end if
     1101
     1102      if (len_trim(this%gas_optics_sw_override_file_name) > 0) then
     1103        if (this%gas_optics_sw_override_file_name(1:1) == '/') then
     1104          this%gas_optics_sw_file_name = trim(this%gas_optics_sw_override_file_name)
     1105        else
     1106          this%gas_optics_sw_file_name = trim(this%directory_name) &
     1107               &  // '/' // trim(this%gas_optics_sw_override_file_name)
     1108        end if
     1109      else
     1110        ! In the IFS, the gas optics files should be specified in
     1111        ! ifs/module/radiation_setup.F90, not here
     1112        this%gas_optics_sw_file_name = trim(this%directory_name) &
     1113             &  // "/ecckd-1.0_sw_climate_rgb-32b_ckd-definition.nc"
     1114      end if
     1115
     1116      if (len_trim(this%gas_optics_lw_override_file_name) > 0) then
     1117        if (this%gas_optics_lw_override_file_name(1:1) == '/') then
     1118          this%gas_optics_lw_file_name = trim(this%gas_optics_lw_override_file_name)
     1119        else
     1120          this%gas_optics_lw_file_name = trim(this%directory_name) &
     1121               &  // '/' // trim(this%gas_optics_lw_override_file_name)
     1122        end if
     1123      else
     1124        ! In the IFS, the gas optics files should be specified in
     1125        ! ifs/module/radiation_setup.F90, not here
     1126        this%gas_optics_lw_file_name = trim(this%directory_name) &
     1127             &  // "/ecckd-1.0_lw_climate_fsck-32b_ckd-definition.nc"
     1128      end if
    10281129
    10291130    end if
    … …  
    10401141      ! In the IFS, the aerosol optics file should be specified in
    10411142      ! ifs/module/radiation_setup.F90, not here
    1042       this%aerosol_optics_file_name &
    1043            &   = trim(this%directory_name) // "/aerosol_ifs_rrtm_45R2.nc"
     1143      if (this%use_general_aerosol_optics) then
     1144         this%aerosol_optics_file_name &
     1145             &   = trim(this%directory_name) // "/aerosol_ifs_48R1.nc"       
     1146      else
     1147        this%aerosol_optics_file_name &
     1148             &   = trim(this%directory_name) // "/aerosol_ifs_rrtm_46R1_with_NI_AM.nc"
     1149      end if
    10441150    end if
    10451151
    … …  
    12291335           &          this%i_gas_model)
    12301336      call print_logical('  Aerosols are', 'use_aerosols', this%use_aerosols)
    1231       call print_logical('  Clouds are', 'do_clouds', this%do_clouds)
     1337      if (this%use_aerosols) then
     1338        call print_logical('  General aerosol optics', &
     1339             &             'use_general_aerosol_optics', this%use_general_aerosol_optics)
     1340      end if
     1341      if (this%do_clouds) then
     1342        write(nulout,'(a)') '  Clouds are ON'
     1343      else
     1344        write(nulout,'(a)') '  Clouds are OFF'
     1345      end if
     1346      if (this%do_sw) then
     1347        call print_logical('  Do cloud/aerosol/surface SW properties per g-point', &
     1348             &  'do_cloud_aerosol_per_sw_g_point', this%do_cloud_aerosol_per_sw_g_point)
     1349      end if
     1350      if (this%do_lw) then
     1351        call print_logical('  Do cloud/aerosol/surface LW properties per g-point', &
     1352             &  'do_cloud_aerosol_per_lw_g_point', this%do_cloud_aerosol_per_lw_g_point)
     1353      end if
    12321354
    12331355      !---------------------------------------------------------------------
    … …  
    12531375             &   'do_nearest_spectral_lw_emiss', this%do_nearest_spectral_lw_emiss)
    12541376      end if
     1377      call print_logical('  Planck-weighted surface albedo/emiss mapping', &
     1378           &   'do_weighted_surface_mapping', this%do_weighted_surface_mapping)
     1379
    12551380      !---------------------------------------------------------------------
    12561381      if (this%do_clouds) then
    … …  
    12601385        call print_real('  Cloud mixing-ratio threshold', &
    12611386             &   'cloud_mixing_ratio_threshold', this%cloud_mixing_ratio_threshold)
    1262         call print_enum('  Liquid optics scheme is', LiquidModelName, &
    1263              &          'i_liq_model',this%i_liq_model)
    1264         call print_enum('  Ice optics scheme is', IceModelName, &
    1265              &          'i_ice_model',this%i_ice_model)
    1266         if (this%i_ice_model == IIceModelFu) then
    1267           call print_logical('  Longwave ice optics bug in Fu scheme is', &
    1268                &   'do_fu_lw_ice_optics_bug',this%do_fu_lw_ice_optics_bug)
     1387        call print_logical('  General cloud optics', &
     1388             &             'use_general_cloud_optics', this%use_general_cloud_optics)
     1389        if (.not. this%use_general_cloud_optics) then
     1390          call print_enum('  Liquid optics scheme is', LiquidModelName, &
     1391               &          'i_liq_model',this%i_liq_model)
     1392          call print_enum('  Ice optics scheme is', IceModelName, &
     1393               &          'i_ice_model',this%i_ice_model)
     1394          if (this%i_ice_model == IIceModelFu) then
     1395            call print_logical('  Longwave ice optics bug in Fu scheme is', &
     1396                 &   'do_fu_lw_ice_optics_bug',this%do_fu_lw_ice_optics_bug)
     1397          end if
    12691398        end if
    12701399        call print_enum('  Cloud overlap scheme is', OverlapName, &
    … …  
    13601489               &   'overhang_factor', this%overhang_factor)
    13611490        end if
     1491
     1492      else if (this%i_solver_sw == ISolverMcICA &
     1493           &  .or. this%i_solver_lw == ISolverMcICA) then
     1494        call print_logical('  Use vectorizable McICA cloud generator', &
     1495             &   'use_vectorizable_generator', this%use_vectorizable_generator)
    13621496      end if
    13631497           
    … …  
    13831517    use parkind1, only : jprb
    13841518    use radiation_io, only : nulout, nulerr, radiation_abort
     1519    use radiation_spectral_definition, only : SolarReferenceTemperature
    13851520
    13861521    class(config_type), intent(in) :: this
    … …  
    13961531    character(len=*), optional, intent(in) :: weighting_name
    13971532
     1533    real(jprb), allocatable   :: mapping(:,:)
     1534
    13981535    ! Internally we deal with wavenumber
    13991536    real(jprb) :: wavenumber1, wavenumber2 ! cm-1
    14001537
     1538    real(jprb) :: wavenumber1_band, wavenumber2_band ! cm-1
     1539
    14011540    integer :: jband ! Loop index for spectral band
    14021541
    14031542    if (this%n_bands_sw <= 0) then
    14041543      write(nulerr,'(a)') '*** Error: get_sw_weights called before number of shortwave bands set'
    1405       call radiation_abort()     
     1544      call radiation_abort('Radiation configuration error')
    14061545    end if
    14071546
    … …  
    14101549    wavenumber2 = 0.01_jprb / wavelength1
    14111550
     1551    call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(SolarReferenceTemperature, &
     1552         &  [wavelength1, wavelength2], [1, 2, 3], mapping, &
     1553         &  use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point), use_fluxes=.true.)
     1554
     1555    ! "mapping" now contains a 3*nband matrix, where mapping(2,:)
     1556    ! contains the weights of interest.  We now find the non-zero weights
    14121557    nweights = 0
    1413 
    1414     do jband = 1,this%n_bands_sw
    1415       if (wavenumber1 < this%wavenumber2_sw(jband) &
    1416            &  .and. wavenumber2 > this%wavenumber1_sw(jband)) then
     1558    do jband = 1,size(mapping,2)
     1559      if (mapping(2,jband) > 0.0_jprb) then
    14171560        nweights = nweights+1
    1418         iband(nweights) = jband
    1419         weight(nweights) = (min(wavenumber2,this%wavenumber2_sw(jband)) &
    1420              &         - max(wavenumber1,this%wavenumber1_sw(jband))) &
    1421              & / (this%wavenumber2_sw(jband) - this%wavenumber1_sw(jband))
     1561        iband(nweights) = jband;
     1562        weight(nweights) = mapping(2,jband)
    14221563      end if
    14231564    end do
    … …  
    14261567      write(nulerr,'(a,e8.4,a,e8.4,a)') '*** Error: wavelength range ', &
    14271568           &  wavelength1, ' to ', wavelength2, ' m is outside shortwave band'
    1428       call radiation_abort()
     1569      call radiation_abort('Radiation configuration error')
    14291570    else if (this%iverbosesetup >= 2 .and. present(weighting_name)) then
    14301571      write(nulout,'(a,a,a,f6.0,a,f6.0,a)') 'Spectral weights for ', &
    14311572           &  weighting_name, ' (', wavenumber1, ' to ', &
    14321573           &  wavenumber2, ' cm-1):'
    1433       do jband = 1, nweights
    1434         write(nulout, '(a,i0,a,f6.0,a,f6.0,a,f8.4)') '  Shortwave band ', &
    1435              &  iband(jband), ' (', this%wavenumber1_sw(iband(jband)), ' to ', &
    1436              &  this%wavenumber2_sw(iband(jband)), ' cm-1): ', weight(jband)
    1437       end do
     1574      if (this%do_cloud_aerosol_per_sw_g_point) then
     1575        do jband = 1, nweights
     1576          write(nulout, '(a,i0,a,f8.4)') '  Shortwave g point ', iband(jband), ': ', weight(jband)
     1577        end do
     1578      else
     1579        do jband = 1, nweights
     1580          wavenumber1_band = this%gas_optics_sw%spectral_def%wavenumber1_band(iband(jband))
     1581          wavenumber2_band = this%gas_optics_sw%spectral_def%wavenumber2_band(iband(jband))
     1582          write(nulout, '(a,i0,a,f6.0,a,f6.0,a,f8.4)') '  Shortwave band ', &
     1583               &  iband(jband), ' (', wavenumber1_band, ' to ', &
     1584               &  wavenumber2_band, ' cm-1): ', weight(jband)
     1585        end do
     1586      end if
    14381587    end if
    14391588
    … …  
    14521601
    14531602    use radiation_io, only : nulerr, radiation_abort
     1603    use radiation_spectral_definition, only : SolarReferenceTemperature
    14541604
    14551605    class(config_type),   intent(inout) :: this
    … …  
    14671617      write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, &
    14681618           &  ' albedo intervals exceeds maximum of ', NMaxAlbedoIntervals
    1469       call radiation_abort();
     1619      call radiation_abort('Radiation configuration error')
    14701620    end if
    14711621
    … …  
    14801630    this%i_sw_albedo_index(ninterval+1:)           = 0
    14811631
     1632    ! If this routine is called before setup_radiation then the
     1633    ! spectral intervals are not yet known
     1634    ! consolidate_sw_albedo_intervals is called later.  Otherwise it
     1635    ! is called immediately and overwrites any existing mapping.
    14821636    if (this%is_consolidated) then
    1483       call this%consolidate_intervals(.true., &
    1484            &  this%do_nearest_spectral_sw_albedo, &
    1485            &  this%sw_albedo_wavelength_bound, this%i_sw_albedo_index, &
    1486            &  this%wavenumber1_sw, this%wavenumber2_sw, &
    1487            &  this%i_albedo_from_band_sw, this%sw_albedo_weights)
     1637      call this%consolidate_sw_albedo_intervals
    14881638    end if
    14891639
    … …  
    14971647
    14981648    use radiation_io, only : nulerr, radiation_abort
     1649    use radiation_spectral_definition, only : TerrestrialReferenceTemperature
    14991650
    15001651    class(config_type),   intent(inout) :: this
    … …  
    15121663      write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, &
    15131664           &  ' emissivity intervals exceeds maximum of ', NMaxAlbedoIntervals
    1514       call radiation_abort();
     1665      call radiation_abort('Radiation configuration error')
    15151666    end if
    15161667
    … …  
    15261677
    15271678    if (this%is_consolidated) then
    1528       call this%consolidate_intervals(.false., &
    1529            &  this%do_nearest_spectral_lw_emiss, &
    1530            &  this%lw_emiss_wavelength_bound, this%i_lw_emiss_index, &
    1531            &  this%wavenumber1_lw, this%wavenumber2_lw, &
    1532            &  this%i_emiss_from_band_lw, this%lw_emiss_weights)
     1679      call this%consolidate_lw_emiss_intervals
    15331680    end if
    15341681
    … …  
    15371684
    15381685  !---------------------------------------------------------------------
    1539   ! This routine consolidates either the input shortwave albedo
    1540   ! intervals with the shortwave bands, or the input longwave
    1541   ! emissivity intervals with the longwave bands, depending on the
    1542   ! arguments provided.
    1543   subroutine consolidate_intervals(this, is_sw, do_nearest, &
    1544        &  wavelength_bound, i_intervals, wavenumber1, wavenumber2, &
    1545        &  i_mapping, weights)
    1546 
    1547     use radiation_io, only : nulout, nulerr, radiation_abort
     1686  ! Set the wavelengths (m) at which monochromatic aerosol properties
     1687  ! are required. This routine must be called before consolidation of
     1688  ! settings.
     1689  subroutine set_aerosol_wavelength_mono(this, wavelength_mono)
     1690
     1691    use radiation_io, only : nulerr, radiation_abort
     1692   
     1693    class(config_type), intent(inout) :: this
     1694    real(jprb),         intent(in)    :: wavelength_mono(:)
     1695
     1696    if (this%is_consolidated) then
     1697      write(nulerr,'(a)') '*** Errror: set_aerosol_wavelength_mono must be called before setup_radiation'
     1698      call radiation_abort('Radiation configuration error')
     1699    end if
     1700   
     1701    if (allocated(this%aerosol_optics%wavelength_mono)) then
     1702      deallocate(this%aerosol_optics%wavelength_mono)
     1703    end if
     1704    allocate(this%aerosol_optics%wavelength_mono(size(wavelength_mono)))
     1705    this%aerosol_optics%wavelength_mono = wavelength_mono
     1706
     1707  end subroutine set_aerosol_wavelength_mono
     1708
     1709
     1710  !---------------------------------------------------------------------
     1711  ! Consolidate the surface shortwave albedo intervals with the
     1712  ! band/g-point intervals
     1713  subroutine consolidate_sw_albedo_intervals(this)
     1714
     1715    use radiation_io, only : nulout
     1716    use radiation_spectral_definition, only : SolarReferenceTemperature
    15481717
    15491718    class(config_type),   intent(inout) :: this
    1550     ! Is this the shortwave?  Otherwise longwave
    1551     logical,    intent(in)    :: is_sw
    1552     ! Do we find the nearest albedo interval to the centre of each
    1553     ! band, or properly weight the contributions? This can be modified
    1554     ! if there is only one albedo intervals.
    1555     logical, intent(inout)    :: do_nearest
    1556     ! Monotonically increasing wavelength bounds between intervals,
    1557     ! not including the outer bounds (which are assumed to be zero and
    1558     ! infinity)
    1559     real(jprb), intent(in)    :: wavelength_bound(NMaxAlbedoIntervals-1)
    1560     ! The albedo band indices corresponding to each interval
    1561     integer,    intent(in)    :: i_intervals(NMaxAlbedoIntervals)
    1562     ! Start and end wavenumber bounds for the ecRad bands (cm-1)
    1563     real(jprb), intent(in)    :: wavenumber1(:), wavenumber2(:)
    1564 
    1565     ! if do_nearest is TRUE then the result is expressed in i_mapping,
    1566     ! which will be allocated to have the same length as wavenumber1,
    1567     ! and contain the index of the albedo interval corresponding to
    1568     ! that band
    1569     integer,    allocatable, intent(inout) :: i_mapping(:)
    1570     ! ...otherwise the result is expressed in "weights", of
    1571     ! size(n_intervals, n_bands) containing how much of each interval
    1572     ! contributes to each band.
    1573     real(jprb), allocatable, intent(inout) :: weights(:,:)
    1574 
    1575     ! Number and loop index of ecRad bands
    1576     integer    :: nband, jband
    1577     ! Number and index of albedo/emissivity intervals
    1578     integer    :: ninterval, iinterval
    1579     ! Sometimes an albedo or emissivity value will be used in more
    1580     ! than one interval, so nvalue indicates how many values will
    1581     ! actually be provided
    1582     integer    :: nvalue
    1583     ! Wavenumber bounds of the albedo/emissivity interval
    1584     real(jprb) :: wavenumber1_albedo, wavenumber2_albedo
    1585     ! Reciprocal of the wavenumber range of the ecRad band
    1586     real(jprb) :: recip_dwavenumber ! cm
    1587     ! Midpoint/bound of wavenumber band
    1588     real(jprb) :: wavenumber_mid, wavenumber_bound ! cm-1
    1589    
    1590     nband = size(wavenumber1)
     1719
     1720    integer :: ninterval, jint, jband
    15911721
    15921722    ! Count the number of albedo/emissivity intervals
    15931723    ninterval = 0
    1594     do iinterval = 1,NMaxAlbedoIntervals
    1595       if (i_intervals(iinterval) > 0) then
    1596         ninterval = iinterval
     1724    do jint = 1,NMaxAlbedoIntervals
     1725      if (this%i_sw_albedo_index(jint) > 0) then
     1726        ninterval = jint
    15971727      else
    15981728        exit
    … …  
    16001730    end do
    16011731
    1602     if (ninterval < 2) then
    1603       ! Zero or one albedo/emissivity intervals found, so we index all
    1604       ! bands to one interval
    1605       if (allocated(i_mapping)) then
    1606         deallocate(i_mapping)
    1607       end if
    1608       allocate(i_mapping(nband))
    1609       i_mapping(:) = 1
    1610       do_nearest = .true.
     1732    if (ninterval < 1) then
     1733      ! The user has not specified shortwave albedo bands - assume
     1734      ! only one
    16111735      ninterval = 1
    1612       nvalue = 1
    1613     else
    1614       ! Check wavelength is monotonically increasing
    1615       do jband = 2,ninterval-1
    1616         if (wavelength_bound(jband) <= wavelength_bound(jband-1)) then
    1617           if (is_sw) then
    1618             write(nulerr, '(a,a)') '*** Error: wavelength bounds for shortwave albedo intervals ', &
    1619                  &  'must be monotonically increasing'
    1620           else
    1621             write(nulerr, '(a,a)') '*** Error: wavelength bounds for longwave emissivity intervals ', &
    1622                  &  'must be monotonically increasing'
    1623           end if
    1624           call radiation_abort()
    1625         end if
    1626       end do
    1627 
    1628       ! What is the maximum index, indicating the number of
    1629       ! albedo/emissivity values to expect?
    1630       nvalue = maxval(i_intervals(1:ninterval))
    1631      
    1632       if (do_nearest) then
    1633         ! Simpler nearest-neighbour mapping from band to
    1634         ! albedo/emissivity interval
    1635         if (allocated(i_mapping)) then
    1636           deallocate(i_mapping)
    1637         end if
    1638         allocate(i_mapping(nband))
    1639 
    1640         ! Loop over bands
    1641         do jband = 1,nband
    1642           ! Compute mid-point of band in wavenumber space (cm-1)
    1643           wavenumber_mid = 0.5_jprb * (wavenumber1(jband) &
    1644                &                     + wavenumber2(jband))
    1645           iinterval = 1
    1646           ! Convert wavelength (m) into wavenumber (cm-1) at the lower
    1647           ! bound of the albedo interval
    1648           wavenumber_bound = 0.01_jprb / wavelength_bound(iinterval)
    1649           ! Find the albedo interval that has the largest overlap with
    1650           ! the band; this approach assumes that the albedo intervals
    1651           ! are larger than the spectral bands
    1652           do while (wavenumber_bound >= wavenumber_mid &
    1653                &    .and. iinterval < ninterval)
    1654             iinterval = iinterval + 1
    1655             if (iinterval < ninterval) then
    1656               wavenumber_bound = 0.01_jprb / wavelength_bound(iinterval)
    1657             else
    1658               ! For the last interval there is no lower bound
    1659               wavenumber_bound = 0.0_jprb
    1660             end if
    1661           end do
    1662           ! Save the index of the band corresponding to the albedo
    1663           ! interval and move onto the next band
    1664           i_mapping(jband) = i_intervals(iinterval)
    1665         end do
    1666       else
    1667         ! More accurate weighting
    1668         if (allocated(weights)) then
    1669           deallocate(weights)
    1670         end if
    1671         allocate(weights(nvalue,nband))
    1672         weights(:,:) = 0.0_jprb
    1673        
    1674         ! Loop over bands
    1675         do jband = 1,nband
    1676           recip_dwavenumber = 1.0_jprb / (wavenumber2(jband) &
    1677                &                        - wavenumber1(jband))
    1678           ! Find the first overlapping albedo band
    1679           iinterval = 1
    1680           ! Convert wavelength (m) into wavenumber (cm-1) at the lower
    1681           ! bound (in wavenumber space) of the albedo/emissivty interval
    1682           wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval)
    1683           do while (wavenumber1_albedo >= wavenumber2(jband) &
    1684                &    .and. iinterval < ninterval)
    1685             iinterval = iinterval + 1
    1686             wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval)
    1687           end do
    1688          
    1689           wavenumber2_albedo = wavenumber2(jband)
    1690          
    1691           ! Add all overlapping bands
    1692           do while (wavenumber2_albedo > wavenumber1(jband) &
    1693                &  .and. iinterval <= ninterval)
    1694             weights(i_intervals(iinterval),jband) &
    1695                  &  = weights(i_intervals(iinterval),jband) &
    1696                  &  + recip_dwavenumber &
    1697                  &  * (min(wavenumber2_albedo,wavenumber2(jband)) &
    1698                  &   - max(wavenumber1_albedo,wavenumber1(jband)))
    1699             wavenumber2_albedo = wavenumber1_albedo
    1700             iinterval = iinterval + 1
    1701             if (iinterval < ninterval) then
    1702               wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval)
    1703             else
    1704               wavenumber1_albedo = 0.0_jprb
    1705             end if
    1706           end do
    1707         end do
    1708       end if
    1709     end if
    1710 
    1711     ! Define how many bands to use for reporting surface downwelling
    1712     ! fluxes for canopy radiation scheme
    1713     if (is_sw) then
     1736      this%i_sw_albedo_index(1) = 1
     1737      this%i_sw_albedo_index(2:) = 0
    17141738      if (this%use_canopy_full_spectrum_sw) then
    17151739        this%n_canopy_bands_sw = this%n_g_sw
    17161740      else
    1717         this%n_canopy_bands_sw = nvalue
     1741        this%n_canopy_bands_sw = 1
     1742      end if
     1743    else
     1744      if (this%use_canopy_full_spectrum_sw) then
     1745        this%n_canopy_bands_sw = this%n_g_sw
     1746      else
     1747        this%n_canopy_bands_sw = maxval(this%i_sw_albedo_index(1:ninterval))
     1748      end if
     1749    end if
     1750   
     1751    if (this%do_weighted_surface_mapping) then
     1752      call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(SolarReferenceTemperature, &
     1753           &  this%sw_albedo_wavelength_bound(1:ninterval-1), this%i_sw_albedo_index(1:ninterval), &
     1754           &  this%sw_albedo_weights, use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point))
     1755    else
     1756      ! Weight each wavenumber equally as in IFS Cycles 48 and earlier
     1757      call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(-1.0_jprb, &
     1758           &  this%sw_albedo_wavelength_bound(1:ninterval-1), this%i_sw_albedo_index(1:ninterval), &
     1759           &  this%sw_albedo_weights, use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point))
     1760    end if
     1761
     1762    ! Legacy method uses input band with largest weight
     1763    if (this%do_nearest_spectral_sw_albedo) then
     1764      allocate(this%i_albedo_from_band_sw(this%n_bands_sw))
     1765      this%i_albedo_from_band_sw = maxloc(this%sw_albedo_weights, dim=1)
     1766    end if
     1767
     1768    if (this%iverbosesetup >= 2) then
     1769      write(nulout, '(a)') 'Surface shortwave albedo'
     1770      if (.not. this%do_nearest_spectral_sw_albedo) then
     1771        call this%gas_optics_sw%spectral_def%print_mapping_from_bands(this%sw_albedo_weights, &
     1772             &       use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point))
     1773      else if (ninterval <= 1) then
     1774        write(nulout, '(a)') 'All shortwave bands will use the same albedo'
     1775      else
     1776        write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_albedo_from_band_sw), &
     1777             &  ' shortwave intervals to albedo intervals:'
     1778        do jband = 1,size(this%i_albedo_from_band_sw)
     1779          write(nulout,'(a,i0)',advance='no') ' ', this%i_albedo_from_band_sw(jband)
     1780        end do
     1781        write(nulout, '()')
     1782      end if
     1783    end if
     1784   
     1785  end subroutine consolidate_sw_albedo_intervals
     1786
     1787
     1788  !---------------------------------------------------------------------
     1789  ! Consolidate the surface longwave emissivity intervals with the
     1790  ! band/g-point intervals
     1791  subroutine consolidate_lw_emiss_intervals(this)
     1792
     1793    use radiation_io, only : nulout
     1794    use radiation_spectral_definition, only : TerrestrialReferenceTemperature
     1795
     1796    class(config_type),   intent(inout) :: this
     1797
     1798    integer :: ninterval, jint, jband
     1799
     1800    ! Count the number of albedo/emissivity intervals
     1801    ninterval = 0
     1802    do jint = 1,NMaxAlbedoIntervals
     1803      if (this%i_lw_emiss_index(jint) > 0) then
     1804        ninterval = jint
     1805      else
     1806        exit
     1807      end if
     1808    end do
     1809
     1810    if (ninterval < 1) then
     1811      ! The user has not specified longwave emissivity bands - assume
     1812      ! only one
     1813      ninterval = 1
     1814      this%i_lw_emiss_index(1) = 1
     1815      this%i_lw_emiss_index(2:) = 0
     1816      if (this%use_canopy_full_spectrum_sw) then
     1817        this%n_canopy_bands_lw = this%n_g_lw
     1818      else
     1819        this%n_canopy_bands_lw = 1
    17181820      end if
    17191821    else
    … …  
    17211823        this%n_canopy_bands_lw = this%n_g_lw
    17221824      else
    1723         this%n_canopy_bands_lw = nvalue
    1724       end if
     1825        this%n_canopy_bands_lw = maxval(this%i_lw_emiss_index(1:ninterval))
     1826      end if
     1827    end if
     1828
     1829    if (this%do_weighted_surface_mapping) then
     1830      call this%gas_optics_lw%spectral_def%calc_mapping_from_bands(TerrestrialReferenceTemperature, &
     1831           &  this%lw_emiss_wavelength_bound(1:ninterval-1), this%i_lw_emiss_index(1:ninterval), &
     1832           &  this%lw_emiss_weights, use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point))
     1833    else
     1834      ! Weight each wavenumber equally as in IFS Cycles 48 and earlier
     1835      call this%gas_optics_lw%spectral_def%calc_mapping_from_bands(-1.0_jprb, &
     1836           &  this%lw_emiss_wavelength_bound(1:ninterval-1), this%i_lw_emiss_index(1:ninterval), &
     1837           &  this%lw_emiss_weights, use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point))
     1838    end if
     1839
     1840    ! Legacy method uses input band with largest weight
     1841    if (this%do_nearest_spectral_lw_emiss) then
     1842      allocate(this%i_emiss_from_band_lw(this%n_bands_lw))
     1843      this%i_emiss_from_band_lw = maxloc(this%lw_emiss_weights, dim=1)
    17251844    end if
    17261845
    17271846    if (this%iverbosesetup >= 2) then
    1728       if (.not. do_nearest) then
    1729         if (is_sw) then
    1730           write(nulout, '(a,i0,a,i0,a)') 'Weighting of ', nvalue, ' albedo values in ', &
    1731              &  nband, ' shortwave bands (wavenumber ranges in cm-1):'
    1732         else
    1733           write(nulout, '(a,i0,a,i0,a)') 'Weighting of ', nvalue, ' emissivity values in ', &
    1734              &  nband, ' longwave bands (wavenumber ranges in cm-1):'
    1735         end if
    1736         do jband = 1,nband
    1737           write(nulout,'(i6,a,i6,a)',advance='no') nint(wavenumber1(jband)), ' to', &
    1738                &                        nint(wavenumber2(jband)), ':'
    1739           do iinterval = 1,nvalue
    1740             write(nulout,'(f5.2)',advance='no') weights(iinterval,jband)
    1741           end do
    1742           write(nulout, '()')
    1743         end do
     1847      write(nulout, '(a)') 'Surface longwave emissivity'
     1848      if (.not. this%do_nearest_spectral_lw_emiss) then
     1849        call this%gas_optics_lw%spectral_def%print_mapping_from_bands(this%lw_emiss_weights, &
     1850             &                          use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point))
    17441851      else if (ninterval <= 1) then
    1745         if (is_sw) then
    1746           write(nulout, '(a)') 'All shortwave bands will use the same albedo'
    1747         else
    1748           write(nulout, '(a)') 'All longwave bands will use the same emissivty'
    1749         end if
     1852        write(nulout, '(a)') 'All longwave bands will use the same emissivty'
    17501853      else
    1751         if (is_sw) then
    1752           write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', nband, &
    1753                &  ' shortwave bands to albedo intervals:'
    1754         else
    1755           write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', nband, &
    1756                &  ' longwave bands to emissivity intervals:'
    1757         end if
    1758         do jband = 1,nband
    1759           write(nulout,'(a,i0)',advance='no') ' ', i_mapping(jband)
     1854        write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_emiss_from_band_lw), &
     1855             &  ' longwave intervals to emissivity intervals:'
     1856        do jband = 1,size(this%i_emiss_from_band_lw)
     1857          write(nulout,'(a,i0)',advance='no') ' ', this%i_emiss_from_band_lw(jband)
    17601858        end do
    17611859        write(nulout, '()')
    … …  
    17631861    end if
    17641862
    1765   end subroutine consolidate_intervals
     1863  end subroutine consolidate_lw_emiss_intervals
    17661864
    17671865
    … …  
    18661964  end subroutine print_enum
    18671965
    1868 
    1869   !---------------------------------------------------------------------
    1870   ! Return .true. if 1D allocatable array "var" is out of physical
    1871   ! range specified by boundmin and boundmax, and issue a warning.
    1872   ! "do_fix" determines whether erroneous values are fixed to lie
    1873   ! within the physical range. To check only a subset of the array,
    1874   ! specify i1 and i2 for the range.
    1875   function out_of_bounds_1d(var, var_name, boundmin, boundmax, do_fix, i1, i2) result (is_bad)
    1876 
    1877     use radiation_io,     only : nulout
    1878 
    1879     real(jprb), allocatable, intent(inout) :: var(:)
    1880     character(len=*),        intent(in) :: var_name
    1881     real(jprb),              intent(in) :: boundmin, boundmax
    1882     logical,                 intent(in) :: do_fix
    1883     integer,       optional, intent(in) :: i1, i2
    1884 
    1885     logical                       :: is_bad
    1886 
    1887     real(jprb) :: varmin, varmax
    1888 
    1889     is_bad = .false.
    1890 
    1891     if (allocated(var)) then
    1892 
    1893       if (present(i1) .and. present(i2)) then
    1894         varmin = minval(var(i1:i2))
    1895         varmax = maxval(var(i1:i2))
    1896       else
    1897         varmin = minval(var)
    1898         varmax = maxval(var)
    1899       end if
    1900 
    1901       if (varmin < boundmin .or. varmax > boundmax) then
    1902         write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') &
    1903              &  '*** Warning: ', var_name, ' range', varmin, ' to', varmax, &
    1904              &  ' is out of physical range', boundmin, 'to', boundmax
    1905         is_bad = .true.
    1906         if (do_fix) then
    1907           if (present(i1) .and. present(i2)) then
    1908             var(i1:i2) = max(boundmin, min(boundmax, var(i1:i2)))
    1909           else
    1910             var = max(boundmin, min(boundmax, var))
    1911           end if
    1912           write(nulout,'(a)') ': corrected'
    1913         else
    1914           write(nulout,'(1x)')
    1915         end if
    1916       end if
    1917 
    1918     end if
    1919    
    1920   end function out_of_bounds_1d
    1921 
    1922 
    1923   !---------------------------------------------------------------------
    1924   ! Return .true. if 2D allocatable array "var" is out of physical
    1925   ! range specified by boundmin and boundmax, and issue a warning.  To
    1926   ! check only a subset of the array, specify i1 and i2 for the range
    1927   ! of the first dimension and j1 and j2 for the range of the second.
    1928   function out_of_bounds_2d(var, var_name, boundmin, boundmax, do_fix, &
    1929        &                    i1, i2, j1, j2) result (is_bad)
    1930 
    1931     use radiation_io,     only : nulout
    1932 
    1933     real(jprb), allocatable, intent(inout) :: var(:,:)
    1934     character(len=*),        intent(in) :: var_name
    1935     real(jprb),              intent(in) :: boundmin, boundmax
    1936     logical,                 intent(in) :: do_fix
    1937     integer,       optional, intent(in) :: i1, i2, j1, j2
    1938 
    1939     ! Local copies of indices
    1940     integer :: ii1, ii2, jj1, jj2
    1941 
    1942     logical                       :: is_bad
    1943 
    1944     real(jprb) :: varmin, varmax
    1945 
    1946     is_bad = .false.
    1947 
    1948     if (allocated(var)) then
    1949 
    1950       if (present(i1) .and. present(i2)) then
    1951         ii1 = i1
    1952         ii2 = i2
    1953       else
    1954         ii1 = lbound(var,1)
    1955         ii2 = ubound(var,1)
    1956       end if
    1957       if (present(j1) .and. present(j2)) then
    1958         jj1 = j1
    1959         jj2 = j2
    1960       else
    1961         jj1 = lbound(var,2)
    1962         jj2 = ubound(var,2)
    1963       end if
    1964       varmin = minval(var(ii1:ii2,jj1:jj2))
    1965       varmax = maxval(var(ii1:ii2,jj1:jj2))
    1966 
    1967       if (varmin < boundmin .or. varmax > boundmax) then
    1968         write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') &
    1969              &  '*** Warning: ', var_name, ' range', varmin, ' to', varmax,&
    1970              &  ' is out of physical range', boundmin, 'to', boundmax
    1971         is_bad = .true.
    1972         if (do_fix) then
    1973           var(ii1:ii2,jj1:jj2) = max(boundmin, min(boundmax, var(ii1:ii2,jj1:jj2)))
    1974           write(nulout,'(a)') ': corrected'
    1975         else
    1976           write(nulout,'(1x)')
    1977         end if
    1978       end if
    1979 
    1980     end if
    1981    
    1982   end function out_of_bounds_2d
    1983 
    1984 
    1985   !---------------------------------------------------------------------
    1986   ! Return .true. if 3D allocatable array "var" is out of physical
    1987   ! range specified by boundmin and boundmax, and issue a warning.  To
    1988   ! check only a subset of the array, specify i1 and i2 for the range
    1989   ! of the first dimension, j1 and j2 for the second and k1 and k2 for
    1990   ! the third.
    1991   function out_of_bounds_3d(var, var_name, boundmin, boundmax, do_fix, &
    1992        &                    i1, i2, j1, j2, k1, k2) result (is_bad)
    1993 
    1994     use radiation_io,     only : nulout
    1995 
    1996     real(jprb), allocatable, intent(inout) :: var(:,:,:)
    1997     character(len=*),        intent(in) :: var_name
    1998     real(jprb),              intent(in) :: boundmin, boundmax
    1999     logical,                 intent(in) :: do_fix
    2000     integer,       optional, intent(in) :: i1, i2, j1, j2, k1, k2
    2001 
    2002     ! Local copies of indices
    2003     integer :: ii1, ii2, jj1, jj2, kk1, kk2
    2004 
    2005     logical                       :: is_bad
    2006 
    2007     real(jprb) :: varmin, varmax
    2008 
    2009     is_bad = .false.
    2010 
    2011     if (allocated(var)) then
    2012 
    2013       if (present(i1) .and. present(i2)) then
    2014         ii1 = i1
    2015         ii2 = i2
    2016       else
    2017         ii1 = lbound(var,1)
    2018         ii2 = ubound(var,1)
    2019       end if
    2020       if (present(j1) .and. present(j2)) then
    2021         jj1 = j1
    2022         jj2 = j2
    2023       else
    2024         jj1 = lbound(var,2)
    2025         jj2 = ubound(var,2)
    2026       end if
    2027       if (present(k1) .and. present(k2)) then
    2028         kk1 = k1
    2029         kk2 = k2
    2030       else
    2031         kk1 = lbound(var,3)
    2032         kk2 = ubound(var,3)
    2033       end if
    2034       varmin = minval(var(ii1:ii2,jj1:jj2,kk1:kk2))
    2035       varmax = maxval(var(ii1:ii2,jj1:jj2,kk1:kk2))
    2036 
    2037       if (varmin < boundmin .or. varmax > boundmax) then
    2038         write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') &
    2039              &  '*** Warning: ', var_name, ' range', varmin, ' to', varmax,&
    2040              &  ' is out of physical range', boundmin, 'to', boundmax
    2041         is_bad = .true.
    2042         if (do_fix) then
    2043           var(ii1:ii2,jj1:jj2,kk1:kk2) = max(boundmin, min(boundmax, &
    2044                &                             var(ii1:ii2,jj1:jj2,kk1:kk2)))
    2045           write(nulout,'(a)') ': corrected'
    2046         else
    2047           write(nulout,'(1x)')
    2048         end if
    2049       end if
    2050 
    2051     end if
    2052    
    2053   end function out_of_bounds_3d
    2054 
    2055 
    20561966end module radiation_config
Note: See TracChangeset for help on using the changeset viewer.