! AI mars 2021
! ====================== Interface between ECRAD and LMDZ ====================
! Depart de la version IFS R.H 
! 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      
!
! Avril 2026 : 
! Module et routines renommees :
! lmdz_ecrad_interface_m
!            lmdz_ecrad_interface / lmdz_ecrad_interface_2call
!             
! ============================================================================
MODULE lmdz_ecrad_interface_m

IMPLICIT NONE        

contains

SUBROUTINE lmdz_ecrad_interface &
! Inputs
     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
     &  namelist_file, ok_2xcall_ecrad, &
     &  debut, ok_volcan, flag_aerosol_strat, &
     &  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, &
     &  ecrad_cloud_cover_sw)

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

! Modules 
USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
USE RADIATION_SETUP
USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant

! 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 mod_phys_lmdz_para
use geometry_mod, only: longitude_deg, latitude_deg
USE read_rsun_ecrad_m,        ONLY : SPSOLARSCAL

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(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
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

! Name of file names specified on command line
character(len=512), INTENT(IN) :: namelist_file
logical, INTENT(IN)            :: ok_2xcall_ecrad, debut, ok_volcan
INTEGER(KIND=JPIM), INTENT(IN) :: flag_aerosol_strat


! 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,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON,KLEV+1)
! 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)

!AI Nov 2023
REAL(KIND=JPRB),  INTENT(OUT) :: ecrad_cloud_cover_sw(KLON)

! 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(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_M
REAL(KIND=JPRB)           :: ZDECORR_LEN_M_1D(KLON)
REAL(KIND=JPRB)           :: ZDECORR_LEN_M_2D(KLON,KLEV)

! 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

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

REAL(KIND=JPRB) ::  inv_cloud_effective_size(KLON,KLEV)
REAL(KIND=JPRB) ::  inv_inhom_effective_size(KLON,KLEV)

integer :: irang


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

  print*,'Entree radiation_scheme, ok_2xcall_ecrad, namelist_file = ', &
          ok_2xcall_ecrad, namelist_file
! A.I juillet 2023 : 
! Initialisation dans radiation_setup au 1er passage dans Ecrad
!$OMP MASTER
  if (debut_ecrad) then
   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
  endif 
!$OMP END MASTER
!$OMP BARRIER 
! Fin partie initialisation et configuration 

!AI print fichiers namelist utilise
!if (is_omp_root) then
! itap_ecrad=itap_ecrad+1
! print*,'Dans radiation_scheme itap_ecrad, mpi_rank, omp_rank, namelist_file : ', &
!         itap_ecrad, mpi_rank, omp_rank, namelist_file
!else
! print*,'mpi_rank omp_rank, namelist_file :', mpi_rank, omp_rank, namelist_file
!endif

! AI 11 23 Allocates depplaces au debut
print*,'*********** ALLOCATES *******************************'
! AI ATTENTION
! Allocate memory in radiation objects
! emissivite avec une seule bande
CALL single_level%allocate(KLON, NSW, 1, &
     &                     use_sw_albedo_direct=.TRUE.)
CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
CALL cloud%allocate(KLON, KLEV)
CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL)
CALL gas%allocate(KLON, KLEV)
CALL flux%allocate(rad_config, 1, KLON, KLEV)

print*,'************* THERMO (input) ************************************'
! 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)
if (rad_config%use_spectral_solar_scaling) then
  single_level%spectral_solar_scaling(:)     = SPSOLARSCAL(:)
endif

! 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
!AI Nov 23
! Simple initialization of the seeds for the Monte Carlo scheme
call single_level%init_seed_simple(kidia, kfdia)

