Index: LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3525)
@@ -113,4 +113,5 @@
   INTEGER :: flag_aerosol
   INTEGER :: flag_aerosol_strat
+  LOGICAL :: flag_aer_feedback
   LOGICAL :: flag_bc_internal_mixture
   LOGICAL :: new_aod
@@ -134,5 +135,7 @@
                    ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan,            &
                    aerosol_couple, chemistry_couple, flag_aerosol,        &
-                   flag_aerosol_strat, new_aod, flag_bc_internal_mixture, &
+                   flag_aerosol_strat,                                    & 
+                   flag_aer_feedback,                                     &
+                   new_aod, flag_bc_internal_mixture,                     &
                    bl95_b0, bl95_b1, read_climoz, alp_offset)
   CALL phys_state_var_init(read_climoz)
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/micphy_tstep.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/micphy_tstep.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/micphy_tstep.F90	(revision 3525)
@@ -1,4 +1,8 @@
+!
+! $Id$
+!
 SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
 
+  USE geometry_mod, ONLY : latitude_deg !NL- latitude corr. to local domain
   USE dimphy, ONLY : klon,klev
   USE aerophys
@@ -9,5 +13,7 @@
   USE sulfate_aer_mod, ONLY : STRAACT
   USE YOMCST, ONLY : RPI, RD, RG
-
+  USE print_control_mod, ONLY: lunout
+  USE strataer_mod
+  
   IMPLICIT NONE
 
@@ -89,5 +95,12 @@
       ! compute nucleation rate in kg(H2SO4)/kgA/s
       CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), &
-             & a_xm,b_xm,c_xm,nucl_rate,ntot,x)
+           & a_xm,b_xm,c_xm,nucl_rate,ntot,x)
+      !NL - add nucleation box (if flag on)
+      IF (flag_nuc_rate_box) THEN
+         IF (latitude_deg(ilon).LE.(nuclat_min) .OR. latitude_deg(ilon).GE.(nuclat_max) &
+              .OR. pplay(ilon,ilev).GE.nucpres_max .AND. pplay(ilon,ilev) .LE. nucpres_min ) THEN
+            nucl_rate=0.0
+         ENDIF
+      ENDIF
       ! compute cond/evap rate in kg(H2SO4)/kgA/s
       CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
@@ -160,5 +173,5 @@
     DO it=1, nbtr
       IF (tr_seri(ilon,ilev,it).LT.0.0) THEN
