source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad.v1.5.1/radiation_scheme.F90

! AI mars 2021
! ====================== Interface between ECRAD and LMDZ ====================
! radiation_scheme.F90 appelee dans radlwsw_m.F90 si iflag_rttm = 2
! revoir toutes les parties avec "AI ATTENTION" 
! Mars 2021 :
!             - Revoir toutes les parties commentees AI ATTENTION
!             1. Traitement des aerosols
!             2. Verifier les parametres times issus de LMDZ (calcul issed)
!             3. Configuration a partir de namelist
!             4. frac_std = 0.75      
! Juillet 2023 :
!             
! ============================================================================

SUBROUTINE RADIATION_SCHEME &
! Inputs
     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
     &  namelist_file, ok_3Deffect, IDAY, TIME, &
     &  PSOLAR_IRRADIANCE, &
     &  PMU0, PTEMPERATURE_SKIN, &
     &  PALBEDO_DIF, PALBEDO_DIR, &
     &  PEMIS, PEMIS_WINDOW, &
     &  PGELAM, PGEMU, &
     &  PPRESSURE_H, PTEMPERATURE_H, PQ, PQSAT, &
     &  PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, &
     &  PCCL4, PO3, PO2, &
     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_SNOW, &
     &  ZRE_LIQUID_UM, ZRE_ICE_UM, &
     &  PAEROSOL_OLD, PAEROSOL, &
! Outputs
     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, &
     &  PFLUX_SW_UP, PFLUX_LW_UP, PFLUX_SW_UP_CLEAR, PFLUX_LW_UP_CLEAR, &
     &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
     &  PEMIS_OUT, PLWDERIVATIVE, &
     &  PSWDIFFUSEBAND, PSWDIRECTBAND)


! RADIATION_SCHEME - Interface to modular radiation scheme

! (C) Copyright 2015- ECMWF.

! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.

! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

! PURPOSE
! -------
!   The modular radiation scheme is contained in a separate
!   library. This routine puts the the IFS arrays into appropriate
!   objects, computing the additional data that is required, and sends
!   it to the radiation scheme.  It returns net fluxes and surface
!   flux components needed by the rest of the model. 

!   Lower case is used for variables and types taken from the
!   radiation library

! INTERFACE
! ---------
!    RADIATION_SCHEME is called from RADLSWR. The
!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
!    should have been run first.

! AUTHOR
! ------
!   Robin Hogan, ECMWF
!   Original: 2015-09-16

! MODIFICATIONS
! -------------

! TO DO
! -----

!-----------------------------------------------------------------------

! Modules from ifs or ifsaux libraries
USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
USE RADIATION_SETUP
USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant
USE RADIATION_SETUP, ONLY : SETUP_RADIATION_SCHEME, &
                         &  config_type, driver_config_type, &
                         &  NWEIGHT_UV,  IBAND_UV,  WEIGHT_UV, &
                         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
                         &  ITYPE_TROP_BG_AER,  TROP_BG_AER_MASS_EXT, &
                         &  ITYPE_STRAT_BG_AER, STRAT_BG_AER_MASS_EXT, &
                         &  ISolverSpartacus

! Modules from radiation library
USE radiation_single_level,   ONLY : single_level_type
USE radiation_thermodynamics, ONLY : thermodynamics_type
USE radiation_gas
USE radiation_cloud,          ONLY : cloud_type
USE radiation_aerosol,        ONLY : aerosol_type
USE radiation_flux,           ONLY : flux_type
USE radiation_interface,      ONLY : radiation, set_gas_units
USE radiation_save,           ONLY : save_inputs

USE lmdz_phys_para

IMPLICIT NONE

! INPUT ARGUMENTS
! *** Array dimensions and ranges
INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
!INTEGER, INTENT(IN) :: KIDIA, KFDIA
INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
!INTEGER, INTENT(IN) :: KLON, KLEV
!INTEGER(KIND=JPIM),INTENT(IN) :: KAEROLMDZ ! Number of aerosol types
INTEGER(KIND=JPIM),INTENT(IN) :: KAEROSOL
INTEGER(KIND=JPIM),INTENT(IN) :: NSW ! Numbe of bands

