Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90	(revision 4977)
+++ 	(revision )
@@ -1,100 +1,0 @@
-!aeropt_5wv_ecrad.F90 2022-09-20 A. Idelkadi et O. Boucher 
-!
-     SUBROUTINE AEROPT_5WV_ECRAD(istartcol,iendcol,istartlev,iendlev, &
-                                 config,thermodynamics,aerosol)
-
-
-      USE DIMPHY
-      USE aero_mod
-      USE phys_local_var_mod, ONLY: od443aer,od550aer,od865aer
-!                                      dryod550aer, 
-!                                      ec550aer,od550lt1aer,abs550aer
-!        USE phys_output_var_mod, ONLY: dryaod_diag
-!        USE YOMCST, ONLY: RD,RG
-      USE phys_local_var_mod, ONLY: rhcl
-
-      use parkind1,                      only : jprb
-      use radiation_config,              only : config_type
-      use radiation_thermodynamics,      only : thermodynamics_type
-      use radiation_aerosol,             only : aerosol_type
-      use radiation_constants,           only : AccelDueToGravity
-      use radiation_aerosol_optics_data, only : aerosol_optics_type, &
-         &  IAerosolClassUndefined,   IAerosolClassIgnored, &
-         &  IAerosolClassHydrophobic, IAerosolClassHydrophilic
-
-      IMPLICIT NONE
-
-      !---ATTENTION n_mono est a redéfinir proprement
-
-      ! Range of levels over which aerosols are present
-      integer, intent(in) :: istartlev, iendlev, istartcol,iendcol
-      type(config_type), intent(in), target :: config
-      type(thermodynamics_type),intent(in)  :: thermodynamics
-      type(aerosol_type),       intent(in)  :: aerosol
-      type(aerosol_optics_type), pointer :: ao
-
-      ! Loop indices for column, level, g point, band and aerosol type, and relative humidity
-      integer :: jcol, jlev, jtype, irh
-      ! indice wavelength in aerosol_optics_lmdz.nc
-      INTEGER, PARAMETER :: la443 = 1
-      INTEGER, PARAMETER :: la550 = 2
-!      INTEGER, PARAMETER :: la670 = 3
-!      INTEGER, PARAMETER :: la765 = 4
-      INTEGER, PARAMETER :: la865 = 5
-      real(jprb) :: factor
-      ! n_bands_sw
-      real(jprb), dimension(config%aerosol_optics%n_mono_wl) :: od_aerosol_mono, local_od_mono
-
-      !--initialization
-      od443aer = 0.0_jprb
-      od550aer = 0.0_jprb
-      od865aer = 0.0_jprb
-
-      ao => config%aerosol_optics
-
-      ! Loop over level
-     do jlev = istartlev,iendlev
-
-        ! Loop over column
-        do jcol = istartcol,iendcol
-
-          ! Compute relative humidity with respect to liquid
-          ! saturation and the index to the relative-humidity index of
-          ! hydrophilic-aerosol data
-          irh = ao%calc_rh_index(rhcl(jcol,jlev))
-
-          factor = ( thermodynamics%pressure_hl(jcol,jlev+1) &
-               &    -thermodynamics%pressure_hl(jcol,jlev  )  ) &
-               &   / AccelDueToGravity
-
-          ! Reset temporary arrays
-          od_aerosol_mono = 0.0_jprb
-
-          do jtype = 1,config%n_aerosol_types
-            ! Add the optical depth for this aerosol type to the total for all aerosols.  
-            ! Note that the following expressions are array-wise
-            if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
-              local_od_mono = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
-                   &  * ao%mass_ext_mono_phobic(:,ao%itype(jtype))
-              od_aerosol_mono = od_aerosol_mono + local_od_mono
-            else if (ao%iclass(jtype) == IAerosolClassHydrophilic) then
-              ! Hydrophilic aerosols require the look-up tables to
-              ! be indexed with irh
-              local_od_mono = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
-                   &  * ao%mass_ext_mono_philic(:,irh,ao%itype(jtype))
-              od_aerosol_mono = od_aerosol_mono + local_od_mono
-            end if
-
-          end do ! Loop over aerosol type
-
-          !--ATTENTION A BIEN FAIRE CORRESPONDRE LES INDICES
-          od443aer(jcol) = od443aer(jcol) + od_aerosol_mono(la443)
-          od550aer(jcol) = od550aer(jcol) + od_aerosol_mono(la550)
-          od865aer(jcol) = od865aer(jcol) + od_aerosol_mono(la865)
-
-        end do ! Loop over column
-
-      end do ! Loop over level
-
-     END SUBROUTINE AEROPT_5WV_ECRAD
-
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.f90	(revision 5268)
@@ -0,0 +1,100 @@
+!aeropt_5wv_ecrad.F90 2022-09-20 A. Idelkadi et O. Boucher 
+!
+     SUBROUTINE AEROPT_5WV_ECRAD(istartcol,iendcol,istartlev,iendlev, &
+                                 config,thermodynamics,aerosol)
+
+
+      USE DIMPHY
+      USE aero_mod
+      USE phys_local_var_mod, ONLY: od443aer,od550aer,od865aer
+!                                      dryod550aer, 
+!                                      ec550aer,od550lt1aer,abs550aer
+!        USE phys_output_var_mod, ONLY: dryaod_diag
+!        USE YOMCST, ONLY: RD,RG
+      USE phys_local_var_mod, ONLY: rhcl
+
+      use parkind1,                      only : jprb
+      use radiation_config,              only : config_type
+      use radiation_thermodynamics,      only : thermodynamics_type
+      use radiation_aerosol,             only : aerosol_type
+      use radiation_constants,           only : AccelDueToGravity
+      use radiation_aerosol_optics_data, only : aerosol_optics_type, &
+         &  IAerosolClassUndefined,   IAerosolClassIgnored, &
+         &  IAerosolClassHydrophobic, IAerosolClassHydrophilic
+
+      IMPLICIT NONE
+
+      !---ATTENTION n_mono est a redéfinir proprement
+
+      ! Range of levels over which aerosols are present
+      integer, intent(in) :: istartlev, iendlev, istartcol,iendcol
+      type(config_type), intent(in), target :: config
+      type(thermodynamics_type),intent(in)  :: thermodynamics
+      type(aerosol_type),       intent(in)  :: aerosol
+      type(aerosol_optics_type), pointer :: ao
+
+      ! Loop indices for column, level, g point, band and aerosol type, and relative humidity
+      integer :: jcol, jlev, jtype, irh
+      ! indice wavelength in aerosol_optics_lmdz.nc
+      INTEGER, PARAMETER :: la443 = 1
+      INTEGER, PARAMETER :: la550 = 2
+!      INTEGER, PARAMETER :: la670 = 3
+!      INTEGER, PARAMETER :: la765 = 4
+      INTEGER, PARAMETER :: la865 = 5
+      real(jprb) :: factor
+      ! n_bands_sw
+      real(jprb), dimension(config%aerosol_optics%n_mono_wl) :: od_aerosol_mono, local_od_mono
+
+      !--initialization
+      od443aer = 0.0_jprb
+      od550aer = 0.0_jprb
+      od865aer = 0.0_jprb
+
+      ao => config%aerosol_optics
+
+      ! Loop over level
+     do jlev = istartlev,iendlev
+
+        ! Loop over column
+        do jcol = istartcol,iendcol
+
+          ! Compute relative humidity with respect to liquid
+          ! saturation and the index to the relative-humidity index of
+          ! hydrophilic-aerosol data
+          irh = ao%calc_rh_index(rhcl(jcol,jlev))
+
+          factor = ( thermodynamics%pressure_hl(jcol,jlev+1) &
+               &    -thermodynamics%pressure_hl(jcol,jlev  )  ) &
+               &   / AccelDueToGravity
+
+          ! Reset temporary arrays
+          od_aerosol_mono = 0.0_jprb
+
+          do jtype = 1,config%n_aerosol_types
+            ! Add the optical depth for this aerosol type to the total for all aerosols.  
+            ! Note that the following expressions are array-wise
+            if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
+              local_od_mono = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
+                   &  * ao%mass_ext_mono_phobic(:,ao%itype(jtype))
+              od_aerosol_mono = od_aerosol_mono + local_od_mono
+            else if (ao%iclass(jtype) == IAerosolClassHydrophilic) then
+              ! Hydrophilic aerosols require the look-up tables to
+              ! be indexed with irh
+              local_od_mono = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
+                   &  * ao%mass_ext_mono_philic(:,irh,ao%itype(jtype))
+              od_aerosol_mono = od_aerosol_mono + local_od_mono
+            end if
+
+          end do ! Loop over aerosol type
+
+          !--ATTENTION A BIEN FAIRE CORRESPONDRE LES INDICES
+          od443aer(jcol) = od443aer(jcol) + od_aerosol_mono(la443)
+          od550aer(jcol) = od550aer(jcol) + od_aerosol_mono(la550)
+          od865aer(jcol) = od865aer(jcol) + od_aerosol_mono(la865)
+
+        end do ! Loop over column
+
+      end do ! Loop over level
+
+     END SUBROUTINE AEROPT_5WV_ECRAD
+
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.F90	(revision 4977)
+++ 	(revision )
@@ -1,155 +1,0 @@
-SUBROUTINE CALCUL_CLOUD_OVERLAP_DECORR_LEN &
-     & (KIDIA, KFDIA, KLON, KLEV, &
-     &  driver_config, &
-     &  pressure_hl, &
-     &  PDECORR_LEN_EDGES_M, PDECORR_LEN_WATER_M, PDECORR_LEN_RATIO)
-
-! A Idelkadi ! -------------
-! 04 2024
-! Depart de la routine de l' IFS
-! But : calcul de la Ld (longueur de decorrelation)
-! Dans le cas overlap_scheme_name=="Exp-Ran"
-! 3 cas en fonction de kdecolat (namelist_ecrad) :
-! 0 => Ld=overlap_decorr_length (namelist_ecrad)
-! 1 => Ld=f(latitude) Shonk et al. (2010) Eq. 13 formula
-! 2 => Ld=f(latitude) Shonk et al. (2010) lissee a l equateur
-! 3 => Ld=f(pres) : 
-!      Ld=low_decorrelation_length si pres >= 680hPa 
-!      Ld=mid_decorrelation_length si pres entre [440,680]
-!      Ld=high_decorrelation_length si pres < 440hPa
-!
-! TO DO 
-! Que faire de PDECORR_LEN_WATER_M ?
-! -------------------------------------------------------------------
-
-USE PARKIND1 , ONLY : JPIM, JPRB
-USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
-USE YOMCST   , ONLY : RPI
-USE geometry_mod, ONLY: latitude_deg
-USE lmdz_cloud_optics_prop_ini , ONLY : prmhc, prlmc
-USE setup_config_from_lmdz,   ONLY : driver_config_type
-USE write_field_phy
-! -------------------------------------------------------------------
-
-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
-! AI 04 2024
-INTEGER(KIND=JPIM),INTENT(IN) :: KLEV
-
-! *** Single-level variables 
-!AI 04 2024
-REAL(KIND=JPRB), intent(in)   :: pressure_hl(KLON,KLEV+1)
-type(driver_config_type),intent(in)   :: driver_config
-
-! OUTPUT ARGUMENTS
-
-! *** Decorrelation lengths for cloud edges and cloud water content, in m
-! A.I 04 2024
-REAL(KIND=JPRB), INTENT(OUT)           :: PDECORR_LEN_EDGES_M(KLON,KLEV)
-REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_WATER_M(KLON,KLEV)
-  
-! Ratio of water-content to cloud-edge decorrelation lengths
-REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_RATIO
-
-! LOCAL VARIABLES
-
-REAL(KIND=JPRB) :: ZRADIANS_TO_DEGREES, ZABS_LAT_DEG, ZCOS_LAT
-
-! AI 04 2024
-REAL(KIND=JPRB) :: PGEMU(KLON) ! Sine of latitude
-REAL(KIND=JPRB) :: RDECORR_CF
-REAL(KIND=JPRB) :: RDECORR_CW = 1.0_JPRB
-
-INTEGER(KIND=JPIM) :: JL, JK
-
-REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
-
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('CLOUD_OVERLAP_DECORR_LEN',0,ZHOOK_HANDLE)
-  
-! -------------------------------------------------------------------
-print*,'driver_config%kdecolat = ', driver_config%kdecolat
-
-IF (driver_config%kdecolat == 0) THEN
-
-  ! Decorrelation lengths are constant values
-  ! overlap_decorr_length lu dans namelist en m
-  RDECORR_CF = driver_config%overlap_decorr_length / 1000.0_JPRB
-  PDECORR_LEN_EDGES_M(KIDIA:KFDIA,:) = driver_config%overlap_decorr_length
-!  IF (PRESENT(PDECORR_LEN_WATER_M)) THEN
-!    PDECORR_LEN_WATER_M = RDECORR_CW * 1000.0_JPRB
-!  ENDIF
-!  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
-!    PDECORR_LEN_RATIO = RDECORR_CW / RDECORR_CF
-!  ENDIF
-
-ELSE
-
-  ZRADIANS_TO_DEGREES = 180.0_JPRB / RPI
-  PGEMU = SIN(latitude_deg)
-
-  IF (driver_config%kdecolat == 1) THEN
-    ! Shonk et al. (2010) Eq. 13 formula
-    DO JL = KIDIA,KFDIA
-      ZABS_LAT_DEG = ABS(ASIN(PGEMU(JL)) * ZRADIANS_TO_DEGREES)
-      PDECORR_LEN_EDGES_M(JL,:) = (2.899_JPRB - 0.02759_JPRB * ZABS_LAT_DEG) * 1000.0_JPRB
-    ENDDO
-  ELSE ! KDECOLAT == 2
-    DO JL = KIDIA,KFDIA
-      ! Shonk et al. (2010) but smoothed over the equator
-      ZCOS_LAT = COS(ASIN(PGEMU(JL)))
-      PDECORR_LEN_EDGES_M(JL,:) = (0.75_JPRB + 2.149_JPRB * ZCOS_LAT*ZCOS_LAT) * 1000.0_JPRB
-    ENDDO
-  ENDIF
-
-
-  ! Both KDECOLAT = 1 and 2 assume that the decorrelation length for
-  ! cloud water content is half that for cloud edges
-!  IF (PRESENT(PDECORR_LEN_WATER_M)) THEN
-!    PDECORR_LEN_WATER_M(KIDIA:KFDIA,:) = PDECORR_LEN_EDGES_M(KIDIA:KFDIA,:) * 0.5_JPRB 
-!  ENDIF
-
-!  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
-!    PDECORR_LEN_RATIO = 0.5_JPRB
-!  ENDIF
-
-ENDIF
-
-!AI 04 2024
-IF (driver_config%kdecolat == 3) THEN
-  DO jk = 1, klev
-      DO JL = 1, klon
-        IF (pressure_hl(JL,jk)<prmhc) THEN
-          PDECORR_LEN_EDGES_M(JL,jk) = driver_config%high_decorrelation_length
-        ELSE IF (pressure_hl(JL,jk)>=prmhc .AND. pressure_hl(JL,jk)<prlmc) THEN
-          PDECORR_LEN_EDGES_M(JL,jk) = driver_config%mid_decorrelation_length
-        ELSE IF (pressure_hl(JL,jk)>=prlmc) THEN
-          PDECORR_LEN_EDGES_M(JL,jk) = driver_config%low_decorrelation_length
-        ENDIF
-      ENDDO
-  ENDDO
-  
-  ! Both KDECOLAT = 1 and 2 assume that the decorrelation length for
-  ! cloud water content is half that for cloud edges
-!  IF (PRESENT(PDECORR_LEN_WATER_M)) THEN
-!    PDECORR_LEN_WATER_M(KIDIA:KFDIA,:) = PDECORR_LEN_EDGES_M(KIDIA:KFDIA,:) * 0.5_JPRB 
-!  ENDIF
-!  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
-!    PDECORR_LEN_RATIO = 0.5_JPRB
-!  ENDIF
-ENDIF
-!CALL writefield_phy('latitude',latitude_deg,1)
-!CALL writefield_phy('pressure_hl',pressure_hl,klev+1)
-!CALL writefield_phy('Ldecorel',PDECORR_LEN_EDGES_M,klev)
-! -------------------------------------------------------------------
-
-IF (LHOOK) CALL DR_HOOK('CLOUD_OVERLAP_DECORR_LEN',1,ZHOOK_HANDLE)
-
-END SUBROUTINE CALCUL_CLOUD_OVERLAP_DECORR_LEN
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.f90	(revision 5268)
@@ -0,0 +1,155 @@
+SUBROUTINE CALCUL_CLOUD_OVERLAP_DECORR_LEN &
+     & (KIDIA, KFDIA, KLON, KLEV, &
+     &  driver_config, &
+     &  pressure_hl, &
+     &  PDECORR_LEN_EDGES_M, PDECORR_LEN_WATER_M, PDECORR_LEN_RATIO)
+
+! A Idelkadi ! -------------
+! 04 2024
+! Depart de la routine de l' IFS
+! But : calcul de la Ld (longueur de decorrelation)
+! Dans le cas overlap_scheme_name=="Exp-Ran"
+! 3 cas en fonction de kdecolat (namelist_ecrad) :
+! 0 => Ld=overlap_decorr_length (namelist_ecrad)
+! 1 => Ld=f(latitude) Shonk et al. (2010) Eq. 13 formula
+! 2 => Ld=f(latitude) Shonk et al. (2010) lissee a l equateur
+! 3 => Ld=f(pres) : 
+!      Ld=low_decorrelation_length si pres >= 680hPa 
+!      Ld=mid_decorrelation_length si pres entre [440,680]
+!      Ld=high_decorrelation_length si pres < 440hPa
+!
+! TO DO 
+! Que faire de PDECORR_LEN_WATER_M ?
+! -------------------------------------------------------------------
+
+USE PARKIND1 , ONLY : JPIM, JPRB
+USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
+USE YOMCST   , ONLY : RPI
+USE geometry_mod, ONLY: latitude_deg
+USE lmdz_cloud_optics_prop_ini , ONLY : prmhc, prlmc
+USE setup_config_from_lmdz,   ONLY : driver_config_type
+USE write_field_phy
+! -------------------------------------------------------------------
+
+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
+! AI 04 2024
+INTEGER(KIND=JPIM),INTENT(IN) :: KLEV
+
+! *** Single-level variables 
+!AI 04 2024
+REAL(KIND=JPRB), intent(in)   :: pressure_hl(KLON,KLEV+1)
+type(driver_config_type),intent(in)   :: driver_config
+
+! OUTPUT ARGUMENTS
+
+! *** Decorrelation lengths for cloud edges and cloud water content, in m
+! A.I 04 2024
+REAL(KIND=JPRB), INTENT(OUT)           :: PDECORR_LEN_EDGES_M(KLON,KLEV)
+REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_WATER_M(KLON,KLEV)
+  
+! Ratio of water-content to cloud-edge decorrelation lengths
+REAL(KIND=JPRB), INTENT(OUT), OPTIONAL :: PDECORR_LEN_RATIO
+
+! LOCAL VARIABLES
+
+REAL(KIND=JPRB) :: ZRADIANS_TO_DEGREES, ZABS_LAT_DEG, ZCOS_LAT
+
+! AI 04 2024
+REAL(KIND=JPRB) :: PGEMU(KLON) ! Sine of latitude
+REAL(KIND=JPRB) :: RDECORR_CF
+REAL(KIND=JPRB) :: RDECORR_CW = 1.0_JPRB
+
+INTEGER(KIND=JPIM) :: JL, JK
+
+REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
+
+! -------------------------------------------------------------------
+
+IF (LHOOK) CALL DR_HOOK('CLOUD_OVERLAP_DECORR_LEN',0,ZHOOK_HANDLE)
+  
+! -------------------------------------------------------------------
+print*,'driver_config%kdecolat = ', driver_config%kdecolat
+
+IF (driver_config%kdecolat == 0) THEN
+
+  ! Decorrelation lengths are constant values
+  ! overlap_decorr_length lu dans namelist en m
+  RDECORR_CF = driver_config%overlap_decorr_length / 1000.0_JPRB
+  PDECORR_LEN_EDGES_M(KIDIA:KFDIA,:) = driver_config%overlap_decorr_length
+!  IF (PRESENT(PDECORR_LEN_WATER_M)) THEN
+!    PDECORR_LEN_WATER_M = RDECORR_CW * 1000.0_JPRB
+!  ENDIF
+!  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
+!    PDECORR_LEN_RATIO = RDECORR_CW / RDECORR_CF
+!  ENDIF
+
+ELSE
+
+  ZRADIANS_TO_DEGREES = 180.0_JPRB / RPI
+  PGEMU = SIN(latitude_deg)
+
+  IF (driver_config%kdecolat == 1) THEN
+    ! Shonk et al. (2010) Eq. 13 formula
+    DO JL = KIDIA,KFDIA
+      ZABS_LAT_DEG = ABS(ASIN(PGEMU(JL)) * ZRADIANS_TO_DEGREES)
+      PDECORR_LEN_EDGES_M(JL,:) = (2.899_JPRB - 0.02759_JPRB * ZABS_LAT_DEG) * 1000.0_JPRB
+    ENDDO
+  ELSE ! KDECOLAT == 2
+    DO JL = KIDIA,KFDIA
+      ! Shonk et al. (2010) but smoothed over the equator
+      ZCOS_LAT = COS(ASIN(PGEMU(JL)))
+      PDECORR_LEN_EDGES_M(JL,:) = (0.75_JPRB + 2.149_JPRB * ZCOS_LAT*ZCOS_LAT) * 1000.0_JPRB
+    ENDDO
+  ENDIF
+
+
+  ! Both KDECOLAT = 1 and 2 assume that the decorrelation length for
+  ! cloud water content is half that for cloud edges
+!  IF (PRESENT(PDECORR_LEN_WATER_M)) THEN
+!    PDECORR_LEN_WATER_M(KIDIA:KFDIA,:) = PDECORR_LEN_EDGES_M(KIDIA:KFDIA,:) * 0.5_JPRB 
+!  ENDIF
+
+!  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
+!    PDECORR_LEN_RATIO = 0.5_JPRB
+!  ENDIF
+
+ENDIF
+
+!AI 04 2024
+IF (driver_config%kdecolat == 3) THEN
+  DO jk = 1, klev
+      DO JL = 1, klon
+        IF (pressure_hl(JL,jk)<prmhc) THEN
+          PDECORR_LEN_EDGES_M(JL,jk) = driver_config%high_decorrelation_length
+        ELSE IF (pressure_hl(JL,jk)>=prmhc .AND. pressure_hl(JL,jk)<prlmc) THEN
+          PDECORR_LEN_EDGES_M(JL,jk) = driver_config%mid_decorrelation_length
+        ELSE IF (pressure_hl(JL,jk)>=prlmc) THEN
+          PDECORR_LEN_EDGES_M(JL,jk) = driver_config%low_decorrelation_length
+        ENDIF
+      ENDDO
+  ENDDO
+  
+  ! Both KDECOLAT = 1 and 2 assume that the decorrelation length for
+  ! cloud water content is half that for cloud edges
+!  IF (PRESENT(PDECORR_LEN_WATER_M)) THEN
+!    PDECORR_LEN_WATER_M(KIDIA:KFDIA,:) = PDECORR_LEN_EDGES_M(KIDIA:KFDIA,:) * 0.5_JPRB 
+!  ENDIF
+!  IF (PRESENT(PDECORR_LEN_RATIO)) THEN
+!    PDECORR_LEN_RATIO = 0.5_JPRB
+!  ENDIF
+ENDIF
+!CALL writefield_phy('latitude',latitude_deg,1)
+!CALL writefield_phy('pressure_hl',pressure_hl,klev+1)
+!CALL writefield_phy('Ldecorel',PDECORR_LEN_EDGES_M,klev)
+! -------------------------------------------------------------------
+
+IF (LHOOK) CALL DR_HOOK('CLOUD_OVERLAP_DECORR_LEN',1,ZHOOK_HANDLE)
+
+END SUBROUTINE CALCUL_CLOUD_OVERLAP_DECORR_LEN
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.F90	(revision 4977)
+++ 	(revision )
@@ -1,1359 +1,0 @@
-! 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      
-! Juillet 2023 :
-!             
-! ============================================================================
-module interface_lmdz_ecrad
-
-IMPLICIT NONE        
-
-contains
-
-SUBROUTINE RADIATION_SCHEME &
-! Inputs
-     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
-     &  namelist_file, ok_3Deffect, &
-     &  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 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
-
-! 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(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_3Deffect, 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)
-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)
-
-!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_3Deffect, namelist_file = ', &
-          ok_3Deffect, 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)
-   debut_ecrad=.false. 
-  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)
-
-! 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_3Deffect) 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) 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
- 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,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)
-!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 RADIATION_SCHEME
-
-SUBROUTINE RADIATION_SCHEME_S2 &
-! Inputs
-     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
-     &  namelist_file, ok_3Deffect, &
-     &  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_3Deffect, 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)
-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)
-
-!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_3Deffect, namelist_file = ', &
-          ok_3Deffect, namelist_file
-! 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 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_3Deffect) 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) 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
- 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,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)
-!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 RADIATION_SCHEME_S2
-
-end module interface_lmdz_ecrad
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.f90	(revision 5268)
@@ -0,0 +1,1359 @@
+! 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      
+! Juillet 2023 :
+!             
+! ============================================================================
+module interface_lmdz_ecrad
+
+IMPLICIT NONE        
+
+contains
+
+SUBROUTINE RADIATION_SCHEME &
+! Inputs
+     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
+     &  namelist_file, ok_3Deffect, &
+     &  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 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
+
+! 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(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_3Deffect, 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)
+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)
+
+!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_3Deffect, namelist_file = ', &
+          ok_3Deffect, 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)
+   debut_ecrad=.false. 
+  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)
+
+! 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_3Deffect) 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) 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
+ 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,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)
+!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 RADIATION_SCHEME
+
+SUBROUTINE RADIATION_SCHEME_S2 &
+! Inputs
+     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
+     &  namelist_file, ok_3Deffect, &
+     &  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_3Deffect, 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)
+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)
+
+!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_3Deffect, namelist_file = ', &
+          ok_3Deffect, namelist_file
+! 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 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_3Deffect) 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) 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
+ 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,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)
+!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 RADIATION_SCHEME_S2
+
+end module interface_lmdz_ecrad
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_setup.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_setup.F90	(revision 4977)
+++ 	(revision )
@@ -1,149 +1,0 @@
-MODULE RADIATION_SETUP
-
-! RADIATION_SETUP - Setting up modular radiation scheme
-!
-! AUTHOR
-! ------
-!   Robin Hogan, ECMWF
-!   Original: 2015-09-16
-!
-! MODIFICATIONS
-! -------------
-!   Abderrahmane Idelkadi LMD, juillet 2023
-!
-!-----------------------------------------------------------------------
-
-  USE PARKIND1,         ONLY : JPRB
-  USE radiation_config, ONLY : config_type, &
-       &                       ISolverMcICA, ISolverSpartacus, &
-       &                       ILiquidModelSlingo, ILiquidModelSOCRATES, &
-       &                       IIceModelFu, IIceModelBaran, IGasModelECCKD, &
-       &                       IGasModelIFSRRTMG, IOverlapExponentialRandom
-
-  USE radiation_interface,      ONLY : setup_radiation
-  USE setup_config_from_lmdz,   ONLY : driver_config_type
-
-  IMPLICIT NONE
-
-
-  ! Ultraviolet weightings
-  INTEGER         :: NWEIGHT_UV
-  INTEGER         :: IBAND_UV(100)
-  REAL(KIND=JPRB) :: WEIGHT_UV(100)
-  ! Photosynthetically active radiation weightings
-  INTEGER         :: NWEIGHT_PAR
-  INTEGER         :: IBAND_PAR(100)
-  REAL(KIND=JPRB) :: WEIGHT_PAR(100)
-
-  ! Background aerosol is specified in an ugly way: using the old
-  ! Tegen fields that are in terms of optical depth, and converted to
-  ! mass mixing ratio via the relevant mass-extinction coefficient
-  INTEGER, PARAMETER :: ITYPE_TROP_BG_AER = 8 ! hydrophobic organic
-  INTEGER, PARAMETER :: ITYPE_STRAT_BG_AER=12 ! non-absorbing sulphate
-  REAL(KIND=JPRB)    :: TROP_BG_AER_MASS_EXT
-  REAL(KIND=JPRB)    :: STRAT_BG_AER_MASS_EXT
-
-CONTAINS
-
-  ! This routine copies information between the LMDZ radiation
-  ! configuration (stored in global variables) and the radiation
-  ! configuration of the modular radiation scheme (stored in
-  ! rad_config).  The optional input logical LOUTPUT controls whether
-  ! to print lots of information during the setup stage (default is
-  ! no).
-! AI At the end of the routine, the parameters are read in namelist
-!    
-  SUBROUTINE SETUP_RADIATION_SCHEME(LOUTPUT,file_name,rad_config,driver_config)
-
-!    USE radiation_config, ONLY : config_type, &
-!        &                       ISolverMcICA, ISolverSpartacus, &
-!        &                       ILiquidModelSlingo, ILiquidModelSOCRATES, &
-!        &                       IIceModelFu, IIceModelBaran, &
-!        &                       IOverlapExponentialRandom
-    USE mod_phys_lmdz_para
-      
-    USE YOMHOOK,  ONLY : LHOOK, DR_HOOK
-    USE YOMLUN,   ONLY : NULOUT, NULERR
-    USE YOESRTWN, ONLY : NMPSRTM
-
- !   USE radiation_interface,      ONLY : setup_radiation
- !   USE setup_config_from_lmdz,   ONLY : driver_config_type
-
-    ! Whether or not to provide information on the radiation scheme
-    ! configuration
-    LOGICAL, INTENT(IN), OPTIONAL :: LOUTPUT
-
-    REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-    character(len=512) :: file_name
-
-    logical :: lprint_setp=.TRUE.
-
-    ! Store configuration information for the radiation scheme in a
-    ! global variable
-    type(config_type) :: rad_config
-    type(driver_config_type) :: driver_config
-
-    IF (LHOOK) CALL DR_HOOK('RADIATION_SETUP:SETUP_RADIATION_SCHEME',0,ZHOOK_HANDLE)
-
-    ! *** GENERAL SETUP ***
-
-    ! Configure verbosity of setup of radiation scheme
- 
-    print*,'********** Dans radiation_setup *****************'
-
-! AI ATTENTION (parameters read in namelist file)
-!    file_name="namelist_ecrad"
-    call rad_config%read(file_name=file_name)
-    call driver_config%read(file_name)
-
-    ! Use configuration data to set-up radiation scheme, including
-    ! reading scattering datafiles
-    CALL setup_radiation(rad_config)
-
-    ! Populate the mapping between the 14 RRTM shortwave bands and the
-    ! 6 albedo inputs. The mapping according to the stated wavelength
-    ! ranges of the 6-band model does not match the hard-wired mapping
-    ! in NMPSRTM, but only the hard-wired values produce sensible
-    ! results...
-    ! Note that NMPSRTM(:)=(/  6, 6, 5, 5, 5, 5, 5, 4, 4, 3, 2, 2, 1, 6 /)
-! AI (6 albedo SW bands) 
-    call rad_config%define_sw_albedo_intervals(6, &
-       &  [0.25e-6_jprb, 0.44e-6_jprb, 0.69e-6_jprb, &
-       &   1.19e-6_jprb, 2.38e-6_jprb], [1,2,3,4,5,6])
-    ! Likewise between the 16 RRTM longwave bands and the 2 emissivity
-    ! inputs (info taken from rrtm_ecrt_140gp_mcica.F90) representing
-    ! outside and inside the window region of the spectrum
-!     rad_config%i_emiss_from_band_lw = (/ 1,1,1,1,1,2,2,2,1,1,1,1,1,1,1,1 /)
-! AI ATTENTION ?????
-!!    call rad_config%define_lw_emiss_intervals(3, &
-!!         &  (/ 8.0e-6_jprb,13.0e-6_jprb /),  (/ 1,2,1 /))
-
-!    ! Get spectral weightings for UV and PAR
-    call rad_config%get_sw_weights(0.2e-6_jprb, 0.4415e-6_jprb, &
-         &  NWEIGHT_UV, IBAND_UV, WEIGHT_UV, 'ultraviolet')
-    call rad_config%get_sw_weights(0.4e-6_jprb, 0.7e-6_jprb, &
-         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
-         &  'photosynthetically active radiation, PAR')
-
-     rad_config%i_aerosol_type_map(1:13) = (/ &
-           &  -1, & 
-           &  -2, & 
-           &  -3, & 
-           &  -4, & 
-           &  -5, & 
-           &  -6, &
-           &  -7, &
-           &   1, &
-           &   2, &
-           &   3, &
-           &  -8, &
-           &  -9, &
-           &   4 /)
-
-
-    IF (LHOOK) CALL DR_HOOK('RADIATION_SETUP:SETUP_RADIATION_SCHEME',1,ZHOOK_HANDLE)
-
-  END SUBROUTINE SETUP_RADIATION_SCHEME
-
-END MODULE RADIATION_SETUP
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_setup.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_setup.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_setup.f90	(revision 5268)
@@ -0,0 +1,149 @@
+MODULE RADIATION_SETUP
+
+! RADIATION_SETUP - Setting up modular radiation scheme
+!
+! AUTHOR
+! ------
+!   Robin Hogan, ECMWF
+!   Original: 2015-09-16
+!
+! MODIFICATIONS
+! -------------
+!   Abderrahmane Idelkadi LMD, juillet 2023
+!
+!-----------------------------------------------------------------------
+
+  USE PARKIND1,         ONLY : JPRB
+  USE radiation_config, ONLY : config_type, &
+       &                       ISolverMcICA, ISolverSpartacus, &
+       &                       ILiquidModelSlingo, ILiquidModelSOCRATES, &
+       &                       IIceModelFu, IIceModelBaran, IGasModelECCKD, &
+       &                       IGasModelIFSRRTMG, IOverlapExponentialRandom
+
+  USE radiation_interface,      ONLY : setup_radiation
+  USE setup_config_from_lmdz,   ONLY : driver_config_type
+
+  IMPLICIT NONE
+
+
+  ! Ultraviolet weightings
+  INTEGER         :: NWEIGHT_UV
+  INTEGER         :: IBAND_UV(100)
+  REAL(KIND=JPRB) :: WEIGHT_UV(100)
+  ! Photosynthetically active radiation weightings
+  INTEGER         :: NWEIGHT_PAR
+  INTEGER         :: IBAND_PAR(100)
+  REAL(KIND=JPRB) :: WEIGHT_PAR(100)
+
+  ! Background aerosol is specified in an ugly way: using the old
+  ! Tegen fields that are in terms of optical depth, and converted to
+  ! mass mixing ratio via the relevant mass-extinction coefficient
+  INTEGER, PARAMETER :: ITYPE_TROP_BG_AER = 8 ! hydrophobic organic
+  INTEGER, PARAMETER :: ITYPE_STRAT_BG_AER=12 ! non-absorbing sulphate
+  REAL(KIND=JPRB)    :: TROP_BG_AER_MASS_EXT
+  REAL(KIND=JPRB)    :: STRAT_BG_AER_MASS_EXT
+
+CONTAINS
+
+  ! This routine copies information between the LMDZ radiation
+  ! configuration (stored in global variables) and the radiation
+  ! configuration of the modular radiation scheme (stored in
+  ! rad_config).  The optional input logical LOUTPUT controls whether
+  ! to print lots of information during the setup stage (default is
+  ! no).
+! AI At the end of the routine, the parameters are read in namelist
+!    
+  SUBROUTINE SETUP_RADIATION_SCHEME(LOUTPUT,file_name,rad_config,driver_config)
+
+!    USE radiation_config, ONLY : config_type, &
+!        &                       ISolverMcICA, ISolverSpartacus, &
+!        &                       ILiquidModelSlingo, ILiquidModelSOCRATES, &
+!        &                       IIceModelFu, IIceModelBaran, &
+!        &                       IOverlapExponentialRandom
+    USE mod_phys_lmdz_para
+      
+    USE YOMHOOK,  ONLY : LHOOK, DR_HOOK
+    USE YOMLUN,   ONLY : NULOUT, NULERR
+    USE YOESRTWN, ONLY : NMPSRTM
+
+ !   USE radiation_interface,      ONLY : setup_radiation
+ !   USE setup_config_from_lmdz,   ONLY : driver_config_type
+
+    ! Whether or not to provide information on the radiation scheme
+    ! configuration
+    LOGICAL, INTENT(IN), OPTIONAL :: LOUTPUT
+
+    REAL(KIND=JPRB) :: ZHOOK_HANDLE
+
+    character(len=512) :: file_name
+
+    logical :: lprint_setp=.TRUE.
+
+    ! Store configuration information for the radiation scheme in a
+    ! global variable
+    type(config_type) :: rad_config
+    type(driver_config_type) :: driver_config
+
+    IF (LHOOK) CALL DR_HOOK('RADIATION_SETUP:SETUP_RADIATION_SCHEME',0,ZHOOK_HANDLE)
+
+    ! *** GENERAL SETUP ***
+
+    ! Configure verbosity of setup of radiation scheme
+ 
+    print*,'********** Dans radiation_setup *****************'
+
+! AI ATTENTION (parameters read in namelist file)
+!    file_name="namelist_ecrad"
+    call rad_config%read(file_name=file_name)
+    call driver_config%read(file_name)
+
+    ! Use configuration data to set-up radiation scheme, including
+    ! reading scattering datafiles
+    CALL setup_radiation(rad_config)
+
+    ! Populate the mapping between the 14 RRTM shortwave bands and the
+    ! 6 albedo inputs. The mapping according to the stated wavelength
+    ! ranges of the 6-band model does not match the hard-wired mapping
+    ! in NMPSRTM, but only the hard-wired values produce sensible
+    ! results...
+    ! Note that NMPSRTM(:)=(/  6, 6, 5, 5, 5, 5, 5, 4, 4, 3, 2, 2, 1, 6 /)
+! AI (6 albedo SW bands) 
+    call rad_config%define_sw_albedo_intervals(6, &
+       &  [0.25e-6_jprb, 0.44e-6_jprb, 0.69e-6_jprb, &
+       &   1.19e-6_jprb, 2.38e-6_jprb], [1,2,3,4,5,6])
+    ! Likewise between the 16 RRTM longwave bands and the 2 emissivity
+    ! inputs (info taken from rrtm_ecrt_140gp_mcica.F90) representing
+    ! outside and inside the window region of the spectrum
+!     rad_config%i_emiss_from_band_lw = (/ 1,1,1,1,1,2,2,2,1,1,1,1,1,1,1,1 /)
+! AI ATTENTION ?????
+!!    call rad_config%define_lw_emiss_intervals(3, &
+!!         &  (/ 8.0e-6_jprb,13.0e-6_jprb /),  (/ 1,2,1 /))
+
+!    ! Get spectral weightings for UV and PAR
+    call rad_config%get_sw_weights(0.2e-6_jprb, 0.4415e-6_jprb, &
+         &  NWEIGHT_UV, IBAND_UV, WEIGHT_UV, 'ultraviolet')
+    call rad_config%get_sw_weights(0.4e-6_jprb, 0.7e-6_jprb, &
+         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
+         &  'photosynthetically active radiation, PAR')
+
+     rad_config%i_aerosol_type_map(1:13) = (/ &
+           &  -1, & 
+           &  -2, & 
+           &  -3, & 
+           &  -4, & 
+           &  -5, & 
+           &  -6, &
+           &  -7, &
+           &   1, &
+           &   2, &
+           &   3, &
+           &  -8, &
+           &  -9, &
+           &   4 /)
+
+
+    IF (LHOOK) CALL DR_HOOK('RADIATION_SETUP:SETUP_RADIATION_SCHEME',1,ZHOOK_HANDLE)
+
+  END SUBROUTINE SETUP_RADIATION_SCHEME
+
+END MODULE RADIATION_SETUP
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.F90	(revision 4977)
+++ 	(revision )
@@ -1,314 +1,0 @@
-! $Id: readaerosol_optic_ecrad.F90
-!
-SUBROUTINE readaerosol_optic_ecrad(debut, aerosol_couple, ok_alw, ok_volcan, &
-     flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, &
-     pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
-     tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) 
-!     tau_aero, piz_aero, cg_aero, &
-!     tausum_aero, drytausum_aero, tau3d_aero )
-
-  ! This routine will :
-  ! 1) recevie the aerosols(already read and interpolated) corresponding to flag_aerosol
-  ! 2) calculate the optical properties for the aerosols
-  !
-
-  USE dimphy
-  USE aero_mod
-  USE phys_local_var_mod, only: sconcso4,sconcno3,sconcoa,sconcbc,sconcss,sconcdust, &
-       concso4,concno3,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
-       loadno3,load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7, & 
-       load_tmp8,load_tmp9,load_tmp10
-
-  USE infotrac_phy, ONLY: tracers, nqtot, nbtr
-  USE YOMCST
-
-  IMPLICIT NONE
-
-  include "clesphys.h"
-
-  ! Input arguments
-  !****************************************************************************************
-  LOGICAL, INTENT(IN)                      :: debut
-  LOGICAL, INTENT(IN)                      :: aerosol_couple
-  LOGICAL, INTENT(IN)                      :: ok_alw
-  LOGICAL, INTENT(IN)                      :: ok_volcan
-  INTEGER, INTENT(IN)                      :: flag_aerosol
-  LOGICAL, INTENT(IN)                      :: flag_bc_internal_mixture
-  INTEGER, INTENT(IN)                      :: itap
-  REAL, INTENT(IN)                         :: rjourvrai
-  REAL, INTENT(IN)                         :: pdtphys
-  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
-  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
-  REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_seri
-  REAL, DIMENSION(klon,klev), INTENT(IN)   :: rhcl   ! humidite relative ciel clair
-  REAL, DIMENSION(klev), INTENT(IN)        :: presnivs
-  REAL, DIMENSION(klon,klev,nbtr), INTENT(IN) :: tr_seri ! concentration tracer
-
-  ! Output arguments
-  !****************************************************************************************
-  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero    ! Total mass for all soluble aerosols
-  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero_pi !     -"-     preindustrial values
-  REAL, DIMENSION(klon,klev,naero_tot), INTENT(OUT) :: m_allaer
-  ! AI a passer par la suite en argument si besoin pour ecrad
-  !REAL, DIMENSION(klon,klev,naero_tot), INTENT(OUT) :: m_allaer_pi !RAF
-
-!  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
-!  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
-!  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: cg_aero     ! asymmetry parameter aerosol
-!  REAL, DIMENSION(klon,nwave,naero_tot), INTENT(OUT)       :: tausum_aero
-!  REAL, DIMENSION(klon,naero_tot), INTENT(OUT)             :: drytausum_aero
-!  REAL, DIMENSION(klon,klev,nwave,naero_tot), INTENT(OUT)  :: tau3d_aero
-
-  ! Local variables
-  !****************************************************************************************
-  REAL, DIMENSION(klon)        :: aerindex      ! POLDER aerosol index 
-  REAL, DIMENSION(klon,klev)   :: sulfacc       ! SO4 accumulation concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: sulfcoarse    ! SO4 coarse concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: bcsol         ! BC soluble concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: bcins         ! BC insoluble concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: pomsol        ! POM soluble concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: pomins        ! POM insoluble concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: cidust        ! DUST aerosol concentration  [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: sscoarse      ! SS Coarse concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: sssupco       ! SS Super Coarse concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: ssacu         ! SS Acumulation concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: nitracc       ! nitrate accumulation concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: nitrcoarse    ! nitrate coarse concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: nitrinscoarse ! nitrate insoluble coarse concentration [ug/m3]
-  REAL, DIMENSION(klon,klev)   :: sulfacc_pi
-  REAL, DIMENSION(klon,klev)   :: sulfcoarse_pi
-  REAL, DIMENSION(klon,klev)   :: bcsol_pi
-  REAL, DIMENSION(klon,klev)   :: bcins_pi
-  REAL, DIMENSION(klon,klev)   :: pomsol_pi
-  REAL, DIMENSION(klon,klev)   :: pomins_pi
-  REAL, DIMENSION(klon,klev)   :: cidust_pi
-  REAL, DIMENSION(klon,klev)   :: sscoarse_pi
-  REAL, DIMENSION(klon,klev)   :: sssupco_pi
-  REAL, DIMENSION(klon,klev)   :: ssacu_pi
-  REAL, DIMENSION(klon,klev)   :: nitracc_pi
-  REAL, DIMENSION(klon,klev)   :: nitrcoarse_pi
-  REAL, DIMENSION(klon,klev)   :: nitrinscoarse_pi
-  REAL, DIMENSION(klon,klev)   :: pdel, zrho
-  REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer_pi
-
-  integer :: id_ASBCM, id_ASPOMM, id_ASSO4M, id_ASMSAM, id_CSSO4M, id_CSMSAM, id_SSSSM
-  integer :: id_CSSSM, id_ASSSM, id_CIDUSTM, id_AIBCM, id_AIPOMM, id_ASNO3M, id_CSNO3M, id_CINO3M
-  INTEGER :: k, i, iq, itr
-
-  !--air density
-  zrho(:,:)=pplay(:,:)/t_seri(:,:)/RD                     !--kg/m3
-
-  !****************************************************************************************
-  ! 1) Get aerosol mass
-  !    
-  !****************************************************************************************
-  !
-  !
-  IF (aerosol_couple) THEN   !--we get aerosols from tr_seri array from INCA
-     !
-     !--copy fields from INCA tr_seri 
-     !--convert to ug m-3 unit for consistency with offline fields
-     !
-     itr = 0
-     DO iq = 1,nqtot
-        IF(.NOT. tracers(iq)%isInPhysics) CYCLE
-        itr = itr+1
-        SELECT CASE(trim(tracers(iq)%name))
-           CASE ("ASBCM");  id_ASBCM  = itr
-           CASE ("ASPOMM"); id_ASPOMM = itr
-           CASE ("ASSO4M"); id_ASSO4M = itr
-           CASE ("ASMSAM"); id_ASMSAM = itr
-           CASE ("CSSO4M"); id_CSSO4M = itr
-           CASE ("CSMSAM"); id_CSMSAM = itr
-           CASE ("SSSSM");  id_SSSSM  = itr
-           CASE ("CSSSM");  id_CSSSM  = itr
-           CASE ("ASSSM");  id_ASSSM  = itr
-           CASE ("CIDUSTM");id_CIDUSTM= itr
-           CASE ("AIBCM");  id_AIBCM  = itr
-           CASE ("AIPOMM"); id_AIPOMM = itr
-           CASE ("ASNO3M"); id_ASNO3M = itr
-           CASE ("CSNO3M"); id_CSNO3M = itr
-           CASE ("CINO3M"); id_CINO3M = itr
-        END SELECT
-     END DO
-
-     bcsol(:,:)        =   tr_seri(:,:,id_ASBCM)                         *zrho(:,:)*1.e9  ! ASBCM
-     pomsol(:,:)       =   tr_seri(:,:,id_ASPOMM)                        *zrho(:,:)*1.e9  ! ASPOMM
-     sulfacc(:,:)      =  (tr_seri(:,:,id_ASSO4M)+tr_seri(:,:,id_ASMSAM))*zrho(:,:)*1.e9  ! ASSO4M (=SO4) + ASMSAM (=MSA)
-     sulfcoarse(:,:)   =  (tr_seri(:,:,id_CSSO4M)+tr_seri(:,:,id_CSMSAM))*zrho(:,:)*1.e9  ! CSSO4M (=SO4) + CSMSAM (=MSA)
-     sssupco(:,:)      =   tr_seri(:,:,id_SSSSM)                         *zrho(:,:)*1.e9  ! SSSSM
-     sscoarse(:,:)     =   tr_seri(:,:,id_CSSSM)                         *zrho(:,:)*1.e9  ! CSSSM
-     ssacu(:,:)        =   tr_seri(:,:,id_ASSSM)                         *zrho(:,:)*1.e9  ! ASSSM
-     cidust(:,:)       =   tr_seri(:,:,id_CIDUSTM)                       *zrho(:,:)*1.e9  ! CIDUSTM
-     bcins(:,:)        =   tr_seri(:,:,id_AIBCM)                         *zrho(:,:)*1.e9  ! AIBCM
-     pomins(:,:)       =   tr_seri(:,:,id_AIPOMM)                        *zrho(:,:)*1.e9  ! AIPOMM
-     nitracc(:,:)      =   tr_seri(:,:,id_ASNO3M)                        *zrho(:,:)*1.e9  ! ASNO3M
-     nitrcoarse(:,:)   =   tr_seri(:,:,id_CSNO3M)                        *zrho(:,:)*1.e9  ! CSNO3M
-     nitrinscoarse(:,:)=   tr_seri(:,:,id_CINO3M)                        *zrho(:,:)*1.e9  ! CINO3M
-     !
-     bcsol_pi(:,:)        =   0.0 ! ASBCM pre-ind
-     pomsol_pi(:,:)       =   0.0 ! ASPOMM pre-ind
-     sulfacc_pi(:,:)      =   0.0 ! ASSO4M (=SO4) + ASMSAM (=MSA) pre-ind
-     sulfcoarse_pi(:,:)   =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) pre-ind
-     sssupco_pi(:,:)      =   0.0 ! SSSSM pre-ind
-     sscoarse_pi(:,:)     =   0.0 ! CSSSM pre-ind
-     ssacu_pi(:,:)        =   0.0 ! ASSSM pre-ind
-     cidust_pi(:,:)       =   0.0 ! CIDUSTM pre-ind
-     bcins_pi(:,:)        =   0.0 ! AIBCM pre-ind
-     pomins_pi(:,:)       =   0.0 ! AIPOMM pre-ind
-     nitracc_pi(:,:)      =   0.0 ! ASNO3M pre-ind
-     nitrcoarse_pi(:,:)   =   0.0 ! CSNO3M pre-ind
-     nitrinscoarse_pi(:,:)=   0.0 ! CINO3M
-     !
-  ELSE !--not aerosol_couple
-     !
-     ! Read and interpolate sulfate
-     IF ( flag_aerosol .EQ. 1 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
-
-        CALL readaerosol_interp(id_ASSO4M_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sulfacc, sulfacc_pi,loadso4)
-     ELSE
-        sulfacc(:,:) = 0. ; sulfacc_pi(:,:) = 0.
-        loadso4=0.
-     ENDIF
-
-     ! Read and interpolate bcsol and bcins
-     IF ( flag_aerosol .EQ. 2 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
-
-        ! Get bc aerosol distribution 
-        CALL readaerosol_interp(id_ASBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi, load_tmp1 )
-        CALL readaerosol_interp(id_AIBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi, load_tmp2 )
-        loadbc(:)=load_tmp1(:)+load_tmp2(:)
-     ELSE
-        bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
-        bcins(:,:) = 0. ; bcins_pi(:,:) = 0.
-        loadbc=0.
-     ENDIF
-
-     ! Read and interpolate pomsol and pomins
-     IF ( flag_aerosol .EQ. 3 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN
-
-        CALL readaerosol_interp(id_ASPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp3)
-        CALL readaerosol_interp(id_AIPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp4)
-        loadoa(:)=load_tmp3(:)+load_tmp4(:)
-     ELSE
-        pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
-        pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
-        loadoa=0.
-     ENDIF
-
-     ! Read and interpolate csssm, ssssm, assssm
-     IF (flag_aerosol .EQ. 4 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
-
-        CALL readaerosol_interp(id_SSSSM_phy ,itap, pdtphys,rjourvrai, &
-        debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp5) 
-        CALL readaerosol_interp(id_CSSSM_phy ,itap, pdtphys,rjourvrai, &
-        debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp6) 
-        CALL readaerosol_interp(id_ASSSM_phy ,itap, pdtphys,rjourvrai, &
-        debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp7) 
-        loadss(:)=load_tmp5(:)+load_tmp6(:)+load_tmp7(:)
-     ELSE
-        sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0. 
-        ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0. 
-        sssupco(:,:)  = 0. ; sssupco_pi = 0. 
-        loadss=0.
-     ENDIF
-
-     ! Read and interpolate cidustm
-     IF (flag_aerosol .EQ. 5 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
-
-        CALL readaerosol_interp(id_CIDUSTM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi, loaddust) 
-
-     ELSE
-        cidust(:,:) = 0. ; cidust_pi(:,:) = 0. 
-        loaddust=0.
-     ENDIF
-     !
-     ! Read and interpolate asno3m, csno3m, cino3m
-     IF (flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
-
-        CALL readaerosol_interp(id_ASNO3M_phy, itap, pdtphys, rjourvrai, & 
-        debut, pplay, paprs, t_seri, nitracc, nitracc_pi, load_tmp8) 
-        CALL readaerosol_interp(id_CSNO3M_phy, itap, pdtphys, rjourvrai, & 
-        debut, pplay, paprs, t_seri, nitrcoarse, nitrcoarse_pi, load_tmp9) 
-        CALL readaerosol_interp(id_CINO3M_phy, itap, pdtphys, rjourvrai, & 
-        debut, pplay, paprs, t_seri, nitrinscoarse, nitrinscoarse_pi, load_tmp10) 
-        loadno3(:)=load_tmp8(:)+load_tmp9(:)+load_tmp10(:)
-
-     ELSE
-        nitracc(:,:)         =   0.0 ; nitracc_pi(:,:)      =   0.0 
-        nitrcoarse(:,:)      =   0.0 ; nitrcoarse_pi(:,:)   =   0.0
-        nitrinscoarse(:,:)   =   0.0 ; nitrinscoarse_pi(:,:)=   0.0
-        loadno3(:)=0.0
-     ENDIF
-     !
-     ! CSSO4M is set to 0 as not reliable
-     sulfcoarse(:,:)      =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) 
-     sulfcoarse_pi(:,:)   =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) pre-ind
-
-  ENDIF !--not aerosol_couple
-
-  !
-  ! Store all aerosols mixing ratios in one variable for radiation scheme (unit kg/kg for ECRAD)
-  ! present-day values
-  m_allaer(:,:,id_ASBCM_phy)  = bcsol(:,:)        /1.e9/zrho(:,:) ! ASBCM
-  m_allaer(:,:,id_ASPOMM_phy) = pomsol(:,:)       /1.e9/zrho(:,:) ! ASPOMM
-  m_allaer(:,:,id_ASSO4M_phy) = sulfacc(:,:)      /1.e9/zrho(:,:) ! ASSO4M (= SO4) 
-  m_allaer(:,:,id_CSSO4M_phy) = sulfcoarse(:,:)   /1.e9/zrho(:,:) ! CSSO4M 
-  m_allaer(:,:,id_SSSSM_phy)  = sssupco(:,:)      /1.e9/zrho(:,:) ! SSSSM
-  m_allaer(:,:,id_CSSSM_phy)  = sscoarse(:,:)     /1.e9/zrho(:,:) ! CSSSM
-  m_allaer(:,:,id_ASSSM_phy)  = ssacu(:,:)        /1.e9/zrho(:,:) ! ASSSM
-  m_allaer(:,:,id_CIDUSTM_phy)= cidust(:,:)       /1.e9/zrho(:,:) ! CIDUSTM
-  m_allaer(:,:,id_AIBCM_phy)  = bcins(:,:)        /1.e9/zrho(:,:) ! AIBCM
-  m_allaer(:,:,id_ASNO3M_phy) = nitracc(:,:)      /1.e9/zrho(:,:) ! ASNO3M
-  m_allaer(:,:,id_CSNO3M_phy) = nitrcoarse(:,:)   /1.e9/zrho(:,:) ! CSNO3M
-  m_allaer(:,:,id_CINO3M_phy) = nitrinscoarse(:,:)/1.e9/zrho(:,:) ! CINO3M
-  m_allaer(:,:,id_AIPOMM_phy) = pomins(:,:)       /1.e9/zrho(:,:) ! AIPOMM
-  m_allaer(:,:,id_STRAT_phy)  = 0.0
-
-  !  pre-industrial (pi) values
-  m_allaer_pi(:,:,id_ASBCM_phy)  = bcsol_pi(:,:)        /1.e9/zrho(:,:) ! ASBCM pre-ind
-  m_allaer_pi(:,:,id_ASPOMM_phy) = pomsol_pi(:,:)       /1.e9/zrho(:,:) ! ASPOMM pre-ind
-  m_allaer_pi(:,:,id_ASSO4M_phy) = sulfacc_pi(:,:)      /1.e9/zrho(:,:) ! ASSO4M (= SO4) pre-ind
-  m_allaer_pi(:,:,id_CSSO4M_phy) = sulfcoarse_pi(:,:)   /1.e9/zrho(:,:) ! CSSO4M pre-ind
-  m_allaer_pi(:,:,id_SSSSM_phy)  = sssupco_pi(:,:)      /1.e9/zrho(:,:) ! SSSSM pre-ind
-  m_allaer_pi(:,:,id_CSSSM_phy)  = sscoarse_pi(:,:)     /1.e9/zrho(:,:) ! CSSSM pre-ind
-  m_allaer_pi(:,:,id_ASSSM_phy)  = ssacu_pi(:,:)        /1.e9/zrho(:,:) ! ASSSM pre-ind
-  m_allaer_pi(:,:,id_CIDUSTM_phy)= cidust_pi(:,:)       /1.e9/zrho(:,:) ! CIDUSTM pre-ind
-  m_allaer_pi(:,:,id_AIBCM_phy)  = bcins_pi(:,:)        /1.e9/zrho(:,:) ! AIBCM pre-ind
-  m_allaer_pi(:,:,id_ASNO3M_phy) = nitracc_pi(:,:)      /1.e9/zrho(:,:) ! ASNO3M pre-ind
-  m_allaer_pi(:,:,id_CSNO3M_phy) = nitrcoarse_pi(:,:)   /1.e9/zrho(:,:) ! CSNO3M pre-ind
-  m_allaer_pi(:,:,id_CINO3M_phy) = nitrinscoarse_pi(:,:)/1.e9/zrho(:,:) ! CINO3M pre-ind
-  m_allaer_pi(:,:,id_AIPOMM_phy) = pomins_pi(:,:)       /1.e9/zrho(:,:) ! AIPOMM pre-ind
-  m_allaer_pi(:,:,id_STRAT_phy)  = 0.0
-
-  !
-  ! Calculate the total mass of all soluble aersosols (in unit ug /m3)
-  ! to be revisited for AR6
-  mass_solu_aero(:,:)    = sulfacc(:,:)    + bcsol(:,:)    + pomsol(:,:)    + nitracc(:,:)    + ssacu(:,:)
-  mass_solu_aero_pi(:,:) = sulfacc_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:) + nitracc_pi(:,:) + ssacu_pi(:,:)
-
-  !****************************************************************************************
-  ! 2) Calculate optical properties for the aerosols
-  !
-  !****************************************************************************************
-  DO k = 1, klev
-     DO i = 1, klon
-        pdel(i,k) = paprs(i,k) - paprs (i,k+1)
-     ENDDO
-  ENDDO
-
-  ! Diagnostics calculation for CMIP5 protocol unit (in unit kg/m3)
-  sconcso4(:)  =m_allaer(:,1,id_ASSO4M_phy)*1.e-9
-  sconcno3(:)  =(m_allaer(:,1,id_ASNO3M_phy)+m_allaer(:,1,id_CSNO3M_phy)+m_allaer(:,1,id_CINO3M_phy))*1.e-9
-  sconcoa(:)   =(m_allaer(:,1,id_ASPOMM_phy)+m_allaer(:,1,id_AIPOMM_phy))*1.e-9
-  sconcbc(:)   =(m_allaer(:,1,id_ASBCM_phy)+m_allaer(:,1,id_AIBCM_phy))*1.e-9
-  sconcss(:)   =(m_allaer(:,1,id_ASSSM_phy)+m_allaer(:,1,id_CSSSM_phy)+m_allaer(:,1,id_SSSSM_phy))*1.e-9
-  sconcdust(:) =m_allaer(:,1,id_CIDUSTM_phy)*1.e-9
-  concso4(:,:) =m_allaer(:,:,id_ASSO4M_phy)*1.e-9
-  concno3(:,:) =(m_allaer(:,:,id_ASNO3M_phy)+m_allaer(:,:,id_CSNO3M_phy)+m_allaer(:,:,id_CINO3M_phy))*1.e-9
-  concoa(:,:)  =(m_allaer(:,:,id_ASPOMM_phy)+m_allaer(:,:,id_AIPOMM_phy))*1.e-9
-  concbc(:,:)  =(m_allaer(:,:,id_ASBCM_phy)+m_allaer(:,:,id_AIBCM_phy))*1.e-9
-  concss(:,:)  =(m_allaer(:,:,id_ASSSM_phy)+m_allaer(:,:,id_CSSSM_phy)+m_allaer(:,:,id_SSSSM_phy))*1.e-9
-  concdust(:,:)=m_allaer(:,:,id_CIDUSTM_phy)*1.e-9
-
-END SUBROUTINE readaerosol_optic_ecrad
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.f90	(revision 5268)
@@ -0,0 +1,314 @@
+! $Id: readaerosol_optic_ecrad.F90
+!
+SUBROUTINE readaerosol_optic_ecrad(debut, aerosol_couple, ok_alw, ok_volcan, &
+     flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, &
+     pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
+     tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) 
+!     tau_aero, piz_aero, cg_aero, &
+!     tausum_aero, drytausum_aero, tau3d_aero )
+
+  ! This routine will :
+  ! 1) recevie the aerosols(already read and interpolated) corresponding to flag_aerosol
+  ! 2) calculate the optical properties for the aerosols
+  !
+
+  USE dimphy
+  USE aero_mod
+  USE phys_local_var_mod, only: sconcso4,sconcno3,sconcoa,sconcbc,sconcss,sconcdust, &
+       concso4,concno3,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
+       loadno3,load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7, & 
+       load_tmp8,load_tmp9,load_tmp10
+
+  USE infotrac_phy, ONLY: tracers, nqtot, nbtr
+  USE YOMCST
+
+  IMPLICIT NONE
+
+  include "clesphys.h"
+
+  ! Input arguments
+  !****************************************************************************************
+  LOGICAL, INTENT(IN)                      :: debut
+  LOGICAL, INTENT(IN)                      :: aerosol_couple
+  LOGICAL, INTENT(IN)                      :: ok_alw
+  LOGICAL, INTENT(IN)                      :: ok_volcan
+  INTEGER, INTENT(IN)                      :: flag_aerosol
+  LOGICAL, INTENT(IN)                      :: flag_bc_internal_mixture
+  INTEGER, INTENT(IN)                      :: itap
+  REAL, INTENT(IN)                         :: rjourvrai
+  REAL, INTENT(IN)                         :: pdtphys
+  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
+  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
+  REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_seri
+  REAL, DIMENSION(klon,klev), INTENT(IN)   :: rhcl   ! humidite relative ciel clair
+  REAL, DIMENSION(klev), INTENT(IN)        :: presnivs
+  REAL, DIMENSION(klon,klev,nbtr), INTENT(IN) :: tr_seri ! concentration tracer
+
+  ! Output arguments
+  !****************************************************************************************
+  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero    ! Total mass for all soluble aerosols
+  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero_pi !     -"-     preindustrial values
+  REAL, DIMENSION(klon,klev,naero_tot), INTENT(OUT) :: m_allaer
+  ! AI a passer par la suite en argument si besoin pour ecrad
+  !REAL, DIMENSION(klon,klev,naero_tot), INTENT(OUT) :: m_allaer_pi !RAF
+
+!  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
+!  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
+!  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: cg_aero     ! asymmetry parameter aerosol
+!  REAL, DIMENSION(klon,nwave,naero_tot), INTENT(OUT)       :: tausum_aero
+!  REAL, DIMENSION(klon,naero_tot), INTENT(OUT)             :: drytausum_aero
+!  REAL, DIMENSION(klon,klev,nwave,naero_tot), INTENT(OUT)  :: tau3d_aero
+
+  ! Local variables
+  !****************************************************************************************
+  REAL, DIMENSION(klon)        :: aerindex      ! POLDER aerosol index 
+  REAL, DIMENSION(klon,klev)   :: sulfacc       ! SO4 accumulation concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sulfcoarse    ! SO4 coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: bcsol         ! BC soluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: bcins         ! BC insoluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: pomsol        ! POM soluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: pomins        ! POM insoluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: cidust        ! DUST aerosol concentration  [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sscoarse      ! SS Coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sssupco       ! SS Super Coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: ssacu         ! SS Acumulation concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: nitracc       ! nitrate accumulation concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: nitrcoarse    ! nitrate coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: nitrinscoarse ! nitrate insoluble coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sulfacc_pi
+  REAL, DIMENSION(klon,klev)   :: sulfcoarse_pi
+  REAL, DIMENSION(klon,klev)   :: bcsol_pi
+  REAL, DIMENSION(klon,klev)   :: bcins_pi
+  REAL, DIMENSION(klon,klev)   :: pomsol_pi
+  REAL, DIMENSION(klon,klev)   :: pomins_pi
+  REAL, DIMENSION(klon,klev)   :: cidust_pi
+  REAL, DIMENSION(klon,klev)   :: sscoarse_pi
+  REAL, DIMENSION(klon,klev)   :: sssupco_pi
+  REAL, DIMENSION(klon,klev)   :: ssacu_pi
+  REAL, DIMENSION(klon,klev)   :: nitracc_pi
+  REAL, DIMENSION(klon,klev)   :: nitrcoarse_pi
+  REAL, DIMENSION(klon,klev)   :: nitrinscoarse_pi
+  REAL, DIMENSION(klon,klev)   :: pdel, zrho
+  REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer_pi
+
+  integer :: id_ASBCM, id_ASPOMM, id_ASSO4M, id_ASMSAM, id_CSSO4M, id_CSMSAM, id_SSSSM
+  integer :: id_CSSSM, id_ASSSM, id_CIDUSTM, id_AIBCM, id_AIPOMM, id_ASNO3M, id_CSNO3M, id_CINO3M
+  INTEGER :: k, i, iq, itr
+
+  !--air density
+  zrho(:,:)=pplay(:,:)/t_seri(:,:)/RD                     !--kg/m3
+
+  !****************************************************************************************
+  ! 1) Get aerosol mass
+  !    
+  !****************************************************************************************
+  !
+  !
+  IF (aerosol_couple) THEN   !--we get aerosols from tr_seri array from INCA
+     !
+     !--copy fields from INCA tr_seri 
+     !--convert to ug m-3 unit for consistency with offline fields
+     !
+     itr = 0
+     DO iq = 1,nqtot
+        IF(.NOT. tracers(iq)%isInPhysics) CYCLE
+        itr = itr+1
+        SELECT CASE(trim(tracers(iq)%name))
+           CASE ("ASBCM");  id_ASBCM  = itr
+           CASE ("ASPOMM"); id_ASPOMM = itr
+           CASE ("ASSO4M"); id_ASSO4M = itr
+           CASE ("ASMSAM"); id_ASMSAM = itr
+           CASE ("CSSO4M"); id_CSSO4M = itr
+           CASE ("CSMSAM"); id_CSMSAM = itr
+           CASE ("SSSSM");  id_SSSSM  = itr
+           CASE ("CSSSM");  id_CSSSM  = itr
+           CASE ("ASSSM");  id_ASSSM  = itr
+           CASE ("CIDUSTM");id_CIDUSTM= itr
+           CASE ("AIBCM");  id_AIBCM  = itr
+           CASE ("AIPOMM"); id_AIPOMM = itr
+           CASE ("ASNO3M"); id_ASNO3M = itr
+           CASE ("CSNO3M"); id_CSNO3M = itr
+           CASE ("CINO3M"); id_CINO3M = itr
+        END SELECT
+     END DO
+
+     bcsol(:,:)        =   tr_seri(:,:,id_ASBCM)                         *zrho(:,:)*1.e9  ! ASBCM
+     pomsol(:,:)       =   tr_seri(:,:,id_ASPOMM)                        *zrho(:,:)*1.e9  ! ASPOMM
+     sulfacc(:,:)      =  (tr_seri(:,:,id_ASSO4M)+tr_seri(:,:,id_ASMSAM))*zrho(:,:)*1.e9  ! ASSO4M (=SO4) + ASMSAM (=MSA)
+     sulfcoarse(:,:)   =  (tr_seri(:,:,id_CSSO4M)+tr_seri(:,:,id_CSMSAM))*zrho(:,:)*1.e9  ! CSSO4M (=SO4) + CSMSAM (=MSA)
+     sssupco(:,:)      =   tr_seri(:,:,id_SSSSM)                         *zrho(:,:)*1.e9  ! SSSSM
+     sscoarse(:,:)     =   tr_seri(:,:,id_CSSSM)                         *zrho(:,:)*1.e9  ! CSSSM
+     ssacu(:,:)        =   tr_seri(:,:,id_ASSSM)                         *zrho(:,:)*1.e9  ! ASSSM
+     cidust(:,:)       =   tr_seri(:,:,id_CIDUSTM)                       *zrho(:,:)*1.e9  ! CIDUSTM
+     bcins(:,:)        =   tr_seri(:,:,id_AIBCM)                         *zrho(:,:)*1.e9  ! AIBCM
+     pomins(:,:)       =   tr_seri(:,:,id_AIPOMM)                        *zrho(:,:)*1.e9  ! AIPOMM
+     nitracc(:,:)      =   tr_seri(:,:,id_ASNO3M)                        *zrho(:,:)*1.e9  ! ASNO3M
+     nitrcoarse(:,:)   =   tr_seri(:,:,id_CSNO3M)                        *zrho(:,:)*1.e9  ! CSNO3M
+     nitrinscoarse(:,:)=   tr_seri(:,:,id_CINO3M)                        *zrho(:,:)*1.e9  ! CINO3M
+     !
+     bcsol_pi(:,:)        =   0.0 ! ASBCM pre-ind
+     pomsol_pi(:,:)       =   0.0 ! ASPOMM pre-ind
+     sulfacc_pi(:,:)      =   0.0 ! ASSO4M (=SO4) + ASMSAM (=MSA) pre-ind
+     sulfcoarse_pi(:,:)   =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) pre-ind
+     sssupco_pi(:,:)      =   0.0 ! SSSSM pre-ind
+     sscoarse_pi(:,:)     =   0.0 ! CSSSM pre-ind
+     ssacu_pi(:,:)        =   0.0 ! ASSSM pre-ind
+     cidust_pi(:,:)       =   0.0 ! CIDUSTM pre-ind
+     bcins_pi(:,:)        =   0.0 ! AIBCM pre-ind
+     pomins_pi(:,:)       =   0.0 ! AIPOMM pre-ind
+     nitracc_pi(:,:)      =   0.0 ! ASNO3M pre-ind
+     nitrcoarse_pi(:,:)   =   0.0 ! CSNO3M pre-ind
+     nitrinscoarse_pi(:,:)=   0.0 ! CINO3M
+     !
+  ELSE !--not aerosol_couple
+     !
+     ! Read and interpolate sulfate
+     IF ( flag_aerosol .EQ. 1 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
+
+        CALL readaerosol_interp(id_ASSO4M_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sulfacc, sulfacc_pi,loadso4)
+     ELSE
+        sulfacc(:,:) = 0. ; sulfacc_pi(:,:) = 0.
+        loadso4=0.
+     ENDIF
+
+     ! Read and interpolate bcsol and bcins
+     IF ( flag_aerosol .EQ. 2 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
+
+        ! Get bc aerosol distribution 
+        CALL readaerosol_interp(id_ASBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi, load_tmp1 )
+        CALL readaerosol_interp(id_AIBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi, load_tmp2 )
+        loadbc(:)=load_tmp1(:)+load_tmp2(:)
+     ELSE
+        bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
+        bcins(:,:) = 0. ; bcins_pi(:,:) = 0.
+        loadbc=0.
+     ENDIF
+
+     ! Read and interpolate pomsol and pomins
+     IF ( flag_aerosol .EQ. 3 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN
+
+        CALL readaerosol_interp(id_ASPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp3)
+        CALL readaerosol_interp(id_AIPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp4)
+        loadoa(:)=load_tmp3(:)+load_tmp4(:)
+     ELSE
+        pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
+        pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
+        loadoa=0.
+     ENDIF
+
+     ! Read and interpolate csssm, ssssm, assssm
+     IF (flag_aerosol .EQ. 4 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
+
+        CALL readaerosol_interp(id_SSSSM_phy ,itap, pdtphys,rjourvrai, &
+        debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp5) 
+        CALL readaerosol_interp(id_CSSSM_phy ,itap, pdtphys,rjourvrai, &
+        debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp6) 
+        CALL readaerosol_interp(id_ASSSM_phy ,itap, pdtphys,rjourvrai, &
+        debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp7) 
+        loadss(:)=load_tmp5(:)+load_tmp6(:)+load_tmp7(:)
+     ELSE
+        sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0. 
+        ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0. 
+        sssupco(:,:)  = 0. ; sssupco_pi = 0. 
+        loadss=0.
+     ENDIF
+
+     ! Read and interpolate cidustm
+     IF (flag_aerosol .EQ. 5 .OR. flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
+
+        CALL readaerosol_interp(id_CIDUSTM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi, loaddust) 
+
+     ELSE
+        cidust(:,:) = 0. ; cidust_pi(:,:) = 0. 
+        loaddust=0.
+     ENDIF
+     !
+     ! Read and interpolate asno3m, csno3m, cino3m
+     IF (flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
+
+        CALL readaerosol_interp(id_ASNO3M_phy, itap, pdtphys, rjourvrai, & 
+        debut, pplay, paprs, t_seri, nitracc, nitracc_pi, load_tmp8) 
+        CALL readaerosol_interp(id_CSNO3M_phy, itap, pdtphys, rjourvrai, & 
+        debut, pplay, paprs, t_seri, nitrcoarse, nitrcoarse_pi, load_tmp9) 
+        CALL readaerosol_interp(id_CINO3M_phy, itap, pdtphys, rjourvrai, & 
+        debut, pplay, paprs, t_seri, nitrinscoarse, nitrinscoarse_pi, load_tmp10) 
+        loadno3(:)=load_tmp8(:)+load_tmp9(:)+load_tmp10(:)
+
+     ELSE
+        nitracc(:,:)         =   0.0 ; nitracc_pi(:,:)      =   0.0 
+        nitrcoarse(:,:)      =   0.0 ; nitrcoarse_pi(:,:)   =   0.0
+        nitrinscoarse(:,:)   =   0.0 ; nitrinscoarse_pi(:,:)=   0.0
+        loadno3(:)=0.0
+     ENDIF
+     !
+     ! CSSO4M is set to 0 as not reliable
+     sulfcoarse(:,:)      =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) 
+     sulfcoarse_pi(:,:)   =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) pre-ind
+
+  ENDIF !--not aerosol_couple
+
+  !
+  ! Store all aerosols mixing ratios in one variable for radiation scheme (unit kg/kg for ECRAD)
+  ! present-day values
+  m_allaer(:,:,id_ASBCM_phy)  = bcsol(:,:)        /1.e9/zrho(:,:) ! ASBCM
+  m_allaer(:,:,id_ASPOMM_phy) = pomsol(:,:)       /1.e9/zrho(:,:) ! ASPOMM
+  m_allaer(:,:,id_ASSO4M_phy) = sulfacc(:,:)      /1.e9/zrho(:,:) ! ASSO4M (= SO4) 
+  m_allaer(:,:,id_CSSO4M_phy) = sulfcoarse(:,:)   /1.e9/zrho(:,:) ! CSSO4M 
+  m_allaer(:,:,id_SSSSM_phy)  = sssupco(:,:)      /1.e9/zrho(:,:) ! SSSSM
+  m_allaer(:,:,id_CSSSM_phy)  = sscoarse(:,:)     /1.e9/zrho(:,:) ! CSSSM
+  m_allaer(:,:,id_ASSSM_phy)  = ssacu(:,:)        /1.e9/zrho(:,:) ! ASSSM
+  m_allaer(:,:,id_CIDUSTM_phy)= cidust(:,:)       /1.e9/zrho(:,:) ! CIDUSTM
+  m_allaer(:,:,id_AIBCM_phy)  = bcins(:,:)        /1.e9/zrho(:,:) ! AIBCM
+  m_allaer(:,:,id_ASNO3M_phy) = nitracc(:,:)      /1.e9/zrho(:,:) ! ASNO3M
+  m_allaer(:,:,id_CSNO3M_phy) = nitrcoarse(:,:)   /1.e9/zrho(:,:) ! CSNO3M
+  m_allaer(:,:,id_CINO3M_phy) = nitrinscoarse(:,:)/1.e9/zrho(:,:) ! CINO3M
+  m_allaer(:,:,id_AIPOMM_phy) = pomins(:,:)       /1.e9/zrho(:,:) ! AIPOMM
+  m_allaer(:,:,id_STRAT_phy)  = 0.0
+
+  !  pre-industrial (pi) values
+  m_allaer_pi(:,:,id_ASBCM_phy)  = bcsol_pi(:,:)        /1.e9/zrho(:,:) ! ASBCM pre-ind
+  m_allaer_pi(:,:,id_ASPOMM_phy) = pomsol_pi(:,:)       /1.e9/zrho(:,:) ! ASPOMM pre-ind
+  m_allaer_pi(:,:,id_ASSO4M_phy) = sulfacc_pi(:,:)      /1.e9/zrho(:,:) ! ASSO4M (= SO4) pre-ind
+  m_allaer_pi(:,:,id_CSSO4M_phy) = sulfcoarse_pi(:,:)   /1.e9/zrho(:,:) ! CSSO4M pre-ind
+  m_allaer_pi(:,:,id_SSSSM_phy)  = sssupco_pi(:,:)      /1.e9/zrho(:,:) ! SSSSM pre-ind
+  m_allaer_pi(:,:,id_CSSSM_phy)  = sscoarse_pi(:,:)     /1.e9/zrho(:,:) ! CSSSM pre-ind
+  m_allaer_pi(:,:,id_ASSSM_phy)  = ssacu_pi(:,:)        /1.e9/zrho(:,:) ! ASSSM pre-ind
+  m_allaer_pi(:,:,id_CIDUSTM_phy)= cidust_pi(:,:)       /1.e9/zrho(:,:) ! CIDUSTM pre-ind
+  m_allaer_pi(:,:,id_AIBCM_phy)  = bcins_pi(:,:)        /1.e9/zrho(:,:) ! AIBCM pre-ind
+  m_allaer_pi(:,:,id_ASNO3M_phy) = nitracc_pi(:,:)      /1.e9/zrho(:,:) ! ASNO3M pre-ind
+  m_allaer_pi(:,:,id_CSNO3M_phy) = nitrcoarse_pi(:,:)   /1.e9/zrho(:,:) ! CSNO3M pre-ind
+  m_allaer_pi(:,:,id_CINO3M_phy) = nitrinscoarse_pi(:,:)/1.e9/zrho(:,:) ! CINO3M pre-ind
+  m_allaer_pi(:,:,id_AIPOMM_phy) = pomins_pi(:,:)       /1.e9/zrho(:,:) ! AIPOMM pre-ind
+  m_allaer_pi(:,:,id_STRAT_phy)  = 0.0
+
+  !
+  ! Calculate the total mass of all soluble aersosols (in unit ug /m3)
+  ! to be revisited for AR6
+  mass_solu_aero(:,:)    = sulfacc(:,:)    + bcsol(:,:)    + pomsol(:,:)    + nitracc(:,:)    + ssacu(:,:)
+  mass_solu_aero_pi(:,:) = sulfacc_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:) + nitracc_pi(:,:) + ssacu_pi(:,:)
+
+  !****************************************************************************************
+  ! 2) Calculate optical properties for the aerosols
+  !
+  !****************************************************************************************
+  DO k = 1, klev
+     DO i = 1, klon
+        pdel(i,k) = paprs(i,k) - paprs (i,k+1)
+     ENDDO
+  ENDDO
+
+  ! Diagnostics calculation for CMIP5 protocol unit (in unit kg/m3)
+  sconcso4(:)  =m_allaer(:,1,id_ASSO4M_phy)*1.e-9
+  sconcno3(:)  =(m_allaer(:,1,id_ASNO3M_phy)+m_allaer(:,1,id_CSNO3M_phy)+m_allaer(:,1,id_CINO3M_phy))*1.e-9
+  sconcoa(:)   =(m_allaer(:,1,id_ASPOMM_phy)+m_allaer(:,1,id_AIPOMM_phy))*1.e-9
+  sconcbc(:)   =(m_allaer(:,1,id_ASBCM_phy)+m_allaer(:,1,id_AIBCM_phy))*1.e-9
+  sconcss(:)   =(m_allaer(:,1,id_ASSSM_phy)+m_allaer(:,1,id_CSSSM_phy)+m_allaer(:,1,id_SSSSM_phy))*1.e-9
+  sconcdust(:) =m_allaer(:,1,id_CIDUSTM_phy)*1.e-9
+  concso4(:,:) =m_allaer(:,:,id_ASSO4M_phy)*1.e-9
+  concno3(:,:) =(m_allaer(:,:,id_ASNO3M_phy)+m_allaer(:,:,id_CSNO3M_phy)+m_allaer(:,:,id_CINO3M_phy))*1.e-9
+  concoa(:,:)  =(m_allaer(:,:,id_ASPOMM_phy)+m_allaer(:,:,id_AIPOMM_phy))*1.e-9
+  concbc(:,:)  =(m_allaer(:,:,id_ASBCM_phy)+m_allaer(:,:,id_AIBCM_phy))*1.e-9
+  concss(:,:)  =(m_allaer(:,:,id_ASSSM_phy)+m_allaer(:,:,id_CSSSM_phy)+m_allaer(:,:,id_SSSSM_phy))*1.e-9
+  concdust(:,:)=m_allaer(:,:,id_CIDUSTM_phy)*1.e-9
+
+END SUBROUTINE readaerosol_optic_ecrad
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_aerosol_optics_lmdz_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_aerosol_optics_lmdz_m.F90	(revision 4977)
+++ 	(revision )
@@ -1,124 +1,0 @@
-module setup_aerosol_optics_lmdz_m
-
-  implicit none
-
-contains
-
-  subroutine setup_aerosol_optics_lmdz(ao, file_name)
-
-    ! Read aerosol optical properties. Note differences with
-    ! "radiation_aerosol_optics_data::setup_aerosol_optics":
-
-    ! -- The input NetCDF file is not flat, it contains NetCDF groups.
-
-    ! -- We do not define ao%ssa_mono_phobic, ao%g_mono_phobic,
-    ! ao%lidar_ratio_mono_phobic, ao%ssa_mono_philic,
-    ! ao%g_mono_philic, ao%lidar_ratio_mono_philic. They are not in
-    ! the input NetCDF file and they are not used by ECRad.
-
-    ! -- We do not define ao%description_phobic_str and
-    ! ao%description_philic_str. We just leave the initialization
-    ! value, which is a blank.
-
-    ! -- We have to cshift the shortwave fields because the the
-    ! shortwave bands are in ascending order in the NetCDF file while
-    ! they are not in ECRad.
-
-    use radiation_aerosol_optics_data, only: aerosol_optics_type, &
-         IAerosolClassUndefined
-    use netcdf95, only: nf95_open, nf95_inq_grp_full_ncid, nf95_close, &
-         nf95_inq_dimid, nf95_inq_varid, nf95_inquire_dimension, &
-         nf95_get_var, nf95_gw_var, nf95_nowrite
-
-    type(aerosol_optics_type), intent(out):: ao
-
-    character(len=*), intent(in):: file_name
-    ! NetCDF file containing the aerosol optics data
-
-    ! Local:
-    integer ncid, grpid, dimid, varid
-
-    !-----------------------------------------------------------------------
-
-    ao%use_hydrophilic = .true.
-    ao%use_monochromatic = .true.
-    call nf95_open(file_name, nf95_nowrite, ncid)
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic", grpid)
-    call nf95_inq_dimid(grpid, "hur", dimid)
-    call nf95_inquire_dimension(grpid, dimid, nclen = ao%nrh)
-    allocate(ao%rh_lower(ao%nrh))
-    call nf95_inq_varid(grpid, "hur_bounds", varid)
-    call nf95_get_var(grpid, varid, ao%rh_lower, count_nc = [1, ao%nrh])
-
-    ! Hydrophilic/LW_bands:
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic/LW_bands", grpid)
-    call nf95_inq_varid(grpid, "asymmetry", varid)
-    call nf95_gw_var(grpid, varid, ao%g_lw_philic)
-    call nf95_inq_varid(grpid, "single_scat_alb", varid)
-    call nf95_gw_var(grpid, varid, ao%ssa_lw_philic)
-    call nf95_inq_varid(grpid, "mass_ext", varid)
-    call nf95_gw_var(grpid, varid, ao%mass_ext_lw_philic)
-
-    ! Hydrophilic/SW_bands:
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic/SW_bands", grpid)
-    call nf95_inq_varid(grpid, "asymmetry", varid)
-    call nf95_gw_var(grpid, varid, ao%g_sw_philic)
-    ao%g_sw_philic = cshift(ao%g_sw_philic, 1)
-    call nf95_inq_varid(grpid, "single_scat_alb", varid)
-    call nf95_gw_var(grpid, varid, ao%ssa_sw_philic)
-    ao%ssa_sw_philic = cshift(ao%ssa_sw_philic, 1)
-    call nf95_inq_varid(grpid, "mass_ext", varid)
-    call nf95_gw_var(grpid, varid, ao%mass_ext_sw_philic)
-    ao%mass_ext_sw_philic = cshift(ao%mass_ext_sw_philic, 1)
-
-    ! Hydrophilic/Monochromatic:
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic/Monochromatic", grpid)
-    call nf95_inq_varid(grpid, "mass_ext", varid)
-    call nf95_gw_var(grpid, varid, ao%mass_ext_mono_philic)
-
-    ! Hydrophobic/LW_bands:
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophobic/LW_bands", grpid)
-    call nf95_inq_varid(grpid, "asymmetry", varid)
-    call nf95_gw_var(grpid, varid, ao%g_lw_phobic)
-    call nf95_inq_varid(grpid, "single_scat_alb", varid)
-    call nf95_gw_var(grpid, varid, ao%ssa_lw_phobic)
-    call nf95_inq_varid(grpid, "mass_ext", varid)
-    call nf95_gw_var(grpid, varid, ao%mass_ext_lw_phobic)
-
-    ! Hydrophobic/SW_bands:
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophobic/SW_bands", grpid)
-    call nf95_inq_varid(grpid, "asymmetry", varid)
-    call nf95_gw_var(grpid, varid, ao%g_sw_phobic)
-    ao%g_sw_phobic = cshift(ao%g_sw_phobic, 1)
-    call nf95_inq_varid(grpid, "single_scat_alb", varid)
-    call nf95_gw_var(grpid, varid, ao%ssa_sw_phobic)
-    ao%ssa_sw_phobic = cshift(ao%ssa_sw_phobic, 1)
-    call nf95_inq_varid(grpid, "mass_ext", varid)
-    call nf95_gw_var(grpid, varid, ao%mass_ext_sw_phobic)
-    ao%mass_ext_sw_phobic = cshift(ao%mass_ext_sw_phobic, 1)
-
-    ! Hydrophobic/Monochromatic:
-    call nf95_inq_grp_full_ncid(ncid, "Hydrophobic/Monochromatic", grpid)
-    call nf95_inq_varid(grpid, "mass_ext", varid)
-    call nf95_gw_var(grpid, varid, ao%mass_ext_mono_phobic)
-
-    call nf95_close(ncid)
-
-    ! Get array sizes
-    ao%n_bands_lw = size(ao%mass_ext_lw_phobic, 1)
-    ao%n_bands_sw = size(ao%mass_ext_sw_phobic, 1)
-    ao%n_mono_wl = size(ao%mass_ext_mono_phobic, 1)
-    ao%n_type_phobic = size(ao%mass_ext_lw_phobic, 2)
-    ao%n_type_philic = size(ao%mass_ext_lw_philic, 3)
-
-    ! Allocate memory for mapping arrays
-    ao%ntype = ao%n_type_phobic + ao%n_type_philic
-    allocate(ao%iclass(ao%ntype))
-    allocate(ao%itype(ao%ntype))
-
-    ao%iclass = IAerosolClassUndefined
-    ao%itype  = 0
-
-  end subroutine setup_aerosol_optics_lmdz
-
-end module setup_aerosol_optics_lmdz_m
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_aerosol_optics_lmdz_m.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_aerosol_optics_lmdz_m.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_aerosol_optics_lmdz_m.f90	(revision 5268)
@@ -0,0 +1,124 @@
+module setup_aerosol_optics_lmdz_m
+
+  implicit none
+
+contains
+
+  subroutine setup_aerosol_optics_lmdz(ao, file_name)
+
+    ! Read aerosol optical properties. Note differences with
+    ! "radiation_aerosol_optics_data::setup_aerosol_optics":
+
+    ! -- The input NetCDF file is not flat, it contains NetCDF groups.
+
+    ! -- We do not define ao%ssa_mono_phobic, ao%g_mono_phobic,
+    ! ao%lidar_ratio_mono_phobic, ao%ssa_mono_philic,
+    ! ao%g_mono_philic, ao%lidar_ratio_mono_philic. They are not in
+    ! the input NetCDF file and they are not used by ECRad.
+
+    ! -- We do not define ao%description_phobic_str and
+    ! ao%description_philic_str. We just leave the initialization
+    ! value, which is a blank.
+
+    ! -- We have to cshift the shortwave fields because the the
+    ! shortwave bands are in ascending order in the NetCDF file while
+    ! they are not in ECRad.
+
+    use radiation_aerosol_optics_data, only: aerosol_optics_type, &
+         IAerosolClassUndefined
+    use netcdf95, only: nf95_open, nf95_inq_grp_full_ncid, nf95_close, &
+         nf95_inq_dimid, nf95_inq_varid, nf95_inquire_dimension, &
+         nf95_get_var, nf95_gw_var, nf95_nowrite
+
+    type(aerosol_optics_type), intent(out):: ao
+
+    character(len=*), intent(in):: file_name
+    ! NetCDF file containing the aerosol optics data
+
+    ! Local:
+    integer ncid, grpid, dimid, varid
+
+    !-----------------------------------------------------------------------
+
+    ao%use_hydrophilic = .true.
+    ao%use_monochromatic = .true.
+    call nf95_open(file_name, nf95_nowrite, ncid)
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic", grpid)
+    call nf95_inq_dimid(grpid, "hur", dimid)
+    call nf95_inquire_dimension(grpid, dimid, nclen = ao%nrh)
+    allocate(ao%rh_lower(ao%nrh))
+    call nf95_inq_varid(grpid, "hur_bounds", varid)
+    call nf95_get_var(grpid, varid, ao%rh_lower, count_nc = [1, ao%nrh])
+
+    ! Hydrophilic/LW_bands:
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic/LW_bands", grpid)
+    call nf95_inq_varid(grpid, "asymmetry", varid)
+    call nf95_gw_var(grpid, varid, ao%g_lw_philic)
+    call nf95_inq_varid(grpid, "single_scat_alb", varid)
+    call nf95_gw_var(grpid, varid, ao%ssa_lw_philic)
+    call nf95_inq_varid(grpid, "mass_ext", varid)
+    call nf95_gw_var(grpid, varid, ao%mass_ext_lw_philic)
+
+    ! Hydrophilic/SW_bands:
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic/SW_bands", grpid)
+    call nf95_inq_varid(grpid, "asymmetry", varid)
+    call nf95_gw_var(grpid, varid, ao%g_sw_philic)
+    ao%g_sw_philic = cshift(ao%g_sw_philic, 1)
+    call nf95_inq_varid(grpid, "single_scat_alb", varid)
+    call nf95_gw_var(grpid, varid, ao%ssa_sw_philic)
+    ao%ssa_sw_philic = cshift(ao%ssa_sw_philic, 1)
+    call nf95_inq_varid(grpid, "mass_ext", varid)
+    call nf95_gw_var(grpid, varid, ao%mass_ext_sw_philic)
+    ao%mass_ext_sw_philic = cshift(ao%mass_ext_sw_philic, 1)
+
+    ! Hydrophilic/Monochromatic:
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophilic/Monochromatic", grpid)
+    call nf95_inq_varid(grpid, "mass_ext", varid)
+    call nf95_gw_var(grpid, varid, ao%mass_ext_mono_philic)
+
+    ! Hydrophobic/LW_bands:
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophobic/LW_bands", grpid)
+    call nf95_inq_varid(grpid, "asymmetry", varid)
+    call nf95_gw_var(grpid, varid, ao%g_lw_phobic)
+    call nf95_inq_varid(grpid, "single_scat_alb", varid)
+    call nf95_gw_var(grpid, varid, ao%ssa_lw_phobic)
+    call nf95_inq_varid(grpid, "mass_ext", varid)
+    call nf95_gw_var(grpid, varid, ao%mass_ext_lw_phobic)
+
+    ! Hydrophobic/SW_bands:
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophobic/SW_bands", grpid)
+    call nf95_inq_varid(grpid, "asymmetry", varid)
+    call nf95_gw_var(grpid, varid, ao%g_sw_phobic)
+    ao%g_sw_phobic = cshift(ao%g_sw_phobic, 1)
+    call nf95_inq_varid(grpid, "single_scat_alb", varid)
+    call nf95_gw_var(grpid, varid, ao%ssa_sw_phobic)
+    ao%ssa_sw_phobic = cshift(ao%ssa_sw_phobic, 1)
+    call nf95_inq_varid(grpid, "mass_ext", varid)
+    call nf95_gw_var(grpid, varid, ao%mass_ext_sw_phobic)
+    ao%mass_ext_sw_phobic = cshift(ao%mass_ext_sw_phobic, 1)
+
+    ! Hydrophobic/Monochromatic:
+    call nf95_inq_grp_full_ncid(ncid, "Hydrophobic/Monochromatic", grpid)
+    call nf95_inq_varid(grpid, "mass_ext", varid)
+    call nf95_gw_var(grpid, varid, ao%mass_ext_mono_phobic)
+
+    call nf95_close(ncid)
+
+    ! Get array sizes
+    ao%n_bands_lw = size(ao%mass_ext_lw_phobic, 1)
+    ao%n_bands_sw = size(ao%mass_ext_sw_phobic, 1)
+    ao%n_mono_wl = size(ao%mass_ext_mono_phobic, 1)
+    ao%n_type_phobic = size(ao%mass_ext_lw_phobic, 2)
+    ao%n_type_philic = size(ao%mass_ext_lw_philic, 3)
+
+    ! Allocate memory for mapping arrays
+    ao%ntype = ao%n_type_phobic + ao%n_type_philic
+    allocate(ao%iclass(ao%ntype))
+    allocate(ao%itype(ao%ntype))
+
+    ao%iclass = IAerosolClassUndefined
+    ao%itype  = 0
+
+  end subroutine setup_aerosol_optics_lmdz
+
+end module setup_aerosol_optics_lmdz_m
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_config_from_lmdz.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_config_from_lmdz.F90	(revision 4977)
+++ 	(revision )
@@ -1,123 +1,0 @@
-module setup_config_from_lmdz
-
-  use parkind1, only : jprb
-
-  implicit none
-
-  public
-
-  type driver_config_type
-     logical    :: ok_effective_size  = .true.
-     logical    :: ok_separation = .false.   
-     real(jprb) :: high_inv_effective_size   = -1.0_jprb ! m-1
-     real(jprb) :: middle_inv_effective_size = -1.0_jprb ! m-1
-     real(jprb) :: low_inv_effective_size    = -1.0_jprb ! m-1
-     real(jprb) :: cloud_inhom_separation_factor  = -1.0_jprb
-     real(jprb) :: cloud_separation_scale_surface = -1.0_jprb
-     real(jprb) :: cloud_separation_scale_toa     = -1.0_jprb
-     real(jprb) :: cloud_separation_scale_power   = -1.0_jprb
-     real(jprb) :: frac_std = 0.75_jprb
-     real(jprb) :: overlap_decorr_length = 2000.0_jprb
-     ! KDECOLAT : 0 > cste, 1 > Shonk et al. (2010) 2 > 1 lisse a l'equateur
-     !            3 : 1+variation sur la verticale
-     !            4 : 2+variation sur la verticale
-     integer    :: kdecolat = 0
-     real(jprb) :: low_decorrelation_length = 2000.0_jprb 
-     real(jprb) :: mid_decorrelation_length = 2000.0_jprb
-     real(jprb) :: high_decorrelation_length = 2000.0_jprb
-
- contains
- procedure :: read => read_config_from_namelist
-
-  end type driver_config_type
-
-contains
-
-  !---------------------------------------------------------------------
-  subroutine read_config_from_namelist(this, file_name, is_success)
-
-    use radiation_io, only : nulerr, radiation_abort       
-
-    class(driver_config_type), intent(inout) :: this
-    character(*), intent(in)          :: file_name
-    logical, intent(out), optional    :: is_success
-    logical    :: ok_effective_size, ok_separation
-    integer    :: iosopen ! Status after calling open
-    real(jprb) :: high_inv_effective_size 
-    real(jprb) :: middle_inv_effective_size
-    real(jprb) :: low_inv_effective_size
-    real(jprb) :: cloud_inhom_separation_factor 
-    real(jprb) :: cloud_separation_scale_surface
-    real(jprb) :: cloud_separation_scale_toa
-    real(jprb) :: cloud_separation_scale_power
-    real(jprb) :: frac_std
-    real(jprb) :: overlap_decorr_length
-    integer    :: kdecolat 
-    real(jprb) :: low_decorrelation_length
-    real(jprb) :: mid_decorrelation_length 
-    real(jprb) :: high_decorrelation_length 
-
-    namelist /radiation_driver/ ok_effective_size, ok_separation, &
-         &  frac_std, overlap_decorr_length, kdecolat, &
-         &  low_decorrelation_length, mid_decorrelation_length, high_decorrelation_length, &
-         &  high_inv_effective_size, middle_inv_effective_size, low_inv_effective_size, &
-         &  cloud_inhom_separation_factor, cloud_separation_scale_surface, &
-         &  cloud_separation_scale_toa, cloud_separation_scale_power 
-
-    ok_effective_size = .false.
-    ok_separation = .false.
-    high_inv_effective_size   = -1.0_jprb 
-    middle_inv_effective_size = -1.0_jprb 
-    low_inv_effective_size    = -1.0_jprb 
-    cloud_inhom_separation_factor  = -1.0_jprb
-    cloud_separation_scale_surface = -1.0_jprb
-    cloud_separation_scale_toa     = -1.0_jprb
-    cloud_separation_scale_power   = -1.0_jprb
-    frac_std = 0.75_jprb
-    overlap_decorr_length = 2000.0_jprb
-    kdecolat = 0
-    low_decorrelation_length = 2000.0_jprb
-    mid_decorrelation_length = 2000.0_jprb
-    high_decorrelation_length = 2000.0_jprb
-
-    ! Open the namelist file and read the radiation_driver namelist
-    open(unit=10, iostat=iosopen, file=trim(file_name))
-    if (iosopen /= 0) then
-      ! An error occurred
-      if (present(is_success)) then
-        is_success = .false.
-        ! We now continue the subroutine so that the default values
-        ! are placed in the config structure
-      else
-        write(nulerr,'(a,a,a)') '*** Error: namelist file "', &
-             &                trim(file_name), '" not found'
-        call radiation_abort('Driver configuration error')
-      end if
-    else
-      ! Read the radiation_driver namelist, noting that it is not an
-      ! error if this namelist is not present, provided all the required
-      ! variables are present in the NetCDF data file instead
-      read(unit=10, nml=radiation_driver)
-      close(unit=10)
-    end if
-
-    ! Copy namelist data into configuration object
-    this%ok_effective_size = ok_effective_size 
-    this%ok_separation = ok_separation
-    this%frac_std = frac_std
-    this%overlap_decorr_length = overlap_decorr_length
-    this%kdecolat = kdecolat
-    this%low_decorrelation_length = low_decorrelation_length
-    this%mid_decorrelation_length = mid_decorrelation_length
-    this%high_decorrelation_length = high_decorrelation_length
-    this%cloud_inhom_separation_factor = cloud_inhom_separation_factor
-    this%cloud_separation_scale_surface = cloud_separation_scale_surface
-    this%cloud_separation_scale_toa = cloud_separation_scale_toa
-    this%cloud_separation_scale_power = cloud_separation_scale_power
-    this%high_inv_effective_size = high_inv_effective_size
-    this%middle_inv_effective_size = middle_inv_effective_size
-    this%low_inv_effective_size = low_inv_effective_size
-
-  end subroutine read_config_from_namelist          
-
-end module setup_config_from_lmdz
Index: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_config_from_lmdz.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_config_from_lmdz.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_config_from_lmdz.f90	(revision 5268)
@@ -0,0 +1,123 @@
+module setup_config_from_lmdz
+
+  use parkind1, only : jprb
+
+  implicit none
+
+  public
+
+  type driver_config_type
+     logical    :: ok_effective_size  = .true.
+     logical    :: ok_separation = .false.   
+     real(jprb) :: high_inv_effective_size   = -1.0_jprb ! m-1
+     real(jprb) :: middle_inv_effective_size = -1.0_jprb ! m-1
+     real(jprb) :: low_inv_effective_size    = -1.0_jprb ! m-1
+     real(jprb) :: cloud_inhom_separation_factor  = -1.0_jprb
+     real(jprb) :: cloud_separation_scale_surface = -1.0_jprb
+     real(jprb) :: cloud_separation_scale_toa     = -1.0_jprb
+     real(jprb) :: cloud_separation_scale_power   = -1.0_jprb
+     real(jprb) :: frac_std = 0.75_jprb
+     real(jprb) :: overlap_decorr_length = 2000.0_jprb
+     ! KDECOLAT : 0 > cste, 1 > Shonk et al. (2010) 2 > 1 lisse a l'equateur
+     !            3 : 1+variation sur la verticale
+     !            4 : 2+variation sur la verticale
+     integer    :: kdecolat = 0
+     real(jprb) :: low_decorrelation_length = 2000.0_jprb 
+     real(jprb) :: mid_decorrelation_length = 2000.0_jprb
+     real(jprb) :: high_decorrelation_length = 2000.0_jprb
+
+ contains
+ procedure :: read => read_config_from_namelist
+
+  end type driver_config_type
+
+contains
+
+  !---------------------------------------------------------------------
+  subroutine read_config_from_namelist(this, file_name, is_success)
+
+    use radiation_io, only : nulerr, radiation_abort       
+
+    class(driver_config_type), intent(inout) :: this
+    character(*), intent(in)          :: file_name
+    logical, intent(out), optional    :: is_success
+    logical    :: ok_effective_size, ok_separation
+    integer    :: iosopen ! Status after calling open
+    real(jprb) :: high_inv_effective_size 
+    real(jprb) :: middle_inv_effective_size
+    real(jprb) :: low_inv_effective_size
+    real(jprb) :: cloud_inhom_separation_factor 
+    real(jprb) :: cloud_separation_scale_surface
+    real(jprb) :: cloud_separation_scale_toa
+    real(jprb) :: cloud_separation_scale_power
+    real(jprb) :: frac_std
+    real(jprb) :: overlap_decorr_length
+    integer    :: kdecolat 
+    real(jprb) :: low_decorrelation_length
+    real(jprb) :: mid_decorrelation_length 
+    real(jprb) :: high_decorrelation_length 
+
+    namelist /radiation_driver/ ok_effective_size, ok_separation, &
+         &  frac_std, overlap_decorr_length, kdecolat, &
+         &  low_decorrelation_length, mid_decorrelation_length, high_decorrelation_length, &
+         &  high_inv_effective_size, middle_inv_effective_size, low_inv_effective_size, &
+         &  cloud_inhom_separation_factor, cloud_separation_scale_surface, &
+         &  cloud_separation_scale_toa, cloud_separation_scale_power 
+
+    ok_effective_size = .false.
+    ok_separation = .false.
+    high_inv_effective_size   = -1.0_jprb 
+    middle_inv_effective_size = -1.0_jprb 
+    low_inv_effective_size    = -1.0_jprb 
+    cloud_inhom_separation_factor  = -1.0_jprb
+    cloud_separation_scale_surface = -1.0_jprb
+    cloud_separation_scale_toa     = -1.0_jprb
+    cloud_separation_scale_power   = -1.0_jprb
+    frac_std = 0.75_jprb
+    overlap_decorr_length = 2000.0_jprb
+    kdecolat = 0
+    low_decorrelation_length = 2000.0_jprb
+    mid_decorrelation_length = 2000.0_jprb
+    high_decorrelation_length = 2000.0_jprb
+
+    ! Open the namelist file and read the radiation_driver namelist
+    open(unit=10, iostat=iosopen, file=trim(file_name))
+    if (iosopen /= 0) then
+      ! An error occurred
+      if (present(is_success)) then
+        is_success = .false.
+        ! We now continue the subroutine so that the default values
+        ! are placed in the config structure
+      else
+        write(nulerr,'(a,a,a)') '*** Error: namelist file "', &
+             &                trim(file_name), '" not found'
+        call radiation_abort('Driver configuration error')
+      end if
+    else
+      ! Read the radiation_driver namelist, noting that it is not an
+      ! error if this namelist is not present, provided all the required
+      ! variables are present in the NetCDF data file instead
+      read(unit=10, nml=radiation_driver)
+      close(unit=10)
+    end if
+
+    ! Copy namelist data into configuration object
+    this%ok_effective_size = ok_effective_size 
+    this%ok_separation = ok_separation
+    this%frac_std = frac_std
+    this%overlap_decorr_length = overlap_decorr_length
+    this%kdecolat = kdecolat
+    this%low_decorrelation_length = low_decorrelation_length
+    this%mid_decorrelation_length = mid_decorrelation_length
+    this%high_decorrelation_length = high_decorrelation_length
+    this%cloud_inhom_separation_factor = cloud_inhom_separation_factor
+    this%cloud_separation_scale_surface = cloud_separation_scale_surface
+    this%cloud_separation_scale_toa = cloud_separation_scale_toa
+    this%cloud_separation_scale_power = cloud_separation_scale_power
+    this%high_inv_effective_size = high_inv_effective_size
+    this%middle_inv_effective_size = middle_inv_effective_size
+    this%low_inv_effective_size = low_inv_effective_size
+
+  end subroutine read_config_from_namelist          
+
+end module setup_config_from_lmdz