-        PRINT *, 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it
+         WRITE(lunout,*) 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it
       ENDIF
     ENDDO
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/strataer_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/strataer_mod.F90	(revision 3525)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/strataer_mod.F90	(revision 3525)
@@ -0,0 +1,233 @@
+! $Id$
+MODULE strataer_mod
+! This module contains information about strato microphysic model parameters
+  
+  IMPLICIT NONE
+
+  ! flag to constraint nucleation rate in a lat/pres box
+  LOGICAL,SAVE :: flag_nuc_rate_box      ! Nucleation rate limit or not to a lat/pres
+  !$OMP THREADPRIVATE(flag_nuc_rate_box)
+  REAL,SAVE    :: nuclat_min             ! min lat to activate nuc rate
+  REAL,SAVE    :: nuclat_max             ! max lat to activate nuc rate
+  REAL,SAVE    :: nucpres_min            ! min pres to activate nuc rate
+  REAL,SAVE    :: nucpres_max            ! max pres to activate nuc rate
+  !$OMP THREADPRIVATE(nuclat_min, nuclat_max, nucpres_min, nucpres_max)
+  
+  ! flag for sulfur emission scenario: (0) background aer ; (1) volcanic eruption ; (2) strato aer injections (SAI)
+  INTEGER,SAVE :: flag_sulf_emit
+  !$OMP THREADPRIVATE(flag_sulf_emit)
+
+  ! flag for sulfur emission altitude distribution: (0) gaussian; (1) uniform
+  INTEGER,SAVE :: flag_sulf_emit_distrib
+  !$OMP THREADPRIVATE(flag_sulf_emit_distrib)
+  
+  !--flag_sulf_emit=1 -- Volcanic eruption(s)
+  INTEGER,SAVE             :: nErupt                    ! number of eruptions specs
+  REAL,SAVE                :: injdur                    ! volcanic injection duration
+  !$OMP THREADPRIVATE(nErupt, injdur)
+  INTEGER,ALLOCATABLE,SAVE :: year_emit_vol(:)          ! year of emission date
+  INTEGER,ALLOCATABLE,SAVE :: mth_emit_vol(:)           ! month of emission date
+  INTEGER,ALLOCATABLE,SAVE :: day_emit_vol(:)           ! day of emission date
+  !$OMP THREADPRIVATE(year_emit_vol, mth_emit_vol, day_emit_vol)
+  REAL,ALLOCATABLE,SAVE    :: m_aer_emiss_vol(:)        ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
+  REAL,ALLOCATABLE,SAVE    :: altemiss_vol(:)           ! emission altitude in m
+  REAL,ALLOCATABLE,SAVE    :: sigma_alt_vol(:)          ! standard deviation of emission altitude in m
+  !$OMP THREADPRIVATE(m_aer_emiss_vol, altemiss_vol, sigma_alt_vol)
+  INTEGER,ALLOCATABLE,SAVE :: ponde_lonlat_vol(:)       ! lon/lat ponderation factor
+  REAL,ALLOCATABLE,SAVE    :: xlat_min_vol(:)           ! min latitude of volcano in degree
+  REAL,ALLOCATABLE,SAVE    :: xlat_max_vol(:)           ! max latitude of volcano in degree
+  REAL,ALLOCATABLE,SAVE    :: xlon_min_vol(:)           ! min longitude of volcano in degree
+  REAL,ALLOCATABLE,SAVE    :: xlon_max_vol(:)           ! max longitude of volcano in degree
+  !$OMP THREADPRIVATE(ponde_lonlat_vol, xlat_min_vol, xlat_max_vol, xlon_min_vol, xlon_max_vol)
+  
+  !--flag_sulf_emit=2 --SAI
+  REAL,SAVE    :: m_aer_emiss_sai        ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
+  REAL,SAVE    :: altemiss_sai           ! emission altitude in m
+  REAL,SAVE    :: sigma_alt_sai          ! standard deviation of emission altitude in m
+  !$OMP THREADPRIVATE(m_aer_emiss_sai, altemiss_sai, sigma_alt_sai)
+  REAL,SAVE    :: xlat_sai               ! latitude of SAI in degree
+  REAL,SAVE    :: xlon_sai               ! longitude of SAI in degree
+  REAL,SAVE    :: dlat, dlon             ! delta latitude and d longitude of grid in degree
+  !$OMP THREADPRIVATE(xlat_sai, xlon_sai, dlat, dlon)
+  
+  !--flag_sulf_emit=3 -- SAI
+  REAL,SAVE    :: xlat_max_sai           ! maximum latitude of SAI in degrees
+  REAL,SAVE    :: xlat_min_sai           ! minimum latitude of SAI in degrees
+  !$OMP THREADPRIVATE(xlat_min_sai,xlat_max_sai)
+  
+CONTAINS
+    
+  SUBROUTINE strataer_init()
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE print_control_mod, ONLY : lunout
+    USE mod_phys_lmdz_para, ONLY : is_master
+
+    ! Local var
+    INTEGER       :: ieru
+ 
+    INTEGER :: i
+    
+    WRITE(lunout,*) 'IN STRATAER INIT WELCOME!'
+   
+    !Config Key  = flag_sulf_emit
+    !Config Desc = aerosol emission mode
+    ! - 0 = background aerosol
+    ! - 1 = volcanic eruption
+    ! - 2 = geo-ingeneering design
+    ! - 3 = geo-engineering between two latitudes
+    !Config Def  = 0
+    !Config Help = Used in physiq.F
+    !
+    flag_sulf_emit = 0
+    nErupt = 0 ! eruption number
+    injdur = 0 ! init injection duration
+    CALL getin_p('flag_sulf_emit',flag_sulf_emit)
+
+    IF (flag_sulf_emit==1) THEN ! Volcano
+       CALL getin_p('nErupt',nErupt)
+       CALL getin_p('injdur',injdur)
+    ELSEIF (flag_sulf_emit == 2) THEN ! SAI
+       CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
+       CALL getin_p('altemiss_sai',altemiss_sai)
+       CALL getin_p('sigma_alt_sai',sigma_alt_sai)
+       CALL getin_p('xlat_sai',xlat_sai)
+       CALL getin_p('xlon_sai',xlon_sai)
+       CALL getin_p('flag_sulf_emit_distrib',flag_sulf_emit_distrib)
+    ELSEIF (flag_sulf_emit == 3) THEN ! SAI between latitudes
+       CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
+       CALL getin_p('altemiss_sai',altemiss_sai)
+       CALL getin_p('sigma_alt_sai',sigma_alt_sai)
+       CALL getin_p('xlon_sai',xlon_sai)
+       CALL getin_p('xlat_max_sai',xlat_max_sai)
+       CALL getin_p('xlat_min_sai',xlat_min_sai)
+       CALL getin_p('flag_sulf_emit_distrib',flag_sulf_emit_distrib)
+    ENDIF
+
+    ALLOCATE(year_emit_vol(nErupt),mth_emit_vol(nErupt),day_emit_vol(nErupt))
+    ALLOCATE(m_aer_emiss_vol(nErupt),altemiss_vol(nErupt),sigma_alt_vol(nErupt))
+    ALLOCATE(xlat_min_vol(nErupt),xlon_min_vol(nErupt))
+    ALLOCATE(xlat_max_vol(nErupt),xlon_max_vol(nErupt))
+    
+    year_emit_vol=0 ; mth_emit_vol=0 ; day_emit_vol=0
+    m_aer_emiss_vol=0. ; altemiss_vol=0. ; sigma_alt_vol=0.
+    xlon_min_vol=0. ; xlon_max_vol=0.
+    xlat_min_vol=0. ; xlat_max_vol=0.
+    
+    CALL getin_p('year_emit_vol',year_emit_vol)
+    CALL getin_p('mth_emit_vol',mth_emit_vol)
+    CALL getin_p('day_emit_vol',day_emit_vol)
+    CALL getin_p('m_aer_emiss_vol',m_aer_emiss_vol)
+    CALL getin_p('altemiss_vol',altemiss_vol)
+    CALL getin_p('sigma_alt_vol',sigma_alt_vol)
+    CALL getin_p('xlon_min_vol',xlon_min_vol)
+    CALL getin_p('xlon_max_vol',xlon_max_vol)
+    CALL getin_p('xlat_min_vol',xlat_min_vol)
+    CALL getin_p('xlat_max_vol',xlat_max_vol)
+    !Config Key  = flag_nuc_rate_box
+    !Config Desc = define or not a box for nucleation rate
+    ! - F = global nucleation
+    ! - T = 2D-box for nucleation need nuclat_min, nuclat_max, nucpres_min and
+    ! nucpres_max
+    !       to define its bounds.
+    !Config Def  = F
+    !Config Help = Used in physiq.F
+    !
+    flag_nuc_rate_box = .FALSE.
+    CALL getin_p('flag_nuc_rate_box',flag_nuc_rate_box)
+    CALL getin_p('nuclat_min',nuclat_min)
+    CALL getin_p('nuclat_max',nuclat_max)
+    CALL getin_p('nucpres_min',nucpres_min)
+    CALL getin_p('nucpres_max',nucpres_max)
+
+    WRITE(lunout,*) 'IN STRATAER INIT2 year_emit_vol = ',year_emit_vol
+    WRITE(lunout,*) 'IN STRATAER INIT2 mth_emit_vol = ',mth_emit_vol
+    WRITE(lunout,*) 'IN STRATAER INIT2 day_emit_vol = ',day_emit_vol
+    
+    !IF (is_master) THEN
+       WRITE(lunout,*) 'IN STRATAER INIT2 year_emit_vol = ',year_emit_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 mth_emit_vol=',mth_emit_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 day_emit_vol=',day_emit_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 =m_aer_emiss_vol',m_aer_emiss_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 =altemiss_vol',altemiss_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 =sigma_alt_vol',sigma_alt_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 xlon_min_vol=',xlon_min_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 xlon_max_vol=',xlon_max_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 xlat_min_vol=',xlat_min_vol
+       WRITE(lunout,*) 'IN STRATAER INIT2 xlat_max_vol=',xlat_max_vol
+       WRITE(lunout,*) 'flag_nuc_rate_box = ',flag_nuc_rate_box
+       WRITE(lunout,*) 'nuclat_min = ',nuclat_min
+       WRITE(lunout,*) 'nuclat_max = ',nuclat_max
+       WRITE(lunout,*) 'nucpres_min = ',nucpres_min
+       WRITE(lunout,*) 'nucpres_max = ',nucpres_max
+       WRITE(lunout,*) 'flag_sulf_emit = ',flag_sulf_emit
+       WRITE(lunout,*) 'injdur = ',injdur
+       WRITE(lunout,*) 'flag_sulf_emit_distrib = ',flag_sulf_emit_distrib
+       WRITE(lunout,*) 'nErupt = ',nErupt
+       WRITE(lunout,*) 'year_emit_vol = ',year_emit_vol
+       WRITE(lunout,*) 'mth_emit_vol = ',mth_emit_vol
+       WRITE(lunout,*) 'day_emit_vol = ',day_emit_vol
+       WRITE(lunout,*) 'm_aer_emiss_vol = ',m_aer_emiss_vol
+       WRITE(lunout,*) 'altemiss_vol = ',altemiss_vol
+       WRITE(lunout,*) 'sigma_alt_vol = ',sigma_alt_vol
+       WRITE(lunout,*) 'xlat_min_vol = ',xlat_min_vol
+       WRITE(lunout,*) 'xlat_max_vol = ',xlat_max_vol
+       WRITE(lunout,*) 'xlon_min_vol = ',xlon_min_vol
+       WRITE(lunout,*) 'xlon_max_vol = ',xlon_max_vol
+       WRITE(lunout,*) 'm_aer_emiss_sai = ',m_aer_emiss_sai
+       WRITE(lunout,*) 'altemiss_sai = ',altemiss_sai
+       WRITE(lunout,*) 'sigma_alt_sai = ',sigma_alt_sai
+       WRITE(lunout,*) 'xlat_sai = ',xlat_sai
+       WRITE(lunout,*) 'xlon_sai = ',xlon_sai
+       WRITE(lunout,*) 'xlat_min_sai = ',xlat_min_sai
+       WRITE(lunout,*) 'xlat_max_sai = ',xlat_max_sai
+    !ENDIF
+
+    CALL strataer_ponde_init
+    WRITE(lunout,*) 'IN STRATAER INT2 END'
+
+  END SUBROUTINE strataer_init
+  
+  ! Compute the ponderation to applicate in each grid point for all eruptions and init
+  ! dlat & dlon variables
+  SUBROUTINE strataer_ponde_init()
+    
+    USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+    USE dimphy, ONLY: klon
+    USE mod_grid_phy_lmdz, ONLY: nbp_lat, nbp_lon
+    USE print_control_mod, ONLY : lunout
+    USE YOMCST, ONLY : RPI
+
+    ! local var
+    REAL                :: pi,lat_reg_deg,lon_reg_deg! latitude and d longitude of grid in degree
+    INTEGER             :: ieru, i, j
+    
+    ALLOCATE(ponde_lonlat_vol(nErupt))
+    
+    !Compute lon/lat ponderation for injection
+    dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
+    dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
+    WRITE(lunout,*) 'IN STRATAER_INIT dlat=',dlat,'dlon=',dlon
+    WRITE(lunout,*) 'IN STRATAER_INIT nErupt=',nErupt
+    WRITE(lunout,*) 'IN STRATAER_INIT xlat_min=',xlat_min_vol,'xlat_max=',xlat_max_vol
+    WRITE(lunout,*) 'IN STRATAER_INIT xlon_min=',xlon_min_vol,'xlon_max=',xlon_max_vol
+    DO ieru=1, nErupt
+       ponde_lonlat_vol(ieru) = 0
+       DO i=1,nbp_lon
+          lon_reg_deg = lon_reg(i)*180./RPI
+          DO j=1,nbp_lat
+             lat_reg_deg = lat_reg(j)*180./RPI
+             IF  ( lat_reg_deg.GE.xlat_min_vol(ieru)-dlat .AND. lat_reg_deg.LT.xlat_max_vol(ieru)+dlat .AND. &
+                  lon_reg_deg.GE.xlon_min_vol(ieru)-dlon .AND. lon_reg_deg.LT.xlon_max_vol(ieru)+dlon ) THEN
+                ponde_lonlat_vol(ieru) = ponde_lonlat_vol(ieru) + 1
+             ENDIF
+          ENDDO
+       ENDDO
+       IF(ponde_lonlat_vol(ieru) == 0) THEN
+          WRITE(lunout,*) 'STRATAER_INIT ERROR: no grid point found for eruption ieru=',ieru
+       ENDIF
+    ENDDO !ieru
+    WRITE(lunout,*) 'IN STRATAER_INIT ponde_lonlat: ',ponde_lonlat_vol
+    
+  END SUBROUTINE strataer_ponde_init
+  
+END MODULE strataer_mod
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/traccoag_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/traccoag_mod.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/traccoag_mod.F90	(revision 3525)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 MODULE traccoag_mod
 !