! AI ATTENTION
!INTEGER, PARAMETER :: KAEROSOL = 12

! *** Single-level fields
REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
! Diffuse and direct components of surface shortwave albedo
!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,YRERAD%NSW)
!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,YRERAD%NSW)
REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,NSW)
REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,NSW)
! Longwave emissivity outside and inside the window region
REAL(KIND=JPRB),   INTENT(IN) :: PEMIS(KLON)
REAL(KIND=JPRB),   INTENT(IN) :: PEMIS_WINDOW(KLON)
! Longitude (radians), sine of latitude
REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
! Land-sea mask
!REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON) 

! *** Variables on half levels
REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)

! *** Gas mass mixing ratios on full levels
REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV) 
! AI
REAL(KIND=JPRB),   INTENT(IN) :: PQSAT(KLON,KLEV)
REAL(KIND=JPRB),   INTENT(IN) :: PCO2
REAL(KIND=JPRB),   INTENT(IN) :: PCH4 
REAL(KIND=JPRB),   INTENT(IN) :: PN2O 
REAL(KIND=JPRB),   INTENT(IN) :: PNO2 
REAL(KIND=JPRB),   INTENT(IN) :: PCFC11
REAL(KIND=JPRB),   INTENT(IN) :: PCFC12
REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22
REAL(KIND=JPRB),   INTENT(IN) :: PCCL4
REAL(KIND=JPRB),   INTENT(IN) :: PO3(KLON,KLEV) ! AI (kg/kg) ATTENTION (Pa*kg/kg) 
REAL(KIND=JPRB),   INTENT(IN) :: PO2

! *** Cloud fraction and hydrometeor mass mixing ratios
REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
!REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)

! *** Aerosol mass mixing ratios
REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)

!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON) 
!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON) 

!AI mars 2021
INTEGER(KIND=JPIM), INTENT(IN) :: IDAY
REAL(KIND=JPRB), INTENT(IN)    :: TIME


! OUTPUT ARGUMENTS

! *** Net fluxes on half-levels (W m-2)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1) 
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1) 
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1) 
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1) 

!*** DN and UP flux on half-levels (W m-2)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1)

! Direct component of surface flux into horizontal plane
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR(KLON)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON)
! As PFLUX_DIR but into a plane perpendicular to the sun
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)

! *** Ultraviolet and photosynthetically active radiation (W m-2)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)

! Diagnosed longwave surface emissivity across the whole spectrum
REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)   

! Partial derivative of total-sky longwave upward flux at each level
! with respect to upward flux at surface, used to correct heating
! rates at gridpoints/timesteps between calls to the full radiation
! scheme.  Note that this version uses the convention of level index
! increasing downwards, unlike the local variable ZLwDerivative that
! is returned from the LW radiation scheme.
REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)

! Surface diffuse and direct downwelling shortwave flux in each
! shortwave albedo band, used in RADINTG to update the surface fluxes
! accounting for high-resolution albedo information
REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSW)
REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,NSW)

! LOCAL VARIABLES
! AI ATTENTION
type(config_type),save         :: rad_config
!!$OMP THREADPRIVATE(rad_config)
type(driver_config_type),save  :: driver_config
!!$OMP THREADPRIVATE(driver_config)
!type(config_type)        :: rad_config
!type(driver_config_type)  :: driver_config
TYPE(single_level_type)   :: single_level
TYPE(thermodynamics_type) :: thermodynamics
TYPE(gas_type)            :: gas
TYPE(cloud_type)          :: cloud
TYPE(aerosol_type)        :: aerosol
TYPE(flux_type)           :: flux

! Mass mixing ratio of ozone (kg/kg)
REAL(KIND=JPRB)           :: ZO3(KLON,KLEV)

! Cloud effective radii in microns
REAL(KIND=JPRB)           :: ZRE_LIQUID_UM(KLON,KLEV)
REAL(KIND=JPRB)           :: ZRE_ICE_UM(KLON,KLEV)

! Cloud overlap decorrelation length for cloud boundaries in km
REAL(KIND=JPRB)           :: ZDECORR_LEN_KM(KLON)