print*,'********** CLOUDS (allocate + input) *******************************************'
!print*,'Appel Allocate clouds'
! 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)
CALL CALCUL_CLOUD_OVERLAP_DECORR_LEN &
     & (KIDIA, KFDIA, KLON, KLEV, &
     &  driver_config, &
     &  thermodynamics%pressure_hl, &
     &  ZDECORR_LEN_M_2D)

 if (driver_config%kdecolat.eq.0) then
  ZDECORR_LEN_M = ZDECORR_LEN_M_2D(1,1)       
  DO JLON = KIDIA,KFDIA
   CALL cloud%set_overlap_param(thermodynamics, &
       &                 ZDECORR_LEN_M, &
       &                 istartcol=JLON, iendcol=JLON)
  ENDDO
 else if (driver_config%kdecolat.eq.1.or.driver_config%kdecolat.eq.2) then
  ZDECORR_LEN_M_1D = ZDECORR_LEN_M_2D(:,1)
  DO JLON = KIDIA,KFDIA
   CALL cloud%set_overlap_param(thermodynamics, &
       &                 ZDECORR_LEN_M_1D, &
       &                 istartcol=JLON, iendcol=JLON)
  ENDDO
 else if (driver_config%kdecolat.eq.3) then
  DO JLON = KIDIA,KFDIA
   CALL cloud%set_overlap_param_var2D(thermodynamics, &
       &                 ZDECORR_LEN_M_2D, KLEV, &
       &                 istartcol=JLON, iendcol=JLON)
  ENDDO
 endif



! 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 (ok_2xcall_ecrad) 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, &
!               &  KIDIA, KFDIA)
!   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, &
!               &  KIDIA, KFDIA)
!   endif     
! else  
 if (rad_config%i_solver_sw == ISolverSPARTACUS &
      & .or.   rad_config%i_solver_lw == ISolverSPARTACUS) then
   ! AI ! Read cloud properties needed by SPARTACUS
   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, &
               &  KIDIA, KFDIA)
   else if (driver_config%ok_separation_eta) 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, &
               &  KIDIA, KFDIA)
   else if (driver_config%ok_separation_tanh) then
     call cloud%param_cloud_effective_separation_tanh(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, &
               &  KIDIA, KFDIA)
   endif
 endif  
!endif

print*,'******** AEROSOLS (input) **************************************'
!IF (NAERMACC > 0) THEN
!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 (input) ************************************************'
!print*,'Appel Allocate gas'
! 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)

! Call radiation scheme
!print*,'*** Appel radiation *** namelist **** omp_rank ****', &
!         omp_rank, namelist_file
!   if (rad_config%i_solver_sw == ISolverSPARTACUS) then 
!    if (driver_config%ok_separation) then
!      print*,'Avant radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
!              mpi_rank, omp_rank, &
!              shape(cloud%inv_cloud_effective_size), &
!              size(cloud%inv_cloud_effective_size)      
!      do jlon=KIDIA, KFDIA      
!        do jlev=1,klev
!         print*,' Avant radiation mpi_rank, omp_rank, jlon, jlev, &
!            &     cloud%inv_cloud_effective_size =', mpi_rank, &
!            &     omp_rank, jlon, jlev, &
!            &     cloud%inv_cloud_effective_size(jlon,jlev)
!        enddo
!      enddo
!       cloud%inv_cloud_effective_size=inv_cloud_effective_size
!       cloud%inv_inhom_effective_size=inv_inhom_effective_size
!    endif
!   endif
if (debut_ecrad .and. driver_config%do_save_inputs) &
     call save_inputs('inputs.nc', rad_config, single_level, thermodynamics, &
     gas, cloud, aerosol, lat = latitude_deg, lon = longitude_deg)

CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
     &  single_level, thermodynamics, gas, cloud, aerosol, flux)

 if (rad_config%use_aerosols) then
    if (rad_config%i_gas_model_sw == IGasModelIFSRRTMG .or. &
     &  rad_config%i_gas_model_lw == IGasModelIFSRRTMG) then     
       CALL aeropt_5wv_ecrad(kidia, kfdia, 1, klev,     &
     &                       rad_config,thermodynamics, aerosol)
    endif                 

    if (flag_aerosol_strat.eq.2) then
       CALL readaerosolstrato_ecrad(rad_config, debut, ok_volcan)
    endif
 endif               