@@ -16,5 +19,5 @@
     USE infotrac
     USE aerophys
-    USE geometry_mod, ONLY : cell_area
+    USE geometry_mod, ONLY : cell_area, boundslat
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
@@ -24,4 +27,7 @@
     USE phys_local_var_mod, ONLY: stratomask
     USE YOMCST
+    USE print_control_mod, ONLY: lunout
+    USE strataer_mod
+    USE phys_cal_mod, ONLY : year_len
 
     IMPLICIT NONE
@@ -52,26 +58,7 @@
 ! Local variables
 !----------------
-! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI)
-    INTEGER,PARAMETER  :: flag_sulf_emit=2 
-!
-!--flag_sulf_emit=1 --example Pinatubo
-    INTEGER,PARAMETER :: year_emit_vol=1991          ! year of emission date
-    INTEGER,PARAMETER :: mth_emit_vol=6              ! month of emission date
-    INTEGER,PARAMETER :: day_emit_vol=15             ! day of emission date 
-    REAL,PARAMETER    :: m_aer_emiss_vol=7.e9   ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
-    REAL,PARAMETER    :: altemiss_vol=17.e3     ! emission altitude in m
-    REAL,PARAMETER    :: sigma_alt_vol=1.e3     ! standard deviation of emission altitude in m
-    REAL,PARAMETER    :: xlat_vol=15.14         ! latitude of volcano in degree
-    REAL,PARAMETER    :: xlon_vol=120.35        ! longitude of volcano in degree
-
-!--flag_sulf_emit=2 --SAI
-    REAL,PARAMETER    :: m_aer_emiss_sai=1.e10  ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
-    REAL,PARAMETER    :: altemiss_sai=17.e3     ! emission altitude in m
-    REAL,PARAMETER    :: sigma_alt_sai=1.e3     ! standard deviation of emission altitude in m
-    REAL,PARAMETER    :: xlat_sai=0.01          ! latitude of SAI in degree
-    REAL,PARAMETER    :: xlon_sai=120.35        ! longitude of SAI in degree
-
-!--other local variables
-    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int
+    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
+    REAL                                   :: sum_emi_so2         ! Test sum of all LON for budg_emi_so2
+    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
     LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
     REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
@@ -90,14 +77,46 @@
     REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
     REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
-    REAL                                   :: dlat, dlon          ! d latitude and d longitude of grid in degree
     REAL                                   :: emission            ! emission
+    REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
+    REAL                                   :: dlat_loc
 
     IF (is_mpi_root) THEN
-      PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
+       WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
+       WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit
+       IF (flag_sulf_emit == 1) THEN
+          WRITE(lunout,*) 'IN traccoag nErupt: ',nErupt
+          WRITE(lunout,*) 'IN traccoag injdur: ',injdur
+          WRITE(lunout,*) 'IN traccoag : year_emit_vol',year_emit_vol
+          WRITE(lunout,*) 'IN traccoag : mth_emit_vol',mth_emit_vol
+          WRITE(lunout,*) 'IN traccoag : day_emit_vol',day_emit_vol
+          WRITE(lunout,*) 'IN traccoag : m_aer_emiss_vol',m_aer_emiss_vol
+          WRITE(lunout,*) 'IN traccoag : altemiss_vol',altemiss_vol
+          WRITE(lunout,*) 'IN traccoag : sigma_alt_vol',sigma_alt_vol
+          WRITE(lunout,*) 'IN traccoag : ponde_lonlat_vol',ponde_lonlat_vol
+          WRITE(lunout,*) 'IN traccoag : xlat_min_vol',xlat_min_vol
+          WRITE(lunout,*) 'IN traccoag : xlat_max_vol',xlat_max_vol
+          WRITE(lunout,*) 'IN traccoag : xlon_min_vol',xlon_min_vol
+          WRITE(lunout,*) 'IN traccoag : xlon_max_vol',xlon_max_vol
+       ELSEIF (flag_sulf_emit == 2) THEN
+          WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai
+          WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai
+          WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai
+          WRITE(lunout,*) 'IN traccoag : xlat_sai',xlat_sai
+          WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai
+       ELSEIF (flag_sulf_emit == 3) THEN
+          WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai
+          WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai
+          WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai
+          WRITE(lunout,*) 'IN traccoag : xlat_min_sai',xlat_min_sai
+          WRITE(lunout,*) 'IN traccoag : xlat_max_sai',xlat_max_sai
+          WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai
+       ENDIF
+       WRITE(lunout,*) 'IN traccoag : flag_nuc_rate_box = ',flag_nuc_rate_box
+       IF (flag_nuc_rate_box) THEN
+          WRITE(lunout,*) 'IN traccoag : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max
+          WRITE(lunout,*) 'IN traccoag : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max
+       ENDIF
     ENDIF
-
-    dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
-    dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
-
+    
     DO it=1, nbtr_bin
       r_bin(it)=mdw(it)/2.
@@ -120,5 +139,5 @@
     IF (debutphy .and. is_mpi_root) THEN
       DO it=1, nbtr_bin
-        PRINT *,'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
+        WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
       ENDDO
     ENDIF
@@ -170,46 +189,82 @@
       !--only emit on day of eruption
       ! stretch emission over one day of Pinatubo eruption
-      IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN 
-!
-        DO i=1,klon
-          !Pinatubo eruption at 15.14N, 120.35E
-          IF  ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
-                xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
-! 
-          PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
-!         compute altLMDz
-            altLMDz(:)=0.0
-            DO k=1, klev
-              zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
-              zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
-              zdz=zdm(k)/zrho                      !thickness of layer in m
-              altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
-            ENDDO
-            !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
-            f_lay_sum=0.0
-            DO k=1, klev
-              f_lay_emiss(k)=0.0
-              DO i_int=1, n_int_alt
-                alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
-                f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* &
-                &              exp(-0.5*((alt-altemiss_vol)/sigma_alt_vol)**2.)*   & 
-                &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 
-              ENDDO
-              f_lay_sum=f_lay_sum+f_lay_emiss(k)
-            ENDDO
-            !correct for step integration error
-            f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
-            !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
-            !vertically distributed emission
-            DO k=1, klev
-              ! stretch emission over one day (minus one timestep) of Pinatubo eruption
-              emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
-              tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 
-              budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
-            ENDDO
-          ENDIF ! emission grid cell
-        ENDDO ! klon loop
-      ENDIF ! emission period
-
+       DO ieru=1, nErupt
+          IF (is_mpi_root) THEN
+             sum_emi_so2 = 0.0 ! Init sum
+          ENDIF
+          IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
+               day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
+             !
+             ! daily injection mass emission - NL
+             m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
+             WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
+                  'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
+                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
+             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
+             DO i=1,klon
+                !Pinatubo eruption at 15.14N, 120.35E
+                dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
+                WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
+                IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
+                     xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
+                   !
+                   WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
+                   WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
+                   !         compute altLMDz
+                   altLMDz(:)=0.0
+                   DO k=1, klev
+                      zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
+                      zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
+                      zdz=zdm(k)/zrho                      !thickness of layer in m
+                      altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
+                   ENDDO
+
+                   SELECT CASE(flag_sulf_emit_distrib)
+                   
+                   CASE(0) ! Gaussian distribution
+                   !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
+                   f_lay_sum=0.0
+                   DO k=1, klev
+                      f_lay_emiss(k)=0.0
+                      DO i_int=1, n_int_alt
+                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* &
+                              &              exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)*   &
+                              &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                      ENDDO
+                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
+                   ENDDO
+                   
+                   CASE(1) ! Uniform distribution
+                   ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
+                   ! height of the injection, centered around altemiss_vol(ieru)
+                   DO k=1, klev
+                      f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
+                      & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
+                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
+                   ENDDO
+
+                   END SELECT        ! End CASE over flag_sulf_emit_distrib)
+
+                   WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
+                   WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
+                   !correct for step integration error
+                   f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
+                   !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
+                   !vertically distributed emission
+                   DO k=1, klev
+                      ! stretch emission over one day of Pinatubo eruption
+                      emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
+                      tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
+                      budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
+                   ENDDO
+                   sum_emi_so2 = sum_emi_so2 + budg_emi_so2(i) ! Sum all LON
+                ENDIF ! emission grid cell
+             ENDDO ! klon loop
+             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") global sum_emi_so2=",sum_emi_so2
+             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
+          ENDIF ! emission period
+       ENDDO ! eruption number
+       
     CASE(2) ! stratospheric aerosol injections (SAI)
 !