! Ratio of cloud overlap decorrelation length for cloud water
! inhomogeneities to that for cloud boundaries (typically 0.5)
!REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO = 0.5_jprb

! The surface net longwave flux if the surface was a black body, used
! to compute the effective broadband surface emissivity
REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)

! Layer mass in kg m-2
REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)

! Time integers
INTEGER :: ITIM

! Loop indices
INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER

REAL(KIND=JPRB) :: ZHOOK_HANDLE

! AI ATTENTION traitement aerosols
INTEGER, PARAMETER :: NAERMACC = 1

! Name of file names specified on command line
character(len=512) :: namelist_file

logical :: loutput=.true.
logical :: lprint_input=.false.
logical :: lprint_config=.true.
logical, save :: debut_ecrad=.true.
!$OMP THREADPRIVATE(debut_ecrad)
integer, save :: itap_ecrad=1
logical :: ok_3Deffect

IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)

! A.I juillet 2023 : 
! Initialisation dans radiation_setup au 1er passage dans Ecrad
!$OMP MASTER
if (.not.ok_3Deffect) then
  if (debut_ecrad) then
   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
   debut_ecrad=.false. 
  endif 
else
   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
endif
!$OMP END MASTER
!$OMP BARRIER 
! Fin partie initialisation et configuration 

!AI juillet 2023 : verif des param de config :
if (lprint_config) then
! IF (is_master) THEN       
   print*,'Parametres de configuration de ecrad, etape ',itap_ecrad     
   print*,'Entree dans radiation_scheme'
   print*,'ok_3Deffect = ',ok_3Deffect
   print*,'Fichier namelist = ',namelist_file

   print*,'do_sw, do_lw, do_sw_direct, do_3d_effects = ', &
           rad_config%do_sw, rad_config%do_lw, rad_config%do_sw_direct, rad_config%do_3d_effects
   print*,'do_lw_side_emissivity, do_clear, do_save_radiative_properties = ', &
           rad_config%do_lw_side_emissivity, rad_config%do_clear, rad_config%do_save_radiative_properties
!   print*,'sw_entrapment_name, sw_encroachment_name = ', &
!           rad_config%sw_entrapment_name, rad_config%sw_encroachment_name
   print*,'do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug = ', &
           rad_config%do_3d_lw_multilayer_effects, rad_config%do_fu_lw_ice_optics_bug
   print*,'do_save_spectral_flux, do_save_gpoint_flux = ', &
           rad_config%do_save_spectral_flux, rad_config%do_save_gpoint_flux
   print*,'do_surface_sw_spectral_flux, do_lw_derivatives = ', &
           rad_config%do_surface_sw_spectral_flux, rad_config%do_lw_derivatives
   print*,'do_lw_aerosol_scattering, do_lw_cloud_scattering = ', &
           rad_config%do_lw_aerosol_scattering, rad_config%do_lw_cloud_scattering
   print*, 'nregions, i_gas_model = ', &
           rad_config%nregions, rad_config%i_gas_model 