print*,'*********** Sortie flux ****************'
! Cloud cover
ecrad_cloud_cover_sw = flux%cloud_cover_sw
! 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)
! Direct component of flux into horizontal plane
PFLUX_DIR(KIDIA:KFDIA,:) = flux%sw_dn_direct(KIDIA:KFDIA,:)
PFLUX_DIR_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_direct_clear(KIDIA:KFDIA,:)
PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
! Direct Surface component of flux into a plane perpendicular to the sun        
  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA,KLEV+1) / 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)
!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

print*,'********** DEALLOCATIONS ************************'
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)
debut_ecrad=.false.

END SUBROUTINE lmdz_ecrad_interface

SUBROUTINE lmdz_ecrad_interface_2call &
! Inputs
     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
     &  namelist_file, ok_2xcall_ecrad, &
     &  debut, ok_volcan, flag_aerosol_strat, &
     &  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, &
     &  ecrad_cloud_cover_sw)

! 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 mod_phys_lmdz_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

! Name of file names specified on command line
character(len=512), INTENT(IN) :: namelist_file
logical, INTENT(IN)            :: ok_2xcall_ecrad, debut, ok_volcan
INTEGER(KIND=JPIM), INTENT(IN) :: flag_aerosol_strat


! 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,KLEV+1)
REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON,KLEV+1)
! 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)

!AI Nov 2023
REAL(KIND=JPRB),  INTENT(OUT) :: ecrad_cloud_cover_sw(KLON)

! 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_M
REAL(KIND=JPRB)           :: ZDECORR_LEN_M_1D(KLON)
REAL(KIND=JPRB)           :: ZDECORR_LEN_M_2D(KLON,KLEV)

! 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

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

REAL(KIND=JPRB) ::  inv_cloud_effective_size(KLON,KLEV)
REAL(KIND=JPRB) ::  inv_inhom_effective_size(KLON,KLEV)

integer :: irang


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

   print*,'Entree radiation_scheme_s2, ok_2xcall_ecrad, namelist_file = ', &
          ok_2xcall_ecrad, namelist_file
! A.I juillet 2023 : 
! Initialisation dans radiation_setup au 1er passage dans Ecrad
!$OMP MASTER
!if (.not.ok_2xcall_ecrad) 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 print fichiers namelist utilise
!if (is_omp_root) then
! itap_ecrad=itap_ecrad+1
! print*,'Dans radiation_scheme itap_ecrad, mpi_rank, omp_rank, namelist_file : ', &
!         itap_ecrad, mpi_rank, omp_rank, namelist_file
!else
! print*,'mpi_rank omp_rank, namelist_file :', mpi_rank, omp_rank, namelist_file
!endif

! AI 11 23 Allocates depplaces au debut
print*,'*********** ALLOCATES *******************************'
! AI ATTENTION
! Allocate memory in radiation objects
! emissivite avec une seule bande
CALL single_level%allocate(KLON, NSW, 1, &
     &                     use_sw_albedo_direct=.TRUE.)
CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
CALL cloud%allocate(KLON, KLEV)
CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL)
CALL gas%allocate(KLON, KLEV)
CALL flux%allocate(rad_config, 1, KLON, KLEV)

print*,'************* THERMO (input) ************************************'
! Set thermodynamic profiles: simply copy over the half-level
! pressure and temperature
! 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
!AI Nov 23
! Simple initialization of the seeds for the Monte Carlo scheme
call single_level%init_seed_simple(kidia, kfdia)