@@ -217,10 +272,12 @@
 !       SAI standard scenario with continuous emission from 1 grid point at the equator
 !       SAI emission on single month
-!       IF  ((mth_cur==4 .AND. &
 !       SAI continuous emission o
-        IF  ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
+        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
+        WRITE(lunout,*) "IN traccoag, dlon=",dlon
+        WRITE(lunout,*) "IN traccoag, dlat=",dlat_loc
+        IF  ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
           &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
 !
-          PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
+          WRITE(lunout,*) 'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
 !         compute altLMDz
           altLMDz(:)=0.0
@@ -231,16 +288,33 @@
             altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
           ENDDO
+
+          SELECT CASE(flag_sulf_emit_distrib)
+
+          CASE(0) ! Gaussian distribution
           !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
           f_lay_sum=0.0
-          DO k=1, klev
-            f_lay_emiss(k)=0.0
-            DO i_int=1, n_int_alt
-              alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
-              f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
-              &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   & 
-              &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 
-            ENDDO
-            f_lay_sum=f_lay_sum+f_lay_emiss(k)
-          ENDDO
+               DO k=1, klev
+                     f_lay_emiss(k)=0.0
+                     DO i_int=1, n_int_alt
+                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
+                         &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   & 
+                         &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 
+                     ENDDO
+                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
+               ENDDO
+
+          CASE(1) ! Uniform distribution
+          f_lay_sum=0.0
+          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
+          ! the height of the injection, centered around altemiss_sai
+               DO k=1, klev
+                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
+                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
+                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
+               ENDDO
+
+          END SELECT ! Gaussian or uniform distribution
+
           !correct for step integration error
           f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
@@ -249,13 +323,89 @@
           DO k=1, klev
             ! stretch emission over whole year (360d)
-            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400.  
+            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400.  
             tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
             budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
           ENDDO
+
 !          !emission as monodisperse particles with 0.1um dry radius (BIN21)
 !          !vertically distributed emission
 !          DO k=1, klev
 !            ! stretch emission over whole year (360d)
-!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400 
+!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 
+!            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 
+!            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
+!          ENDDO
+        ENDIF ! emission grid cell
+      ENDDO ! klon loop
+
+    CASE(3) ! --- SAI injection over a single band of longitude and between
+            !     lat_min and lat_max
+
+    WRITE(lunout,*) 'IN traccoag, dlon=',dlon
+    DO i=1,klon
+!       SAI scenario with continuous emission
+        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
+        WRITE(lunout,*) 'IN traccoag, dlat = ',dlat_loc
+        theta_min = max(xlat(i)-dlat_loc,xlat_min_sai)
+        theta_max = min(xlat(i)+dlat_loc,xlat_max_sai)
+        IF  ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. &
+          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
+!
+!         compute altLMDz
+          altLMDz(:)=0.0
+          DO k=1, klev
+            zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
+            zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
+            zdz=zdm(k)/zrho                      !thickness of layer in m
+            altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
+          ENDDO
+
+          SELECT CASE(flag_sulf_emit_distrib)
+
+          CASE(0) ! Gaussian distribution
+          !compute distribution of emission to vertical model layers (based on
+          !Gaussian peak in altitude)
+          f_lay_sum=0.0
+               DO k=1, klev
+                     f_lay_emiss(k)=0.0
+                     DO i_int=1, n_int_alt
+                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
+                         & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
+                         & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                     ENDDO
+                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
+               ENDDO
+
+          CASE(1) ! Uniform distribution
+          f_lay_sum=0.0
+          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
+          ! the height of the injection, centered around altemiss_sai
+               DO k=1, klev
+                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
+                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
+                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
+               ENDDO
+
+          END SELECT ! Gaussian or uniform distribution
+
+          !correct for step integration error
+          f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
+          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
+          !vertically distributed emission
+          DO k=1, klev
+            ! stretch emission over whole year (360d)
+            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ &
+                      & year_len/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & 
+                      & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
+            tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
+            budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
+          ENDDO
+
+!          !emission as monodisperse particles with 0.1um dry radius (BIN21)
+!          !vertically distributed emission
+!          DO k=1, klev
+!            ! stretch emission over whole year (360d)
+!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 
 !            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 
 !            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
@@ -291,5 +441,6 @@
         IF (mdw(it) .LT. 2.5e-6) THEN
           !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
-          !assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity
+          !assume that particles consist of ammonium sulfate at the surface (132g/mol) 
+          !and are dry at T = 20 deg. C and 50 perc. humidity
           surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
                            & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90	(revision 3525)
@@ -18,5 +18,6 @@
        iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
        ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
-       chemistry_couple, flag_aerosol, flag_aerosol_strat, new_aod, &
+       chemistry_couple, flag_aerosol, flag_aerosol_strat,         &
+       flag_aer_feedback, new_aod, &
        flag_bc_internal_mixture, bl95_b0, bl95_b1,&
        read_climoz, &
@@ -30,18 +31,16 @@
     USE print_control_mod, ONLY: lunout
 
-    include "conema3.h"
-    include "fisrtilp.h"
-    include "nuage.h"
-    include "YOMCST.h"
-    include "YOMCST2.h"
-
-    include "thermcell.h"
-
+    INCLUDE "conema3.h"
+    INCLUDE "fisrtilp.h"
+    INCLUDE "nuage.h"
+    INCLUDE "YOMCST.h"
+    INCLUDE "YOMCST2.h"
+    INCLUDE "thermcell.h"
 
     !IM : on inclut/initialise les taux de CH4, N2O, CFC11 et CFC12
-    include "clesphys.h"
-    include "compbl.h"
-    include "comsoil.h"
-    include "YOEGWD.h"
+    INCLUDE "clesphys.h"
+    INCLUDE "compbl.h"
+    INCLUDE "comsoil.h"
+    INCLUDE "YOEGWD.h"
     !
     ! Configuration de la "physique" de LMDZ a l'aide de la fonction
@@ -49,6 +48,4 @@
     !
     ! LF 05/2001
-    !
-
     !
     ! type_ocean:      type d'ocean (force, slab, couple)
@@ -80,4 +77,5 @@
     INTEGER              :: flag_aerosol
     INTEGER              :: flag_aerosol_strat
+    LOGICAL              :: flag_aer_feedback
     LOGICAL              :: flag_bc_internal_mixture
     LOGICAL              :: new_aod
@@ -101,4 +99,5 @@
     INTEGER, SAVE       :: flag_aerosol_omp
     INTEGER, SAVE       :: flag_aerosol_strat_omp
+    LOGICAL, SAVE       :: flag_aer_feedback_omp
     LOGICAL, SAVE       :: flag_bc_internal_mixture_omp
     LOGICAL, SAVE       :: new_aod_omp
@@ -1051,5 +1050,16 @@
     ENDIF
 
-    !
+    !Config Key  = flag_aer_feedback
+    !Config Desc = (des)activate aerosol radiative feedback
+    ! - F = no aerosol radiative feedback
+    ! - T = aerosol radiative feedback
+    !Config Def  = T
+    !Config Help = Used in physiq.F
+    !
+    flag_aer_feedback_omp = .TRUE. 
+    IF (iflag_rrtm_omp==1) THEN
+       CALL getin('flag_aer_feedback',flag_aer_feedback_omp)
+    ENDIF
+
     !Config Key  = iflag_cld_th 
     !Config Desc =  
@@ -2298,4 +2308,5 @@
     flag_aerosol=flag_aerosol_omp
     flag_aerosol_strat=flag_aerosol_strat_omp
+    flag_aer_feedback=flag_aer_feedback_omp
     flag_bc_internal_mixture=flag_bc_internal_mixture_omp
     new_aod=new_aod_omp
@@ -2522,206 +2533,207 @@
     !$OMP MASTER
 
-    write(lunout,*)' ##############################################'
-    write(lunout,*)' Configuration des parametres de la physique: '
-    write(lunout,*)' Type ocean = ', type_ocean
-    write(lunout,*)' Version ocean = ', version_ocean
-    write(lunout,*)' Config veget = ', ok_veget,type_veget
-    write(lunout,*)' Snow model SISVAT : ok_snow = ', ok_snow
-    write(lunout,*)' Config xml pour XIOS : ok_all_xml = ', ok_all_xml
-    write(lunout,*)' Sortie journaliere = ', ok_journe
-    write(lunout,*)' Sortie haute frequence = ', ok_hf
-    write(lunout,*)' Sortie mensuelle = ', ok_mensuel
-    write(lunout,*)' Sortie instantanee = ', ok_instan
-    write(lunout,*)' Frequence appel simulateur ISCCP, freq_ISCCP =', freq_ISCCP
-    write(lunout,*)' Frequence appel simulateur ISCCP, ecrit_ISCCP =', ecrit_ISCCP
-    write(lunout,*)' Frequence appel simulateur COSP, freq_COSP =', freq_COSP
-    write(lunout,*)' Frequence appel simulateur AIRS, freq_AIRS =', freq_AIRS
-    write(lunout,*)' Sortie bilan d''energie, ip_ebil_phy =', ip_ebil_phy
-    write(lunout,*)' Excentricite = ',R_ecc
-    write(lunout,*)' Equinoxe = ',R_peri
-    write(lunout,*)' Inclinaison =',R_incl
-    write(lunout,*)' Constante solaire =',solaire
-    write(lunout,*)' ok_suntime_rrtm =',ok_suntime_rrtm
-    write(lunout,*)' co2_ppm =',co2_ppm
-    write(lunout,*)' RCO2_act = ',RCO2_act
-    write(lunout,*)' CH4_ppb =',CH4_ppb,' RCH4_act = ',RCH4_act
-    write(lunout,*)' N2O_ppb =',N2O_ppb,' RN2O_act=  ',RN2O_act
-    write(lunout,*)' CFC11_ppt=',CFC11_ppt,' RCFC11_act=  ',RCFC11_act
-    write(lunout,*)' CFC12_ppt=',CFC12_ppt,' RCFC12_act=  ',RCFC12_act
-    write(lunout,*)' RCO2_per = ',RCO2_per,' RCH4_per = ', RCH4_per
-    write(lunout,*)' RN2O_per = ',RN2O_per,' RCFC11_per = ', RCFC11_per
-    write(lunout,*)' RCFC12_per = ',RCFC12_per
-    write(lunout,*)' cvl_comp_threshold=', cvl_comp_threshold
-    write(lunout,*)' cvl_sig2feed=', cvl_sig2feed
-    write(lunout,*)' cvl_corr=', cvl_corr
-    write(lunout,*)'ok_lic_melt=', ok_lic_melt
-    write(lunout,*)'ok_lic_cond=', ok_lic_cond
-    write(lunout,*)'iflag_cycle_diurne=',iflag_cycle_diurne
-    write(lunout,*)'soil_model=',soil_model
-    write(lunout,*)'new_oliq=',new_oliq
-    write(lunout,*)'ok_orodr=',ok_orodr
-    write(lunout,*)'ok_orolf=',ok_orolf
-    write(lunout,*)'ok_limitvrai=',ok_limitvrai
-    write(lunout,*)'nbapp_rad=',nbapp_rad
-    write(lunout,*)'iflag_con=',iflag_con
-    write(lunout,*)'nbapp_cv=',nbapp_cv
-    write(lunout,*)'nbapp_wk=',nbapp_wk
-    write(lunout,*)'iflag_ener_conserv=',iflag_ener_conserv
-    write(lunout,*)'ok_conserv_q=',ok_conserv_q
-    write(lunout,*)'iflag_fisrtilp_qsat=',iflag_fisrtilp_qsat
-    write(lunout,*)'iflag_bergeron=',iflag_bergeron
-    write(lunout,*)' epmax = ', epmax
-    write(lunout,*)' coef_epmax_cape = ', coef_epmax_cape
-    write(lunout,*)' ok_adj_ema = ', ok_adj_ema
-    write(lunout,*)' iflag_clw = ', iflag_clw
-    write(lunout,*)' cld_lc_lsc = ', cld_lc_lsc
-    write(lunout,*)' cld_lc_con = ', cld_lc_con
-    write(lunout,*)' cld_tau_lsc = ', cld_tau_lsc
-    write(lunout,*)' cld_tau_con = ', cld_tau_con
-    write(lunout,*)' ffallv_lsc = ', ffallv_lsc
-    write(lunout,*)' ffallv_con = ', ffallv_con
-    write(lunout,*)' coef_eva = ', coef_eva
-    write(lunout,*)' reevap_ice = ', reevap_ice
-    write(lunout,*)' iflag_pdf = ', iflag_pdf
-    write(lunout,*)' iflag_cld_th = ', iflag_cld_th
-    write(lunout,*)' iflag_cld_cv = ', iflag_cld_cv
-    write(lunout,*)' tau_cld_cv = ', tau_cld_cv
-    write(lunout,*)' coefw_cld_cv = ', coefw_cld_cv
-    write(lunout,*)' iflag_radia = ', iflag_radia
-    write(lunout,*)' iflag_rrtm = ', iflag_rrtm
-    write(lunout,*)' NSW = ', NSW
-    write(lunout,*)' iflag_albedo = ', iflag_albedo !albedo SB
-    write(lunout,*)' ok_chlorophyll =',ok_chlorophyll ! albedo SB
-    write(lunout,*)' iflag_ratqs = ', iflag_ratqs
-    write(lunout,*)' seuil_inversion = ', seuil_inversion
-    write(lunout,*)' fact_cldcon = ', fact_cldcon
-    write(lunout,*)' facttemps = ', facttemps
-    write(lunout,*)' ok_newmicro = ',ok_newmicro 
-    write(lunout,*)' ratqsbas = ',ratqsbas 
-    write(lunout,*)' ratqshaut = ',ratqshaut 
-    write(lunout,*)' tau_ratqs = ',tau_ratqs 
-    write(lunout,*)' top_height = ',top_height 
-    write(lunout,*)' rad_froid = ',rad_froid
-    write(lunout,*)' rad_chau1 = ',rad_chau1
-    write(lunout,*)' rad_chau2 = ',rad_chau2
-    write(lunout,*)' t_glace_min = ',t_glace_min
-    write(lunout,*)' t_glace_max = ',t_glace_max
-    write(lunout,*)' exposant_glace = ',exposant_glace
-    write(lunout,*)' iflag_t_glace = ',iflag_t_glace
-    write(lunout,*)' iflag_cloudth_vert = ',iflag_cloudth_vert
-    write(lunout,*)' iflag_rain_incloud_vol = ',iflag_rain_incloud_vol
-    write(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo
-    write(lunout,*)' rei_min = ',rei_min
-    write(lunout,*)' rei_max = ',rei_max
-    write(lunout,*)' overlap = ',overlap 
-    write(lunout,*)' cdmmax = ',cdmmax 
-    write(lunout,*)' cdhmax = ',cdhmax 
-    write(lunout,*)' ksta = ',ksta 
-    write(lunout,*)' ksta_ter = ',ksta_ter 
-    write(lunout,*)' f_ri_cd_min = ',f_ri_cd_min 
-    write(lunout,*)' ok_kzmin = ',ok_kzmin 
-    write(lunout,*)' pbl_lmixmin_alpha = ',pbl_lmixmin_alpha
-    write(lunout,*)' fmagic = ',fmagic
-    write(lunout,*)' pmagic = ',pmagic
-    write(lunout,*)' ok_ade = ',ok_ade
-    write(lunout,*)' ok_volcan = ',ok_volcan
-    write(lunout,*)' ok_aie = ',ok_aie
-    write(lunout,*)' ok_alw = ',ok_alw
-    write(lunout,*)' aerosol_couple = ', aerosol_couple
-    write(lunout,*)' chemistry_couple = ', chemistry_couple
-    write(lunout,*)' flag_aerosol = ', flag_aerosol
-    write(lunout,*)' flag_aerosol_strat= ', flag_aerosol_strat
-    write(lunout,*)' new_aod = ', new_aod
-    write(lunout,*)' aer_type = ',aer_type
-    write(lunout,*)' bl95_b0 = ',bl95_b0
-    write(lunout,*)' bl95_b1 = ',bl95_b1
-    write(lunout,*)' lev_histhf = ',lev_histhf 
-    write(lunout,*)' lev_histday = ',lev_histday 
-    write(lunout,*)' lev_histmth = ',lev_histmth 
-    write(lunout,*)' lev_histins = ',lev_histins
-    write(lunout,*)' lev_histLES = ',lev_histLES
-    write(lunout,*)' lev_histdayNMC = ',lev_histdayNMC
-    write(lunout,*)' levout_histNMC = ',levout_histNMC
-    write(lunout,*)' ok_histNMC = ',ok_histNMC
-    write(lunout,*)' freq_outNMC = ',freq_outNMC
-    write(lunout,*)' freq_calNMC = ',freq_calNMC
-    write(lunout,*)' iflag_pbl = ', iflag_pbl
+    WRITE(lunout,*)' ##############################################'
+    WRITE(lunout,*)' Configuration des parametres de la physique: '
+    WRITE(lunout,*)' Type ocean = ', type_ocean
+    WRITE(lunout,*)' Version ocean = ', version_ocean
+    WRITE(lunout,*)' Config veget = ', ok_veget,type_veget
+    WRITE(lunout,*)' Snow model SISVAT : ok_snow = ', ok_snow
+    WRITE(lunout,*)' Config xml pour XIOS : ok_all_xml = ', ok_all_xml
+    WRITE(lunout,*)' Sortie journaliere = ', ok_journe
+    WRITE(lunout,*)' Sortie haute frequence = ', ok_hf
+    WRITE(lunout,*)' Sortie mensuelle = ', ok_mensuel
+    WRITE(lunout,*)' Sortie instantanee = ', ok_instan
+    WRITE(lunout,*)' Frequence appel simulateur ISCCP, freq_ISCCP =', freq_ISCCP
+    WRITE(lunout,*)' Frequence appel simulateur ISCCP, ecrit_ISCCP =', ecrit_ISCCP
+    WRITE(lunout,*)' Frequence appel simulateur COSP, freq_COSP =', freq_COSP
+    WRITE(lunout,*)' Frequence appel simulateur AIRS, freq_AIRS =', freq_AIRS
+    WRITE(lunout,*)' Sortie bilan d''energie, ip_ebil_phy =', ip_ebil_phy
+    WRITE(lunout,*)' Excentricite = ',R_ecc
+    WRITE(lunout,*)' Equinoxe = ',R_peri
+    WRITE(lunout,*)' Inclinaison =',R_incl
+    WRITE(lunout,*)' Constante solaire =',solaire
+    WRITE(lunout,*)' ok_suntime_rrtm =',ok_suntime_rrtm
+    WRITE(lunout,*)' co2_ppm =',co2_ppm
+    WRITE(lunout,*)' RCO2_act = ',RCO2_act
+    WRITE(lunout,*)' CH4_ppb =',CH4_ppb,' RCH4_act = ',RCH4_act
+    WRITE(lunout,*)' N2O_ppb =',N2O_ppb,' RN2O_act=  ',RN2O_act
+    WRITE(lunout,*)' CFC11_ppt=',CFC11_ppt,' RCFC11_act=  ',RCFC11_act
+    WRITE(lunout,*)' CFC12_ppt=',CFC12_ppt,' RCFC12_act=  ',RCFC12_act
+    WRITE(lunout,*)' RCO2_per = ',RCO2_per,' RCH4_per = ', RCH4_per
+    WRITE(lunout,*)' RN2O_per = ',RN2O_per,' RCFC11_per = ', RCFC11_per
+    WRITE(lunout,*)' RCFC12_per = ',RCFC12_per
+    WRITE(lunout,*)' cvl_comp_threshold=', cvl_comp_threshold
+    WRITE(lunout,*)' cvl_sig2feed=', cvl_sig2feed
+    WRITE(lunout,*)' cvl_corr=', cvl_corr
+    WRITE(lunout,*)'ok_lic_melt=', ok_lic_melt
+    WRITE(lunout,*)'ok_lic_cond=', ok_lic_cond
+    WRITE(lunout,*)'iflag_cycle_diurne=',iflag_cycle_diurne
+    WRITE(lunout,*)'soil_model=',soil_model
+    WRITE(lunout,*)'new_oliq=',new_oliq
+    WRITE(lunout,*)'ok_orodr=',ok_orodr
+    WRITE(lunout,*)'ok_orolf=',ok_orolf
+    WRITE(lunout,*)'ok_limitvrai=',ok_limitvrai
+    WRITE(lunout,*)'nbapp_rad=',nbapp_rad
+    WRITE(lunout,*)'iflag_con=',iflag_con
+    WRITE(lunout,*)'nbapp_cv=',nbapp_cv
+    WRITE(lunout,*)'nbapp_wk=',nbapp_wk
+    WRITE(lunout,*)'iflag_ener_conserv=',iflag_ener_conserv
+    WRITE(lunout,*)'ok_conserv_q=',ok_conserv_q
+    WRITE(lunout,*)'iflag_fisrtilp_qsat=',iflag_fisrtilp_qsat
+    WRITE(lunout,*)'iflag_bergeron=',iflag_bergeron
+    WRITE(lunout,*)' epmax = ', epmax
+    WRITE(lunout,*)' coef_epmax_cape = ', coef_epmax_cape
+    WRITE(lunout,*)' ok_adj_ema = ', ok_adj_ema
+    WRITE(lunout,*)' iflag_clw = ', iflag_clw
+    WRITE(lunout,*)' cld_lc_lsc = ', cld_lc_lsc
+    WRITE(lunout,*)' cld_lc_con = ', cld_lc_con
+    WRITE(lunout,*)' cld_tau_lsc = ', cld_tau_lsc
+    WRITE(lunout,*)' cld_tau_con = ', cld_tau_con
+    WRITE(lunout,*)' ffallv_lsc = ', ffallv_lsc
+    WRITE(lunout,*)' ffallv_con = ', ffallv_con
+    WRITE(lunout,*)' coef_eva = ', coef_eva
+    WRITE(lunout,*)' reevap_ice = ', reevap_ice
+    WRITE(lunout,*)' iflag_pdf = ', iflag_pdf
+    WRITE(lunout,*)' iflag_cld_th = ', iflag_cld_th
+    WRITE(lunout,*)' iflag_cld_cv = ', iflag_cld_cv
+    WRITE(lunout,*)' tau_cld_cv = ', tau_cld_cv
+    WRITE(lunout,*)' coefw_cld_cv = ', coefw_cld_cv
+    WRITE(lunout,*)' iflag_radia = ', iflag_radia
+    WRITE(lunout,*)' iflag_rrtm = ', iflag_rrtm
+    WRITE(lunout,*)' NSW = ', NSW
+    WRITE(lunout,*)' iflag_albedo = ', iflag_albedo !albedo SB
+    WRITE(lunout,*)' ok_chlorophyll =',ok_chlorophyll ! albedo SB
+    WRITE(lunout,*)' iflag_ratqs = ', iflag_ratqs
+    WRITE(lunout,*)' seuil_inversion = ', seuil_inversion
+    WRITE(lunout,*)' fact_cldcon = ', fact_cldcon
+    WRITE(lunout,*)' facttemps = ', facttemps
+    WRITE(lunout,*)' ok_newmicro = ',ok_newmicro 
+    WRITE(lunout,*)' ratqsbas = ',ratqsbas 
+    WRITE(lunout,*)' ratqshaut = ',ratqshaut 
+    WRITE(lunout,*)' tau_ratqs = ',tau_ratqs 
+    WRITE(lunout,*)' top_height = ',top_height 
+    WRITE(lunout,*)' rad_froid = ',rad_froid
+    WRITE(lunout,*)' rad_chau1 = ',rad_chau1
+    WRITE(lunout,*)' rad_chau2 = ',rad_chau2
+    WRITE(lunout,*)' t_glace_min = ',t_glace_min
+    WRITE(lunout,*)' t_glace_max = ',t_glace_max
+    WRITE(lunout,*)' exposant_glace = ',exposant_glace
+    WRITE(lunout,*)' iflag_t_glace = ',iflag_t_glace
+    WRITE(lunout,*)' iflag_cloudth_vert = ',iflag_cloudth_vert
+    WRITE(lunout,*)' iflag_rain_incloud_vol = ',iflag_rain_incloud_vol
+    WRITE(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo
+    WRITE(lunout,*)' rei_min = ',rei_min
+    WRITE(lunout,*)' rei_max = ',rei_max
+    WRITE(lunout,*)' overlap = ',overlap 
+    WRITE(lunout,*)' cdmmax = ',cdmmax 
+    WRITE(lunout,*)' cdhmax = ',cdhmax 
+    WRITE(lunout,*)' ksta = ',ksta 
+    WRITE(lunout,*)' ksta_ter = ',ksta_ter 
+    WRITE(lunout,*)' f_ri_cd_min = ',f_ri_cd_min 
+    WRITE(lunout,*)' ok_kzmin = ',ok_kzmin 
+    WRITE(lunout,*)' pbl_lmixmin_alpha = ',pbl_lmixmin_alpha
+    WRITE(lunout,*)' fmagic = ',fmagic
+    WRITE(lunout,*)' pmagic = ',pmagic
+    WRITE(lunout,*)' ok_ade = ',ok_ade
+    WRITE(lunout,*)' ok_volcan = ',ok_volcan
+    WRITE(lunout,*)' ok_aie = ',ok_aie
+    WRITE(lunout,*)' ok_alw = ',ok_alw
+    WRITE(lunout,*)' aerosol_couple = ', aerosol_couple
+    WRITE(lunout,*)' chemistry_couple = ', chemistry_couple
+    WRITE(lunout,*)' flag_aerosol = ', flag_aerosol
+    WRITE(lunout,*)' flag_aerosol_strat= ', flag_aerosol_strat
+    WRITE(lunout,*) ' flag_aer_feedback= ', flag_aer_feedback
+    WRITE(lunout,*)' new_aod = ', new_aod
+    WRITE(lunout,*)' aer_type = ',aer_type
+    WRITE(lunout,*)' bl95_b0 = ',bl95_b0
+    WRITE(lunout,*)' bl95_b1 = ',bl95_b1
+    WRITE(lunout,*)' lev_histhf = ',lev_histhf 
+    WRITE(lunout,*)' lev_histday = ',lev_histday 
+    WRITE(lunout,*)' lev_histmth = ',lev_histmth 
+    WRITE(lunout,*)' lev_histins = ',lev_histins
+    WRITE(lunout,*)' lev_histLES = ',lev_histLES
+    WRITE(lunout,*)' lev_histdayNMC = ',lev_histdayNMC
+    WRITE(lunout,*)' levout_histNMC = ',levout_histNMC
+    WRITE(lunout,*)' ok_histNMC = ',ok_histNMC
+    WRITE(lunout,*)' freq_outNMC = ',freq_outNMC
+    WRITE(lunout,*)' freq_calNMC = ',freq_calNMC
+    WRITE(lunout,*)' iflag_pbl = ', iflag_pbl
 !FC
-    write(lunout,*)' ifl_pbltree = ', ifl_pbltree
-    write(lunout,*)' Cd_frein = ', Cd_frein
-    write(lunout,*)' iflag_pbl_split = ', iflag_pbl_split
-    write(lunout,*)' iflag_order2_sollw = ', iflag_order2_sollw
-    write(lunout,*)' iflag_thermals = ', iflag_thermals
-    write(lunout,*)' iflag_thermals_ed = ', iflag_thermals_ed
-    write(lunout,*)' fact_thermals_ed_dz = ', fact_thermals_ed_dz
-    write(lunout,*)' iflag_thermals_optflux = ', iflag_thermals_optflux
-    write(lunout,*)' iflag_thermals_closure = ', iflag_thermals_closure
-    write(lunout,*)' iflag_clos = ', iflag_clos
-    write(lunout,*)' coef_clos_ls = ', coef_clos_ls
-    write(lunout,*)' type_run = ',type_run 
-    write(lunout,*)' ok_cosp = ',ok_cosp
-    write(lunout,*)' ok_airs = ',ok_airs
-
-    write(lunout,*)' ok_mensuelCOSP = ',ok_mensuelCOSP
-    write(lunout,*)' ok_journeCOSP = ',ok_journeCOSP
-    write(lunout,*)' ok_hfCOSP =',ok_hfCOSP
-    write(lunout,*)' solarlong0 = ', solarlong0
-    write(lunout,*)' qsol0 = ', qsol0
-    write(lunout,*)' evap0 = ', evap0
-    write(lunout,*)' albsno0 = ', albsno0
-    write(lunout,*)' iflag_sic = ', iflag_sic
-    write(lunout,*)' inertie_sol = ', inertie_sol
-    write(lunout,*)' inertie_sic = ', inertie_sic
-    write(lunout,*)' inertie_lic = ', inertie_lic
-    write(lunout,*)' inertie_sno = ', inertie_sno
-    write(lunout,*)' f_cdrag_ter = ',f_cdrag_ter
-    write(lunout,*)' f_cdrag_oce = ',f_cdrag_oce
-    write(lunout,*)' f_rugoro = ',f_rugoro
-    write(lunout,*)' z0min = ',z0min
-    write(lunout,*)' supcrit1 = ', supcrit1
-    write(lunout,*)' supcrit2 = ', supcrit2
-    write(lunout,*)' iflag_mix = ', iflag_mix
-    write(lunout,*)' iflag_mix_adiab = ', iflag_mix_adiab
-    write(lunout,*)' scut = ', scut
-    write(lunout,*)' qqa1 = ', qqa1
-    write(lunout,*)' qqa2 = ', qqa2
-    write(lunout,*)' gammas = ', gammas
-    write(lunout,*)' Fmax = ', Fmax
-    write(lunout,*)' tmax_fonte_cv = ', tmax_fonte_cv
-    write(lunout,*)' alphas = ', alphas
-    write(lunout,*)' iflag_wake = ', iflag_wake
-    write(lunout,*)' alp_offset = ', alp_offset
+    WRITE(lunout,*)' ifl_pbltree = ', ifl_pbltree
+    WRITE(lunout,*)' Cd_frein = ', Cd_frein
+    WRITE(lunout,*)' iflag_pbl_split = ', iflag_pbl_split
+    WRITE(lunout,*)' iflag_order2_sollw = ', iflag_order2_sollw
+    WRITE(lunout,*)' iflag_thermals = ', iflag_thermals
+    WRITE(lunout,*)' iflag_thermals_ed = ', iflag_thermals_ed
+    WRITE(lunout,*)' fact_thermals_ed_dz = ', fact_thermals_ed_dz
+    WRITE(lunout,*)' iflag_thermals_optflux = ', iflag_thermals_optflux
+    WRITE(lunout,*)' iflag_thermals_closure = ', iflag_thermals_closure
+    WRITE(lunout,*)' iflag_clos = ', iflag_clos
+    WRITE(lunout,*)' coef_clos_ls = ', coef_clos_ls
+    WRITE(lunout,*)' type_run = ',type_run 
+    WRITE(lunout,*)' ok_cosp = ',ok_cosp
+    WRITE(lunout,*)' ok_airs = ',ok_airs
+
+    WRITE(lunout,*)' ok_mensuelCOSP = ',ok_mensuelCOSP
+    WRITE(lunout,*)' ok_journeCOSP = ',ok_journeCOSP
+    WRITE(lunout,*)' ok_hfCOSP =',ok_hfCOSP
+    WRITE(lunout,*)' solarlong0 = ', solarlong0
+    WRITE(lunout,*)' qsol0 = ', qsol0
+    WRITE(lunout,*)' evap0 = ', evap0
+    WRITE(lunout,*)' albsno0 = ', albsno0
+    WRITE(lunout,*)' iflag_sic = ', iflag_sic
+    WRITE(lunout,*)' inertie_sol = ', inertie_sol
+    WRITE(lunout,*)' inertie_sic = ', inertie_sic
+    WRITE(lunout,*)' inertie_lic = ', inertie_lic
+    WRITE(lunout,*)' inertie_sno = ', inertie_sno
+    WRITE(lunout,*)' f_cdrag_ter = ',f_cdrag_ter
+    WRITE(lunout,*)' f_cdrag_oce = ',f_cdrag_oce
+    WRITE(lunout,*)' f_rugoro = ',f_rugoro
+    WRITE(lunout,*)' z0min = ',z0min
+    WRITE(lunout,*)' supcrit1 = ', supcrit1
+    WRITE(lunout,*)' supcrit2 = ', supcrit2
+    WRITE(lunout,*)' iflag_mix = ', iflag_mix
+    WRITE(lunout,*)' iflag_mix_adiab = ', iflag_mix_adiab
+    WRITE(lunout,*)' scut = ', scut
+    WRITE(lunout,*)' qqa1 = ', qqa1
+    WRITE(lunout,*)' qqa2 = ', qqa2
+    WRITE(lunout,*)' gammas = ', gammas
+    WRITE(lunout,*)' Fmax = ', Fmax
+    WRITE(lunout,*)' tmax_fonte_cv = ', tmax_fonte_cv
+    WRITE(lunout,*)' alphas = ', alphas
+    WRITE(lunout,*)' iflag_wake = ', iflag_wake
+    WRITE(lunout,*)' alp_offset = ', alp_offset
     ! nrlmd le 10/04/2012
-    write(lunout,*)' iflag_trig_bl = ', iflag_trig_bl
-    write(lunout,*)' s_trig = ', s_trig
-    write(lunout,*)' tau_trig_shallow = ', tau_trig_shallow
-    write(lunout,*)' tau_trig_deep = ', tau_trig_deep
-    write(lunout,*)' iflag_clos_bl = ', iflag_clos_bl
+    WRITE(lunout,*) ' iflag_trig_bl = ', iflag_trig_bl
+    WRITE(lunout,*) ' s_trig = ', s_trig
+    WRITE(lunout,*) ' tau_trig_shallow = ', tau_trig_shallow
+    WRITE(lunout,*) ' tau_trig_deep = ', tau_trig_deep
+    WRITE(lunout,*) ' iflag_clos_bl = ', iflag_clos_bl
     ! fin nrlmd le 10/04/2012
 
-    write(lunout,*)' lonmin lonmax latmin latmax bilKP_ins =',&
+    WRITE(lunout,*) ' lonmin lonmax latmin latmax bilKP_ins =',&
          lonmin_ins, lonmax_ins, latmin_ins, latmax_ins
-    write(lunout,*)' ecrit_ hf, ins, day, mth, reg, tra, ISCCP, LES',&
+    WRITE(lunout,*) ' ecrit_ hf, ins, day, mth, reg, tra, ISCCP, LES',&
          ecrit_hf, ecrit_ins, ecrit_day, ecrit_mth, ecrit_reg, ecrit_tra, ecrit_ISCCP, ecrit_LES
 
-    write(lunout,*) 'ok_strato = ', ok_strato
-    write(lunout,*) 'ok_hines = ',  ok_hines
-    write(lunout,*) 'ok_gwd_rando = ',  ok_gwd_rando
-    write(lunout,*) 'ok_qch4 = ',  ok_qch4
-    write(lunout,*) 'gwd_rando_ruwmax = ', gwd_rando_ruwmax
-    write(lunout,*) 'gwd_rando_sat = ', gwd_rando_sat
-    write(lunout,*) 'gwd_front_ruwmax = ', gwd_front_ruwmax
-    write(lunout,*) 'gwd_front_sat = ', gwd_front_sat
-    write(lunout,*) 'SSO gkdrag =',gkdrag
-    write(lunout,*) 'SSO grahilo=',grahilo
-    write(lunout,*) 'SSO grcrit=',grcrit
-    write(lunout,*) 'SSO gfrcrit=',gfrcrit
-    write(lunout,*) 'SSO gkwake=',gkwake
-    write(lunout,*) 'SSO gklift=',gklift
-    write(lunout,*) 'adjust_tropopause = ', adjust_tropopause
-    write(lunout,*) 'ok_daily_climoz = ',ok_daily_climoz
-    write(lunout,*) 'read_climoz = ', read_climoz
-    write(lunout,*) 'carbon_cycle_tr = ', carbon_cycle_tr
-    write(lunout,*) 'carbon_cycle_cpl = ', carbon_cycle_cpl
+    WRITE(lunout,*) ' ok_strato = ', ok_strato
+    WRITE(lunout,*) ' ok_hines = ',  ok_hines
+    WRITE(lunout,*) ' ok_gwd_rando = ',  ok_gwd_rando
+    WRITE(lunout,*) ' ok_qch4 = ',  ok_qch4
+    WRITE(lunout,*) ' gwd_rando_ruwmax = ', gwd_rando_ruwmax
+    WRITE(lunout,*) ' gwd_rando_sat = ', gwd_rando_sat
+    WRITE(lunout,*) ' gwd_front_ruwmax = ', gwd_front_ruwmax
+    WRITE(lunout,*) ' gwd_front_sat = ', gwd_front_sat
+    WRITE(lunout,*) ' SSO gkdrag =',gkdrag
+    WRITE(lunout,*) ' SSO grahilo=',grahilo
+    WRITE(lunout,*) ' SSO grcrit=',grcrit
+    WRITE(lunout,*) ' SSO gfrcrit=',gfrcrit
+    WRITE(lunout,*) ' SSO gkwake=',gkwake
+    WRITE(lunout,*) ' SSO gklift=',gklift
+    WRITE(lunout,*) ' adjust_tropopause = ', adjust_tropopause
+    WRITE(lunout,*) ' ok_daily_climoz = ',ok_daily_climoz
+    WRITE(lunout,*) ' read_climoz = ', read_climoz
+    WRITE(lunout,*) ' carbon_cycle_tr = ', carbon_cycle_tr
+    WRITE(lunout,*) ' carbon_cycle_cpl = ', carbon_cycle_cpl
 
     !$OMP END MASTER
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3525)
@@ -1067,4 +1067,9 @@
     INTEGER, SAVE :: flag_aerosol_strat
     !$OMP THREADPRIVATE(flag_aerosol_strat)
+    !
+    !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
+    LOGICAL, SAVE :: flag_aer_feedback
+    !$OMP THREADPRIVATE(flag_aer_feedback)
+
     !c-fin STRAT AEROSOL
     !
@@ -1234,5 +1239,6 @@
             iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
             ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
-            chemistry_couple, flag_aerosol, flag_aerosol_strat, new_aod, &
+            chemistry_couple, flag_aerosol, flag_aerosol_strat,         &
+            flag_aer_feedback, new_aod, &
             flag_bc_internal_mixture, bl95_b0, bl95_b1, &
                                 ! nv flags pour la convection et les
@@ -3908,5 +3914,5 @@
                cldfrarad, cldemirad, cldtaurad, &
                ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
-               flag_aerosol, flag_aerosol_strat, &
+               flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
                tau_aero, piz_aero, cg_aero, &
                tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & 
@@ -3994,5 +4000,5 @@
                      cldfrarad, cldemirad, cldtaurad, &
                      ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
-                     flag_aerosol, flag_aerosol_strat, &
+                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
                      tau_aero, piz_aero, cg_aero, &
                      tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
@@ -4873,4 +4879,5 @@
       dryaod_diag=.FALSE.
       ok_4xCO2atm= .FALSE.
+!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
 
       IF (is_master) then
@@ -4911,4 +4918,5 @@
       call bcast(dryaod_diag)
       call bcast(ok_4xCO2atm)
+!      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
 #endif
     endif
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3525)
@@ -17,5 +17,5 @@
    cldfra, cldemi, cldtaupd,&
    ok_ade, ok_aie, ok_volcan, flag_aerosol,&
-   flag_aerosol_strat,&
+   flag_aerosol_strat, flag_aer_feedback, &
    tau_aero, piz_aero, cg_aero,&
    tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB pour RRTM
@@ -104,4 +104,5 @@
   ! flag_aerosol-input-I- aerosol flag from 0 to 6
   ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (0, 1, 2)
+  ! flag_aer_feedback-input-I- activate aerosol radiative feedback (T, F)
   ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
   ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
@@ -201,4 +202,5 @@
   INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
   INTEGER, INTENT(in)  :: flag_aerosol_strat                             ! use stratospheric aerosols
+  LOGICAL, INTENT(in)  :: flag_aer_feedback                              ! activate aerosol radiative feedback
   REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
   REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
@@ -873,5 +875,6 @@
          ZTOPLWAIAERO,ZSOLLWAIAERO, &
          ZLWADAERO, & !--NL
-         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat) ! flags aerosols
+         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, &
+         flag_aer_feedback) ! flags aerosols
             
 !        print *,'RADLWSW: apres RECMWF'
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/recmwf_aero.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/recmwf_aero.F90	(revision 3521)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/recmwf_aero.F90	(revision 3525)
@@ -37,5 +37,5 @@
  & PLWADAERO,& !--NL
 !..end
- & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat)
+ & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,flag_aer_feedback)
 !--fin
 
@@ -87,4 +87,5 @@
 ! flag_aerosol-input-I- aerosol flag from 0 to 7
 ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (T/F)
+! flag_aer_feedback-input-I- use aerosols radiative effect flag (T/F)
 ! PPIZA_NAT  : (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosol
 ! PCGA_NAT   : (KPROMA,KLEV,NSW); Assymetry factor for natural aerosol
@@ -218,5 +219,6 @@
 INTEGER, INTENT(in)  :: flag_aerosol           ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
 LOGICAL, INTENT(in)  :: flag_aerosol_strat     ! use stratospheric aerosols
-REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWADAERO(KPROMA), PSOLSWADAERO(KPROMA)       ! Aerosol direct forcing at TOA and surface
+LOGICAL, INTENT(in)  :: flag_aer_feedback      ! use aerosols radiative feedback
+REAL(KIND=JPRB)   ,INTENT(out)   :: PTOPSWADAERO(KPROMA), PSOLSWADAERO(KPROMA)       ! Aerosol direct forcing at TOA and surface
 REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWAD0AERO(KPROMA), PSOLSWAD0AERO(KPROMA)     ! Aerosol direct forcing at TOA and surface
 REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWAIAERO(KPROMA), PSOLSWAIAERO(KPROMA)       ! ditto, indirect
@@ -305,5 +307,5 @@
 ! FALSE: fluxes use no aerosols (case 1) 
 ! to be used only for maintaining bit reproducibility with aerosol diagnostics activated
-LOGICAL :: AEROSOLFEEDBACK_ACTIVE = .TRUE.
+ LOGICAL :: AEROSOLFEEDBACK_ACTIVE ! now externalized from .def files
 
 !OB - Fluxes including aerosol effects
@@ -342,4 +344,7 @@
 IBEG=KST
 IEND=KEND
+
+AEROSOLFEEDBACK_ACTIVE = flag_aer_feedback !NL: externalize aer feedback
+
 
 !*       1.    PREPARATORY WORK