!   print*, 'ice_optics_override_file_name, liq_optics_override_file_name = ', &
!           rad_config%ice_optics_override_file_name, rad_config%liq_optics_override_file_name
!   print*, 'aerosol_optics_override_file_name, cloud_pdf_override_file_name = ', &
!           rad_config%aerosol_optics_override_file_name, rad_config%cloud_pdf_override_file_name
!   print*, 'gas_optics_sw_override_file_name, gas_optics_lw_override_file_name = ', &
!           rad_config%gas_optics_sw_override_file_name, rad_config%gas_optics_lw_override_file_name
   print*, 'i_liq_model, i_ice_model, max_3d_transfer_rate = ', &
           rad_config%i_liq_model, rad_config%i_ice_model, rad_config%max_3d_transfer_rate
   print*, 'min_cloud_effective_size, overhang_factor = ', &
           rad_config%min_cloud_effective_size, rad_config%overhang_factor
   print*, 'use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw = ', &
           rad_config%use_canopy_full_spectrum_sw, rad_config%use_canopy_full_spectrum_lw
   print*, 'do_canopy_fluxes_sw, do_canopy_fluxes_lw = ', &
           rad_config%do_canopy_fluxes_sw, rad_config%do_canopy_fluxes_lw
   print*, 'do_canopy_gases_sw, do_canopy_gases_lw = ', &
           rad_config%do_canopy_gases_sw, rad_config%do_canopy_gases_lw 
   print*, 'use_general_cloud_optics, use_general_aerosol_optics = ', &
           rad_config%use_general_cloud_optics, rad_config%use_general_aerosol_optics
   print*, 'do_sw_delta_scaling_with_gases, i_overlap_scheme = ', &
           rad_config%do_sw_delta_scaling_with_gases, rad_config%i_overlap_scheme
   print*, 'i_solver_sw, i_solver_sw, use_beta_overlap, use_vectorizable_generator = ', &
           rad_config%i_solver_sw, rad_config%i_solver_lw, &
           rad_config%use_beta_overlap, rad_config%use_vectorizable_generator
   print*, 'use_expm_everywhere, iverbose, iverbosesetup = ', &
           rad_config%use_expm_everywhere, rad_config%iverbose, rad_config%iverbosesetup
   print*, 'cloud_inhom_decorr_scaling, cloud_fraction_threshold = ', &
           rad_config%cloud_inhom_decorr_scaling, rad_config%cloud_fraction_threshold
   print*, 'clear_to_thick_fraction, max_gas_od_3d, max_cloud_od = ', &
           rad_config%clear_to_thick_fraction, rad_config%max_gas_od_3d, rad_config%max_cloud_od
   print*, 'cloud_mixing_ratio_threshold, overhead_sun_factor =', &
           rad_config%cloud_mixing_ratio_threshold, rad_config%overhead_sun_factor
   print*, 'n_aerosol_types, i_aerosol_type_map, use_aerosols = ', &
           rad_config%n_aerosol_types, rad_config%i_aerosol_type_map, rad_config%use_aerosols
   print*, 'mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od = ', &
           rad_config%mono_lw_wavelength, rad_config%mono_lw_total_od,rad_config% mono_sw_total_od
   print*, 'mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo = ', &
           rad_config%mono_lw_single_scattering_albedo, rad_config%mono_sw_single_scattering_albedo
   print*, 'mono_lw_asymmetry_factor, mono_sw_asymmetry_factor = ', &
           rad_config%mono_lw_asymmetry_factor, rad_config%mono_sw_asymmetry_factor
   print*, 'i_cloud_pdf_shape = ', & 
           rad_config%i_cloud_pdf_shape
!           cloud_type_name, use_thick_cloud_spectral_averaging = ', &
!           rad_config%i_cloud_pdf_shape, rad_config%cloud_type_name, &
!           rad_config%use_thick_cloud_spectral_averaging
   print*, 'do_nearest_spectral_sw_albedo, do_nearest_spectral_lw_emiss = ', &
           rad_config%do_nearest_spectral_sw_albedo, rad_config%do_nearest_spectral_lw_emiss
   print*, 'sw_albedo_wavelength_bound, lw_emiss_wavelength_bound = ', &
           rad_config%sw_albedo_wavelength_bound, rad_config%lw_emiss_wavelength_bound
   print*, 'i_sw_albedo_index, i_lw_emiss_index = ', &
           rad_config%i_sw_albedo_index, rad_config%i_lw_emiss_index
   print*, 'do_cloud_aerosol_per_lw_g_point = ', &
           rad_config%do_cloud_aerosol_per_lw_g_point
   print*, 'do_cloud_aerosol_per_sw_g_point, do_weighted_surface_mapping = ', &
           rad_config%do_cloud_aerosol_per_sw_g_point, rad_config%do_weighted_surface_mapping
   print*, 'n_bands_sw, n_bands_lw, n_g_sw, n_g_lw = ', &
           rad_config%n_bands_sw, rad_config%n_bands_lw, rad_config%n_g_sw, rad_config%n_g_lw

           itap_ecrad=itap_ecrad+1
!   ENDIF          
endif           

! AI ATTENTION
! Allocate memory in radiation objects
! emissivite avec une seule bande
CALL single_level%allocate(KLON, NSW, 1, &
     &                     use_sw_albedo_direct=.TRUE.)