print*,'********** CLOUDS (allocate + input) *******************************************'
!print*,'Appel Allocate clouds'
! 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
   CALL CALCUL_CLOUD_OVERLAP_DECORR_LEN &
     & (KIDIA, KFDIA, KLON, KLEV, &
     &  driver_config, &
     &  thermodynamics%pressure_hl , &
     &  ZDECORR_LEN_M_2D)

 if (driver_config%kdecolat.eq.0) then
  ZDECORR_LEN_M = ZDECORR_LEN_M_2D(1,1)
  DO JLON = KIDIA,KFDIA
   CALL cloud%set_overlap_param(thermodynamics, &
       &                 ZDECORR_LEN_M, &
       &                 istartcol=JLON, iendcol=JLON)
  ENDDO
 else if (driver_config%kdecolat.eq.1.or.driver_config%kdecolat.eq.2) then
  ZDECORR_LEN_M_1D = ZDECORR_LEN_M_2D(:,1)
  DO JLON = KIDIA,KFDIA
   CALL cloud%set_overlap_param(thermodynamics, &
       &                 ZDECORR_LEN_M_1D, &
       &                 istartcol=JLON, iendcol=JLON)
  ENDDO
 else if (driver_config%kdecolat.eq.3) then
  DO JLON = KIDIA,KFDIA
   CALL cloud%set_overlap_param_var2D(thermodynamics, &
       &                 ZDECORR_LEN_M_2D, KLEV, &
       &                 istartcol=JLON, iendcol=JLON)
  ENDDO
 endif
   
! 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 (ok_2xcall_ecrad) 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, &
!               &  KIDIA, KFDIA)
!   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, &
!               &  KIDIA, KFDIA)
!   endif     
! else  
 if (rad_config%i_solver_sw == ISolverSPARTACUS &
      & .or.   rad_config%i_solver_lw == ISolverSPARTACUS) then
   ! AI ! Read cloud properties needed by SPARTACUS
   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, &
               &  KIDIA, KFDIA)
   else if (driver_config%ok_separation_eta) 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, &
               &  KIDIA, KFDIA)
   else if (driver_config%ok_separation_tanh) then
     call cloud%param_cloud_effective_separation_tanh(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, &
               &  KIDIA, KFDIA)
   endif
 endif  
!endif

print*,'******** AEROSOLS (input) **************************************'
!IF (NAERMACC > 0) THEN
!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 (input) ************************************************'
!print*,'Appel Allocate gas'
! 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)

! Call radiation scheme
!print*,'*** Appel radiation *** namelist **** omp_rank ****', &
!         omp_rank, namelist_file
!   if (rad_config%i_solver_sw == ISolverSPARTACUS) then 
!    if (driver_config%ok_separation) then
!      print*,'Avant radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
!              mpi_rank, omp_rank, &
!              shape(cloud%inv_cloud_effective_size), &
!              size(cloud%inv_cloud_effective_size)      
!      do jlon=KIDIA, KFDIA      
!        do jlev=1,klev
!         print*,' Avant radiation mpi_rank, omp_rank, jlon, jlev, &
!            &     cloud%inv_cloud_effective_size =', mpi_rank, &
!            &     omp_rank, jlon, jlev, &
!            &     cloud%inv_cloud_effective_size(jlon,jlev)
!        enddo
!      enddo
!       cloud%inv_cloud_effective_size=inv_cloud_effective_size
!       cloud%inv_inhom_effective_size=inv_inhom_effective_size
!    endif
!   endif
CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
     &  single_level, thermodynamics, gas, cloud, aerosol, flux)

 if (rad_config%use_aerosols) then
    if (rad_config%i_gas_model_sw == IGasModelIFSRRTMG .or. &
     &    rad_config%i_gas_model_lw == IGasModelIFSRRTMG) then     
       CALL aeropt_5wv_ecrad(kidia, kfdia, 1, klev, &
                             rad_config,thermodynamics,aerosol)
    endif                 
    if (flag_aerosol_strat.eq.2) then
       CALL readaerosolstrato_ecrad(rad_config, debut, ok_volcan)
    endif
 endif               

print*,'*********** Sortie flux ****************'
! Cloud cover
ecrad_cloud_cover_sw = flux%cloud_cover_sw
! 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,:)
PFLUX_DIR_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_direct_clear(KIDIA:KFDIA,:)
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,KLEV+1) / 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)
!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

print*,'********** DEALLOCATIONS ************************'
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 lmdz_ecrad_interface_2call
end module lmdz_ecrad_interface_m
