Changeset 4853


Ignore:
Timestamp:
Mar 19, 2024, 3:34:21 PM (7 weeks ago)
Author:
idelkadi
Message:

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

  • version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
  • interface routines between lmdz and ecrad grouped in a new "lmdz" directory
  • double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
  • clean-up in the read.F routine: delete unitules arguments
  • modification of default gas model in namelist (default: ECCKD)
Location:
LMDZ6/trunk
Files:
8 added
7 deleted
37 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/DefLists/namelist_ecrad

    r4571 r4853  
    5252 use_aerosols           = false,          ! Include aerosols in radiation calculations?
    5353 n_aerosol_types         = 13,             
    54  !aerosol_optics_override_file_name = "aerosol_optics_lmdz.nc"
     54 !aerosol_optics_override_file_name = "aer_opt_LMDZ_ECCKD.nc",
     55 !aerosol_optics_override_file_name = "aer_opt_LMDZ_RRTMG.nc",
    5556 do_save_spectral_flux   = false,          ! Save spectral fluxes in output file?
    5657 do_save_gpoint_flux     = false,          ! Save fluxes per g-point in output file?
    57  gas_model_name          = "RRTMG-IFS",    ! Gas model: "Monochromatic", *"RRTMG-IFS"*, "RRTMG-PSRAD"
     58 gas_model_name          = "ECCKD",    ! Gas model: "Monochromatic", "RRTMG-IFS", "ECCKD"
    5859/

  • LMDZ6/trunk/bld.cfg.ecrad

    r4773 r4853  
    3535src::ext_src %EXT_SRC
    3636src::Ocean_skin %SRC_PATH/%PHYS/Ocean_skin
     37src::lmdz %RAD/lmdz
    3738src::radiation %RAD/radiation
    3839src::ifsrrtm %RAD/ifsrrtm

  • LMDZ6/trunk/libf/phylmd/ecrad/CHANGELOG

    r4773 r4853  
    1515        - Increased security value in single-precision SW
    1616          reflectance-transmittance calculation from 1e-12 to 1e-6
     17        - Added radiation/ecrad_config.h file to provide site-specific
     18          optimization options as preprocessor parameters
     19        - Replaced slow "sum" intrinsic function in McICA and Tripleclouds
     20          solvers with optimized versions for x86-64 and DWD's NEC
     21        - Enable different gas optics models for shortwave and longwave:
     22          namelist gas_model_name can still be used, or use
     23          sw_gas_model_name and lw_gas_model_name to specify separate models
     24        - Fix the general aerosol optics + RRTMG combination, so that
     25          optical properties correctly weighted by solar or terrestrial
     26          Planck function, rather than unweighted. This changes fluxes by up
     27          to 0.3 W m-2 in SW and 0.03 in LW, heating rates by up to 0.015
     28          K/d in SW and 0.0015 in LW. Reverts bug introduced between git
     29          commits a405cca (5 Dec 2022) and 7182230 (8 Dec 2022).
     30        - Added ifs/satur.F90, used in the IFS to compute relative
     31          humidity for aerosol hydration, and the ability to call it from
     32          ecrad_driver.F90, but note that this routine computes saturation
     33          with respect to ice at colder temperatures. If used, the impact on
     34          fluxes is up to around 1 W m-2 in the SW and 0.05 in the LW, and
     35          for heating rates 0.002 K/d in the SW and 0.003 in the LW.
    1736
    1837version 1.6.0 (27 April 2023)

  • LMDZ6/trunk/libf/phylmd/ecrad/driver/Makefile

    r4773 r4853  
    2020test_programs: $(TEST_PROGRAMS)
    2121
     22# Link ecrad executable; add "-lifs" if you want to use the "satur"
     23# routine in ecrad_driver.F90
    2224$(EXECUTABLE): $(OBJECTS) ../lib/*.a ecrad_driver.o
    2325        $(FC) $(FCFLAGS) ecrad_driver.o $(OBJECTS) $(LIBS) -o $(EXECUTABLE)

  • LMDZ6/trunk/libf/phylmd/ecrad/driver/ecrad_driver.F90

    r4773 r4853  
    5555  implicit none
    5656
     57  ! Uncomment this if you want to use the "satur" routine below
     58!#include "satur.intfb.h"
     59 
    5760  ! The NetCDF file containing the input profiles
    5861  type(netcdf_file)         :: file
    … …  
    269272
    270273  ! Compute saturation with respect to liquid (needed for aerosol
    271   ! hydration) call
     274  ! hydration) call...
    272275  call thermodynamics%calc_saturation_wrt_liquid(driver_config%istartcol,driver_config%iendcol)
    273276
     277  ! ...or alternatively use the "satur" function in the IFS (requires
     278  ! adding -lifs to the linker command line) but note that this
     279  ! computes saturation with respect to ice at colder temperatures,
     280  ! which is almost certainly incorrect
     281  !allocate(thermodynamics%h2o_sat_liq(ncol,nlev))
     282  !call satur(driver_config%istartcol, driver_config%iendcol, ncol, 1, nlev, .false., &
     283  !     0.5_jprb * (thermodynamics.pressure_hl(:,1:nlev)+thermodynamics.pressure_hl(:,2:nlev)), &
     284  !     0.5_jprb * (thermodynamics.temperature_hl(:,1:nlev)+thermodynamics.temperature_hl(:,2:nlev)), &
     285  !     thermodynamics%h2o_sat_liq, 2)
     286 
    274287  ! Check inputs are within physical bounds, printing message if not
    275288  is_out_of_bounds =     gas%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, &

  • LMDZ6/trunk/libf/phylmd/ecrad/ifs/Makefile

    r4773 r4853  
    22        radiation_scheme.F90 radiation_setup.F90 yoerdu.F90 \
    33        yomrip.F90 yoephy.F90 yoecld.F90 yoe_spectral_planck.F90 \
    4         cloud_overlap_decorr_len.F90 yoerad.F90
     4        cloud_overlap_decorr_len.F90 yoerad.F90 yoethf.F90 satur.F90
    55
    66OBJECTS := $(SOURCES:.F90=.o)
    … …  
    3131cos_sza.o ice_effective_radius.o liquid_effective_radius.o radiation_scheme.o radiation_setup.o: yoerad.o
    3232yoerad.o: yoe_spectral_planck.o
     33satur.o: yoethf.o

  • LMDZ6/trunk/libf/phylmd/ecrad/ifs/yoephy.F90

    r4773 r4853  
    1111
    1212USE PARKIND1, ONLY : JPRB, JPIM
    13 USE ISO_C_BINDING
     13!USE ISO_C_BINDING
    1414
    1515IMPLICIT NONE

  • LMDZ6/trunk/libf/phylmd/ecrad/ifsaux/yomcst.F90

    r4773 r4853  
    7171REAL(KIND=JPRB), PARAMETER :: RMCCL4 = 153.823_JPRB
    7272
     73REAL(KIND=JPRB), PARAMETER :: RCPD  = 3.5_JPRB*RD
     74REAL(KIND=JPRB), PARAMETER :: RLMLT = RLSTT-RLVTT
     75
    7376END MODULE YOMCST

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/Makefile

    r4773 r4853  
    111111
    112112radiation_general_cloud_optics.o: radiation_config.o radiation_cloud.o radiation_thermodynamics.o radiation_constants.o
     113
     114*.o: ecrad_config.h

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90

    r4773 r4853  
    1919!   2022-11-22  P. Ukkonen / R. Hogan  Optimizations to enhance vectorization
    2020
     21#include "ecrad_config.h"
     22
    2123module radiation_aerosol_optics
    2224
    … …  
    9698    use parkind1,                      only : jprb
    9799    use yomhook,                       only : lhook, dr_hook, jphook
     100#ifdef EASY_NETCDF_READ_MPI
     101    use easy_netcdf_read_mpi,          only : netcdf_file
     102#else
    98103    use easy_netcdf,                   only : netcdf_file
     104#endif
    99105    use radiation_config,              only : config_type
    100106    use radiation_aerosol_optics_data, only : aerosol_optics_type
    … …  
    341347    use parkind1,                      only : jprb
    342348    use yomhook,                       only : lhook, dr_hook, jphook
     349#ifdef EASY_NETCDF_READ_MPI
     350    use easy_netcdf_read_mpi,          only : netcdf_file
     351#else
    343352    use easy_netcdf,                   only : netcdf_file
     353#endif
    344354    use radiation_config,              only : config_type
    345355    use radiation_aerosol_optics_data, only : aerosol_optics_type

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90

    r4773 r4853  
    1717!   2018-04-20  A. Bozzo  Read optical properties at selected wavelengths
    1818
     19#include "ecrad_config.h"
    1920
    2021module radiation_aerosol_optics_data
    … …  
    157158
    158159    use yomhook,              only : lhook, dr_hook, jphook
     160#ifdef EASY_NETCDF_READ_MPI
     161    use easy_netcdf_read_mpi, only : netcdf_file
     162#else
    159163    use easy_netcdf,          only : netcdf_file
     164#endif
    160165    use radiation_io,         only : nulerr, radiation_abort
    161166
    … …  
    401406  subroutine save_aerosol_optics(this, file_name, iverbose)
    402407
    403     use yomhook,     only : lhook, dr_hook, jphook
    404     use easy_netcdf, only : netcdf_file
     408    use yomhook,              only : lhook, dr_hook, jphook
     409    use easy_netcdf,          only : netcdf_file
    405410
    406411    class(aerosol_optics_type), intent(inout) :: this

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_description.F90

    r4773 r4853  
    1313! Email:   r.j.hogan@ecmwf.int
    1414!
     15
     16#include "ecrad_config.h"
    1517
    1618module radiation_aerosol_optics_description
    … …  
    6062    logical, allocatable :: is_preferred_philic(:)
    6163
     64    ! Verbosity level
     65    integer :: iverbose
     66   
    6267  contains
    6368    procedure :: read
    … …  
    7580
    7681    use yomhook,              only : lhook, dr_hook, jphook
     82#ifdef EASY_NETCDF_READ_MPI
     83    use easy_netcdf_read_mpi, only : netcdf_file
     84#else
    7785    use easy_netcdf,          only : netcdf_file
     86#endif
    7887
    7988    class(aerosol_optics_description_type), intent(inout) :: this
    … …  
    8493    type(netcdf_file)  :: file
    8594
    86     real(jphook) :: hook_handle
     95    real(jphook)       :: hook_handle
    8796
    8897    if (lhook) call dr_hook('radiation_aerosol_optics_description:load',0,hook_handle)
    … …  
    108117    call file%close()
    109118
     119    if (present(iverbose)) then
     120      this%iverbose = iverbose
     121    else
     122      this%iverbose = 3
     123    end if
     124   
    110125    if (lhook) call dr_hook('radiation_aerosol_optics_description:load',1,hook_handle)
    111126
    … …  
    124139
    125140    use yomhook,              only : lhook, dr_hook, jphook
    126 
     141    use radiation_io,         only : nulout, nulerr, radiation_abort
     142   
    127143    class(aerosol_optics_description_type), intent(inout) :: this
    128144    character(len=2), intent(in) :: code_str
    … …  
    132148    integer :: ja
    133149
    134     real(jphook) :: hook_handle
     150    logical :: is_found, is_philic, is_phobic
     151
     152    real(jphook)         :: hook_handle
    135153
    136154    if (lhook) call dr_hook('radiation_aerosol_optics_description:preferred_optical_model',0,hook_handle)
    … …  
    138156    ! Check for empty string
    139157    if (optical_model_str == ' ') then
     158      if (lhook) call dr_hook('radiation_aerosol_optics_description:preferred_optical_model',1,hook_handle)
    140159      return
    141160    end if
     161
     162    is_found  = .false.
     163    is_philic = .false.
     164    is_phobic = .false.
    142165   
    143166    ! Loop over hydrophilic types
    … …  
    145168      ! Check if we have a match
    146169      if (to_string(this%code_philic(:,ja)) == code_str &
    147            &  .and. to_string(this%optical_model_philic(1:len(optical_model_str),ja)) &
     170           &  .and. trim(to_string(this%optical_model_philic(:,ja))) &
    148171           &          == optical_model_str) then
    149172        this%is_preferred_philic(ja) = .true.
     173        is_found  = .true.
     174        is_philic = .true.
    150175      end if
    151176    end do
    … …  
    153178    do ja = 1,size(this%bin_phobic)
    154179      if (to_string(this%code_phobic(:,ja)) == code_str &
    155            &  .and. to_string(this%optical_model_phobic(1:len(optical_model_str),ja)) &
     180           &  .and. trim(to_string(this%optical_model_phobic(:,ja))) &
    156181           &          == optical_model_str) then
    157182        this%is_preferred_phobic(ja) = .true.
     183        is_found  = .true.
     184        is_phobic = .true.
    158185      end if
    159186    end do
    160187
     188    if (.not. is_found) then
     189      write(nulerr,'(a,a2,a,a,a)') '*** Error: Preferred "', code_str ,'" aerosol optical model "', &
     190           &  trim(optical_model_str), '" not found in file'
     191      call radiation_abort()
     192    else if (this%iverbose > 2) then
     193      write(nulout,'(a,a2,a,a,a)',advance='no') 'Preferred "', code_str, '" aerosol optical model set to "', &
     194           &  trim(optical_model_str), '" ('
     195      if (is_phobic) then
     196        write(nulout,'(a)',advance='no') ' hydrophobic'
     197      end if
     198      if (is_philic) then
     199        write(nulout,'(a)',advance='no') ' hydrophilic'
     200      end if
     201      write(nulout,'(a)') ' )'
     202    end if
     203   
    161204    if (lhook) call dr_hook('radiation_aerosol_optics_description:preferred_optical_model',1,hook_handle)
    162205
    163206  end subroutine preferred_optical_model
    164207
     208 
    165209  !---------------------------------------------------------------------
    166210  ! Return the index to the aerosol optical properties corresponding
    … …  
    179223   
    180224    use yomhook,              only : lhook, dr_hook, jphook
    181     use easy_netcdf,          only : netcdf_file
    182225    use radiation_io,         only : nulout
    183226
    … …  
    234277          end if
    235278          if (present(optical_model_str)) then
    236             if (to_string(this%optical_model_philic(1:len(optical_model_str),ja)) &
     279            if (trim(to_string(this%optical_model_philic(:,ja))) &
    237280                 &  == optical_model_str) then
    238281              ! Requested optical model matches
    … …  
    285328          end if
    286329          if (present(optical_model_str)) then
    287             if (to_string(this%optical_model_phobic(1:len(optical_model_str),ja)) &
     330            if (trim(to_string(this%optical_model_phobic(:,ja))) &
    288331                 &  == optical_model_str) then
    289332              ! Requested optical model matches
    … …  
    315358
    316359    if (is_ambiguous) then
    317       write(nulout,'(a,a2,a,l1,a)') 'Warning: get_index("', code_str, '",', lhydrophilic, &
     360      write(nulout,'(a,a2,a,l,a)') 'Warning: radiation_aerosol_optics_description:get_index("', &
     361           &  code_str, '",', lhydrophilic, &
    318362           &  ',...) does not unambiguously identify an aerosol optical property index'
    319363    end if
    320 
     364   
    321365    if (lhook) call dr_hook('radiation_aerosol_optics_description:get_index',1,hook_handle)
    322366
    … …  
    325369  !---------------------------------------------------------------------
    326370  ! Utility function to convert an array of single characters to a
    327   ! character string (yes Fortran's string handling is a bit rubbish)
     371  ! character string (yes Fortran's string handling is a bit
     372  ! rubbish). We set NULL characters (ASCII code 0) returned from the
     373  ! NetCDF library to spaces, so that TRIM can remove them.
    328374  pure function to_string(arr) result(str)
    329375    character, intent(in)  :: arr(:)
    … …  
    331377    integer :: jc
    332378    do jc = 1,size(arr)
    333       str(jc:jc) = arr(jc)
     379      if (ichar(arr(jc)) == 0) then
     380        ! Replace NULL character with a space
     381        str(jc:jc) = ' '
     382      else
     383        str(jc:jc) = arr(jc)
     384      end if
    334385    end do
    335386  end function to_string

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_cloud.F90

    r4773 r4853  
    6464    ! gridbox area for use in representing 3D effects. This variable
    6565    ! is dimensioned (ncol,nlev).
    66     real(jprb), allocatable,  dimension(:,:) :: inv_cloud_effective_size ! m-1
     66    real(jprb), allocatable, dimension(:,:) :: inv_cloud_effective_size ! m-1
    6767
    6868    ! Similarly for the in-cloud heterogeneities, used to compute the
    … …  
    606606
    607607    use yomhook,                  only : lhook, dr_hook, jphook
    608 !    USE mod_phys_lmdz_para
    609608
    610609    class(cloud_type), intent(inout) :: this

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_cloud_cover.F90

    r4773 r4853  
    1919!   2020-10-07  R. Hogan  Ensure iobj1 initialized in case of alpha_obj==0
    2020
     21#include "ecrad_config.h"
     22
    2123module radiation_cloud_cover
    2224
    … …  
    254256           &  * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1))
    255257! Added for DWD (2020)
    256 #ifdef __SX__
     258#ifdef DWD_VECTOR_OPTIMIZATIONS
    257259    end do
    258260    do jlev = 1,nlev-1

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_cloud_generator.F90

    r4773 r4853  
    540540    use radiation_pdf_sampler, only : pdf_sampler_type
    541541    implicit none
    542 !#if defined(__GFORTRAN__) || defined(__PGI) || defined(__NEC__)
    543 !#else
    544 !    !$omp declare simd(sample_from_pdf_simd) uniform(this) &
    545 !    !$omp linear(ref(fsd)) linear(ref(cdf))
    546 !#endif
     542#if defined(__GFORTRAN__) || defined(__PGI) || defined(__NEC__) || defined(__INTEL_LLVM_COMPILER)
     543#else
     544    !!$omp declare simd(sample_from_pdf_simd) uniform(this) &
     545    !!$omp linear(ref(fsd)) linear(ref(cdf))
     546#endif
    547547    type(pdf_sampler_type), intent(in)  :: this
    548548

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_cloud_optics_data.F90

    r4773 r4853  
    1313! Email:   r.j.hogan@ecmwf.int
    1414!
     15
     16#include "ecrad_config.h"
    1517
    1618module radiation_cloud_optics_data
    … …  
    4850  subroutine setup_cloud_optics(this, liq_file_name, ice_file_name, iverbose)
    4951   
    50     use yomhook,  only : lhook, dr_hook, jphook
    51     use easy_netcdf, only : netcdf_file
     52    use yomhook,              only : lhook, dr_hook, jphook
     53#ifdef EASY_NETCDF_READ_MPI
     54    use easy_netcdf_read_mpi, only : netcdf_file
     55#else
     56    use easy_netcdf,          only : netcdf_file
     57#endif
    5258
    5359    class(cloud_optics_type), intent(inout) :: this

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

    r4773 r4853  
    3434!
    3535
     36#include "ecrad_config.h"
     37
    3638module radiation_config
    3739
    … …  
    202204    ! more random numbers being generated?  This is the default on NEC
    203205    ! SX.
    204 #ifdef __SX__
     206#ifdef DWD_VECTOR_OPTIMIZATIONS
    205207    logical :: use_vectorizable_generator = .true.
    206208#else
    … …  
    276278
    277279    ! Codes describing the gas model
    278     integer :: i_gas_model = IGasModelIFSRRTMG
     280    integer :: i_gas_model_sw = IGasModelIFSRRTMG
     281    integer :: i_gas_model_lw = IGasModelIFSRRTMG
    279282
    280283    ! Optics if i_gas_model==IGasModelMonochromatic.
    … …  
    702705    character(511) :: ssi_override_file_name
    703706    character(63)  :: liquid_model_name, ice_model_name, gas_model_name
     707    character(63)  :: sw_gas_model_name, lw_gas_model_name
    704708    character(63)  :: sw_solver_name, lw_solver_name, overlap_scheme_name
    705709    character(63)  :: sw_entrapment_name, sw_encroachment_name, cloud_pdf_shape_name
    … …  
    709713         &     .false.,.false.,.false.,.false.,.false.,.false.]
    710714    integer :: i_aerosol_type_map(NMaxAerosolTypes) ! More than 256 is an error
    711 
     715   
    712716    logical :: do_nearest_spectral_sw_albedo
    713717    logical :: do_nearest_spectral_lw_emiss
    … …  
    716720    integer :: i_sw_albedo_index(NMaxAlbedoIntervals)
    717721    integer :: i_lw_emiss_index (NMaxAlbedoIntervals)
     722    integer :: i_gas_model
    718723
    719724    integer :: iunit ! Unit number of namelist file
    … …  
    726731         &  do_surface_sw_spectral_flux, do_lw_derivatives, do_toa_spectral_flux, &
    727732         &  do_lw_aerosol_scattering, do_lw_cloud_scattering, &
    728          &  n_regions, directory_name, gas_model_name, &
     733         &  n_regions, directory_name, gas_model_name, sw_gas_model_name, lw_gas_model_name, &
    729734         &  ice_optics_override_file_name, liq_optics_override_file_name, &
    730735         &  aerosol_optics_override_file_name, cloud_pdf_override_file_name, &
    … …  
    815820    encroachment_scaling = -1.0_jprb
    816821    gas_model_name = '' !DefaultGasModelName
     822    sw_gas_model_name = '' !DefaultGasModelName
     823    lw_gas_model_name = '' !DefaultGasModelName
    817824    liquid_model_name = '' !DefaultLiquidModelName
    818825    ice_model_name = '' !DefaultIceModelName
    … …  
    10141021         &            'ice_model_name', this%i_ice_model)
    10151022
    1016     ! Determine gas optics model
     1023    ! Determine gas optics model(s) - firstly try the generic gas_model_name
     1024    i_gas_model = -1
    10171025    call get_enum_code(gas_model_name, GasModelName, &
    1018          &            'gas_model_name', this%i_gas_model)
    1019 
     1026         &            'gas_model_name', i_gas_model)
     1027    if (i_gas_model > -1) then
     1028      this%i_gas_model_sw = i_gas_model
     1029      this%i_gas_model_lw = i_gas_model
     1030    end if
     1031    ! ...then the band-specific values
     1032    call get_enum_code(sw_gas_model_name, GasModelName, &
     1033         &            'sw_gas_model_name', this%i_gas_model_sw)
     1034    call get_enum_code(lw_gas_model_name, GasModelName, &
     1035         &            'lw_gas_model_name', this%i_gas_model_lw)
     1036   
    10201037    ! Determine solvers
    10211038    call get_enum_code(sw_solver_name, SolverName, &
    … …  
    10551072    end if
    10561073
    1057     if (this%i_gas_model == IGasModelIFSRRTMG &
    1058          & .and. (this%use_general_cloud_optics &
    1059          &        .or. this%use_general_aerosol_optics)) then
    1060       if (this%do_sw .and. this%do_cloud_aerosol_per_sw_g_point) then
     1074    if (this%use_general_cloud_optics .or. this%use_general_aerosol_optics) then
     1075      if (this%do_sw .and. this%do_cloud_aerosol_per_sw_g_point &
     1076           &  .and. this%i_gas_model_sw == IGasModelIFSRRTMG) then
    10611077        write(nulout,'(a)') 'Warning: RRTMG SW only supports cloud/aerosol/surface optical properties per band, not per g-point'
    10621078        this%do_cloud_aerosol_per_sw_g_point = .false.
    10631079      end if
    1064       if (this%do_lw .and. this%do_cloud_aerosol_per_lw_g_point) then
     1080      if (this%do_lw .and. this%do_cloud_aerosol_per_lw_g_point &
     1081           &  .and. this%i_gas_model_lw == IGasModelIFSRRTMG) then
    10651082        write(nulout,'(a)') 'Warning: RRTMG LW only supports cloud/aerosol/surface optical properties per band, not per g-point'
    10661083        this%do_cloud_aerosol_per_lw_g_point = .false.
    … …  
    11281145
    11291146    ! If ecCKD gas optics model is being used set relevant file names
    1130     if (this%i_gas_model == IGasModelECCKD) then
     1147    if (this%i_gas_model_sw == IGasModelECCKD .or. this%i_gas_model_lw == IGasModelECCKD) then
    11311148
    11321149      ! This gas optics model usually used with general cloud and
    … …  
    11381155        write(nulout,'(a)') 'Warning: ecCKD gas optics model usually used with general aerosol optics'
    11391156      end if
     1157
     1158    end if
     1159
     1160    if (this%i_gas_model_sw == IGasModelECCKD) then
    11401161
    11411162      if (len_trim(this%gas_optics_sw_override_file_name) > 0) then
    … …  
    11531174      end if
    11541175
     1176    end if
     1177
     1178    if (this%i_gas_model_lw == IGasModelECCKD) then
     1179
    11551180      if (len_trim(this%gas_optics_lw_override_file_name) > 0) then
    11561181        if (this%gas_optics_lw_override_file_name(1:1) == '/') then
    … …  
    11701195
    11711196    if (this%use_spectral_solar_cycle) then
    1172       if (this%i_gas_model /= IGasModelECCKD) then
     1197      if (this%i_gas_model_sw /= IGasModelECCKD) then
    11731198        write(nulerr,'(a)') '*** Error: solar cycle only available with ecCKD gas optics model'
    11741199        call radiation_abort('Radiation configuration error')
    … …  
    12781303    end if
    12791304
    1280     ! In the monochromatic case we need to override the liquid, ice
    1281     ! and aerosol models to ensure compatibility
    1282     if (this%i_gas_model == IGasModelMonochromatic) then
     1305    if (this%i_gas_model_sw == IGasModelMonochromatic .or. this%i_gas_model_lw == IGasModelMonochromatic) then
     1306
     1307      if (this%i_gas_model_sw /= this%i_gas_model_lw) then
     1308        write(nulerr,'(a,i0)') '*** Error: Monochromatic gas optics model must be used in shortwave and longwave'
     1309        call radiation_abort('Radiation configuration error')
     1310      end if
     1311   
     1312      ! In the monochromatic case we need to override the liquid, ice
     1313      ! and aerosol models to ensure compatibility
    12831314      this%i_liq_model = ILiquidModelMonochromatic
    12841315      this%i_ice_model = IIceModelMonochromatic
    12851316      this%use_aerosols = .false.
     1317     
    12861318    end if
    12871319
    … …  
    13921424      call print_logical('  Saving spectral flux profiles', &
    13931425           &   'do_save_spectral_flux', this%do_save_spectral_flux)
    1394       call print_enum('  Gas model is', GasModelName, 'i_gas_model', &
    1395            &          this%i_gas_model)
     1426      call print_enum('  Shortwave gas model is', GasModelName, 'i_gas_model_sw', &
     1427           &          this%i_gas_model_sw)
     1428      call print_enum('  Longwave gas model is', GasModelName, 'i_gas_model_lw', &
     1429           &          this%i_gas_model_lw)
    13961430      call print_logical('  Aerosols are', 'use_aerosols', this%use_aerosols)
    13971431      if (this%use_aerosols) then
    … …  
    14811515             &          'i_solver_sw', this%i_solver_sw)
    14821516       
    1483         if (this%i_gas_model == IGasModelMonochromatic) then
     1517        if (this%i_gas_model_sw == IGasModelMonochromatic) then
    14841518          call print_real('  Shortwave atmospheric optical depth', &
    14851519               &   'mono_sw_total_od', this%mono_sw_total_od)
    … …  
    15021536             &          this%i_solver_lw)
    15031537
    1504         if (this%i_gas_model == IGasModelMonochromatic) then
     1538        if (this%i_gas_model_lw == IGasModelMonochromatic) then
    15051539          if (this%mono_lw_wavelength > 0.0_jprb) then
    15061540            call print_real('  Longwave effective wavelength (m)', &

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ecckd.F90

    r4773 r4853  
    1515!
    1616
     17#include "ecrad_config.h"
     18
    1719module radiation_ecckd
    1820
    … …  
    105107! Vectorized version of the optical depth look-up performs better on
    106108! NEC, but slower on x86
    107 #ifdef __SX__
     109#ifdef DWD_VECTOR_OPTIMIZATIONS
    108110    procedure :: calc_optical_depth => calc_optical_depth_ckd_model_vec
    109111#else
    … …  
    125127  subroutine read_ckd_model(this, filename, iverbose)
    126128
    127     use easy_netcdf,  only : netcdf_file
     129#ifdef EASY_NETCDF_READ_MPI
     130    use easy_netcdf_read_mpi, only : netcdf_file
     131#else
     132    use easy_netcdf,          only : netcdf_file
     133#endif
    128134    !use radiation_io, only : nulerr, radiation_abort
    129     use yomhook,      only : lhook, dr_hook, jphook
     135    use yomhook,              only : lhook, dr_hook, jphook
    130136
    131137    class(ckd_model_type), intent(inout) :: this
    … …  
    166172    this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1
    167173    call file%get('temperature', temperature_full)
    168     ! AI oct 2023 ajout pour le double appel de ecrad   
    169     if (allocated(this%temperature1)) deallocate(this%temperature1)
    170174    allocate(this%temperature1(this%npress));
    171175    this%temperature1 = temperature_full(:,1)
    … …  
    200204    ! Read gases
    201205    call file%get('n_gases', this%ngas)
    202     if (allocated(this%single_gas)) deallocate(this%single_gas)
    203206    allocate(this%single_gas(this%ngas))
    204207    call file%get_global_attribute('constituent_id', constituent_id)
    … …  
    292295  subroutine read_spectral_solar_cycle(this, filename, iverbose, use_updated_solar_spectrum)
    293296
    294     use easy_netcdf,  only : netcdf_file
    295     use radiation_io, only : nulout, nulerr, radiation_abort
    296     use yomhook,      only : lhook, dr_hook, jphook
     297#ifdef EASY_NETCDF_READ_MPI
     298    use easy_netcdf_read_mpi, only : netcdf_file
     299#else
     300    use easy_netcdf,          only : netcdf_file
     301#endif
     302    use radiation_io,         only : nulout, nulerr, radiation_abort
     303    use yomhook,              only : lhook, dr_hook, jphook
    297304
    298305    ! Reference total solar irradiance (W m-2)
    … …  
    450457  subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
    451458       &  pressure_hl, temperature_fl, mole_fraction_fl, &
    452        &  optical_depth_fl, rayleigh_od_fl)
     459       &  optical_depth_fl, rayleigh_od_fl, concentration_scaling)
    453460
    454461    use yomhook,             only : lhook, dr_hook, jphook
    … …  
    466473    ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
    467474    real(jprb),            intent(in)  :: mole_fraction_fl(ncol,nlev,nmaxgas)
     475    ! Optional concentration scaling of each gas
     476    real(jprb), optional,  intent(in)  :: concentration_scaling(nmaxgas)
    468477   
    469478    ! Output variables
    … …  
    486495
    487496    real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1
     497    real(jprb) :: scaling
    488498
    489499    ! Indices and weights in temperature, pressure and concentration interpolation
    … …  
    547557            multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode)
    548558
     559            if (present(concentration_scaling)) then
     560              multiplier = multiplier * concentration_scaling(igascode)
     561            end if
     562           
    549563            do jlev = 1,nlev
    550564              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    … …  
    557571          case (IConcDependenceRelativeLinear)
    558572            molar_abs => this%single_gas(jgas)%molar_abs
    559             multiplier = simple_multiplier  * (mole_fraction_fl(jcol,:,igascode) &
    560                  &                            - single_gas%reference_mole_frac)
     573
     574            if (present(concentration_scaling)) then
     575              multiplier = simple_multiplier &
     576                   &  * (mole_fraction_fl(jcol,:,igascode)*concentration_scaling(igascode) &
     577                   &     - single_gas%reference_mole_frac)
     578            else
     579              multiplier = simple_multiplier  * (mole_fraction_fl(jcol,:,igascode) &
     580                   &                            - single_gas%reference_mole_frac)
     581            end if
     582           
    561583            do jlev = 1,nlev
    562584              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    … …  
    579601
    580602          case (IConcDependenceLUT)
     603
     604            if (present(concentration_scaling)) then
     605              scaling = concentration_scaling(igascode)
     606            else
     607              scaling = 1.0_jprb
     608            end if
     609           
    581610            ! Logarithmic interpolation in concentration space
    582611            molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
    … …  
    584613            do jlev = 1,nlev
    585614              ! Take care of mole_fraction == 0
    586               log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
     615              log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode)*scaling, mole_frac1))
    587616              cindex1  = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
    588617              cindex1  = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
    … …  
    601630            do jlev = 1,nlev
    602631              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    603                    &  + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( &
     632                   &  + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode) * scaling) * ( &
    604633                   &      (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) &
    605634                   &     +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) &

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ecckd_gas.F90

    r4773 r4853  
    1414! License: see the COPYING file for details
    1515!
     16
     17#include "ecrad_config.h"
    1618
    1719module radiation_ecckd_gas
    … …  
    8284  subroutine read_ckd_gas(this, file, gas_name, i_gas_code)
    8385
    84     use easy_netcdf, only : netcdf_file
     86#ifdef EASY_NETCDF_READ_MPI
     87    use easy_netcdf_read_mpi, only : netcdf_file
     88#else
     89    use easy_netcdf,          only : netcdf_file
     90#endif
    8591
    8692    class(ckd_gas_type), intent(inout) :: this

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90

    r4773 r4853  
    3939    if (lhook) call dr_hook('radiation_ecckd_interface:setup_gas_optics',0,hook_handle)
    4040
    41     if (config%do_sw) then
     41    if (config%do_sw .and. config%i_gas_model_sw == IGasModelECCKD) then
    4242
    4343      ! Read shortwave ecCKD gas optics NetCDF file
    … …  
    5757      end if
    5858
    59 !      if (allocated(config%i_band_from_g_sw)) deallocate(config%i_band_from_g_sw)
    60 !      allocate(config%i_band_from_g_sw(config%n_g_sw))
    61 !      if (allocated(config%i_band_from_reordered_g_sw)) deallocate(config%i_band_from_reordered_g_sw)
    62 !      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
    63 !      if (allocated(config%i_g_from_reordered_g_sw)) deallocate(config%i_g_from_reordered_g_sw)
    64 !      allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
    65       if (.not.allocated(config%i_band_from_g_sw)) &
    66               allocate(config%i_band_from_g_sw(config%n_g_sw))
    67       if (.not.allocated(config%i_band_from_reordered_g_sw)) &
    68               allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
    69       if (.not.allocated(config%i_g_from_reordered_g_sw)) &
    70               allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
     59      allocate(config%i_band_from_g_sw          (config%n_g_sw))
     60      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
     61      allocate(config%i_g_from_reordered_g_sw   (config%n_g_sw))
    7162       
    7263      if (config%do_cloud_aerosol_per_sw_g_point) then
    … …  
    9384    end if
    9485
    95     if (config%do_lw) then
     86    if (config%do_lw .and. config%i_gas_model_lw == IGasModelECCKD) then
    9687
    9788      ! Read longwave ecCKD gas optics NetCDF file
    … …  
    111102      end if
    112103
    113 !      if (allocated(config%i_band_from_g_lw)) deallocate(config%i_band_from_g_lw)
    114 !      allocate(config%i_band_from_g_lw          (config%n_g_lw))
    115 !      if (allocated(config%i_band_from_reordered_g_lw)) deallocate(config%i_band_from_reordered_g_lw)
    116 !      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
    117 !      if (allocated(config%i_g_from_reordered_g_lw)) deallocate(config%i_g_from_reordered_g_lw)
    118 !      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
    119       if (.not.allocated(config%i_band_from_g_lw)) &
    120          allocate(config%i_band_from_g_lw          (config%n_g_lw))
    121       if (.not.allocated(config%i_band_from_reordered_g_lw)) &
    122          allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
    123       if (.not.allocated(config%i_g_from_reordered_g_lw)) &
    124          allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
    125 
     104      allocate(config%i_band_from_g_lw          (config%n_g_lw))
     105      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
     106      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
    126107
    127108      if (config%do_cloud_aerosol_per_lw_g_point) then
    … …  
    199180    use yomhook,  only : lhook, dr_hook, jphook
    200181
    201     use radiation_config,         only : config_type
     182    use radiation_config,         only : config_type, IGasModelECCKD
    202183    use radiation_thermodynamics, only : thermodynamics_type
    203184    use radiation_single_level,   only : single_level_type
    … …  
    242223    real(jprb) :: temperature_fl(istartcol:iendcol,nlev)
    243224
    244     integer :: jcol
     225    real(jprb) :: concentration_scaling(NMaxGases)
     226   
     227    logical :: is_volume_mixing_ratio
     228   
     229    integer :: jcol, jlev, jg
    245230
    246231    real(jphook) :: hook_handle
    … …  
    259244         &  / (thermodynamics%pressure_hl(istartcol:iendcol,1:nlev) &
    260245         &    +thermodynamics%pressure_hl(istartcol:iendcol,2:nlev+1))
    261  
    262     if (config%do_sw) then
    263 
    264       call config%gas_optics_sw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
    265            &  NMaxGases, thermodynamics%pressure_hl, &
    266            &  temperature_fl, &
    267            &  gas%mixing_ratio, &
    268 !           &  reshape(gas%mixing_ratio(istartcol:iendcol,:,:), &
    269 !           &          [nlev,iendcol-istartcol+1,NMaxGases],order=[2,1,3]), &
    270            &  od_sw, rayleigh_od_fl=ssa_sw)
     246
     247    ! Check that the gas concentrations are stored in volume mixing
     248    ! ratio with no scaling; if not, return a vector of scalings
     249    call gas%assert_units(IVolumeMixingRatio, scale_factor=1.0_jprb, &
     250         &                istatus=is_volume_mixing_ratio)
     251    if (.not. is_volume_mixing_ratio) then
     252      call gas%get_scaling(IVolumeMixingRatio, concentration_scaling)
     253    else
     254      concentration_scaling = 1.0_jprb
     255    end if
     256   
     257    if (config%do_sw .and. config%i_gas_model_sw == IGasModelECCKD) then
     258
     259      if (is_volume_mixing_ratio) then
     260        call config%gas_optics_sw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
     261             &  NMaxGases, thermodynamics%pressure_hl, &
     262             &  temperature_fl, gas%mixing_ratio, &
     263             &  od_sw, rayleigh_od_fl=ssa_sw)
     264      else
     265        call config%gas_optics_sw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
     266             &  NMaxGases, thermodynamics%pressure_hl, &
     267             &  temperature_fl, gas%mixing_ratio, &
     268             &  od_sw, rayleigh_od_fl=ssa_sw, concentration_scaling=concentration_scaling)
     269      end if
     270     
    271271      ! At this point od_sw = absorption optical depth and ssa_sw =
    272272      ! rayleigh optical depth: convert to total optical depth and
    273273      ! single-scattering albedo
    274       od_sw = od_sw + ssa_sw
    275       ssa_sw = ssa_sw / od_sw
     274      do jcol = istartcol,iendcol
     275        do jlev = 1, nlev
     276          do jg = 1, config%n_g_sw
     277            od_sw(jg,jlev,jcol)  = od_sw(jg,jlev,jcol) + ssa_sw(jg,jlev,jcol)
     278            ssa_sw(jg,jlev,jcol) = ssa_sw(jg,jlev,jcol) / od_sw(jg,jlev,jcol)
     279          end do
     280        end do
     281      end do
    276282
    277283      if (present(incoming_sw)) then
    … …  
    287293    end if
    288294
    289     if (config%do_lw) then
    290 
    291       call config%gas_optics_lw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
    292            &  NMaxGases, thermodynamics%pressure_hl, &
    293            &  temperature_fl, &
    294            &  gas%mixing_ratio, &
    295 !           &  reshape(gas%mixing_ratio(istartcol:iendcol,:,:), &
    296 !           &          [nlev,iendcol-istartcol+1,NMaxGases],order=[2,1,3]), &
    297            &  od_lw)
     295    if (config%do_lw .and. config%i_gas_model_lw == IGasModelECCKD) then
     296
     297      if (is_volume_mixing_ratio) then
     298        call config%gas_optics_lw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
     299             &  NMaxGases, thermodynamics%pressure_hl, &
     300             &  temperature_fl, gas%mixing_ratio, &
     301             &  od_lw)
     302      else
     303        call config%gas_optics_lw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
     304             &  NMaxGases, thermodynamics%pressure_hl, &
     305             &  temperature_fl, gas%mixing_ratio, &
     306             &  od_lw, concentration_scaling=concentration_scaling)
     307      end if
    298308
    299309      ! Calculate the Planck function for each g point
    … …  
    305315           &  single_level%skin_temperature(istartcol:iendcol), &
    306316           &  lw_emission(:,:))
     317!NEC$ forced_collapse
    307318      lw_emission = lw_emission * (1.0_jprb - lw_albedo)
    308319

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_flux.F90

    r4773 r4853  
    2020!   2021-01-20  R. Hogan  Added heating_rate_out_of_physical_bounds function
    2121!   2022-12-07  R. Hogan  Added top-of-atmosphere spectral output
     22
     23#include "ecrad_config.h"
    2224
    2325module radiation_flux
    … …  
    117119
    118120! Added for DWD (2020)
    119 #ifdef __SX__
     121#ifdef DWD_VECTOR_OPTIMIZATIONS
    120122      logical, parameter :: use_indexed_sum_vec = .true.
    121123#else

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_gas.F90

    r4773 r4853  
    7272     procedure :: assert_units => assert_units_gas
    7373     procedure :: get        => get_gas
     74     procedure :: get_scaling
    7475     procedure :: reverse    => reverse_gas
    7576     procedure :: out_of_physical_bounds
    … …  
    355356    real(jprb), optional, intent(in)    :: scale_factor
    356357
    357     integer :: ig
     358    integer :: jg
    358359
    359360    ! Scaling factor to convert from old to new
    … …  
    396397      end if
    397398    else
    398       do ig = 1,this%ntype
    399         call this%set_units(iunits, igas=this%icode(ig), scale_factor=new_sf)
     399      do jg = 1,this%ntype
     400        call this%set_units(iunits, igas=this%icode(jg), scale_factor=new_sf)
    400401      end do
    401402    end if
    … …  
    403404  end subroutine set_units_gas
    404405
    405 
     406 
     407  !---------------------------------------------------------------------
     408  ! Return a vector indicating the scaling that one would need to
     409  ! apply to each gas in order to obtain the dimension units in
     410  ! "iunits" (which can be IVolumeMixingRatio or IMassMixingRatio)
     411  subroutine get_scaling(this, iunits, scaling)
     412    class(gas_type), intent(in)  :: this
     413    integer,         intent(in)  :: iunits
     414    real(jprb),      intent(out) :: scaling(NMaxGases)
     415    integer :: jg
     416   
     417    scaling = this%scale_factor
     418    do jg = 1,NMaxGases
     419      if (iunits == IMassMixingRatio .and. this%iunits(jg) == IVolumeMixingRatio) then
     420        scaling(jg) = scaling(jg) * GasMolarMass(jg) / AirMolarMass
     421      else if (iunits == IVolumeMixingRatio .and. this%iunits(jg) == IMassMixingRatio) then
     422        scaling(jg) = scaling(jg) * AirMolarMass / GasMolarMass(jg)
     423      end if
     424    end do
     425   
     426  end subroutine get_scaling
     427
     428 
    406429  !---------------------------------------------------------------------
    407430  ! Assert that gas mixing ratio units are "iunits", applying to gas
    408431  ! with ID "igas" if present, otherwise to all gases. Otherwise the
    409   ! program will exit. Otional argument scale factor specifies any
    410   ! subsequent multiplication to apply; for PPMV one would use
    411   ! iunits=IVolumeMixingRatio and scale_factor=1.0e6.
    412   recursive subroutine assert_units_gas(this, iunits, igas, scale_factor)
     432  ! program will exit, except if the optional argument "istatus" is
     433  ! provided in which case it will return true if the units are
     434  ! correct and false if they are not. Optional argument scale factor
     435  ! specifies any subsequent multiplication to apply; for PPMV one
     436  ! would use iunits=IVolumeMixingRatio and scale_factor=1.0e6.
     437  recursive subroutine assert_units_gas(this, iunits, igas, scale_factor, istatus)
    413438
    414439    use radiation_io,   only : nulerr, radiation_abort
    415440
    416     class(gas_type),      intent(in) :: this
    417     integer,              intent(in) :: iunits
    418     integer,    optional, intent(in) :: igas
    419     real(jprb), optional, intent(in) :: scale_factor
    420 
    421     integer :: ig
     441    class(gas_type),      intent(in)  :: this
     442    integer,              intent(in)  :: iunits
     443    integer,    optional, intent(in)  :: igas
     444    real(jprb), optional, intent(in)  :: scale_factor
     445    logical,    optional, intent(out) :: istatus
     446
     447    integer :: jg
    422448
    423449    real(jprb) :: sf
    … …  
    429455    end if
    430456
     457    if (present(istatus)) then
     458      istatus = .true.
     459    end if
     460   
    431461    if (present(igas)) then
    432462      if (this%is_present(igas)) then
    433463        if (iunits /= this%iunits(igas)) then
    434           write(nulerr,'(a,a,a)') '*** Error: ', trim(GasName(igas)), &
    435                &  ' is not in the required units'
    436           call radiation_abort()
     464          if (present(istatus)) then
     465            istatus = .false.
     466          else
     467            write(nulerr,'(a,a,a)') '*** Error: ', trim(GasName(igas)), &
     468                 &  ' is not in the required units'
     469            call radiation_abort()
     470          end if
    437471        else if (sf /= this%scale_factor(igas)) then
    438           write(nulerr,'(a,a,a,e12.4,a,e12.4)') '*** Error: ', GasName(igas), &
    439                &  ' scaling of ', this%scale_factor(igas), &
    440                &  ' does not match required ', sf
    441           call radiation_abort()
     472          if (present(istatus)) then
     473            istatus = .false.
     474          else
     475            write(nulerr,'(a,a,a,e12.4,a,e12.4)') '*** Error: ', GasName(igas), &
     476                 &  ' scaling of ', this%scale_factor(igas), &
     477                 &  ' does not match required ', sf
     478            call radiation_abort()
     479          end if
    442480        end if
    443481      end if
    444482    else
    445       do ig = 1,this%ntype
    446         call this%assert_units(iunits, igas=this%icode(ig), scale_factor=sf)
     483      do jg = 1,this%ntype
     484        call this%assert_units(iunits, igas=this%icode(jg), scale_factor=sf, istatus=istatus)
    447485      end do
    448486    end if

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90

    r4773 r4853  
    7575    ! Allocate structures
    7676    if (config%do_sw) then
    77       if (allocated(config%cloud_optics_sw)) deallocate(config%cloud_optics_sw)     
    7877      allocate(config%cloud_optics_sw(config%n_cloud_types))
    7978    end if
    8079
    8180    if (config%do_lw) then
    82       if (allocated(config%cloud_optics_lw)) deallocate(config%cloud_optics_lw)     
    8381      allocate(config%cloud_optics_lw(config%n_cloud_types))
    8482    end if
    … …  
    172170    real(jprb), dimension(istartcol:iendcol,nlev) :: water_path
    173171
    174     integer :: jtype, jcol, jlev
     172    integer :: jtype, jcol, jlev, jg
    175173
    176174    real(jphook) :: hook_handle
    … …  
    275273          if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    276274            ! Scale to get asymmetry factor and single scattering albedo
    277             g_sw_cloud(:,jlev,jcol) = g_sw_cloud(:,jlev,jcol) &
    278                  &  / max(ssa_sw_cloud(:,jlev,jcol), 1.0e-15_jprb)
    279             ssa_sw_cloud(:,jlev,jcol) = ssa_sw_cloud(:,jlev,jcol) &
    280                  &  / max(od_sw_cloud(:,jlev,jcol),  1.0e-15_jprb)
     275            do jg = 1, config%n_bands_sw
     276              g_sw_cloud(jg,jlev,jcol) = g_sw_cloud(jg,jlev,jcol) &
     277                 &  / max(ssa_sw_cloud(jg,jlev,jcol), 1.0e-15_jprb)
     278              ssa_sw_cloud(jg,jlev,jcol) = ssa_sw_cloud(jg,jlev,jcol) &
     279                 &  / max(od_sw_cloud(jg,jlev,jcol),  1.0e-15_jprb)
     280            end do
    281281          end if
    282282        end do

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics_data.F90

    r4802 r4853  
    1515!
    1616
     17#include "ecrad_config.h"
     18
    1719module radiation_general_cloud_optics_data
    1820
    … …  
    7375
    7476    use yomhook,                       only : lhook, dr_hook, jphook
    75     use easy_netcdf,                   only : netcdf_file
     77#ifdef EASY_NETCDF_READ_MPI
     78    use easy_netcdf_read_mpi, only : netcdf_file
     79#else
     80    use easy_netcdf,          only : netcdf_file
     81#endif
    7682    use radiation_spectral_definition, only : spectral_definition_type
    7783    use radiation_io,                  only : nulout, nulerr, radiation_abort
    … …  
    181187
    182188    ! Thin averaging
    183 
    184     ! Circumvent ifort bug (see https://github.com/ecmwf-ifs/ecrad/issues/13):
    185     allocate(this%mass_ext(size(mapping, 1), nre))
    186 
    187189    this%mass_ext  = matmul(mapping, mass_ext)
    188190    this%ssa       = matmul(mapping, mass_ext*ssa) / this%mass_ext
    … …  
    267269         &     scat_asymmetry ! Scattering optical depth x asymmetry factor
    268270
    269     real(jprb) :: od_local(ng)
     271    real(jprb) :: od_local
    270272
    271273    real(jprb) :: re_index, weight1, weight2
    272274    integer :: ire
    273275
    274     integer :: jcol, jlev
     276    integer :: jcol, jlev, jg
    275277
    276278    real(jphook) :: hook_handle
    … …  
    287289            weight2 = re_index - ire
    288290            weight1 = 1.0_jprb - weight2
    289             od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) &
    290                  &                              +weight2*this%mass_ext(:,ire+1))
    291             od(:,jlev,jcol) = od(:,jlev,jcol) + od_local
    292             od_local = od_local * (weight1*this%ssa(:,ire) &
    293                  &                +weight2*this%ssa(:,ire+1))
    294             scat_od(:,jlev,jcol) = scat_od(:,jlev,jcol) + od_local
    295             scat_asymmetry(:,jlev,jcol) = scat_asymmetry(:,jlev,jcol) &
    296                  & + od_local * (weight1*this%asymmetry(:,ire) &
    297                  &              +weight2*this%asymmetry(:,ire+1))
     291            do jg = 1, ng
     292              od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(jg,ire) &
     293                 &                                +weight2*this%mass_ext(jg,ire+1))
     294              od(jg,jlev,jcol) = od(jg,jlev,jcol) + od_local
     295              od_local = od_local * (weight1*this%ssa(jg,ire) &
     296                 &                  +weight2*this%ssa(jg,ire+1))
     297              scat_od(jg,jlev,jcol) = scat_od(jg,jlev,jcol) + od_local
     298              scat_asymmetry(jg,jlev,jcol) = scat_asymmetry(jg,jlev,jcol) &
     299                 & + od_local * (weight1*this%asymmetry(jg,ire) &
     300                 &              +weight2*this%asymmetry(jg,ire+1))
     301            end do
     302
    298303          end if
    299304        end do

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90

    r4773 r4853  
    4141
    4242    use radiation_config
     43    use radiation_spectral_definition, only &
     44         &  : SolarReferenceTemperature, TerrestrialReferenceTemperature
    4345
    4446    type(config_type), intent(inout), target :: config
    … …  
    6567          &   44, 107, 94, 14, 108, 15, 16, 109, 17, 18, 110, 111, 112 &
    6668          & /)
     69
     70    logical :: do_sw, do_lw
    6771   
    6872    real(jphook) :: hook_handle
    … …  
    7781    if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',0,hook_handle)
    7882
     83    do_sw = (config%do_sw .and. config%i_gas_model_sw == IGasModelIFSRRTMG)
     84    do_lw = (config%do_lw .and. config%i_gas_model_lw == IGasModelIFSRRTMG)
     85   
    7986    ! The IFS implementation of RRTMG uses many global variables.  In
    8087    ! the IFS these will have been set up already; otherwise set them
    … …  
    8592      call SURRTPK
    8693      call SURRTRF
    87       call RRTM_INIT_140GP(directory)
    88       call SRTM_INIT(directory)
     94      if (do_lw) then
     95        call RRTM_INIT_140GP(directory)
     96      end if
     97      if (do_sw) then
     98        call SRTM_INIT(directory)
     99      end if
    89100    end if
    90101
    91     ! Cloud and aerosol properties can only be defined per band
    92     config%do_cloud_aerosol_per_sw_g_point = .false.
    93     config%do_cloud_aerosol_per_lw_g_point = .false.
    94 
    95     config%n_g_sw = jpgsw
    96     config%n_g_lw = jpglw
    97     config%n_bands_sw = 14
    98     config%n_bands_lw = 16
    99 
    100     ! Wavenumber ranges of each band may be needed so that the user
    101     ! can compute UV and photosynthetically active radiation for a
    102     ! particular wavelength range
    103     call config%gas_optics_sw%spectral_def%allocate_bands_only( &
    104          &  [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
    105          &   8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
    106          &  [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
    107          &   12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
    108     call config%gas_optics_lw%spectral_def%allocate_bands_only( &
    109          &  [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
    110          &   1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
    111          &  [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
    112          &   1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
    113 
    114     allocate(config%i_band_from_g_sw          (config%n_g_sw))
    115     allocate(config%i_band_from_g_lw          (config%n_g_lw))
    116     allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
    117     allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
    118     allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
    119     allocate(config%i_g_from_reordered_g_lw(config%n_g_lw))
    120 
    121     ! Shortwave starts at 16: need to start at 1
    122     config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
    123     config%i_band_from_g_lw = ngb_lw
    124 
    125     if (config%i_solver_sw == ISolverSpartacus) then
    126       ! SPARTACUS requires g points ordered in approximately
    127       ! increasing order of optical depth
    128       config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
    129     else
    130       ! Implied-do for no reordering
    131 !      config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
    132       config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
     102    if (do_sw) then
     103     
     104      ! Cloud and aerosol properties can only be defined per band
     105      config%do_cloud_aerosol_per_sw_g_point = .false.
     106      config%n_g_sw = jpgsw
     107      config%n_bands_sw = 14
     108      ! Wavenumber ranges of each band may be needed so that the user
     109      ! can compute UV and photosynthetically active radiation for a
     110      ! particular wavelength range
     111      call config%gas_optics_sw%spectral_def%allocate_bands_only(SolarReferenceTemperature, &
     112           &  [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
     113           &   8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
     114           &  [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
     115           &   12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
     116      allocate(config%i_band_from_g_sw          (config%n_g_sw))
     117      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
     118      allocate(config%i_g_from_reordered_g_sw   (config%n_g_sw))
     119      ! Shortwave starts at 16: need to start at 1
     120      config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
     121
     122      if (config%i_solver_sw == ISolverSpartacus) then
     123        ! SPARTACUS requires g points ordered in approximately
     124        ! increasing order of optical depth
     125        config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
     126      else
     127        ! Implied-do for no reordering
     128        !      config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
     129        config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
     130      end if
     131
     132      config%i_band_from_reordered_g_sw &
     133           = config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
     134
     135      ! The i_spec_* variables are used solely for storing spectral
     136      ! data, and this can either be by band or by g-point
     137      if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
     138        if (config%do_save_gpoint_flux) then
     139          config%n_spec_sw = config%n_g_sw
     140          config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
     141        else
     142          config%n_spec_sw = config%n_bands_sw
     143          config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
     144        end if
     145      else
     146        config%n_spec_sw = 0
     147        nullify(config%i_spec_from_reordered_g_sw)
     148      end if
     149     
    133150    end if
    134151
    135     if (config%i_solver_lw == ISolverSpartacus) then
    136       ! SPARTACUS requires g points ordered in approximately
    137       ! increasing order of optical depth
    138       config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
    139     else
    140       ! Implied-do for no reordering
    141       config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
     152    if (do_lw) then
     153      ! Cloud and aerosol properties can only be defined per band
     154      config%do_cloud_aerosol_per_lw_g_point = .false.
     155      config%n_g_lw = jpglw
     156      config%n_bands_lw = 16
     157      call config%gas_optics_lw%spectral_def%allocate_bands_only(TerrestrialReferenceTemperature, &
     158           &  [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
     159           &   1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
     160           &  [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
     161           &   1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
     162      allocate(config%i_band_from_g_lw          (config%n_g_lw))
     163      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
     164      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
     165      config%i_band_from_g_lw = ngb_lw
     166
     167      if (config%i_solver_lw == ISolverSpartacus) then
     168        ! SPARTACUS requires g points ordered in approximately
     169        ! increasing order of optical depth
     170        config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
     171      else
     172        ! Implied-do for no reordering
     173        config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
     174      end if
     175
     176      config%i_band_from_reordered_g_lw &
     177           = config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
     178
     179      ! The i_spec_* variables are used solely for storing spectral
     180      ! data, and this can either be by band or by g-point
     181      if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
     182        if (config%do_save_gpoint_flux) then
     183          config%n_spec_lw = config%n_g_lw
     184          config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
     185        else
     186          config%n_spec_lw = config%n_bands_lw
     187          config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
     188        end if
     189      else
     190        config%n_spec_lw = 0
     191        nullify(config%i_spec_from_reordered_g_lw)
     192      end if
     193
    142194    end if
    143 
    144     config%i_band_from_reordered_g_sw &
    145          = config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
    146 
    147     config%i_band_from_reordered_g_lw &
    148          = config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
    149 
    150     ! The i_spec_* variables are used solely for storing spectral
    151     ! data, and this can either be by band or by g-point
    152     if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
    153       if (config%do_save_gpoint_flux) then
    154         config%n_spec_sw = config%n_g_sw
    155         config%n_spec_lw = config%n_g_lw
    156         config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
    157         config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
    158       else
    159         config%n_spec_sw = config%n_bands_sw
    160         config%n_spec_lw = config%n_bands_lw
    161         config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
    162         config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
    163       end if
    164     else
    165       config%n_spec_sw = 0
    166       config%n_spec_lw = 0
    167       nullify(config%i_spec_from_reordered_g_sw)
    168       nullify(config%i_spec_from_reordered_g_lw)
    169     end if
    170 
     195   
    171196    if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',1,hook_handle)
    172197
    … …  
    201226    use yomhook  , only : lhook, dr_hook, jphook
    202227
    203     use radiation_config,         only : config_type, ISolverSpartacus
     228    use radiation_config,         only : config_type, ISolverSpartacus, IGasModelIFSRRTMG
    204229    use radiation_thermodynamics, only : thermodynamics_type
    205230    use radiation_single_level,   only : single_level_type
    … …  
    334359    integer :: istartlev, iendlev
    335360
     361    logical :: do_sw, do_lw
     362   
    336363    integer :: jlev, jgreorder, jg, ig, iband, jcol
    337364
    … …  
    345372
    346373    if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',0,hook_handle)
     374
     375    do_sw = (config%do_sw .and. config%i_gas_model_sw == IGasModelIFSRRTMG)
     376    do_lw = (config%do_lw .and. config%i_gas_model_lw == IGasModelIFSRRTMG)
    347377
    348378    ! Compute start and end levels for indexing the gas mixing ratio
    … …  
    392422         &  ZPAVEL , ZTAVEL , ZPZ , ZTZ, IREFLECT) 
    393423
    394     CALL RRTM_SETCOEF_140GP &
    395          &( istartcol, iendcol, nlev , ZCOLDRY  , ZWBRODL , ZWKL , &
    396          &  ZFAC00 , ZFAC01   , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
    397          &  ZCOLH2O, ZCOLCO2  , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
    398          &  ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
    399          &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
    400          &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
    401          &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
    402          &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)   
    403 
    404     ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
    405 
    406     CALL RRTM_GAS_OPTICAL_DEPTH &
    407          &( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
    408          &  ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
    409          &  JP, JT, JT1, ZONEMINUS ,&
    410          &  ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
    411          &  ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
    412          &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
    413          &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
    414          &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
    415          &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)     
    416 
    417     if (present(lw_albedo)) then
    418    
    419       call planck_function_atmos(nlev, istartcol, iendcol, config, &
    420            &                     thermodynamics, ZPFRAC, planck_hl)
    421 
    422       if (single_level%is_simple_surface) then
    423         call planck_function_surf(istartcol, iendcol, config, &
    424              &                    single_level%skin_temperature, ZPFRAC(:,:,1), &
    425              &                    lw_emission)
     424    if (do_lw) then
     425   
     426      CALL RRTM_SETCOEF_140GP &
     427           &( istartcol, iendcol, nlev , ZCOLDRY  , ZWBRODL , ZWKL , &
     428           &  ZFAC00 , ZFAC01   , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
     429           &  ZCOLH2O, ZCOLCO2  , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
     430           &  ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
     431           &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
     432           &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
     433           &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
     434           &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)   
     435
     436      ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
     437
     438      CALL RRTM_GAS_OPTICAL_DEPTH &
     439           &( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
     440           &  ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
     441           &  JP, JT, JT1, ZONEMINUS ,&
     442           &  ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
     443           &  ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
     444           &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
     445           &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
     446           &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
     447           &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)     
     448   
     449      if (present(lw_albedo)) then
     450   
     451        call planck_function_atmos(nlev, istartcol, iendcol, config, &
     452             &                     thermodynamics, ZPFRAC, planck_hl)
     453
     454        if (single_level%is_simple_surface) then
     455          call planck_function_surf(istartcol, iendcol, config, &
     456               &                    single_level%skin_temperature, ZPFRAC(:,:,1), &
     457               &                    lw_emission)
     458         
     459          ! The following can be used to extract the parameters defined at
     460          ! the top of the planck_function routine below:
     461          !write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
     462          !&  sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
    426463       
    427         ! The following can be used to extract the parameters defined at
    428         ! the top of the planck_function routine below:
    429         !write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
    430         !&  sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
    431        
    432         ! lw_emission at this point is actually the planck function of
    433         ! the surface
    434         lw_emission = lw_emission * (1.0_jprb - lw_albedo)
     464          ! lw_emission at this point is actually the planck function of
     465          ! the surface
     466          lw_emission = lw_emission * (1.0_jprb - lw_albedo)
     467        else
     468          ! Longwave emission has already been computed
     469          if (config%use_canopy_full_spectrum_lw) then
     470            lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
     471          else
     472            lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
     473                 & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
     474          end if
     475        end if
     476
     477      end if
     478
     479      if (config%i_solver_lw == ISolverSpartacus) then
     480        ! We need to rearrange the gas optics info in memory: reordering
     481        ! the g points in order of approximately increasing optical
     482        ! depth (for efficient 3D processing on only the regions of the
     483        ! spectrum that are optically thin for gases) and reorder in
     484        ! pressure since the the functions above treat pressure
     485        ! decreasing with increasing index.  Note that the output gas
     486        ! arrays have dimensions in a different order to the inputs,
     487        ! so there is some inefficiency here.
     488        do jgreorder = 1,config%n_g_lw
     489          iband = config%i_band_from_reordered_g_lw(jgreorder)
     490          ig = config%i_g_from_reordered_g_lw(jgreorder)
     491         
     492          ! Top-of-atmosphere half level
     493          do jlev = 1,nlev
     494            do jcol = istartcol,iendcol
     495              ! Some g points can return negative optical depths;
     496              ! specifically original g points 54-56 which causes
     497              ! unphysical single-scattering albedo when combined with
     498              ! aerosol
     499              od_lw(jgreorder,jlev,jcol) &
     500                   &   = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
     501            end do
     502          end do
     503        end do
    435504      else
    436         ! Longwave emission has already been computed
    437         if (config%use_canopy_full_spectrum_lw) then
    438           lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
    439         else
    440           lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
    441                & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
    442         end if
     505        ! G points have not been reordered
     506        do jcol = istartcol,iendcol
     507          do jlev = 1,nlev
     508            ! Check for negative optical depth
     509            od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
     510          end do
     511        end do
    443512      end if
    444513
    445514    end if
    446515
    447     if (config%i_solver_lw == ISolverSpartacus) then
    448       ! We need to rearrange the gas optics info in memory: reordering
    449       ! the g points in order of approximately increasing optical
    450       ! depth (for efficient 3D processing on only the regions of the
    451       ! spectrum that are optically thin for gases) and reorder in
    452       ! pressure since the the functions above treat pressure
    453       ! decreasing with increasing index.  Note that the output gas
    454       ! arrays have dimensions in a different order to the inputs,
    455       ! so there is some inefficiency here.
    456       do jgreorder = 1,config%n_g_lw
    457         iband = config%i_band_from_reordered_g_lw(jgreorder)
    458         ig = config%i_g_from_reordered_g_lw(jgreorder)
    459        
    460         ! Top-of-atmosphere half level
    461         do jlev = 1,nlev
    462           do jcol = istartcol,iendcol
    463             ! Some g points can return negative optical depths;
    464             ! specifically original g points 54-56 which causes
    465             ! unphysical single-scattering albedo when combined with
    466             ! aerosol
    467             od_lw(jgreorder,jlev,jcol) &
    468                  &   = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
     516    if (do_sw) then
     517   
     518      CALL SRTM_SETCOEF &
     519           & ( istartcol, iendcol, nlev,&
     520           & ZPAVEL  , ZTAVEL,&
     521           & ZCOLDRY , ZWKL,&
     522           & ILAYTROP,&
     523           & ZCOLCH4  , ZCOLCO2 , ZCOLH2O , ZCOLMOL  , ZCOLO2 , ZCOLO3,&
     524           & ZFORFAC , ZFORFRAC , INDFOR  , ZSELFFAC, ZSELFFRAC, INDSELF, &
     525           & ZFAC00  , ZFAC01   , ZFAC10  , ZFAC11,&
     526           & JP      , JT       , JT1     , single_level%cos_sza(istartcol:iendcol)  &
     527           & ) 
     528   
     529      ! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
     530      ! is below the horizon, so we do it here
     531      ZOD_SW(istartcol:iendcol,:,:)  = 0.0_jprb
     532      ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
     533      ZINCSOL(istartcol:iendcol,:)   = 0.0_jprb
     534
     535      CALL SRTM_GAS_OPTICAL_DEPTH &
     536           &( istartcol, iendcol , nlev  , ZONEMINUS_ARRAY,&
     537           & single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
     538           & ZCOLCH4 , ZCOLCO2  , ZCOLH2O, ZCOLMOL , ZCOLO2   , ZCOLO3,&
     539           & ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
     540           & ZFAC00  , ZFAC01   , ZFAC10 , ZFAC11  ,&
     541           & JP      , JT       , JT1    ,&
     542           & ZOD_SW  , ZSSA_SW  , ZINCSOL )
     543   
     544      ! Scale the incoming solar per band, if requested
     545      if (config%use_spectral_solar_scaling) then
     546        do jg = 1,JPGPT_SW
     547          do jcol = istartcol,iendcol
     548            ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
     549                 &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
    469550          end do
    470551        end do
    471       end do
    472     else
    473       ! G points have not been reordered
    474       do jcol = istartcol,iendcol
    475         do jlev = 1,nlev
    476           ! Check for negative optical depth
    477           od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
    478         end do
    479       end do
    480     end if
    481    
    482     CALL SRTM_SETCOEF &
    483          & ( istartcol, iendcol, nlev,&
    484          & ZPAVEL  , ZTAVEL,&
    485          & ZCOLDRY , ZWKL,&
    486          & ILAYTROP,&
    487          & ZCOLCH4  , ZCOLCO2 , ZCOLH2O , ZCOLMOL  , ZCOLO2 , ZCOLO3,&
    488          & ZFORFAC , ZFORFRAC , INDFOR  , ZSELFFAC, ZSELFFRAC, INDSELF, &
    489          & ZFAC00  , ZFAC01   , ZFAC10  , ZFAC11,&
    490          & JP      , JT       , JT1     , single_level%cos_sza(istartcol:iendcol)  &
    491          & ) 
    492    
    493     ! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
    494     ! is below the horizon, so we do it here
    495     ZOD_SW(istartcol:iendcol,:,:)  = 0.0_jprb
    496     ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
    497     ZINCSOL(istartcol:iendcol,:)   = 0.0_jprb
    498 
    499     CALL SRTM_GAS_OPTICAL_DEPTH &
    500          &( istartcol, iendcol , nlev  , ZONEMINUS_ARRAY,&
    501          & single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
    502          & ZCOLCH4 , ZCOLCO2  , ZCOLH2O, ZCOLMOL , ZCOLO2   , ZCOLO3,&
    503          & ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
    504          & ZFAC00  , ZFAC01   , ZFAC10 , ZFAC11  ,&
    505          & JP      , JT       , JT1    ,&
    506          & ZOD_SW  , ZSSA_SW  , ZINCSOL )
    507    
    508     ! Scale the incoming solar per band, if requested
    509     if (config%use_spectral_solar_scaling) then
    510       do jg = 1,JPGPT_SW
    511         do jcol = istartcol,iendcol
    512           ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
    513             &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
    514         end do
    515       end do
    516     end if
    517 
    518     ! Scaling factor to ensure that the total solar irradiance is as
    519     ! requested.  Note that if the sun is below the horizon then
    520     ! ZINCSOL will be zero.
    521     if (present(incoming_sw)) then
    522       do jcol = istartcol,iendcol
    523         if (single_level%cos_sza(jcol) > 0.0_jprb) then
     552      end if
     553
     554      ! Scaling factor to ensure that the total solar irradiance is as
     555      ! requested.  Note that if the sun is below the horizon then
     556      ! ZINCSOL will be zero.
     557      if (present(incoming_sw)) then
     558        do jcol = istartcol,iendcol
     559          if (single_level%cos_sza(jcol) > 0.0_jprb) then
    524560! Added for DWD (2020)
    525561!NEC$ nounroll
    526           incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
    527         else
    528           incoming_sw_scale(jcol) = 1.0_jprb
    529         end if
    530       end do
    531     end if
    532 
    533     if (config%i_solver_sw == ISolverSpartacus) then
    534 !    if (.true.) then
    535       ! Account for reordered g points
    536       do jgreorder = 1,config%n_g_sw
    537         ig = config%i_g_from_reordered_g_sw(jgreorder)
    538         do jlev = 1,nlev
    539           do jcol = istartcol,iendcol
    540             ! Check for negative optical depth
    541             od_sw (jgreorder,nlev+1-jlev,jcol) &
    542                  &  = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
    543             ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
    544            end do
     562            incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
     563          else
     564            incoming_sw_scale(jcol) = 1.0_jprb
     565          end if
    545566        end do
    546         if (present(incoming_sw)) then
    547           incoming_sw(jgreorder,:) &
    548                &  = incoming_sw_scale(:) * ZINCSOL(:,ig)
    549         end if
    550       end do
    551     else
    552       ! G points have not been reordered
    553       do jcol = istartcol,iendcol
    554         do jlev = 1,nlev
    555           do jg = 1,config%n_g_sw
    556             ! Check for negative optical depth
    557             od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
    558             ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
    559           end do
     567      end if
     568
     569      if (config%i_solver_sw == ISolverSpartacus) then
     570!      if (.true.) then
     571        ! Account for reordered g points
     572        do jgreorder = 1,config%n_g_sw
     573          ig = config%i_g_from_reordered_g_sw(jgreorder)
     574          do jlev = 1,nlev
     575            do jcol = istartcol,iendcol
     576              ! Check for negative optical depth
     577              od_sw (jgreorder,nlev+1-jlev,jcol) &
     578                   &  = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
     579              ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
     580            end do
     581          end do
     582          if (present(incoming_sw)) then
     583            incoming_sw(jgreorder,:) &
     584                 &  = incoming_sw_scale(:) * ZINCSOL(:,ig)
     585          end if
    560586        end do
    561         if (present(incoming_sw)) then
    562           do jg = 1,config%n_g_sw
    563             incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
    564           end do
    565         end if
    566       end do
     587      else
     588        ! G points have not been reordered
     589        do jcol = istartcol,iendcol
     590          do jlev = 1,nlev
     591            do jg = 1,config%n_g_sw
     592              ! Check for negative optical depth
     593              od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
     594              ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
     595            end do
     596          end do
     597          if (present(incoming_sw)) then
     598            do jg = 1,config%n_g_sw
     599              incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
     600            end do
     601          end if
     602        end do
     603      end if
     604
    567605    end if
    568606   

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_interface.F90

    r4773 r4853  
    6868
    6969    ! Load the look-up tables from files in the specified directory
    70     if (config%i_gas_model == IGasModelMonochromatic) then
     70    if (config%i_gas_model_sw == IGasModelMonochromatic) then
    7171      call setup_gas_optics_mono(config, trim(config%directory_name))
    72     else if (config%i_gas_model == IGasModelIFSRRTMG) then
    73       call setup_gas_optics_rrtmg(config, trim(config%directory_name))
    74     else if (config%i_gas_model == IGasModelECCKD) then
    75       call setup_gas_optics_ecckd(config)
     72    else
     73      ! Note that we can run RRTMG and ECCKD for different parts of
     74      ! the spectrum: the setup routines only configure the relevant
     75      ! part.
     76      if (config%i_gas_model_sw == IGasModelIFSRRTMG .or. config%i_gas_model_lw == IGasModelIFSRRTMG) then
     77        call setup_gas_optics_rrtmg(config, trim(config%directory_name))
     78      end if
     79      if (config%i_gas_model_sw == IGasModelECCKD .or. config%i_gas_model_lw == IGasModelECCKD) then
     80        call setup_gas_optics_ecckd(config)
     81      end if
    7682    end if
    7783
    … …  
    122128
    123129    if (config%do_clouds) then
    124       if (config%i_gas_model == IGasModelMonochromatic) then
     130      if (config%i_gas_model_sw == IGasModelMonochromatic) then
    125131        !      call setup_cloud_optics_mono(config)
    126132      elseif (config%use_general_cloud_optics) then
    … …  
    132138
    133139    if (config%use_aerosols) then
    134       if (config%i_gas_model == IGasModelMonochromatic) then
     140      if (config%i_gas_model_sw == IGasModelMonochromatic) then
    135141!        call setup_aerosol_optics_mono(config)
    136142      else
    … …  
    167173    type(gas_type),    intent(inout) :: gas
    168174
    169     if (config%i_gas_model == IGasModelMonochromatic) then
     175    if (config%i_gas_model_sw == IGasModelMonochromatic) then
    170176      call set_gas_units_mono(gas)
    171     elseif (config%i_gas_model == IGasModelECCKD) then
     177    elseif (config%i_gas_model_sw == IGasModelIFSRRTMG &
     178         &  .or. config%i_gas_model_lw == IGasModelIFSRRTMG) then
     179      ! Convert to mass-mixing ratio for RRTMG; note that ecCKD can
     180      ! work with this but performs an internal scaling
     181      call set_gas_units_ifs(gas)
     182    else
     183      ! Use volume mixing ratio preferred by ecCKD
    172184      call set_gas_units_ecckd(gas)
    173     else
    174       call set_gas_units_ifs(gas)
    175185    end if
    176186
    … …  
    196206    use radiation_io,             only : nulout
    197207    use radiation_config,         only : config_type, &
    198          &   IGasModelMonochromatic, IGasModelIFSRRTMG, &
     208         &   IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD, &
    199209         &   ISolverMcICA, ISolverSpartacus, ISolverHomogeneous, &
    200210         &   ISolverTripleclouds
    … …  
    227237    use radiation_general_cloud_optics, only : general_cloud_optics
    228238    use radiation_aerosol_optics, only : add_aerosol_optics
    229     USE mod_phys_lmdz_para
    230 
    231239
    232240    ! Inputs
    … …  
    238246    type(thermodynamics_type),intent(in)   :: thermodynamics
    239247    type(gas_type),           intent(in)   :: gas
    240     type(cloud_type),        intent(inout):: cloud
     248    type(cloud_type),         intent(inout):: cloud
    241249    type(aerosol_type),       intent(in)   :: aerosol
    242250    ! Output
    … …  
    297305
    298306    real(jphook) :: hook_handle
    299     integer :: jcol, jlev
    300307
    301308    if (lhook) call dr_hook('radiation_interface:radiation',0,hook_handle)
    302 
    303     if (config%i_solver_sw == ISolverSPARTACUS) then
    304          print*,'Dans radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
    305               mpi_rank, omp_rank, &
    306               shape(cloud%inv_cloud_effective_size), &
    307               size(cloud%inv_cloud_effective_size)   
    308 !         do jcol=istartcol, iendcol
    309 !              do jlev=1,nlev
    310 !              print*,'Entree radiation_interf, mpi_rank, omp_rank, jcol, jlev &
    311 !           &  cloud%inv_cloud_effective_size =',mpi_rank, omp_rank, jcol, jlev, &
    312 !           &  cloud%inv_cloud_effective_size(jcol,jlev)
    313 !              enddo
    314 !          enddo
    315     endif
    316 !    cloud%inv_cloud_effective_size=0.05_jprb
    317309
    318310    if (thermodynamics%pressure_hl(istartcol,2) &
    … …  
    339331      ! incoming shortwave flux at each g-point, for the specified
    340332      ! range of atmospheric columns
    341       if (config%i_gas_model == IGasModelMonochromatic) then
     333      if (config%i_gas_model_sw == IGasModelMonochromatic) then
    342334        call gas_optics_mono(ncol,nlev,istartcol,iendcol, config, &
    343335             &  single_level, thermodynamics, gas, lw_albedo, &
    344336             &  od_lw, od_sw, ssa_sw, &
    345337             &  planck_hl, lw_emission, incoming_sw)
    346       else if (config%i_gas_model == IGasModelIFSRRTMG) then
    347         call gas_optics_rrtmg(ncol,nlev,istartcol,iendcol, config, &
    348              &  single_level, thermodynamics, gas, &
    349              &  od_lw, od_sw, ssa_sw, lw_albedo=lw_albedo, &
    350              &  planck_hl=planck_hl, lw_emission=lw_emission, &
    351              &  incoming_sw=incoming_sw)
    352338      else
    353         call gas_optics_ecckd(ncol,nlev,istartcol,iendcol, config, &
    354              &  single_level, thermodynamics, gas, &
    355              &  od_lw, od_sw, ssa_sw, lw_albedo=lw_albedo, &
    356              &  planck_hl=planck_hl, lw_emission=lw_emission, &
    357              &  incoming_sw=incoming_sw)
     339        if (config%i_gas_model_sw == IGasModelIFSRRTMG &
     340             &   .or. config%i_gas_model_lw == IGasModelIFSRRTMG) then
     341          call gas_optics_rrtmg(ncol,nlev,istartcol,iendcol, config, &
     342               &  single_level, thermodynamics, gas, &
     343               &  od_lw, od_sw, ssa_sw, lw_albedo=lw_albedo, &
     344               &  planck_hl=planck_hl, lw_emission=lw_emission, &
     345               &  incoming_sw=incoming_sw)
     346        end if
     347        if (config%i_gas_model_sw == IGasModelECCKD &
     348             &   .or. config%i_gas_model_lw == IGasModelECCKD) then
     349          call gas_optics_ecckd(ncol,nlev,istartcol,iendcol, config, &
     350               &  single_level, thermodynamics, gas, &
     351               &  od_lw, od_sw, ssa_sw, lw_albedo=lw_albedo, &
     352               &  planck_hl=planck_hl, lw_emission=lw_emission, &
     353               &  incoming_sw=incoming_sw)
     354        end if
    358355      end if
    359356
    … …  
    368365        ! Compute hydrometeor absorption/scattering properties in each
    369366        ! shortwave and longwave band
    370         if (config%i_gas_model == IGasModelMonochromatic) then
     367        if (config%i_gas_model_sw == IGasModelMonochromatic) then
    371368          call cloud_optics_mono(nlev, istartcol, iendcol, &
    372369               &  config, thermodynamics, cloud, &
    … …  
    387384
    388385      if (config%use_aerosols) then
    389         if (config%i_gas_model == IGasModelMonochromatic) then
     386        if (config%i_gas_model_sw == IGasModelMonochromatic) then
    390387!          call add_aerosol_optics_mono(nlev,istartcol,iendcol, &
    391388!               &  config, thermodynamics, gas, aerosol, &
    … …  
    474471        else if (config%i_solver_sw == ISolverSPARTACUS) then
    475472          ! Compute fluxes using the SPARTACUS shortwave solver
    476 !             cloud%inv_cloud_effective_size=0.05_jprb
    477 !             do jcol=istartcol, iendcol
    478 !              do jlev=1,nlev
    479 !              print*,'jcol, jlev, dans radiation_interf i &
    480 !               &       cloud%inv_cloud_effective_size =',jcol, jlev, &
    481 !               cloud%inv_cloud_effective_size(jcol,jlev)
    482 !              enddo
    483 !             enddo
    484473          call solver_spartacus_sw(nlev,istartcol,iendcol, &
    485474               &  config, single_level, thermodynamics, cloud, &

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90

    r4773 r4853  
    1919!   2017-10-23  R. Hogan  Renamed single-character variables
    2020
     21#include "ecrad_config.h"
     22
    2123module radiation_mcica_lw
    2224
    … …  
    124126    ! Identify clear-sky layers
    125127    logical :: is_clear_sky_layer(nlev)
     128
     129    ! Temporary storage for more efficient summation
     130#ifdef DWD_REDUCTION_OPTIMIZATIONS
     131    real(jprb), dimension(nlev+1,2) :: sum_aux
     132#else
     133    real(jprb) :: sum_up, sum_dn
     134#endif
    126135
    127136    ! Index of the highest cloudy layer
    … …  
    179188
    180189      ! Sum over g-points to compute broadband fluxes
    181       flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
    182       flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
     190#ifdef DWD_REDUCTION_OPTIMIZATIONS
     191      sum_aux(:,:) = 0.0_jprb
     192      do jg = 1,ng
     193        do jlev = 1,nlev+1
     194          sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up_clear(jg,jlev)
     195          sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_clear(jg,jlev)
     196        end do
     197      end do
     198      flux%lw_up_clear(jcol,:) = sum_aux(:,1)
     199      flux%lw_dn_clear(jcol,:) = sum_aux(:,2)
     200#else
     201      do jlev = 1,nlev+1
     202        sum_up = 0.0_jprb
     203        sum_dn = 0.0_jprb
     204        !$omp simd reduction(+:sum_up, sum_dn)
     205        do jg = 1,ng
     206          sum_up = sum_up + flux_up_clear(jg,jlev)
     207          sum_dn = sum_dn + flux_dn_clear(jg,jlev)
     208        end do
     209        flux%lw_up_clear(jcol,jlev) = sum_up
     210        flux%lw_dn_clear(jcol,jlev) = sum_dn
     211      end do
     212#endif
     213
    183214      ! Store surface spectral downwelling fluxes
    184215      flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
    … …  
    279310          else
    280311            ! Clear-sky layer: copy over clear-sky values
    281             reflectance(:,jlev) = ref_clear(:,jlev)
    282             transmittance(:,jlev) = trans_clear(:,jlev)
    283             source_up(:,jlev) = source_up_clear(:,jlev)
    284             source_dn(:,jlev) = source_dn_clear(:,jlev)
     312            do jg = 1,ng
     313              reflectance(jg,jlev) = ref_clear(jg,jlev)
     314              transmittance(jg,jlev) = trans_clear(jg,jlev)
     315              source_up(jg,jlev) = source_up_clear(jg,jlev)
     316              source_dn(jg,jlev) = source_dn_clear(jg,jlev)
     317            end do
    285318          end if
    286319        end do
    … …  
    307340       
    308341        ! Store overcast broadband fluxes
    309         flux%lw_up(jcol,:) = sum(flux_up,1)
    310         flux%lw_dn(jcol,:) = sum(flux_dn,1)
     342#ifdef DWD_REDUCTION_OPTIMIZATIONS
     343        sum_aux(:,:) = 0._jprb
     344        do jg = 1, ng
     345          do jlev = 1, nlev+1
     346            sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
     347            sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn(jg,jlev)
     348          end do
     349        end do
     350        flux%lw_up(jcol,:) = sum_aux(:,1)
     351        flux%lw_dn(jcol,:) = sum_aux(:,2)
     352#else
     353        do jlev = 1,nlev+1
     354          sum_up = 0.0_jprb
     355          sum_dn = 0.0_jprb
     356          !$omp simd reduction(+:sum_up, sum_dn)
     357          do jg = 1,ng
     358            sum_up = sum_up + flux_up(jg,jlev)
     359            sum_dn = sum_dn + flux_dn(jg,jlev)
     360          end do
     361          flux%lw_up(jcol,jlev) = sum_up
     362          flux%lw_dn(jcol,jlev) = sum_dn
     363        end do
     364#endif
    311365
    312366        ! Cloudy flux profiles currently assume completely overcast
    313367        ! skies; perform weighted average with clear-sky profile
    314         flux%lw_up(jcol,:) =  total_cloud_cover *flux%lw_up(jcol,:) &
    315              &  + (1.0_jprb - total_cloud_cover)*flux%lw_up_clear(jcol,:)
    316         flux%lw_dn(jcol,:) =  total_cloud_cover *flux%lw_dn(jcol,:) &
    317              &  + (1.0_jprb - total_cloud_cover)*flux%lw_dn_clear(jcol,:)
     368        do jlev = 1,nlev+1
     369          flux%lw_up(jcol,jlev) =  total_cloud_cover *flux%lw_up(jcol,jlev) &
     370             &       + (1.0_jprb - total_cloud_cover)*flux%lw_up_clear(jcol,jlev)
     371          flux%lw_dn(jcol,jlev) =  total_cloud_cover *flux%lw_dn(jcol,jlev) &
     372             &       + (1.0_jprb - total_cloud_cover)*flux%lw_dn_clear(jcol,jlev)
     373        end do
    318374        ! Store surface spectral downwelling fluxes
    319375        flux%lw_dn_surf_g(:,jcol) = total_cloud_cover*flux_dn(:,nlev+1) &
    … …  
    335391        ! No cloud in profile and clear-sky fluxes already
    336392        ! calculated: copy them over
    337         flux%lw_up(jcol,:) = flux%lw_up_clear(jcol,:)
    338         flux%lw_dn(jcol,:) = flux%lw_dn_clear(jcol,:)
     393        do jlev = 1,nlev+1
     394          flux%lw_up(jcol,jlev) = flux%lw_up_clear(jcol,jlev)
     395          flux%lw_dn(jcol,jlev) = flux%lw_dn_clear(jcol,jlev)
     396        end do
    339397        flux%lw_dn_surf_g(:,jcol) = flux%lw_dn_surf_clear_g(:,jcol)
    340398        if (config%do_lw_derivatives) then

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90

    r4773 r4853  
    1818!   2017-10-23  R. Hogan  Renamed single-character variables
    1919
     20#include "ecrad_config.h"
     21
    2022module radiation_mcica_sw
    2123
    … …  
    119121    ! Total cloud cover output from the cloud generator
    120122    real(jprb) :: total_cloud_cover
     123
     124    ! Temporary storage for more efficient summation
     125#ifdef DWD_REDUCTION_OPTIMIZATIONS
     126    real(jprb), dimension(nlev+1,3) :: sum_aux
     127#else
     128    real(jprb) :: sum_up, sum_dn_diff, sum_dn_dir
     129#endif
    121130
    122131    ! Number of g points
    … …  
    175184       
    176185        ! Sum over g-points to compute and save clear-sky broadband
    177         ! fluxes
    178         flux%sw_up_clear(jcol,:) = sum(flux_up,1)
     186        ! fluxes. Note that the built-in "sum" function is very slow,
     187        ! and before being replaced by the alternatives below
     188        ! accounted for around 40% of the total cost of this routine.
     189#ifdef DWD_REDUCTION_OPTIMIZATIONS
     190        ! Optimized summation for the NEC architecture
     191        sum_aux(:,:) = 0.0_jprb
     192        do jg = 1,ng
     193          do jlev = 1,nlev+1
     194            sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
     195            sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_direct(jg,jlev)
     196            sum_aux(jlev,3) = sum_aux(jlev,3) + flux_dn_diffuse(jg,jlev)
     197          end do
     198        end do
     199        flux%sw_up_clear(jcol,:) = sum_aux(:,1)
     200        flux%sw_dn_clear(jcol,:) = sum_aux(:,2) + sum_aux(:,3)
    179201        if (allocated(flux%sw_dn_direct_clear)) then
    180           flux%sw_dn_direct_clear(jcol,:) &
    181                &  = sum(flux_dn_direct,1)
    182           flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
    183                &  + flux%sw_dn_direct_clear(jcol,:)
    184         else
    185           flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
    186                &  + sum(flux_dn_direct,1)
     202          flux%sw_dn_direct_clear(jcol,:) = sum_aux(:,2)
    187203        end if
     204#else
     205        ! Optimized summation for the x86-64 architecture
     206        do jlev = 1,nlev+1
     207          sum_up      = 0.0_jprb
     208          sum_dn_diff = 0.0_jprb
     209          sum_dn_dir  = 0.0_jprb
     210          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
     211          do jg = 1,ng
     212            sum_up      = sum_up      + flux_up(jg,jlev)
     213            sum_dn_diff = sum_dn_diff + flux_dn_diffuse(jg,jlev)
     214            sum_dn_dir  = sum_dn_dir  + flux_dn_direct(jg,jlev)
     215          end do
     216          flux%sw_up_clear(jcol,jlev) = sum_up
     217          flux%sw_dn_clear(jcol,jlev) = sum_dn_diff + sum_dn_dir
     218          if (allocated(flux%sw_dn_direct_clear)) then
     219            flux%sw_dn_direct_clear(jcol,jlev) = sum_dn_dir
     220          end if
     221        end do
     222#endif
     223       
    188224        ! Store spectral downwelling fluxes at surface
    189         flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
    190         flux%sw_dn_direct_surf_clear_g(:,jcol)  = flux_dn_direct(:,nlev+1)
     225        do jg = 1,ng
     226          flux%sw_dn_diffuse_surf_clear_g(jg,jcol) = flux_dn_diffuse(jg,nlev+1)
     227          flux%sw_dn_direct_surf_clear_g(jg,jcol)  = flux_dn_direct(jg,nlev+1)
     228        end do
    191229
    192230        ! Do cloudy-sky calculation
    … …  
    249287            else
    250288              ! Clear-sky layer: copy over clear-sky values
    251               reflectance(:,jlev) = ref_clear(:,jlev)
    252               transmittance(:,jlev) = trans_clear(:,jlev)
    253               ref_dir(:,jlev) = ref_dir_clear(:,jlev)
    254               trans_dir_diff(:,jlev) = trans_dir_diff_clear(:,jlev)
    255               trans_dir_dir(:,jlev) = trans_dir_dir_clear(:,jlev)
     289              do jg = 1,ng
     290                reflectance(jg,jlev) = ref_clear(jg,jlev)
     291                transmittance(jg,jlev) = trans_clear(jg,jlev)
     292                ref_dir(jg,jlev) = ref_dir_clear(jg,jlev)
     293                trans_dir_diff(jg,jlev) = trans_dir_diff_clear(jg,jlev)
     294                trans_dir_dir(jg,jlev) = trans_dir_dir_clear(jg,jlev)
     295              end do
    256296            end if
    257297          end do
    … …  
    264304         
    265305          ! Store overcast broadband fluxes
    266           flux%sw_up(jcol,:) = sum(flux_up,1)
     306#ifdef DWD_REDUCTION_OPTIMIZATIONS
     307          sum_aux(:,:) = 0.0_jprb
     308          do jg = 1,ng
     309            do jlev = 1,nlev+1
     310              sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
     311              sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_direct(jg,jlev)
     312              sum_aux(jlev,3) = sum_aux(jlev,3) + flux_dn_diffuse(jg,jlev)
     313            end do
     314          end do
     315          flux%sw_up(jcol,:) = sum_aux(:,1)
     316          flux%sw_dn(jcol,:) = sum_aux(:,2) + sum_aux(:,3)
    267317          if (allocated(flux%sw_dn_direct)) then
    268             flux%sw_dn_direct(jcol,:) = sum(flux_dn_direct,1)
    269             flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
    270                  &  + flux%sw_dn_direct(jcol,:)
    271           else
    272             flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
    273                  &  + sum(flux_dn_direct,1)
     318            flux%sw_dn_direct(jcol,:) = sum_aux(:,2)
    274319          end if
    275 
     320#else
     321          do jlev = 1,nlev+1
     322            sum_up      = 0.0_jprb
     323            sum_dn_diff = 0.0_jprb
     324            sum_dn_dir  = 0.0_jprb
     325            !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
     326            do jg = 1,ng
     327              sum_up      = sum_up      + flux_up(jg,jlev)
     328              sum_dn_diff = sum_dn_diff + flux_dn_diffuse(jg,jlev)
     329              sum_dn_dir  = sum_dn_dir  + flux_dn_direct(jg,jlev)
     330            end do
     331            flux%sw_up(jcol,jlev) = sum_up
     332            flux%sw_dn(jcol,jlev) = sum_dn_diff + sum_dn_dir
     333            if (allocated(flux%sw_dn_direct)) then
     334              flux%sw_dn_direct(jcol,jlev) = sum_dn_dir
     335            end if
     336          end do
     337#endif
     338         
    276339          ! Cloudy flux profiles currently assume completely overcast
    277340          ! skies; perform weighted average with clear-sky profile
    278           flux%sw_up(jcol,:) =  total_cloud_cover *flux%sw_up(jcol,:) &
    279                &  + (1.0_jprb - total_cloud_cover)*flux%sw_up_clear(jcol,:)
    280           flux%sw_dn(jcol,:) =  total_cloud_cover *flux%sw_dn(jcol,:) &
    281                &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_clear(jcol,:)
    282           if (allocated(flux%sw_dn_direct)) then
    283             flux%sw_dn_direct(jcol,:) = total_cloud_cover *flux%sw_dn_direct(jcol,:) &
    284                  &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_clear(jcol,:)
    285           end if
     341          do jlev = 1, nlev+1
     342            flux%sw_up(jcol,jlev) =  total_cloud_cover *flux%sw_up(jcol,jlev) &
     343                 &     + (1.0_jprb - total_cloud_cover)*flux%sw_up_clear(jcol,jlev)
     344            flux%sw_dn(jcol,jlev) =  total_cloud_cover *flux%sw_dn(jcol,jlev) &
     345                 &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_clear(jcol,jlev)
     346            if (allocated(flux%sw_dn_direct)) then
     347              flux%sw_dn_direct(jcol,jlev) = total_cloud_cover *flux%sw_dn_direct(jcol,jlev) &
     348                   &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_clear(jcol,jlev)
     349            end if
     350          end do
    286351          ! Likewise for surface spectral fluxes
    287           flux%sw_dn_diffuse_surf_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
    288           flux%sw_dn_direct_surf_g(:,jcol)  = flux_dn_direct(:,nlev+1)
    289           flux%sw_dn_diffuse_surf_g(:,jcol) = total_cloud_cover *flux%sw_dn_diffuse_surf_g(:,jcol) &
    290                &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_diffuse_surf_clear_g(:,jcol)
    291           flux%sw_dn_direct_surf_g(:,jcol) = total_cloud_cover *flux%sw_dn_direct_surf_g(:,jcol) &
    292                &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_surf_clear_g(:,jcol)
    293          
     352          do jg = 1,ng
     353            flux%sw_dn_diffuse_surf_g(jg,jcol) = flux_dn_diffuse(jg,nlev+1)
     354            flux%sw_dn_direct_surf_g(jg,jcol)  = flux_dn_direct(jg,nlev+1)
     355            flux%sw_dn_diffuse_surf_g(jg,jcol) = total_cloud_cover *flux%sw_dn_diffuse_surf_g(jg,jcol) &
     356                 &                 + (1.0_jprb - total_cloud_cover)*flux%sw_dn_diffuse_surf_clear_g(jg,jcol)
     357            flux%sw_dn_direct_surf_g(jg,jcol)  = total_cloud_cover *flux%sw_dn_direct_surf_g(jg,jcol) &
     358                 &                 + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_surf_clear_g(jg,jcol)
     359          end do
     360
    294361        else
    295362          ! No cloud in profile and clear-sky fluxes already
    296363          ! calculated: copy them over
    297           flux%sw_up(jcol,:) = flux%sw_up_clear(jcol,:)
    298           flux%sw_dn(jcol,:) = flux%sw_dn_clear(jcol,:)
    299           if (allocated(flux%sw_dn_direct)) then
    300             flux%sw_dn_direct(jcol,:) = flux%sw_dn_direct_clear(jcol,:)
    301           end if
    302           flux%sw_dn_diffuse_surf_g(:,jcol) = flux%sw_dn_diffuse_surf_clear_g(:,jcol)
    303           flux%sw_dn_direct_surf_g(:,jcol)  = flux%sw_dn_direct_surf_clear_g(:,jcol)
     364          do jlev = 1, nlev+1
     365            flux%sw_up(jcol,jlev) = flux%sw_up_clear(jcol,jlev)
     366            flux%sw_dn(jcol,jlev) = flux%sw_dn_clear(jcol,jlev)
     367            if (allocated(flux%sw_dn_direct)) then
     368              flux%sw_dn_direct(jcol,jlev) = flux%sw_dn_direct_clear(jcol,jlev)
     369            end if
     370          end do
     371          do jg = 1,ng
     372            flux%sw_dn_diffuse_surf_g(jg,jcol) = flux%sw_dn_diffuse_surf_clear_g(jg,jcol)
     373            flux%sw_dn_direct_surf_g(jg,jcol)  = flux%sw_dn_direct_surf_clear_g(jg,jcol)
     374          end do
    304375
    305376        end if ! Cloud is present in profile
    … …  
    307378      else
    308379        ! Set fluxes to zero if sun is below the horizon
    309         flux%sw_up(jcol,:) = 0.0_jprb
    310         flux%sw_dn(jcol,:) = 0.0_jprb
    311         if (allocated(flux%sw_dn_direct)) then
    312           flux%sw_dn_direct(jcol,:) = 0.0_jprb
    313         end if
    314         flux%sw_up_clear(jcol,:) = 0.0_jprb
    315         flux%sw_dn_clear(jcol,:) = 0.0_jprb
    316         if (allocated(flux%sw_dn_direct_clear)) then
    317           flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
    318         end if
    319         flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
    320         flux%sw_dn_direct_surf_g(:,jcol)  = 0.0_jprb
    321         flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
    322         flux%sw_dn_direct_surf_clear_g(:,jcol)  = 0.0_jprb
     380        do jlev = 1, nlev+1
     381          flux%sw_up(jcol,jlev) = 0.0_jprb
     382          flux%sw_dn(jcol,jlev) = 0.0_jprb
     383          if (allocated(flux%sw_dn_direct)) then
     384            flux%sw_dn_direct(jcol,jlev) = 0.0_jprb
     385          end if
     386          flux%sw_up_clear(jcol,jlev) = 0.0_jprb
     387          flux%sw_dn_clear(jcol,jlev) = 0.0_jprb
     388          if (allocated(flux%sw_dn_direct_clear)) then
     389            flux%sw_dn_direct_clear(jcol,jlev) = 0.0_jprb
     390          end if
     391        end do
     392        do jg = 1,ng
     393          flux%sw_dn_diffuse_surf_g(jg,jcol) = 0.0_jprb
     394          flux%sw_dn_direct_surf_g(jg,jcol)  = 0.0_jprb
     395          flux%sw_dn_diffuse_surf_clear_g(jg,jcol) = 0.0_jprb
     396          flux%sw_dn_direct_surf_clear_g(jg,jcol)  = 0.0_jprb
     397        end do
    323398      end if ! Sun above horizon
    324399

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_save.F90

    r4773 r4853  
    8686    end if
    8787
    88     if (config%i_gas_model == IGasModelMonochromatic &
     88    if (config%i_gas_model_lw == IGasModelMonochromatic &
    8989         .and. config%mono_lw_wavelength > 0.0_jprb) then
    9090      lw_units_str = 'W m-3'
    … …  
    515515    end if
    516516
    517     if (config%i_gas_model == IGasModelMonochromatic &
     517    if (config%i_gas_model_lw == IGasModelMonochromatic &
    518518         .and. config%mono_lw_wavelength > 0.0_jprb) then
    519519      lw_units_str = 'W m-3'

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90

    r4773 r4853  
    8787    use radiation_constants, only      : Pi, GasConstantDryAir, &
    8888         &                               AccelDueToGravity
    89     USE mod_phys_lmdz_para
    9089
    9190    implicit none
    … …  
    327326      write(nulout,'(a)',advance='no') '  Processing columns'
    328327    end if
    329 
    330       print*,'Dans radiation_spartacus, mpi_rank, omp_rank, &
    331               size, chape inv_cloud = ',&
    332               mpi_rank, omp_rank, &
    333               shape(cloud%inv_cloud_effective_size), &
    334               size(cloud%inv_cloud_effective_size)
    335328
    336329    ! Main loop over columns
    … …  
    504497               &  .not. (nreg == 2 .and. cloud%fraction(jcol,jlev) &
    505498               &  > 1.0-config%cloud_fraction_threshold)) then
    506 !           print*,' Dans radiation_spartacus mpi_rank, omp_rank, jcol, jlev, &
    507 !            &     cloud%inv_cloud_effective_size =', mpi_rank, &
    508 !            &     omp_rank, jcol, jlev, &
    509 !            &     cloud%inv_cloud_effective_size(jcol,jlev)
    510499            if (cloud%inv_cloud_effective_size(jcol,jlev) > 0.0_jprb) then
    511500              ! 3D effects are only simulated if

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_spectral_definition.F90

    r4773 r4853  
    1515!
    1616
     17#include "ecrad_config.h"
     18
    1719module radiation_spectral_definition
    1820
    … …  
    2325  public
    2426
    25   real(jprb), parameter :: SolarReferenceTemperature = 5777.0_jprb ! K
     27  real(jprb), parameter :: SolarReferenceTemperature       = 5777.0_jprb ! K
    2628  real(jprb), parameter :: TerrestrialReferenceTemperature = 273.15_jprb ! K
    2729
    … …  
    8587  subroutine read_spectral_definition(this, file)
    8688
    87     use easy_netcdf, only : netcdf_file
     89#ifdef EASY_NETCDF_READ_MPI
     90    use easy_netcdf_read_mpi, only : netcdf_file
     91#else
     92    use easy_netcdf,          only : netcdf_file
     93#endif
    8894    use yomhook,     only : lhook, dr_hook, jphook
    8995
    … …  
    132138
    133139  !---------------------------------------------------------------------
    134   ! Store a simple band description by copying over the lower and
    135   ! upper wavenumbers of each band
    136   subroutine allocate_bands_only(this, wavenumber1, wavenumber2)
     140  ! Store a simple band description by copying over the reference
     141  ! temperature and the lower and upper wavenumbers of each band
     142  subroutine allocate_bands_only(this, reference_temperature, wavenumber1, wavenumber2)
    137143
    138144    use yomhook,     only : lhook, dr_hook, jphook
    139145
    140146    class(spectral_definition_type), intent(inout) :: this
    141     real(jprb),        dimension(:), intent(in)    :: wavenumber1, wavenumber2
     147    real(jprb),                      intent(in)    :: reference_temperature    ! K
     148    real(jprb),        dimension(:), intent(in)    :: wavenumber1, wavenumber2 ! cm-1
    142149
    143150    real(jphook) :: hook_handle
    … …  
    152159    this%wavenumber1_band = wavenumber1
    153160    this%wavenumber2_band = wavenumber2
    154 
     161    this%reference_temperature = reference_temperature
     162   
    155163    if (lhook) call dr_hook('radiation_spectral_definition:allocate_bands_only',1,hook_handle)
    156164
    … …  
    167175    this%ng    = 0
    168176    this%nband = 0
     177    this%reference_temperature = -1.0_jprb
    169178
    170179    if (allocated(this%wavenumber1))      deallocate(this%wavenumber1)

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90

    r4773 r4853  
    170170    logical :: is_clear_sky_layer(0:nlev+1)
    171171
     172    ! Temporaries to speed up summations
     173    real(jprb) :: sum_dn, sum_up
     174   
    172175    ! Index of the highest cloudy layer
    173176    integer :: i_cloud_top
    … …  
    261264      if (config%do_clear) then
    262265        ! Sum over g-points to compute broadband fluxes
    263         flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
    264         flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
     266        do jlev = 1,nlev+1
     267          sum_up = 0.0_jprb
     268          sum_dn = 0.0_jprb
     269          !$omp simd reduction(+:sum_up, sum_dn)
     270          do jg = 1,ng
     271            sum_up = sum_up + flux_up_clear(jg,jlev)
     272            sum_dn = sum_dn + flux_dn_clear(jg,jlev)
     273          end do
     274          flux%lw_up_clear(jcol,jlev) = sum_up
     275          flux%lw_dn_clear(jcol,jlev) = sum_dn
     276        end do
     277
    265278        ! Store surface spectral downwelling fluxes / TOA upwelling
    266         flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
    267         flux%lw_up_toa_clear_g (:,jcol) = flux_up_clear(:,1)
     279        do jg = 1,ng
     280          flux%lw_dn_surf_clear_g(jg,jcol) = flux_dn_clear(jg,nlev+1)
     281          flux%lw_up_toa_clear_g (jg,jcol) = flux_up_clear(jg,1)
     282        end do
    268283        ! Save the spectral fluxes if required
    269284        if (config%do_save_spectral_flux) then
    … …  
    453468          end if
    454469        else
    455           flux%lw_dn(jcol,:) = sum(flux_dn_clear(:,jlev))
     470          sum_dn = 0.0_jprb
     471          !$omp simd reduction(+:sum_dn)
     472          do jg = 1,ng
     473            sum_dn = sum_dn + flux_dn_clear(jg,jlev)
     474          end do
     475          flux%lw_dn(jcol,jlev) = sum_dn
    456476          if (config%do_save_spectral_flux) then
    457477            call indexed_sum(flux_dn_clear(:,jlev), &
    … …  
    470490           &  + total_albedo(:,1,i_cloud_top)*flux_dn_clear(:,i_cloud_top)
    471491      flux_up(:,2:) = 0.0_jprb
    472       flux%lw_up(jcol,i_cloud_top) = sum(flux_up(:,1))
     492
     493      sum_up = 0.0_jprb
     494      !$omp simd reduction(+:sum_up)
     495      do jg = 1,ng
     496        sum_up = sum_up + flux_up(jg,1)
     497      end do
     498      flux%lw_up(jcol,i_cloud_top) = sum_up
     499
    473500      if (config%do_save_spectral_flux) then
    474501        call indexed_sum(flux_up(:,1), &
    … …  
    478505      do jlev = i_cloud_top-1,1,-1
    479506        flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev)
    480         flux%lw_up(jcol,jlev) = sum(flux_up(:,1))
     507        sum_up = 0.0_jprb
     508        !$omp simd reduction(+:sum_up)
     509        do jg = 1,ng
     510          sum_up = sum_up + flux_up(jg,1)
     511        end do
     512        flux%lw_up(jcol,jlev) = sum_up
    481513        if (config%do_save_spectral_flux) then
    482514          call indexed_sum(flux_up(:,1), &
    … …  
    528560
    529561        ! Store the broadband fluxes
    530         flux%lw_up(jcol,jlev+1) = sum(sum(flux_up,1))
    531         flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn,1))
     562        sum_up = 0.0_jprb
     563        sum_dn = 0.0_jprb
     564        do jreg = 1,nregions
     565          !$omp simd reduction(+:sum_up, sum_dn)
     566          do jg = 1,ng
     567            sum_up = sum_up + flux_up(jg,jreg)
     568            sum_dn = sum_dn + flux_dn(jg,jreg)
     569          end do
     570        end do
     571        flux%lw_up(jcol,jlev+1) = sum_up
     572        flux%lw_dn(jcol,jlev+1) = sum_dn
    532573
    533574        ! Save the spectral fluxes if required

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90

    r4773 r4853  
    7474    ! Gas and aerosol optical depth, single-scattering albedo and
    7575    ! asymmetry factor at each shortwave g-point
    76 !    real(jprb), intent(in), dimension(istartcol:iendcol,nlev,config%n_g_sw) :: &
    77     real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: &
    78          &  od, ssa, g
     76    real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) &
     77         &  :: od, ssa, g
    7978
    8079    ! Cloud and precipitation optical depth, single-scattering albedo and
    8180    ! asymmetry factor in each shortwave band
    82     real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
    83          &  od_cloud, ssa_cloud, g_cloud
     81    real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) &
     82         &  :: od_cloud, ssa_cloud, g_cloud
    8483
    8584    ! Optical depth, single scattering albedo and asymmetry factor in
    … …  
    9291    ! flux into a plane perpendicular to the incoming radiation at
    9392    ! top-of-atmosphere in each of the shortwave g points
    94     real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
    95          &  albedo_direct, albedo_diffuse, incoming_sw
     93    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) &
     94         &  :: albedo_direct, albedo_diffuse, incoming_sw
    9695
    9796    ! Output
    … …  
    166165    real(jprb) :: scat_od, scat_od_cloud
    167166
     167    ! Temporaries to speed up summations
     168    real(jprb) :: sum_dn_diff, sum_dn_dir, sum_up
     169
     170    ! Local cosine of solar zenith angle
    168171    real(jprb) :: mu0
    169172
    … …  
    444447      end if
    445448     
    446       ! Store the TOA broadband fluxes
    447       flux%sw_up(jcol,1) = sum(sum(flux_up,1))
    448       flux%sw_dn(jcol,1) = mu0 * sum(sum(direct_dn,1))
     449      ! Store the TOA broadband fluxes, noting that there is no
     450      ! diffuse downwelling at TOA. The intrinsic "sum" command has
     451      ! been found to be very slow; better performance is found on
     452      ! x86-64 architecture with explicit loops and the "omp simd
     453      ! reduction" directive.
     454      sum_up     = 0.0_jprb
     455      sum_dn_dir = 0.0_jprb
     456      do jreg = 1,nregions
     457        !$omp simd reduction(+:sum_up, sum_dn_dir)
     458        do jg = 1,ng
     459          sum_up     = sum_up     + flux_up(jg,jreg)
     460          sum_dn_dir = sum_dn_dir + direct_dn(jg,jreg)
     461        end do
     462      end do
     463      flux%sw_up(jcol,1) = sum_up
     464      flux%sw_dn(jcol,1) = mu0 * sum_dn_dir
    449465      if (allocated(flux%sw_dn_direct)) then
    450466        flux%sw_dn_direct(jcol,1) = flux%sw_dn(jcol,1)
    451467      end if
    452468      if (config%do_clear) then
    453         flux%sw_up_clear(jcol,1) = sum(flux_up_clear)
    454         flux%sw_dn_clear(jcol,1) = mu0 * sum(direct_dn_clear)
     469        sum_up     = 0.0_jprb
     470        sum_dn_dir = 0.0_jprb
     471        !$omp simd reduction(+:sum_up, sum_dn_dir)
     472        do jg = 1,ng
     473          sum_up     = sum_up     + flux_up_clear(jg)
     474          sum_dn_dir = sum_dn_dir + direct_dn_clear(jg)
     475        end do
     476        flux%sw_up_clear(jcol,1) = sum_up
     477        flux%sw_dn_clear(jcol,1) = mu0 * sum_dn_dir
    455478        if (allocated(flux%sw_dn_direct_clear)) then
    456479          flux%sw_dn_direct_clear(jcol,1) = flux%sw_dn_clear(jcol,1)
    … …  
    467490             &           config%i_spec_from_reordered_g_sw, &
    468491             &           flux%sw_dn_band(:,jcol,1))
    469         flux%sw_dn_band(:,jcol,1) = &
    470              &  mu0 * flux%sw_dn_band(:,jcol,1)
     492        flux%sw_dn_band(:,jcol,1) = mu0 * flux%sw_dn_band(:,jcol,1)
    471493        if (allocated(flux%sw_dn_direct_band)) then
    472494          flux%sw_dn_direct_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1)
    … …  
    549571               ! nothing to do
    550572
    551         ! Store the broadband fluxes
    552         flux%sw_up(jcol,jlev+1) = sum(sum(flux_up,1))
     573        ! Store the broadband fluxes. The intrinsic "sum" command has
     574        ! been found to be very slow; better performance is found on
     575        ! x86-64 architecture with explicit loops and the "omp simd
     576        ! reduction" directive.
     577        sum_up      = 0.0_jprb
     578        sum_dn_dir  = 0.0_jprb
     579        sum_dn_diff = 0.0_jprb
     580        do jreg = 1,nregions
     581          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
     582          do jg = 1,ng
     583            sum_up      = sum_up      + flux_up(jg,jreg)
     584            sum_dn_diff = sum_dn_diff + flux_dn(jg,jreg)
     585            sum_dn_dir  = sum_dn_dir  + direct_dn(jg,jreg)
     586          end do
     587        end do
     588        flux%sw_up(jcol,jlev+1) = sum_up
     589        flux%sw_dn(jcol,jlev+1) = mu0 * sum_dn_dir + sum_dn_diff
    553590        if (allocated(flux%sw_dn_direct)) then
    554           flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1))
    555           flux%sw_dn(jcol,jlev+1) &
    556                &  = flux%sw_dn_direct(jcol,jlev+1) + sum(sum(flux_dn,1))
    557         else
    558           flux%sw_dn(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1)) + sum(sum(flux_dn,1))   
     591          flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum_dn_dir
    559592        end if
    560593        if (config%do_clear) then
    561           flux%sw_up_clear(jcol,jlev+1) = sum(flux_up_clear)
     594          sum_up      = 0.0_jprb
     595          sum_dn_dir  = 0.0_jprb
     596          sum_dn_diff = 0.0_jprb
     597          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
     598          do jg = 1,ng
     599            sum_up      = sum_up      + flux_up_clear(jg)
     600            sum_dn_diff = sum_dn_diff + flux_dn_clear(jg)
     601            sum_dn_dir  = sum_dn_dir  + direct_dn_clear(jg)
     602          end do
     603          flux%sw_up_clear(jcol,jlev+1) = sum_up
     604          flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum_dn_dir + sum_dn_diff
    562605          if (allocated(flux%sw_dn_direct_clear)) then
    563             flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear)
    564             flux%sw_dn_clear(jcol,jlev+1) &
    565                  &  = flux%sw_dn_direct_clear(jcol,jlev+1) + sum(flux_dn_clear)
    566           else
    567             flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear) &
    568                  &  + sum(flux_dn_clear)
     606            flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum_dn_dir
    569607          end if
    570608        end if
    … …  
    605643          end if
    606644        end if
    607 
    608645      end do ! Final loop over levels
    609646     

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_two_stream.F90

    r4773 r4853  
    2222!   2023-09-28  R Hogan  Increased security for single-precision SW "k"
    2323
     24#include "ecrad_config.h"
     25
    2426module radiation_two_stream
    2527
    … …  
    4345
    4446contains
    45 
    46 #ifdef FAST_EXPONENTIAL
    47   !---------------------------------------------------------------------
    48   ! Fast exponential for negative arguments: a Pade approximant that
    49   ! doesn't go negative for negative arguments, applied to arg/8, and
    50   ! the result is then squared three times
    51   elemental function exp_fast(arg) result(ex)
    52     real(jprd) :: arg, ex
    53     ex = 1.0_jprd / (1.0_jprd + arg*(-0.125_jprd &
    54          + arg*(0.0078125_jprd - 0.000325520833333333_jprd * arg)))
    55     ex = ex*ex
    56     ex = ex*ex
    57     ex = ex*ex
    58   end function exp_fast
    59 #else
    60 #define exp_fast exp
    61 #endif
    6247
    6348  !---------------------------------------------------------------------
    … …  
    214199        k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
    215200             1.0e-12_jprd)) ! Eq 18 of Meador & Weaver (1980)
    216         exponential = exp_fast(-k_exponent*od(jg))
     201        exponential = exp(-k_exponent*od(jg))
    217202        exponential2 = exponential*exponential
    218203        reftrans_factor = 1.0 / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2)
    … …  
    306291#endif
    307292
    308 ! Added for DWD (2020)
    309 !NEC$ shortloop
    310293    do jg = 1, ng
    311294      factor = (LwDiffusivityWP * 0.5_jprb) * ssa(jg)
    … …  
    315298           1.0e-12_jprb)) ! Eq 18 of Meador & Weaver (1980)
    316299      if (od(jg) > 1.0e-3_jprb) then
    317         exponential = exp_fast(-k_exponent*od(jg))
     300        exponential = exp(-k_exponent*od(jg))
    318301        exponential2 = exponential*exponential
    319302        reftrans_factor = 1.0_jprb / (k_exponent + gamma1 + (k_exponent - gamma1)*exponential2)
    … …  
    392375#endif
    393376
    394     transmittance = exp_fast(-LwDiffusivityWP*od)
    395 
    396 ! Added for DWD (2020)
    397 !NEC$ shortloop
     377#ifndef DWD_TWO_STREAM_OPTIMIZATIONS
     378    transmittance = exp(-LwDiffusivityWP*od)
     379#endif
     380
    398381    do jg = 1, ng
    399382      ! Compute upward and downward emission assuming the Planck
    … …  
    401384      ! (e.g. Wiscombe , JQSRT 1976).
    402385      coeff = LwDiffusivityWP*od(jg)
     386#ifdef DWD_TWO_STREAM_OPTIMIZATIONS
     387      transmittance(jg) = exp(-coeff)
     388#endif
    403389      if (od(jg) > 1.0e-3_jprb) then
    404390        ! Simplified from calc_reflectance_transmittance_lw above
    … …  
    516502        k_gamma4 = k_exponent*gamma4
    517503        ! Check for mu0 <= 0!
    518         exponential0 = exp_fast(-od_over_mu0)
     504        exponential0 = exp(-od_over_mu0)
    519505        trans_dir_dir(jg) = exponential0
    520         exponential = exp_fast(-k_exponent*od(jg))
     506        exponential = exp(-k_exponent*od(jg))
    521507       
    522508        exponential2 = exponential*exponential
    … …  
    609595    ! The three transfer coefficients from the two-stream
    610596    ! differentiatial equations
     597#ifndef DWD_TWO_STREAM_OPTIMIZATIONS
    611598    real(jprb), dimension(ng) :: gamma1, gamma2, gamma3, gamma4
    612599    real(jprb), dimension(ng) :: alpha1, alpha2, k_exponent
    613600    real(jprb), dimension(ng) :: exponential ! = exp(-k_exponent*od)
     601#else
     602    real(jprb) :: gamma1, gamma2, gamma3, gamma4
     603    real(jprb) :: alpha1, alpha2, k_exponent
     604    real(jprb) :: exponential ! = exp(-k_exponent*od)
     605#endif
    614606   
    615607    real(jprb) :: reftrans_factor, factor
    … …  
    625617#endif
    626618
     619#ifndef DWD_TWO_STREAM_OPTIMIZATIONS
    627620    ! GCC 9.3 strange error: intermediate values of ~ -8000 cause a
    628621    ! FPE when vectorizing exp(), but not in non-vectorized loop, nor
    629622    ! with larger negative values!
    630623    trans_dir_dir = max(-max(od * (1.0_jprb/mu0), 0.0_jprb),-1000.0_jprb)
    631     trans_dir_dir = exp_fast(trans_dir_dir)
    632 
    633 ! Added for DWD (2020)
    634 !NEC$ shortloop
     624    trans_dir_dir = exp(trans_dir_dir)
     625
    635626    do jg = 1, ng
    636627
    … …  
    659650    end do
    660651
    661     exponential = exp_fast(-k_exponent*od)
    662 
    663 !NEC$ shortloop
     652    exponential = exp(-k_exponent*od)
     653
    664654    do jg = 1, ng
    665655      k_mu0 = k_exponent(jg)*mu0
    … …  
    705695      trans_dir_diff(jg) = max(0.0_jprb, min(trans_dir_diff(jg), mu0*(1.0_jprb-trans_dir_dir(jg))-ref_dir(jg)))
    706696    end do
    707    
     697
     698#else
     699    ! GPU-capable and vector-optimized version for ICON
     700    do jg = 1, ng
     701
     702      trans_dir_dir(jg) = max(-max(od(jg) * (1.0_jprb/mu0),0.0_jprb),-1000.0_jprb)
     703      trans_dir_dir(jg) = exp(trans_dir_dir(jg))
     704
     705      ! Zdunkowski "PIFM" (Zdunkowski et al., 1980; Contributions to
     706      ! Atmospheric Physics 53, 147-66)
     707      factor = 0.75_jprb*asymmetry(jg)
     708
     709      gamma1 = 2.0_jprb  - ssa(jg) * (1.25_jprb + factor)
     710      gamma2 = ssa(jg) * (0.75_jprb - factor)
     711      gamma3 = 0.5_jprb  - mu0*factor
     712      gamma4 = 1.0_jprb - gamma3
     713
     714      alpha1 = gamma1*gamma4 + gamma2*gamma3 ! Eq. 16
     715      alpha2 = gamma1*gamma3 + gamma2*gamma4 ! Eq. 17
     716#ifdef PARKIND1_SINGLE
     717      k_exponent = sqrt(max((gamma1 - gamma2) * (gamma1 + gamma2), 1.0e-6_jprb))  ! Eq 18
     718#else
     719      k_exponent = sqrt(max((gamma1 - gamma2) * (gamma1 + gamma2), 1.0e-12_jprb)) ! Eq 18
     720#endif
     721
     722      exponential = exp(-k_exponent*od(jg))
     723
     724      k_mu0 = k_exponent*mu0
     725      one_minus_kmu0_sqr = 1.0_jprb - k_mu0*k_mu0
     726      k_gamma3 = k_exponent*gamma3
     727      k_gamma4 = k_exponent*gamma4
     728      exponential2 = exponential*exponential
     729      k_2_exponential = 2.0_jprb * k_exponent * exponential
     730      reftrans_factor = 1.0_jprb / (k_exponent + gamma1 + (k_exponent - gamma1)*exponential2)
     731       
     732      ! Meador & Weaver (1980) Eq. 25
     733      ref_diff(jg) = gamma2 * (1.0_jprb - exponential2) * reftrans_factor
     734       
     735      ! Meador & Weaver (1980) Eq. 26
     736      trans_diff(jg) = k_2_exponential * reftrans_factor
     737       
     738      ! Here we need mu0 even though it wasn't in Meador and Weaver
     739      ! because we are assuming the incoming direct flux is defined to
     740      ! be the flux into a plane perpendicular to the direction of the
     741      ! sun, not into a horizontal plane
     742      reftrans_factor = mu0 * ssa(jg) * reftrans_factor &
     743            &  / merge(one_minus_kmu0_sqr, epsilon(1.0_jprb), abs(one_minus_kmu0_sqr) > epsilon(1.0_jprb))
     744     
     745      ! Meador & Weaver (1980) Eq. 14, multiplying top & bottom by
     746      ! exp(-k_exponent*od) in case of very high optical depths
     747      ref_dir(jg) = reftrans_factor &
     748           &  * ( (1.0_jprb - k_mu0) * (alpha2 + k_gamma3) &
     749           &     -(1.0_jprb + k_mu0) * (alpha2 - k_gamma3)*exponential2 &
     750           &     -k_2_exponential*(gamma3 - alpha2*mu0)*trans_dir_dir(jg) )
     751       
     752      ! Meador & Weaver (1980) Eq. 15, multiplying top & bottom by
     753      ! exp(-k_exponent*od), minus the 1*exp(-od/mu0) term
     754      ! representing direct unscattered transmittance.
     755      trans_dir_diff(jg) = reftrans_factor * ( k_2_exponential*(gamma4 + alpha1*mu0) &
     756           & - trans_dir_dir(jg) &
     757           & * ( (1.0_jprb + k_mu0) * (alpha1 + k_gamma4) &
     758           &    -(1.0_jprb - k_mu0) * (alpha1 - k_gamma4) * exponential2) )
     759
     760      ! Final check that ref_dir + trans_dir_diff <= 1
     761      ref_dir(jg)        = max(0.0_jprb, min(ref_dir(jg), mu0*(1.0_jprb-trans_dir_dir(jg))))
     762      trans_dir_diff(jg) = max(0.0_jprb, min(trans_dir_diff(jg), mu0*(1.0_jprb-trans_dir_dir(jg))-ref_dir(jg)))
     763
     764    end do
     765#endif
     766
    708767#ifdef DO_DR_HOOK_TWO_STREAM
    709768    if (lhook) call dr_hook('radiation_two_stream:calc_ref_trans_sw',1,hook_handle)
    … …  
    757816      k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
    758817           &       1.0e-12_jprd)) ! Eq 18
    759       exponential = exp_fast(-k_exponent*od(jg))
     818      exponential = exp(-k_exponent*od(jg))
    760819      exponential2 = exponential*exponential
    761820      k_2_exponential = 2.0_jprd * k_exponent * exponential
    … …  
    766825      ! Until 1.1.8, used LwDiffusivity instead of 2.0, although the
    767826      ! effect is very small
    768       !      frac_scat_diffuse(jg) = 1.0_jprb - min(1.0_jprb,exp_fast(-LwDiffusivity*od(jg)) &
     827      !      frac_scat_diffuse(jg) = 1.0_jprb - min(1.0_jprb,exp(-LwDiffusivity*od(jg)) &
    769828      !           &  / max(1.0e-8_jprb, k_2_exponential * reftrans_factor))
    770829      frac_scat_diffuse(jg) = 1.0_jprb &
    771            &  - min(1.0_jprb,exp_fast(-2.0_jprb*od(jg)) &
     830           &  - min(1.0_jprb,exp(-2.0_jprb*od(jg)) &
    772831           &  / max(1.0e-8_jprb, k_2_exponential * reftrans_factor))
    773832    end do

  • LMDZ6/trunk/libf/phylmd/ecrad/test/ifs/Makefile

    r4773 r4853  
    1515CONFIG = configCY49R1.nam
    1616CONFIG_ECCKD = configCY49R1_ecckd.nam
     17CONFIG_MIXED = configCY49R1_mixed.nam
    1718
    1819# Typing "make" will run radiation scheme on IFS profiles
    … …  
    111112test_ecckd_tc:
    112113        $(DRIVER) $(CONFIG_ECCKD) $(INPUT).nc $(INPUT)_ecckd_tc_out.nc
     114
     115test_mixed_gas:
     116        $(DRIVER) $(CONFIG_MIXED) $(INPUT).nc $(INPUT)_sw_ecckd_lw_ecckd_out.nc
     117        $(CHANGENAM) $(CONFIG_MIXED) config_mix.nam lw_gas_model_name='"RRTMG-IFS"' do_cloud_aerosol_per_lw_g_point=false
     118        $(DRIVER) config_mix.nam $(INPUT).nc $(INPUT)_sw_ecckd_lw_rrtmg_out.nc
     119        $(CHANGENAM) $(CONFIG_MIXED) config_mix.nam sw_gas_model_name='"RRTMG-IFS"' do_cloud_aerosol_per_sw_g_point=false
     120        $(DRIVER) config_mix.nam $(INPUT).nc $(INPUT)_sw_rrtmg_lw_ecckd_out.nc
     121        $(CHANGENAM) $(CONFIG_MIXED) config_mix.nam sw_gas_model_name='"RRTMG-IFS"' lw_gas_model_name='"RRTMG-IFS"' do_cloud_aerosol_per_lw_g_point=false do_cloud_aerosol_per_sw_g_point=false
     122        $(DRIVER) config_mix.nam $(INPUT).nc $(INPUT)_sw_rrtmg_lw_rrtmg_out.nc
    113123
    114124# ecCKD with no aerosols

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r4843 r4853  
    43644364                  flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
    43654365                  pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
    4366                   tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
    4367                   tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
    4368                   tausum_aero, drytausum_aero, tau3d_aero)
     4366                  tr_seri, mass_solu_aero, mass_solu_aero_pi) 
    43694367#else
    43704368                abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'

  • LMDZ6/trunk/libf/phylmd/radlwsw_m.F90

    r4790 r4853  
    13541354      endif
    13551355! lldebug_for_offline
    1356  
     1356
     1357  if (namelist_ecrad_file.eq.'namelist_ecrad') then
     1358      print*,' 1er apell Ecrad : ok_3Deffect, namelist_ecrad_file = ', &
     1359          ok_3Deffect, namelist_ecrad_file   
    13571360      CALL RADIATION_SCHEME &
    13581361      & (ist, iend, klon, klev, naero_grp, NSW, &
    … …  
    13961399      & PSFSWDIF, PSFSWDIR, &
    13971400      & cloud_cover_sw)
     1401  else
     1402   print*,' 2e apell Ecrad : ok_3Deffect, namelist_ecrad_file = ', &
     1403          ok_3Deffect, namelist_ecrad_file       
     1404   CALL RADIATION_SCHEME_S2 &
     1405      & (ist, iend, klon, klev, naero_grp, NSW, &
     1406      & namelist_ecrad_file, ok_3Deffect, &
     1407      & debut, ok_volcan, flag_aerosol_strat, &
     1408      & day_cur, current_time, &
     1409!       Cste solaire/(d_Terre-Soleil)**2
     1410      & SOLARIRAD, &
     1411!       Cos(angle zin), temp sol             
     1412      & rmu0, tsol, &
     1413!       Albedo diffuse et directe
     1414      & PALBD_NEW,PALBP_NEW, &
     1415!       Emessivite : PEMIS_WINDOW (???), &
     1416      & ZEMIS, ZEMISW, &
     1417!       longitude(rad), sin(latitude), PMASQ_ ???
     1418      & ZGELAM, ZGEMU, &
     1419!       Temp et pres aux interf, vapeur eau, Satur spec humid
     1420      & paprs_i, ZTH_i, q_i, qsat_i, &
     1421!       Gas
     1422       & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, &
     1423       & ZCCL4, POZON_i(:,:,1), ZO2, &
     1424!       nuages :
     1425      & cldfra_i, flwc_i, fiwc_i, ZQ_SNOW, &
     1426!       rayons effectifs des gouttelettes             
     1427      & ref_liq_i, ref_ice_i, &
     1428!       aerosols
     1429     & ZAEROSOL_OLD, ZAEROSOL, &
     1430! Outputs
     1431!       Net flux :
     1432      & ZSWFT_i, ZLWFT_i, ZSWFT0_ii, ZLWFT0_ii, &
     1433!       DWN flux :
     1434      & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
     1435!       UP flux :
     1436      & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
     1437!       Surf Direct flux : ATTENTION
     1438      & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
     1439!       UV and para flux
     1440      & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
     1441!      & ZFLUX_SW_DN_TOA,
     1442      & ZEMIS_OUT, ZLWDERIVATIVE, &
     1443      & PSFSWDIF, PSFSWDIR, &
     1444      & cloud_cover_sw)
     1445  endif
     1446
    13981447
    13991448      print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
Note: See TracChangeset for help on using the changeset viewer.