print*,'************* THERMO (allocate + input) ************************************'
! Set thermodynamic profiles: simply copy over the half-level
! pressure and temperature
!print*,'Appel allocate thermo'
CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
!print*,'Definir les champs thermo'
! AI
! pressure_hl > paprs
! temperature_hl calculee dans radlsw de la meme facon que pour RRTM
thermodynamics%pressure_hl   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
thermodynamics%temperature_hl(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)

!print*,'Compute saturation specific humidity'
! Compute saturation specific humidity, used to hydrate aerosols. The
! "2" for the last argument indicates that the routine is not being
! called from within the convection scheme.
!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, &
!     &  PPRESSURE, PTEMPERATURE, thermodynamics%h2o_sat_liq, 2)
! Alternative approximate version using temperature and pressure from
! the thermodynamics structure
!CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
!AI ATTENTION
thermodynamics%h2o_sat_liq = PQSAT

print*,'********** SINGLE LEVEL VARS **********************************'
!AI ATTENTION
! Set single-level fileds
single_level%solar_irradiance              = PSOLAR_IRRADIANCE
single_level%cos_sza(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
single_level%sw_albedo(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
single_level%lw_emissivity(KIDIA:KFDIA,1)  = PEMIS(KIDIA:KFDIA)
!single_level%lw_emissivity(KIDIA:KFDIA,2)  = PEMIS_WINDOW(KIDIA:KFDIA)

! Create the relevant seed from date and time get the starting day
! and number of minutes since start
!IDAY = NDD(NINDAT)
!cur_day
!ITIM = NINT(NSTEP * YRRIP%TSTEP / 60.0_JPRB)
ITIM = NINT(TIME / 60.0_JPRB)
!current_time
!allocate(single_level%iseed(KIDIA:KFDIA))
DO JLON = KIDIA, KFDIA
  ! This method gives a unique value for roughly every 1-km square
  ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
  ! latitude in degrees, which we multiply by 100 to give a unique
  ! value for roughly every km. PGELAM*60*100 gives a unique number
  ! for roughly every km of longitude around the equator, which we
  ! multiply by 180*100 so there is no overlap with the latitude
  ! values.  The result can be contained in a 32-byte integer (but
  ! since random numbers are generated with the help of integer
  ! overflow, it should not matter if the number did overflow).
  single_level%iseed(JLON) = ITIM + IDAY & 
       &  +  NINT(PGELAM(JLON)*108000000.0_JPRB &
       &          + ASIN(PGEMU(JLON))*6000.0_JPRB)
ENDDO

print*,'********** CLOUDS (allocate + input) *******************************************'
!print*,'Appel Allocate clouds'
CALL cloud%allocate(KLON, KLEV)
! Set cloud fields
cloud%q_liq(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
cloud%q_ice(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
cloud%fraction(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)

!!! ok AI ATTENTION a voir avec JL
! Compute effective radi and convert to metres
! AI. : on passe directement les champs de LMDZ
cloud%re_liq(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:)
cloud%re_ice(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:)

! Get the cloud overlap decorrelation length (for cloud boundaries),
! in km, according to the parameterization specified by NDECOLAT,
! and insert into the "cloud" object. Also get the ratio of
! decorrelation lengths for cloud water content inhomogeneities and
! cloud boundaries, and set it in the "rad_config" object.
! IFS :
!CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA, KFDIA, KLON, PGEMU, YRERAD%NDECOLAT, &
!     &    ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
! AI valeur dans namelist
! rad_config%cloud_inhom_decorr_scaling = ZDECORR_LEN_RATIO

!AI ATTENTION meme valeur que dans offline
! A mettre dans namelist
ZDECORR_LEN_KM = driver_config%overlap_decorr_length
DO JLON = KIDIA,KFDIA
  CALL cloud%set_overlap_param(thermodynamics, &
       &                 ZDECORR_LEN_KM(JLON), &
       &                 istartcol=JLON, iendcol=JLON)
ENDDO

! IFS :
! Cloud water content fractional standard deviation is configurable
! from namelist NAERAD but must be globally constant. Before it was
! hard coded at 1.0.
!CALL cloud%create_fractional_std(KLON, KLEV, YRERAD%RCLOUD_FRAC_STD)
! AI ATTENTION frac_std=0.75 meme valeur que dans la version offline
CALL cloud%create_fractional_std(KLON, KLEV, driver_config%frac_std)

if (rad_config%i_solver_sw == ISolverSPARTACUS &
      & .or.   rad_config%i_solver_lw == ISolverSPARTACUS) then
! AI ! Read cloud properties needed by SPARTACUS
!AI ATTENTION meme traitement dans le version offline

! By default mid and high cloud effective size is 10 km
!CALL cloud%create_inv_cloud_effective_size(KLON,KLEV,1.0_JPRB/10000.0_JPRB)

!  if (driver_config%low_inv_effective_size >= 0.0_jprb &
!     & .or. driver_config%middle_inv_effective_size >= 0.0_jprb &
!     & .or. driver_config%high_inv_effective_size >= 0.0_jprb) then
  if (driver_config%ok_effective_size) then
     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
               &  thermodynamics%pressure_hl, &
               &  driver_config%low_inv_effective_size, &
               &  driver_config%middle_inv_effective_size, &
               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb)
!  else if (driver_config%cloud_separation_scale_surface > 0.0_jprb &
!         .AND. driver_config%cloud_separation_scale_toa > 0.0_jprb) then
  else if (driver_config%ok_separation) then
      call cloud%param_cloud_effective_separation_eta(klon, klev, &
               &  thermodynamics%pressure_hl, &
               &  driver_config%cloud_separation_scale_surface, &
               &  driver_config%cloud_separation_scale_toa, &
               &  driver_config%cloud_separation_scale_power, &
               &  driver_config%cloud_inhom_separation_factor)
  endif
endif  

print*,'******** AEROSOLS (allocate + input) **************************************'
!IF (NAERMACC > 0) THEN
  CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL) ! MACC climatology
!ELSE
!  CALL aerosol%allocate(KLON, 1, KLEV, 6) ! Tegen climatology
!ENDIF
! Compute the dry mass of each layer neglecting humidity effects, in
! kg m-2, needed to scale some of the aerosol inputs
! AI commente ATTENTION
!CALL thermodynamics%get_layer_mass(ZLAYER_MASS)

! Copy over aerosol mass mixing ratio
!IF (NAERMACC > 0) THEN

  ! MACC aerosol climatology - this is already in mass mixing ratio
  ! units with the required array orientation so we can copy it over
  ! directly
  aerosol%mixing_ratio(KIDIA:KFDIA,:,:) = PAEROSOL(KIDIA:KFDIA,:,:)

  ! Add the tropospheric and stratospheric backgrounds contained in the
  ! old Tegen arrays - this is very ugly!
! AI ATTENTION
!  IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
!         &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:) &
!         &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
!  ENDIF
!  IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
!         &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:) &
!         &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
!  ENDIF

!ELSE

  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
  ! 550-nm optical depth in each layer. The optics data file
  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
  ! coefficient, but a scaling factor that the 550-nm optical depth
  ! should be multiplied by to obtain the optical depth in each
  ! spectral band.  Therefore, in order for the units to work out, we
  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
  ! need to permute the array.
!  DO JLEV = 1,KLEV
!    DO JAER = 1,6
!      aerosol%mixing_ratio(KIDIA:KFDIA,JLEV,JAER) &
!         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV) &
!         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
!    ENDDO
!  ENDDO

!ENDIF

print*,'********** GAS (allocate + input) ************************************************'
!print*,'Appel Allocate gas'
CALL gas%allocate(KLON, KLEV)

! Convert ozone Pa*kg/kg to kg/kg
! AI ATTENTION
!DO JLEV = 1,KLEV
!  DO JLON = KIDIA,KFDIA
!    ZO3(JLON,JLEV) = PO3_DP(JLON,JLEV) &
!         &         / (PPRESSURE_H(JLON,JLEV+1)-PPRESSURE_H(JLON,JLEV))
!  ENDDO
!ENDDO

!  Insert gas mixing ratios
!print*,'Insert gas mixing ratios'
CALL gas%put(IH2O,    IMassMixingRatio, PQ)
CALL gas%put(IO3,     IMassMixingRatio, PO3)
CALL gas%put_well_mixed(ICO2,    IMAssMixingRatio, PCO2)
CALL gas%put_well_mixed(ICH4,    IMassMixingRatio, PCH4)
CALL gas%put_well_mixed(IN2O,    IMassMixingRatio, PN2O)
CALL gas%put_well_mixed(ICFC11,  IMassMixingRatio, PCFC11)
CALL gas%put_well_mixed(ICFC12,  IMassMixingRatio, PCFC12)
CALL gas%put_well_mixed(IHCFC22, IMassMixingRatio, PHCFC22)
CALL gas%put_well_mixed(ICCL4,   IMassMixingRatio, PCCL4)
CALL gas%put_well_mixed(IO2,     IMassMixingRatio, PO2)
! Ensure the units of the gas mixing ratios are what is required by
! the gas absorption model
call set_gas_units(rad_config, gas)

print*,'************** FLUX (allocate) ***********************'
CALL flux%allocate(rad_config, 1, KLON, KLEV)

! Call radiation scheme
print*,'******** Appel radiation scheme **************************'
CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
     &  single_level, thermodynamics, gas, cloud, aerosol, flux)

! Compute required output fluxes
! DN and UP flux
PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)

! First the net fluxes
PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
     &  = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
PFLUX_LW_CLEAR(KIDIA:KFDIA,:) &
     &  = flux%lw_dn_clear(KIDIA:KFDIA,:) - flux%lw_up_clear(KIDIA:KFDIA,:)

! Now the surface fluxes
!PFLUX_SW_DN_SURF(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,KLEV+1)
!PFLUX_LW_DN_SURF(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,KLEV+1)
!PFLUX_SW_UP_SURF(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,KLEV+1)
!PFLUX_LW_UP_SURF(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,KLEV+1)
!PFLUX_SW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_dn_clear(KIDIA:KFDIA,KLEV+1)
!PFLUX_LW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_dn_clear(KIDIA:KFDIA,KLEV+1)
!PFLUX_SW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_up_clear(KIDIA:KFDIA,KLEV+1)
!PFLUX_LW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_up_clear(KIDIA:KFDIA,KLEV+1)
PFLUX_DIR(KIDIA:KFDIA) = flux%sw_dn_direct(KIDIA:KFDIA,KLEV+1)
PFLUX_DIR_CLEAR(KIDIA:KFDIA) = flux%sw_dn_direct_clear(KIDIA:KFDIA,KLEV+1)
PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA) / PMU0(KIDIA:KFDIA)
END WHERE

! Top-of-atmosphere downwelling flux
!PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
!PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
!PFLUX_LW_DN_TOA(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,1)
!PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)

! Compute UV fluxes as weighted sum of appropriate shortwave bands
!AI ATTENTION
if (0.eq.1) then
PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
DO JBAND = 1,NWEIGHT_UV
  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
ENDDO

! Compute photosynthetically active radiation similarly
PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
DO JBAND = 1,NWEIGHT_PAR
  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
       &  + WEIGHT_PAR(JBAND) &
       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
ENDDO
endif
! Compute effective broadband emissivity
ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5) 
  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
END WHERE

! Copy longwave derivatives
! AI ATTENTION
!IF (YRERAD%LAPPROXLWUPDATE) THEN
IF (rad_config%do_lw_derivatives) THEN
  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
END IF

! Store the shortwave downwelling fluxes in each albedo band
!AI ATTENTION
!IF (YRERAD%LAPPROXSWUPDATE) THEN
if (0.eq.1) then
IF (rad_config%do_surface_sw_spectral_flux) THEN
  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
  DO JBAND = 1,rad_config%n_bands_sw
    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
    DO JLON = KIDIA,KFDIA
      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
           &  + flux%sw_dn_surf_band(JBAND,JLON) &
           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
    ENDDO
  ENDDO
ENDIF
endif
CALL single_level%deallocate
CALL thermodynamics%deallocate
CALL gas%deallocate
CALL cloud%deallocate
CALL aerosol%deallocate
CALL flux%deallocate

IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)

END SUBROUTINE RADIATION_SCHEME