Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/aerophys.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/aerophys.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/aerophys.F90	(revision 4743)
@@ -17,4 +17,11 @@
   REAL, PARAMETER                        :: mSatom=32.06*1.66E-27    ! Mass of a S atom [kg]
   REAL, PARAMETER                        :: mOCSmol=60.07*1.66E-27   ! Mass of an OCS molecule [kg]
+  REAL, PARAMETER                        :: mClatom=35.45*1.66E-27   ! Mass of an Cl atom [kg]
+  REAL, PARAMETER                        :: mHClmol=36.46*1.66E-27   ! Mass of an HCl molecule [kg]
+  REAL, PARAMETER                        :: mBratom=79.90*1.66E-27   ! Mass of an Br atom [kg]
+  REAL, PARAMETER                        :: mHBrmol=80.92*1.66E-27   ! Mass of an HBr molecule [kg]
+  REAL, PARAMETER                        :: mNOmol=30.01*1.66E-27    ! Mass of an NO molecule [kg]
+  REAL, PARAMETER                        :: mNO2mol=46.01*1.66E-27   ! Mass of an NO2 molecule [kg]
+  REAL, PARAMETER                        :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg] 
 !
 END MODULE aerophys
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/interp_sulf_input.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/interp_sulf_input.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/interp_sulf_input.F90	(revision 4743)
@@ -12,5 +12,5 @@
   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
   USE phys_local_var_mod, ONLY : budg_3D_backgr_ocs, budg_3D_backgr_so2
-  USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime
+  USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime, H2SO4_lifetime, O3_clim
   USE mod_phys_lmdz_para 
   USE dimphy
@@ -19,4 +19,5 @@
   USE aerophys
   USE YOMCST
+  USE strataer_local_var_mod, ONLY : flag_newclim_file,flag_verbose_strataer
 
   IMPLICIT NONE
@@ -42,14 +43,14 @@
   REAL paprs_glo(klon_glo,klev+1)
 
-  REAL, POINTER:: latitude(:)
+  REAL, ALLOCATABLE :: latitude(:)
 ! (of input data sorted in strictly ascending order)
 
-  REAL, POINTER:: longitude(:)
+  REAL, ALLOCATABLE :: longitude(:)
 ! (of input data sorted in strictly ascending order)
 
-  REAL, POINTER:: time(:)
+  REAL, ALLOCATABLE :: time(:)
 ! (of input data sorted in strictly ascending order)
 
-  REAL, POINTER:: lev(:)
+  REAL, ALLOCATABLE :: lev(:)
 ! levels of input data
 
@@ -62,4 +63,11 @@
   REAL OCS_clim_glo(klon_glo,klev)
   REAL SO2_clim_glo(klon_glo,klev)
+  
+  REAL, ALLOCATABLE :: O3_clim_in(:, :, :, :)
+  REAL, ALLOCATABLE :: O3_clim_mth(:, :, :)
+  REAL, ALLOCATABLE :: O3_clim_tmp(:, :)
+  REAL O3_clim_glo(klon_glo,klev)
+  
+  ! fixed climatos
   REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :)
   REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :)
@@ -70,4 +78,10 @@
   REAL OCS_lifetime_glo(klon_glo,klev)
   REAL SO2_lifetime_glo(klon_glo,klev)
+  
+  REAL, ALLOCATABLE :: H2SO4_lifetime_in(:, :, :, :)
+  REAL, ALLOCATABLE :: H2SO4_lifetime_mth(:, :, :)
+  REAL, ALLOCATABLE :: H2SO4_lifetime_tmp(:, :)
+
+  REAL H2SO4_lifetime_glo(klon_glo,klev)
 !
   REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:)
@@ -79,6 +93,8 @@
 
 ! For NetCDF:
-  INTEGER ncid_in  ! IDs for input files
-  INTEGER varid, ncerr
+  INTEGER             :: ncid_in  ! IDs for input files
+  INTEGER             :: varid, ncerr
+  CHARACTER (LEN=3)   :: nc_lat, nc_lon
+  CHARACTER (LEN=256) :: nc_fname
     
   INTEGER, PARAMETER :: lev_input=17
@@ -103,44 +119,70 @@
   IF (is_mpi_root.AND.is_omp_root) THEN
 
-!--reading emission files
-    CALL nf95_open("ocs_so2_annual_lmdz.nc", nf90_nowrite, ncid_in)
-
+    !--init ncdf variables
+    IF(flag_newclim_file) THEN
+      nc_fname = "ocs_so2_h2so4_annual_lmdz.nc"
+      nc_lat   = "LAT"
+      nc_lon   = "LON"
+    ELSE
+      ! old file for retro compatibility
+      nc_fname = "ocs_so2_annual_lmdz.nc"
+      nc_lat   = "lat"
+      nc_lon   = "lon"
+    ENDIF
+
+    !--reading emission files
+    CALL nf95_open(nc_fname, nf90_nowrite, ncid_in)
+    IF(flag_verbose_strataer) print *, "Open file ", nc_fname
+    
     CALL nf95_inq_varid(ncid_in, "LEV", varid)
     CALL nf95_gw_var(ncid_in, varid, lev)
     n_lev = size(lev)
-
-    CALL nf95_inq_varid(ncid_in, "lat", varid)
+    
+    CALL nf95_inq_varid(ncid_in, nc_lat, varid)
     CALL nf95_gw_var(ncid_in, varid, latitude)
     n_lat = size(latitude)
-
-    CALL nf95_inq_varid(ncid_in, "lon", varid)
+    
+    CALL nf95_inq_varid(ncid_in, nc_lon, varid)
     CALL nf95_gw_var(ncid_in, varid, longitude)
     n_lon = size(longitude)
-
+    
     CALL nf95_inq_varid(ncid_in, "TIME", varid)
     CALL nf95_gw_var(ncid_in, varid, time)
     n_mth = size(time)
-
+    
     IF (.NOT.ALLOCATED(OCS_clim_in))     ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth))
     IF (.NOT.ALLOCATED(SO2_clim_in))     ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth))
+    IF (.NOT.ALLOCATED(O3_clim_in))      ALLOCATE(O3_clim_in(n_lon, n_lat, n_lev, n_mth))
+    
     IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
     IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
-
+    IF (.NOT.ALLOCATED(H2SO4_lifetime_in)) ALLOCATE(H2SO4_lifetime_in(n_lon, n_lat, n_lev, n_mth))
+    
     CALL nf95_inq_varid(ncid_in, "OCS", varid)
     ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in)
-    print *,'code erreur OCS=', ncerr, varid
+    IF(flag_verbose_strataer) print *,'code erreur OCS=', ncerr, varid
 
     CALL nf95_inq_varid(ncid_in, "SO2", varid)
     ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in)
-    print *,'code erreur SO2=', ncerr, varid
+    IF(flag_verbose_strataer) print *,'code erreur SO2=', ncerr, varid
 
     CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid)
     ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in)
-    print *,'code erreur OCS_lifetime_in=', ncerr, varid
+    IF(flag_verbose_strataer) print *,'code erreur OCS_lifetime_in=', ncerr, varid
 
     CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid)
     ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in)
-    print *,'code erreur SO2_lifetime_in=', ncerr, varid
-
+    IF(flag_verbose_strataer) print *,'code erreur SO2_lifetime_in=', ncerr, varid
+    
+    IF(flag_newclim_file) THEN
+       CALL nf95_inq_varid(ncid_in, "O3", varid)
+       ncerr = nf90_get_var(ncid_in, varid, O3_clim_in)
+       IF(flag_verbose_strataer) print *,'code erreur O3=', ncerr, varid
+       
+       CALL nf95_inq_varid(ncid_in, "H2SO4_LIFET", varid)
+       ncerr = nf90_get_var(ncid_in, varid, H2SO4_lifetime_in)
+       IF(flag_verbose_strataer) print *,'code erreur H2SO4_lifetime_in=', ncerr, varid
+    ENDIF
+    
     CALL nf95_close(ncid_in)
 
@@ -149,9 +191,14 @@
     IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev))
     IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev))
+    IF (.NOT.ALLOCATED(O3_clim_mth))  ALLOCATE(O3_clim_mth(n_lon, n_lat, n_lev))
+    IF (.NOT.ALLOCATED(O3_clim_tmp))  ALLOCATE(O3_clim_tmp(klon_glo, n_lev))
+    
     IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev))
     IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev))
     IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev))
     IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev))
-
+    IF (.NOT.ALLOCATED(H2SO4_lifetime_mth)) ALLOCATE(H2SO4_lifetime_mth(n_lon, n_lat, n_lev))
+    IF (.NOT.ALLOCATED(H2SO4_lifetime_tmp)) ALLOCATE(H2SO4_lifetime_tmp(klon_glo, n_lev))
+    
 !---select the correct month, undo multiplication with 1.e12 (precision reasons)
 !---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio
@@ -161,4 +208,14 @@
       SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur)
       OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur)
+      
+      ! O3 from 2d model is not tracer, in VMR
+      IF(flag_newclim_file) THEN
+         H2SO4_lifetime_mth(:,j,:) = H2SO4_lifetime_in(:,n_lat+1-j,:,mth_cur)
+         ! new input files
+         O3_clim_mth(:,j,:) = 1.e-6*O3_clim_in(:,n_lat+1-j,:,mth_cur)
+      ELSE
+         H2SO4_lifetime_mth(:,j,:) = 1.e-6
+         O3_clim_mth(:,j,:) = 1.e-6
+      ENDIF
     ENDDO
 
@@ -166,10 +223,25 @@
     CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp)
     CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp)
+    CALL grid2dTo1d_glo(O3_clim_mth,O3_clim_tmp)
+    
     CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp)
     CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp)
+    CALL grid2dTo1d_glo(H2SO4_lifetime_mth,H2SO4_lifetime_tmp)
 
 !--set lifetime to very high value in uninsolated areas
     DO i=1, klon_glo
       DO kk=1, n_lev
+!--added tests on VMR of species for realism
+        IF (OCS_clim_tmp(i,kk).LT.1.e-20) THEN
+          OCS_clim_tmp(i,kk)=1.0e-20
+        ENDIF
+        IF (SO2_clim_tmp(i,kk).LT.1.e-20) THEN
+          SO2_clim_tmp(i,kk)=1.0e-20
+        ENDIF
+        !       50 ppbv min for O3
+        IF (O3_clim_tmp(i,kk).LT.50.0e-9) THEN
+           O3_clim_tmp(i,kk)=50.0e-9
+        ENDIF
+        
         IF (OCS_lifetime_tmp(i,kk)==0.0) THEN
           OCS_lifetime_tmp(i,kk)=1.0e12
@@ -178,4 +250,7 @@
           SO2_lifetime_tmp(i,kk)=1.0e12
         ENDIF
+       IF (H2SO4_lifetime_tmp(i,kk)==0.0) THEN
+          H2SO4_lifetime_tmp(i,kk)=1.0e12
+       ENDIF
       ENDDO
     ENDDO
@@ -186,6 +261,10 @@
      OCS_lifetime_glo(i,k)=0.0
      SO2_lifetime_glo(i,k)=0.0
+     O3_clim_glo(i,k)=0.0
+     
      OCS_clim_glo(i,k)=0.0
      SO2_clim_glo(i,k)=0.0
+     H2SO4_lifetime_glo(i,k)=0.0
+     
      DO kk=1, n_lev
       OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
@@ -195,4 +274,11 @@
            MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
            *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      IF(flag_newclim_file) THEN
+         H2SO4_lifetime_glo(i,k)=H2SO4_lifetime_glo(i,k)+ &
+              MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
+              -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
+              *H2SO4_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      ENDIF
+      
       OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ &
            MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
@@ -201,4 +287,10 @@
            MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
            *SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      IF(flag_newclim_file) THEN
+         O3_clim_glo(i,k)=O3_clim_glo(i,k)+ &
+              MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
+              -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
+              *O3_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      ENDIF
       ENDDO
     ENDDO
@@ -213,15 +305,17 @@
   CALL scatter(OCS_clim_glo, OCS_clim)
   CALL scatter(SO2_clim_glo, SO2_clim)
+  CALL scatter(O3_clim_glo, O3_clim)
   CALL scatter(OCS_lifetime_glo, OCS_lifetime)
   CALL scatter(SO2_lifetime_glo, SO2_lifetime)
-
+  CALL scatter(H2SO4_lifetime_glo, H2SO4_lifetime)
+  
   IF (is_mpi_root.AND.is_omp_root) THEN
 !
-    DEALLOCATE(OCS_clim_in,SO2_clim_in)
-    DEALLOCATE(OCS_clim_mth,SO2_clim_mth)
-    DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp)
-    DEALLOCATE(OCS_lifetime_in,SO2_lifetime_in)
-    DEALLOCATE(OCS_lifetime_mth,SO2_lifetime_mth)
-    DEALLOCATE(OCS_lifetime_tmp,SO2_lifetime_tmp)
+    DEALLOCATE(OCS_clim_in,SO2_clim_in,O3_clim_in)
+    DEALLOCATE(OCS_clim_mth,SO2_clim_mth,O3_clim_mth)
+    DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp,O3_clim_tmp)
+    DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in, H2SO4_lifetime_in)
+    DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth, H2SO4_lifetime_mth)
+    DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp, H2SO4_lifetime_tmp)
 !
   ENDIF !--is_mpi_root and is_omp_root
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/micphy_tstep.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/micphy_tstep.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/micphy_tstep.F90	(revision 4743)
@@ -14,5 +14,5 @@
   USE YOMCST, ONLY : RPI, RD, RG
   USE print_control_mod, ONLY: lunout
-  USE strataer_mod
+  USE strataer_local_var_mod
   
   IMPLICIT NONE
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/nucleation_tstep_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/nucleation_tstep_mod.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/nucleation_tstep_mod.F90	(revision 4743)
@@ -11,9 +11,9 @@
   USE infotrac_phy
   USE YOMCST, ONLY : RPI, RD, RMD, RKBOL, RNAVO
-
+  USE strataer_local_var_mod, ONLY : flag_new_nucl
+  
   IMPLICIT NONE
 
   ! input variables
-  LOGICAL, PARAMETER :: flag_new_nucl=.TRUE.
   REAL rhoa    ! H2SO4 number density [molecules/cm3]
   REAL t_seri  ! temperature (K)
@@ -170,31 +170,24 @@
 
   IF (t < 190.15) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)'
      t=190.15
   ENDIF
 
   IF (t > 300.15) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature > 300.15 K, using 300.15 K'
      t=300.15
   ENDIF
 
   IF (rh < 0.0001) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water < 0.0001, using 0.0001'
      rh=0.0001
   ENDIF
 
   IF (rh > 1.0) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water > 1 using 1'
-!     print *, 'rh=',rh
      rh=1.0
   ENDIF
 
   IF (rhoa < 1.e4) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'
      rhoa=1.e4
   ENDIF
 
   IF (rhoa > 1.e11) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3'
      rhoa=1.e11
   ENDIF
@@ -281,11 +274,8 @@
 
   IF (jnuc < 1.E-7) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev'): nucleation rate < 1E-7/cm3s, using 0.0/cm3s,'
      jnuc=0.0
   ENDIF
 
   IF (jnuc > 1.e10) THEN
-!     print *,'Warning (ilon=',ilon,'ilev=',ilev, &
-!        & '): nucleation rate > 1e10/cm3s, using 1e10/cm3s'
      jnuc=1.e10
   ENDIF
@@ -425,13 +415,10 @@
   !Temperature bounds
   IF (t.LE.165.) THEN
-!     print *,'Warning: temperature < 165.0 K, using 165.0 K in neutral nucleation calculation'
      tln=165.0
   ENDIF
   IF (t.LE.195.) THEN
-!     print *,'Warning: temperature < 195.0 K, using 195.0 K in ion-induced nucleation calculation'
      tli=195.0
   ENDIF
   IF (t.GE.400.) THEN
-!     print *,'Warning: temperature > 400. K, using 400. K in nucleation calculations'
      tln=400.
      tli=400.
@@ -440,17 +427,13 @@
   ! Saturation ratio bounds
   IF (satrat.LT.1.E-7) THEN
-!     print *,'Warning: saturation ratio of water < 1.E-7, using 1.E-7 in ion-induced nucleation calculation'
      satratli=1.E-7
   ENDIF
   IF (satrat.LT.1.E-5) THEN
-!     print *,'Warning: saturation ratio of water < 1.E-5, using 1.E-5 in neutral nucleation calculation'
      satratln=1.E-5
   ENDIF
   IF (satrat.GT.0.95) THEN
-!     print *,'Warning: saturation ratio of water > 0.95, using 0.95 in ion-induced nucleation calculation'
      satratli=0.95
   ENDIF
   IF (satrat.GT.1.0) THEN
-!     print *,'Warning: saturation ratio of water > 1 using 1 in neutral nucleation calculation'
      satratln=1.0
   ENDIF
@@ -458,14 +441,11 @@
   ! Sulfuric acid concentration bounds
   IF (rhoa.LE.1.E4) THEN
-!     print *,'Warning: sulfuric acid < 1e4 1/cm3, using 1e4 1/cm3 in nucleation calculation'
      rhoaln=1.E4
      rhoali=1.E4
   ENDIF
   IF (rhoa.GT.1.E13) THEN
-!     print *,'Warning: sulfuric acid > 1e13 1/cm3, using 1e13 1/cm3 in neutral nucleation calculation'
      rhoaln=1.E13
   ENDIF
   IF (rhoa.GT.1.E16) THEN
-!     print *,'Warning: sulfuric acid concentration > 1e16 1/cm3, using 1e16 1/cm3 in ion-induced nucleation calculation'
      rhoali=1.E16
   ENDIF
@@ -521,5 +501,4 @@
      ! Dimer formation rate
      jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))/2.*SQRT(t)*rhoa**2.
-!     jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))/2.*SQRT(t)*rhoa**2.
      ntot_n=1. !set to 1 
      na_n=1.   ! The critical cluster contains one molecules but the produced cluster contains 2 molecules
@@ -604,5 +583,4 @@
      na_n=x_n*ntot_n
      IF (na_n .LT. 1.) THEN
-!        print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'
         na_n=1.0
      ENDIF
@@ -664,5 +642,4 @@
      
      IF (kinetic_i) THEN    
-!        jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))*  &
         jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))*  &
              &  SQRT(tli)*rhoali !1/cm3s  
@@ -816,5 +793,4 @@
         na_i=x_i*ntot_i
         IF (na_i .LT. 1.) THEN
-!           print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'
            na_n=1.0
         ENDIF
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/ocs_to_so2.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/ocs_to_so2.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/ocs_to_so2.F90	(revision 4743)
@@ -8,5 +8,6 @@
   USE infotrac_phy
   USE YOMCST, ONLY : RG
-  USE phys_local_var_mod, ONLY : OCS_lifetime, budg_3D_ocs_to_so2, budg_ocs_to_so2 
+  USE phys_local_var_mod, ONLY : OCS_lifetime, budg_3D_ocs_to_so2, budg_ocs_to_so2
+  USE strataer_local_var_mod, ONLY : flag_min_rreduce
 
   IMPLICIT NONE
@@ -23,5 +24,6 @@
   ! local variables
   INTEGER                                       :: i,j,k,nb,ilon,ilev
-
+  REAL                                          :: rreduce
+  
 !--convert OCS to SO2
   budg_3D_ocs_to_so2(:,:)=0.0
@@ -29,17 +31,24 @@
 
   DO ilon=1, klon
-  DO ilev=1, klev
-  !only in the stratosphere
-  IF (is_strato(ilon,ilev)) THEN
-    IF (OCS_lifetime(ilon,ilev).GT.0.0) THEN
-      budg_3D_ocs_to_so2(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/OCS_lifetime(ilon,ilev)))
-    ENDIF
-    tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - budg_3D_ocs_to_so2(ilon,ilev)
-    tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*budg_3D_ocs_to_so2(ilon,ilev)
-    !convert budget from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic
-    budg_3D_ocs_to_so2(ilon,ilev)=budg_3D_ocs_to_so2(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
-    budg_ocs_to_so2(ilon)=budg_ocs_to_so2(ilon)+budg_3D_ocs_to_so2(ilon,ilev) 
-  ENDIF
-  ENDDO
+     DO ilev=1, klev
+        !only in the stratosphere
+        IF (is_strato(ilon,ilev)) THEN
+          IF (OCS_lifetime(ilon,ilev).GT.0.0.AND.OCS_lifetime(ilon,ilev).LT.1.E10) THEN
+            rreduce = OCS_lifetime(ilon,ilev)
+            ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72
+            IF(flag_min_rreduce) THEN
+               IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys
+            ENDIF
+            budg_3D_ocs_to_so2(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/rreduce))
+            tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - budg_3D_ocs_to_so2(ilon,ilev)
+            tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*budg_3D_ocs_to_so2(ilon,ilev)
+            !convert budget from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic
+            budg_3D_ocs_to_so2(ilon,ilev)=budg_3D_ocs_to_so2(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
+            budg_ocs_to_so2(ilon)=budg_ocs_to_so2(ilon)+budg_3D_ocs_to_so2(ilon,ilev)
+         ENDIF
+         ! END IF(OCS_lifetime(ilon,ilev).GT.0.0.AND.OCS_lifetime(ilon,ilev).LT.1.E10)
+      ENDIF
+      ! END IF(is_strato(ilon,ilev))
+   ENDDO
   ENDDO
 
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/so2_to_h2so4.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/so2_to_h2so4.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/so2_to_h2so4.F90	(revision 4743)
@@ -7,7 +7,9 @@
   USE aerophys
   USE infotrac_phy
-  USE YOMCST, ONLY : RG
-  USE phys_local_var_mod, ONLY : SO2_lifetime, budg_3D_so2_to_h2so4, budg_so2_to_h2so4
-
+  USE YOMCST, ONLY : RG, RD
+  ! lifetime (sec) et O3_clim (VMR)
+  USE phys_local_var_mod, ONLY : SO2_lifetime, H2SO4_lifetime, O3_clim, budg_3D_so2_to_h2so4, budg_so2_to_h2so4
+  USE strataer_local_var_mod, ONLY : flag_OH_reduced, flag_H2SO4_photolysis, flag_min_rreduce
+  
   IMPLICIT NONE
 
@@ -24,5 +26,7 @@
   ! local variables in coagulation routine
   INTEGER                                       :: i,j,k,nb,ilon,ilev
-
+  REAL                                          :: dummyso4toso2,rkho2o3,rkoho3,rkso2oh
+  REAL                                          :: rmairdens,rrak0,rrak1,rrate,rreduce
+  
 !--convert SO2 to H2SO4
   budg_3D_so2_to_h2so4(:,:)=0.0
@@ -30,18 +34,119 @@
 
   DO ilon=1, klon
-  DO ilev=1, klev
-  !only in the stratosphere
-  IF (is_strato(ilon,ilev)) THEN
-    IF (SO2_lifetime(ilon,ilev).GT.0.0) THEN
-      budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/SO2_lifetime(ilon,ilev)))
-    ENDIF
-    tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev)
-    tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev)
-    !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
-    budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
-    budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev)
-  ENDIF
-  ENDDO
-  ENDDO
-
+     DO ilev=1, klev
+     !only in the stratosphere
+        IF (is_strato(ilon,ilev)) THEN
+           IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN
+              IF (flag_OH_reduced) THEN
+                 !--convert SO2 to H2SO4 (slimane)
+                 ! Reduce OH because of SO2
+                 ! d[OH]/dt = k1[HO2][O] +k2[HO2][O3] +k3[HO2][NO]
+                 ! − (k4[O] +k5[O3] +k6[CO] +ks[SO2])[OH]
+                 ! In steady-state: d[OH]/dt = 0, so the ratio [OH]/[HO2] =
+                 ! (k1[O] +k2[O3] +k3[NO]) / (k4[O] +k5[O3] +k6[CO] +ks[SO2])
+                 ! In the lower and middle stratosphere,
+                 ! [OH]/[HO2] ~(k2[O3] +k3[NO]) / (k5[O3] +k6[CO] +ks[SO2])
+                 ! CO is only important near the tropopause. NO should not dominate over O3.
+                 ! So when consideringthe big terms: [OH]/[HO2] ~ ~ k2[O3] / (k5[O3] +ks[SO2])
+                 ! For c: control run c (no effect of SO2), [OHc]/[HO2c] ~ ~ k2[O3] / k5[O3] = 1/ Rc
+                 ! For p: perturbed run (high SO2), [OHp]/[HO2p] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) =1/Rp
+                 ! If HOx = OH + HO2 = constante,
+                 ! OHc + HO2c = OHp + HO2p
+                 ! OHc (1+Rc) = OHp (1+Rp)
+                 ! OHp = OHc (1+Rc) / (1+Rp)
+                 ! 1+Rc = (k2[O3] +k5[O3] ) / k2[O3]
+                 ! 1+Rp = (k2[O3] +k5[O3] +ks[SO2]) / k2[O3]
+                 ! (1+Rc) / (1+Rp)= (k2[O3] +k5[O3] )/ (k2[O3] +k5[O3] +ks[SO2])
+                 ! OHp = OHc * (k(HO2+OH)*[O3] +k(OH+O3)[O3] ) /
+                 ! (k(HO2+OH)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2])
+                 ! Final: OHp = OHc * (k(HO2+O3)*[O3] +k(OH+O3)[O3] ) /
+                 !  (k(HO2+O3)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2])
+                 
+                 !   latest jpl 2020
+                 rkho2o3 = 1.0e-14*exp( -490./t_seri(ilon,ilev) )
+                 rkoho3 = 1.7e-12*exp( -940./t_seri(ilon,ilev) )
+                 
+                 ! rkso2oh= so2 + oh + m  ->hso3 + m-> h2so4 + m
+                 ! air density M(molecule/cm3) = Pa/(K*1.38e-19)
+                 rmairdens =  pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19)
+                 ! JPL 2020
+                 rrak0 = 3.3e-31*( (300./t_seri(ilon,ilev))**4.3 )
+                 rrak1 = 1.6e-12
+                 rrate = (rrak0*rmairdens) / ( 1. + rrak0*rmairdens/rrak1 )
+                 rreduce = 1./( 1. + log10((rrak0*rmairdens)/rrak1)**2 )
+                 rkso2oh = rrate*(0.6**rreduce)
+                 
+                 ! O3 (molec/cm3): O3_clim in vmr
+                 rrak0 = O3_clim(ilon,ilev)*rmairdens
+                 ! SO2 (molec/cm3): convert from kg/kgA
+                 rrak1 = tr_seri(ilon,ilev,id_SO2_strat) &
+                      & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mSO2mol
+                 
+                 IF (rrak1 .GE. 0.0) THEN
+                    rreduce =( (rkho2o3+rkoho3)*rrak0 + rkso2oh*rrak1 ) / &
+                         (  (rkho2o3+rkoho3)*rrak0 )
+                    ! reduce OH, so extend SO2 lifetime
+                    rreduce = rreduce*SO2_lifetime(ilon,ilev)
+                 ELSE
+                    rreduce =  SO2_lifetime(ilon,ilev)
+                 ENDIF
+                 ! end of IF (rrak1) THEN
+              ELSE
+                 rreduce =  SO2_lifetime(ilon,ilev)
+              ENDIF
+              ! end of IF (flag_OH_reduced) THEN
+              
+              ! Check lifetime rreduce < timestep*1.5 (such as SO2 loss > 0.5*SO2) with exp(-1/1.5)=0.52
+              ! Check lifetime rreduce < timestep*3 (such as SO2 loss > 0.28*SO2) with exp(-1/3)=0.72
+              IF(flag_min_rreduce) THEN
+                 IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys
+              ENDIF
+              budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce))
+              tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev)
+              tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev)
+           ENDIF
+           ! IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN
+           
+           
+           IF (flag_H2SO4_photolysis) THEN
+              !--convert H2SO4 to SO2 using HCl photolysis (Rinsland et al, 1995)
+              ! 4.6e-8 sec-1: J strato taken from Feierabend, et al., Chem. Phys. Lett., 2006.
+              ! error before:  ...exp(-pdtphys/4.6e-8) so negligible photolysis,
+              !                 should have been exp(-pdtphys*4.6e-8)
+              ! Also: Lane and Kjaergaard: ,J. Phys. Chem. A, 2008
+              ! We find that the dominant photodissociation mechanism of H2SO4
+              ! below 70 km is absorption in the visible region by OH stretching overtone
+              ! transitions, whereas above 70 km, absorption of Lyman-R ...
+              ! In 2D model, 0.3 Rinsland et al, 1995 but Burkholder, Mills and McKeen say 0.016
+              ! JH2SO4=JHCL*0.3   (*0.016 would be too slow)
+              ! test: quick sequential integration for SO2 to H2SO4 and reverse
+              
+              IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN
+                 
+                 rreduce = H2SO4_lifetime(ilon,ilev)
+                 dummyso4toso2 = 0.0
+                 
+                 ! Check lifetime rreduce < timestep*1.5 (such as H2SO4 loss > 0.5*H2SO4) with exp(-1/1.5)=0.52
+                 ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72
+                 IF(flag_min_rreduce) THEN
+                    IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys
+                 ENDIF
+                 dummyso4toso2 = (mSO2mol/mH2SO4mol)*tr_seri(ilon,ilev,id_H2SO4_strat)*(1.0-exp(-pdtphys/rreduce))
+                 budg_3D_so2_to_h2so4(ilon,ilev) = budg_3D_so2_to_h2so4(ilon,ilev) + dummyso4toso2
+                 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + dummyso4toso2
+                 tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) - (mH2SO4mol/mSO2mol)*dummyso4toso2
+              ENDIF
+              ! IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN
+           ENDIF
+           ! end of IF (flag_H2SO4_photolysis) THEN
+           
+           !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
+           budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol* &
+                (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
+           budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev)
+        ENDIF
+        ! IF (is_strato(ilon,ilev)) THEN
+     ENDDO ! DO (klev)
+  ENDDO ! DO (klon)
+  
 END SUBROUTINE SO2_TO_H2SO4
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratH2O_methox.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratH2O_methox.F90	(revision 4743)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratH2O_methox.F90	(revision 4743)
@@ -0,0 +1,162 @@
+!
+! $Id: stratH2O_methox.F90 3677 2020-05-06 15:18:32Z oboucher $
+!
+SUBROUTINE stratH2O_methox(debutphy,paprs,dq_ch4mmr)
+!
+! output: CH4VMR in mmr (mass mixing ratio/sec: kg H2O/kg air)
+
+  USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, & 
+                      nf95_inq_varid, nf95_inquire_dimension, nf95_open
+  USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
+
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+  USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+
+  USE mod_phys_lmdz_para 
+  USE dimphy
+  USE phys_cal_mod, ONLY : mth_cur
+  USE infotrac_phy
+  USE aerophys
+  USE YOMCST
+  USE strataer_local_var_mod, ONLY : flag_newclim_file
+  
+  IMPLICIT NONE
+  
+  include "dimensions.h"
+  
+! Variable input
+  REAL paprs(klon,klev+1)
+  LOGICAL, INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
+! Variable output
+! tendance buffer pour appel de add_phys_tend
+  REAL, INTENT(INOUT), DIMENSION(klon,klev)  :: dq_ch4mmr
+  
+! Variables locales
+  INTEGER n_lat   ! number of latitudes in the input data
+  INTEGER n_lon   ! number of longitudes in the input data
+  INTEGER, SAVE :: n_lev   ! number of levels in the input data
+!$OMP THREADPRIVATE(n_lev)
+  INTEGER n_mth   ! number of months in the input data
+  INTEGER, SAVE :: mth_pre
+!$OMP THREADPRIVATE(mth_pre)
+  
+! Champs reconstitues
+  REAL paprs_glo(klon_glo,klev+1)
+  
+  REAL, ALLOCATABLE:: latitude(:)
+! (of input data sorted in strictly ascending order)
+  REAL, ALLOCATABLE:: longitude(:)
+! (of input data sorted in strictly ascending order)
+  REAL, ALLOCATABLE:: time(:)
+! (of input data sorted in strictly ascending order)
+  REAL, ALLOCATABLE:: lev(:)
+! levels of input data
+  
+!stratospheric H2O source from CH4 oxidation
+! fixed climatos
+! H2O production in VMR/sec)
+  REAL, ALLOCATABLE :: CH4RVMR_in(:, :, :, :)
+  REAL, ALLOCATABLE :: CH4RVMR_mth(:, :, :)
+  REAL, ALLOCATABLE :: CH4RVMR_tmp(:, :)
+  REAL CH4RVMR_glo(klon_glo,klev)
+!
+  INTEGER i, k, kk, j
+  
+! For NetCDF:
+  INTEGER ncid_in  ! IDs for input files
+  INTEGER varid, ncerr
+    
+  INTEGER, PARAMETER :: lev_input=17
+!--pressure at interfaces of input data (in Pa)
+  REAL, DIMENSION(lev_input+1), PARAMETER ::          & 
+                    paprs_input=(/                    &
+  1.00000002e+05,   6.06530673e+04,   3.67879449e+04, &
+  2.23130165e+04,   1.35335286e+04,   8.20850004e+03, &
+  4.97870695e+03,   3.01973841e+03,   1.83156393e+03, &
+  1.11089968e+03,   6.73794715e+02,   4.08677153e+02, &
+  2.47875223e+02,   1.50343923e+02,   9.11881985e+01, &
+  5.53084382e+01,   3.35462635e+01,   0.0           /)
+!
+  
+  IF (debutphy .OR. mth_cur .NE. mth_pre) THEN
+     
+!--preparation of global fields
+     CALL gather(paprs, paprs_glo)
+     
+     IF (is_mpi_root.AND.is_omp_root) THEN
+        
+        print *,'In stratH2O_methox read file: ch4r_annual_lmdz.nc'
+        
+        !--reading strat. H2O source (CH4 oxidation) files
+        CALL nf95_open("ch4r_annual_lmdz.nc", nf90_nowrite, ncid_in)
+        
+        CALL  nf95_inq_varid(ncid_in, "LEV", varid)
+        CALL nf95_gw_var(ncid_in, varid, lev)
+        n_lev = size(lev)
+        
+        CALL nf95_inq_varid(ncid_in, "LAT", varid)
+        CALL nf95_gw_var(ncid_in, varid, latitude)
+        n_lat = size(latitude)
+        
+        CALL nf95_inq_varid(ncid_in, "LON", varid)
+        CALL nf95_gw_var(ncid_in, varid, longitude)
+        n_lon = size(longitude)
+        
+        CALL nf95_inq_varid(ncid_in, "TIME", varid)
+        CALL nf95_gw_var(ncid_in, varid, time)
+        n_mth = size(time)
+        
+        IF (.NOT.ALLOCATED(CH4RVMR_in)) ALLOCATE(CH4RVMR_in(n_lon, n_lat, n_lev, n_mth))
+        
+        CALL nf95_inq_varid(ncid_in, "CH4RVMR", varid)
+        ncerr = nf90_get_var(ncid_in, varid, CH4RVMR_in)
+        print *,'code erreur CH4RVMR=', ncerr, varid
+        
+        CALL nf95_close(ncid_in)
+        
+        IF (.NOT.ALLOCATED(CH4RVMR_mth)) ALLOCATE(CH4RVMR_mth(n_lon, n_lat, n_lev))
+        IF (.NOT.ALLOCATED(CH4RVMR_tmp)) ALLOCATE(CH4RVMR_tmp(klon_glo, n_lev))
+        
+        !---select the correct month,
+        !---correct latitudinal order,convert input from volume mixing ratio to mass mixing ratio
+        DO j=1,n_lat
+           ! convert VMR/sec in mmr (mass mixing ratio/sec: kg H2O/kg air)
+           ! x2 because CH4->2*H2O
+           CH4RVMR_mth(:,j,:) = 2*CH4RVMR_in(:,n_lat+1-j,:,mth_cur)*mH2Omol/mAIRmol
+        ENDDO
+        
+        !---reduce to a klon_glo grid but keep the levels
+        CALL grid2dTo1d_glo(CH4RVMR_mth,CH4RVMR_tmp)
+        
+        !---regrid weighted climatologies
+        DO i=1, klon_glo
+           DO k=1, klev
+              
+              CH4RVMR_glo(i,k)=0.0
+              
+              DO kk=1, n_lev
+                 !
+                 CH4RVMR_glo(i,k)=CH4RVMR_glo(i,k)+ &
+                      MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
+                      -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
+                      *CH4RVMR_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+              ENDDO ! kk loop
+           ENDDO ! k loop
+        ENDDO ! i loop
+     ENDIF !--is_mpi_root and is_omp_root
+     
+     !--keep memory of previous month
+     mth_pre=mth_cur
+     
+     !--scatter global fields around
+     CALL scatter(CH4RVMR_glo, dq_ch4mmr)
+     
+     IF (is_mpi_root.AND.is_omp_root) THEN
+        DEALLOCATE(CH4RVMR_in,CH4RVMR_mth,CH4RVMR_tmp)
+     ENDIF !--is_mpi_root and is_omp_root
+  ENDIF ! debutphy.OR.new month
+  
+  RETURN
+  
+END SUBROUTINE stratH2O_methox
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_emiss_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_emiss_mod.F90	(revision 4743)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_emiss_mod.F90	(revision 4743)
@@ -0,0 +1,248 @@
+MODULE strataer_emiss_mod
+! This module contains information about strato microphysic model emission parameters
+  
+CONTAINS
+    
+  SUBROUTINE strataer_emiss_init()
+    
+    USE strataer_local_var_mod
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE print_control_mod, ONLY : lunout
+    USE mod_phys_lmdz_para, ONLY : is_master
+    
+    ! Local variables
+    INTEGER                  :: ispec
+    
+    WRITE(lunout,*) 'IN STRATAER_EMISS INIT WELCOME!'
+    
+    IF (flag_emit.EQ.1 .OR. flag_emit.EQ.4) THEN ! Volcano
+       CALL getin_p('nErupt',nErupt) !eruption nb
+       CALL getin_p('injdur',injdur) !injection duration
+       
+       IF (flag_emit==1) THEN
+          CALL getin_p('nAerErupt', nAerErupt) !sulfur aer nb
+       ELSEIF (flag_emit==4) THEN
+          CALL getin_p('nSpeciesErupt', nSpeciesErupt) !chimical species nb
+       ENDIF
+       
+       IF (nErupt.GT.0) THEN
+          ALLOCATE(year_emit_vol(nErupt),mth_emit_vol(nErupt),day_emit_vol(nErupt))
+          ALLOCATE(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))
+          IF (flag_emit==1) THEN
+             ALLOCATE(m_sulf_emiss_vol(nErupt))
+             ALLOCATE(m_aer_emiss_vol(nErupt,nAerErupt))
+             IF (ok_qemiss) then
+                ALLOCATE(m_H2O_emiss_vol(nErupt))
+                ALLOCATE(m_H2O_emiss_vol_daily(nErupt))
+!                ALLOCATE(d_q_emiss(klon,klev))
+                ALLOCATE(budg_emi(klon,nAerErupt+1))
+                m_H2O_emiss_vol(:)=0.
+                m_H2O_emiss_vol_daily(:)=0.
+!                d_q_emiss(:,:)=0.
+             ELSE
+                ALLOCATE(budg_emi(klon,nAerErupt))
+             ENDIF
+             budg_emi(:,:)=0.
+          ELSEIF (flag_emit==4) THEN
+             ALLOCATE(m_Chlore_emiss_vol(nErupt))
+             ALLOCATE(id_HCl) 
+             ALLOCATE(m_Brome_emiss_vol(nErupt))
+             ALLOCATE(id_HBr)
+             ALLOCATE(id_species(nSpeciesErupt))
+             ALLOCATE(m_species_emiss_vol(nErupt,nSpeciesErupt))
+             ALLOCATE(m_NOx_emiss_vol(nErupt))
+             ALLOCATE(m_H2O_emiss_vol(nErupt))
+             ALLOCATE(id_NOx)
+             ALLOCATE(id_H2O)
+             ALLOCATE(budg_emi(klon,nSpeciesErupt))
+             budg_emi(:,:)=0.
+          ENDIF
+       ELSE
+          WRITE(lunout,*) 'ERROR : Using flag_emit=1 or 4 (ie Volcanic eruption) but nErupt (',nErupt,') <=0 !'
+          CALL abort_physic('strataer_emiss_mod', &
+               'No eruption define in physiq_def (nErupt=0). Add at one eruption or use background condition.',1)
+       ENDIF ! fin if nerupt
+       
+       ! injection params (dates, loc, injections params)
+       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('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)
+       IF (flag_emit==1) THEN
+          CALL getin_p('m_sulf_emiss_vol',m_sulf_emiss_vol)
+        if (ok_qemiss) then
+         CALL getin_p('m_H2O_emiss_vol',m_H2O_emiss_vol)
+        endif
+       ELSEIF (flag_emit==4) THEN
+          CALL getin_p('id_species',id_species)
+          CALL getin_p('m_Chlore_emiss_vol',m_Chlore_emiss_vol)  
+          CALL getin_p('id_HCl',id_HCl)
+          CALL getin_p('m_Brome_emiss_vol',m_Brome_emiss_vol)
+          CALL getin_p('id_HBr',id_HBr)
+          CALL getin_p('m_NOx_emiss_vol',m_NOx_emiss_vol)
+          CALL getin_p('id_NOx',id_NOx)
+          CALL getin_p('m_H2O_emiss_vol',m_H2O_emiss_vol)
+          CALL getin_p('id_H2O',id_H2O)
+       ENDIF
+       
+    ELSEIF (flag_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('year_emit_sai_start',year_emit_sai_start)
+       CALL getin_p('year_emit_sai_end',year_emit_sai_end)
+       CALL getin_p('mth_emit_sai_start',mth_emit_sai_start)
+       CALL getin_p('mth_emit_sai_end',mth_emit_sai_end)
+       CALL getin_p('day_emit_sai_start',day_emit_sai_start)
+       CALL getin_p('day_emit_sai_end',day_emit_sai_end)
+       
+    ELSEIF (flag_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)
+    ENDIF
+    
+    IF (flag_emit == 1) THEN
+       DO ispec=1,nAerErupt
+          m_aer_emiss_vol(:,ispec) = m_sulf_emiss_vol(:)
+       ENDDO
+    ELSEIF (flag_emit== 4) THEN
+       DO ispec=1,nSpeciesErupt
+          IF(id_species(ispec) == id_HCl) THEN
+             m_species_emiss_vol(:,ispec) = m_Chlore_emiss_vol(:)
+          ENDIF
+          IF (id_species(ispec) == id_HBr) THEN 
+             m_species_emiss_vol(:,ispec) = m_Brome_emiss_vol(:)
+          ENDIF
+          IF (id_species(ispec) == id_NOx) THEN
+             m_species_emiss_vol(:,ispec) = m_NOx_emiss_vol(:)
+          ENDIF
+          IF (id_species(ispec) == id_H2O) THEN
+             m_species_emiss_vol(:,ispec) = m_H2O_emiss_vol(:)
+          ENDIF
+       ENDDO
+    ENDIF
+    
+    !============= Injection ponderation =============
+    IF (flag_emit > 0) THEN
+       CALL strataer_ponde_init
+       WRITE(lunout,*) 'IN STRATAER INIT : ponde_lonlat_vol',ponde_lonlat_vol
+    ENDIF
+    
+    !============= Print params =============
+    IF (is_master) THEN
+       IF (flag_emit == 1 .OR. flag_emit == 4) THEN
+          WRITE(lunout,*) 'IN STRATAER nErupt: ',nErupt
+          WRITE(lunout,*) 'IN STRATAER injdur: ',injdur
+          WRITE(lunout,*) 'IN STRATAER nAerErupt: ',nAerErupt
+          
+          WRITE(lunout,*) 'IN STRATAER : year_emit_vol',year_emit_vol
+          WRITE(lunout,*) 'IN STRATAER : mth_emit_vol',mth_emit_vol
+          WRITE(lunout,*) 'IN STRATAER : day_emit_vol',day_emit_vol
+          WRITE(lunout,*) 'IN STRATAER : altemiss_vol',altemiss_vol
+          WRITE(lunout,*) 'IN STRATAER : sigma_alt_vol',sigma_alt_vol
+          WRITE(lunout,*) 'IN STRATAER : xlat_min_vol',xlat_min_vol
+          WRITE(lunout,*) 'IN STRATAER : xlat_max_vol',xlat_max_vol
+          WRITE(lunout,*) 'IN STRATAER : xlon_min_vol',xlon_min_vol
+          WRITE(lunout,*) 'IN STRATAER : xlon_max_vol',xlon_max_vol
+          IF (flag_emit==1) THEN
+             WRITE(lunout,*) 'IN STRATAEREMISS : m_sulf_emiss_vol',m_sulf_emiss_vol
+             WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_vol',m_aer_emiss_vol
+             IF (ok_qemiss) then
+                WRITE(lunout,*) 'IN STRATAEREMISS : m_H2O_emiss_vol',m_H2O_emiss_vol
+             ENDIF
+          ENDIF
+       ELSEIF (flag_emit == 2) THEN
+          WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
+          WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
+          WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
+          WRITE(lunout,*) 'IN STRATAER : xlat_sai',xlat_sai
+          WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
+       ELSEIF (flag_emit == 3) THEN
+          WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
+          WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
+          WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
+          WRITE(lunout,*) 'IN STRATAER : year_emit_sai start/end',year_emit_sai_start,year_emit_sai_end
+          WRITE(lunout,*) 'IN STRATAER : mth_emit_sai start/end',mth_emit_sai_start,mth_emit_sai_end
+          WRITE(lunout,*) 'IN STRATAER : day_emit_sai start/end',day_emit_sai_start,day_emit_sai_end
+          WRITE(lunout,*) 'IN STRATAER : xlat_min_sai',xlat_min_sai
+          WRITE(lunout,*) 'IN STRATAER : xlat_max_sai',xlat_max_sai
+          WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
+       ENDIF
+       IF(flag_emit == 4) THEN
+          WRITE(lunout,*) 'IN STRATAER : nSpeciesErupt: ',nSpeciesErupt
+          WRITE(lunout,*) 'IN STRATAER : id_species = ',id_species
+          WRITE(lunout,*) 'IN STRATAER : id_HCl = ',id_HCl
+          WRITE(lunout,*) 'IN STRATAER : id_HBr = ',id_HBr
+          WRITE(lunout,*) 'IN STRATAER : id_NOx = ',id_NOx
+          WRITE(lunout,*) 'IN STRATAER : id_H2O = ',id_H2O
+          WRITE(lunout,*) 'IN STRATAER : m_Chlore_emiss_vol = ',m_Chlore_emiss_vol
+          WRITE(lunout,*) 'IN STRATAER : m_Brome_emiss_vol = ',m_Brome_emiss_vol
+          WRITE(lunout,*) 'IN STRATAER : m_NOx_emiss_vol = ',m_NOx_emiss_vol
+          WRITE(lunout,*) 'IN STRATAER : m_H2O_emiss_vol = ',m_H2O_emiss_vol
+       ENDIF
+    ENDIF ! if master
+    
+    WRITE(lunout,*) 'IN STRATAER_EMISS END'
+  END SUBROUTINE strataer_emiss_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 strataer_local_var_mod
+    
+    INCLUDE "YOMCST.h"  !--RPI
+    
+    ! local var
+    REAL                :: lat_reg_deg,lon_reg_deg ! lat and lon of grid points 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_PONDE_INIT ponde_lonlat: ', ponde_lonlat_vol
+    
+  END SUBROUTINE strataer_ponde_init
+  
+END MODULE strataer_emiss_mod
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_local_var_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 4743)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 4743)
@@ -0,0 +1,245 @@
+! $Id: strataer_local_var_mod.F90 2022-10-17 19:39:09Z nlebas $
+MODULE strataer_local_var_mod
+! This module contains strato microphysic model parameters & variables
+  
+  IMPLICIT NONE
+  
+  !============= GENERAL PARAMS =============
+  !flag for type emission scenario: (0) background aer ; (1) volcanic eruption with Sulfur ;
+  !(2,3) strato aer injections (SAI) ; (4) volcanic eruption chimistry; (5) rocket
+  INTEGER,SAVE :: flag_emit
+  !$OMP THREADPRIVATE(flag_emit)
+  
+  ! flag for emission altitude distribution: (0) gaussian; (1) uniform
+  INTEGER,SAVE :: flag_emit_distrib
+  !$OMP THREADPRIVATE(flag_emit_distrib)
+  
+  ! flag to choose nucleation nucleation method
+  LOGICAL,SAVE :: flag_new_nucl   ! T=new routine from A. Maattanen (LATMOS), F=older routine from H. Vehkamäki (FMI)
+  !$OMP THREADPRIVATE(flag_new_nucl)
+  
+  ! Use relative humidity from 2D model stratospheric H2O because LMDz is too dry in the stratosphere
+  ! (no CH4 oxidation)
+  LOGICAL,SAVE :: flag_H2O2d_nucleation
+  !$OMP THREADPRIVATE(flag_H2O2d_nucleation)
+  
+  ! OH reduction from SO2. OH is reduced when its reaction with SO2 competes sufficiently with its reaction
+  ! with O3 (Bekki, 1995). As a result, SO2 lifetime is extended. 2D model O3 climatologies are needed.
+  LOGICAL,SAVE :: flag_OH_reduced
+  !$OMP THREADPRIVATE(flag_OH_reduced)
+  
+  ! H2SO4 photolysis: H2SO4 is converted into SO2 by complex photolytic processes. Here simplified approach
+  ! by setting H2SO4 cross-sections = 0.3*HCl cross-sections (Rinsland et al., 1995)
+  LOGICAL,SAVE :: flag_H2SO4_photolysis
+  !$OMP THREADPRIVATE(flag_H2SO4_photolysis)
+  
+  ! flag for minimum lifetime (=1.5 pdt phys)
+  LOGICAL,SAVE :: flag_min_rreduce
+  !$OMP THREADPRIVATE(flag_min_rreduce)
+  
+  ! flag to read new climato (O3, H2O & H2SO4_LIFET)
+  LOGICAL,SAVE :: flag_newclim_file
+  !$OMP THREADPRIVATE(flag_newclim_file)
+  
+  ! Verbose mode to get more print info
+  LOGICAL, SAVE :: flag_verbose_strataer
+  !$OMP THREADPRIVATE(flag_verbose_strataer)
+  
+  
+  !============= NUCLEATION VARS =============
+  ! 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)
+
+  LOGICAL,SAVE :: ok_qemiss
+  !$OMP THREADPRIVATE(ok_qemiss)
+  INTEGER,SAVE :: flh2o  ! ds stratemit : flh2o =0 (tr_seri), flh2o=1 (dq) 
+  !$OMP THREADPRIVATE(flh2o)
+!  REAL,ALLOCATABLE,SAVE    :: d_q_emiss(:,:)
+!  !$OMP THREADPRIVATE(d_q_emiss)
+  
+  REAL,ALLOCATABLE,SAVE    :: budg_emi(:,:)            !DIMENSION(klon,n)
+  !$OMP THREADPRIVATE(budg_emi)
+  
+  
+  !============= EMISSION VARS =============
+  !--flag_emit=1 OR == 4 -- 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    :: altemiss_vol(:)           ! emission altitude in m
+  REAL,ALLOCATABLE,SAVE    :: sigma_alt_vol(:)          ! standard deviation of emission altitude in m
+  !$OMP THREADPRIVATE(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_emit=1
+  INTEGER,SAVE             :: nAerErupt                 ! number Aerosol 
+  !$OMP THREADPRIVATE(nAerErupt)
+  REAL,ALLOCATABLE,SAVE    :: m_sulf_emiss_vol(:)        ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
+  REAL,ALLOCATABLE,SAVE    :: m_aer_emiss_vol(:,:)
+  !$OMP THREADPRIVATE(m_sulf_emiss_vol,m_aer_emiss_vol)
+  
+  !--flag_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)
+  INTEGER,SAVE    :: year_emit_sai_start
+  INTEGER,SAVE    :: year_emit_sai_end
+  INTEGER,SAVE    :: mth_emit_sai_start
+  INTEGER,SAVE    :: mth_emit_sai_end
+  INTEGER,SAVE    :: day_emit_sai_start
+  INTEGER,SAVE    :: day_emit_sai_end
+  !$OMP THREADPRIVATE(year_emit_sai_start, year_emit_sai_end)
+  !$OMP THREADPRIVATE(mth_emit_sai_start, mth_emit_sai_end)
+  !$OMP THREADPRIVATE(day_emit_sai_start, day_emit_sai_end)
+  REAL,SAVE    :: xlat_sai               ! latitude of SAI in degree
+  REAL,SAVE    :: xlon_sai               ! longitude of SAI in degree
+  !$OMP THREADPRIVATE(xlat_sai, xlon_sai)
+  
+  !--flag_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)
+  
+  !--flag_emit=4 -- volc species
+  INTEGER,SAVE             :: nSpeciesErupt            ! number of species Repr
+  INTEGER,ALLOCATABLE,SAVE :: id_species(:)            ! indice species Repr
+  REAL,ALLOCATABLE,SAVE    :: m_species_emiss_vol(:,:) ! emitted species 
+  !$OMP THREADPRIVATE(nSpeciesErupt,id_species,m_species_emiss_vol)
+  INTEGER,ALLOCATABLE,SAVE :: id_HCl
+  INTEGER,ALLOCATABLE,SAVE :: id_HBr
+  INTEGER,ALLOCATABLE,SAVE :: id_NOx
+  INTEGER,ALLOCATABLE,SAVE :: id_H2O
+  !$OMP THREADPRIVATE(id_HCl,id_HBr,id_NOx,id_H2O)
+  REAL,ALLOCATABLE,SAVE    :: m_Chlore_emiss_vol(:)   ! emitted Chlore mass 
+  REAL,ALLOCATABLE,SAVE    :: m_Brome_emiss_vol(:)    ! emitted Brome mass
+  REAL,ALLOCATABLE,SAVE    :: m_NOx_emiss_vol(:)      ! emitted NOx mass
+  REAL,ALLOCATABLE,SAVE    :: m_H2O_emiss_vol(:)      ! emitted H2O mass
+  REAL,ALLOCATABLE,SAVE    :: m_H2O_emiss_vol_daily(:)
+  !$OMP THREADPRIVATE(m_Chlore_emiss_vol,m_Brome_emiss_vol,m_NOx_emiss_vol,m_H2O_emiss_vol)
+  !$OMP THREADPRIVATE(m_H2O_emiss_vol_daily)
+  
+  !--flag_emit=5 -- Rockets Emitted
+  INTEGER, SAVE            :: ifreqroc        ! frequence (=2 ex: tous les 2 mois)
+  !$OMP THREADPRIVATE(ifreqroc)
+  INTEGER,ALLOCATABLE,SAVE :: day_emit_roc(:) ! day of emission date
+  !$OMP THREADPRIVATE(day_emit_roc)
+  
+  REAL,SAVE    :: dlat, dlon             ! delta latitude and d longitude of grid in degree
+  !$OMP THREADPRIVATE(dlat, dlon)
+  
+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
+    USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas
+    
+    WRITE(lunout,*) 'IN STRATAER_LOCAL_VAR INIT WELCOME!'
+    
+    !============= Check Sulfur aerosols ID =============
+    WRITE(lunout,*) 'STRATAER_LOCAL_VAR INIT: id_OCS_strat=',id_OCS_strat,' id_SO2_strat=',id_SO2_strat,' id_H2SO4_strat=',id_H2SO4_strat
+    
+    IF(id_OCS_strat < 0 .OR. id_OCS_strat > nbtr_sulgas) THEN
+       WRITE(lunout,*) 'ERROR : OCS index id_OCS_strat=',id_OCS_strat,' is negative or superior than the total sulfur gases !'
+       CALL abort_physic('strataer_local_var_mod','Wrong OCS index, check your tracer.def file.',1)
+    ELSEIF(id_SO2_strat < 0 .OR. id_SO2_strat > nbtr_sulgas) THEN
+       WRITE(lunout,*) 'ERROR : SO2 index id_SO2_strat=',id_SO2_strat,' is negative or superior than the total sulfur gases !'
+       CALL abort_physic('strataer_local_var_mod','Wrong SO2 index, check your tracer.def file.',1)
+    ELSEIF(id_H2SO4_strat < 0 .OR. id_H2SO4_strat > nbtr_sulgas) THEN
+       WRITE(lunout,*) 'ERROR : H2SO4 index id_H2SO4_strat=',id_H2SO4_strat,' is negative or superior than the total sulfur gases !'
+       CALL abort_physic('strataer_local_var_mod','Wrong H2SO4 index, check your tracer.def file.',1)
+    ENDIF
+    
+    !============= Init params =============
+    flag_emit = 0                   ! Background (default)
+    flag_emit_distrib = 0           ! Gaussian (default)
+    flag_new_nucl = .TRUE.          ! Define nucleation routine (default: A. Maattanen - LATMOS)
+    flag_verbose_strataer = .FALSE. ! verbose mode
+    flag_newclim_file = .TRUE.      ! Define input climato file (default: all climato)
+    flag_H2O2d_nucleation = .FALSE. ! Use H2O 2D climato (default: No)
+    flag_OH_reduced = .FALSE.       ! OH reduce (default: No)
+    flag_H2SO4_photolysis = .FALSE. ! H2SO4 photolysis (default: No)
+    flag_min_rreduce = .TRUE.       ! Minimum lifetime=1.5 pdt phys (default: Yes)
+    ok_qemiss = .FALSE.             ! H2O emission flag
+    
+    ! nuc init
+    flag_nuc_rate_box = .FALSE.
+    nuclat_min=0  ; nuclat_max=0
+    nucpres_min=0 ; nucpres_max=0
+    
+    ! emiss init
+    nErupt = 0 ! eruption number
+    injdur = 0 ! init injection duration
+    nAerErupt = 1 ; nSpeciesErupt = 1 ; id_species = 3
+    year_emit_vol=0 ; mth_emit_vol=0 ; day_emit_vol=0
+    m_species_emiss_vol=0. ; m_aer_emiss_vol=0. ; m_sulf_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.
+    m_Chlore_emiss_vol=0. ; m_Brome_emiss_vol=0.
+    m_NOx_emiss_vol=0. ; m_H2O_emiss_vol=0.
+    day_emit_roc=0 ; ifreqroc=2 ; flh2o=0
+    
+    !============= Read params =============
+    CALL getin_p('flag_emit',flag_emit)
+    CALL getin_p('flag_emit_distrib',flag_emit_distrib)
+    CALL getin_p('flag_verbose_strataer',flag_verbose_strataer)
+    CALL getin_p('flag_new_nucl',flag_new_nucl)
+    CALL getin_p('flag_newclim_file',flag_newclim_file)
+    CALL getin_p('flag_H2O2d_nucleation',flag_H2O2d_nucleation)
+    CALL getin_p('flag_OH_reduced',flag_OH_reduced)
+    CALL getin_p('flag_H2SO4_photolysis',flag_H2SO4_photolysis)
+    CALL getin_p('flag_min_rreduce',flag_min_rreduce)
+    CALL getin_p('ok_qemiss',ok_qemiss)
+    
+    !============= Test flag coherence =============
+    IF (.NOT. flag_newclim_file) THEN
+       IF (flag_H2SO4_photolysis .OR. flag_OH_reduced .OR. flag_H2O2d_nucleation) THEN
+          WRITE(lunout,*) 'ERROR : flag_newclim_file=',flag_newclim_file, &
+               ' whereas flag_H2SO4_photolysis=',flag_H2SO4_photolysis,', flag_OH_reduced=',flag_OH_reduced, &
+               ' and flag_H2O2d_nucleation=',flag_H2O2d_nucleation
+          CALL abort_physic('strataer_local_var_mod','Incompatible options in physiq_def file !',1)
+       ENDIF
+       IF(flag_min_rreduce) THEN
+          WRITE(lunout,*) 'Warning : flag_min_rreduce will be ignored with old climato file !'
+       ENDIF
+    ENDIF
+    
+    !============= Print params =============
+    IF (is_master) THEN
+       WRITE(lunout,*) 'flag_emit = ',flag_emit
+       WRITE(lunout,*) 'IN STRATAER : flag_new_nucl = ',flag_new_nucl
+       WRITE(lunout,*) 'IN STRATAER : flag_newclim_file = ',flag_newclim_file
+       WRITE(lunout,*) 'IN STRATAER : flag_emit_distrib = ',flag_emit_distrib
+       WRITE(lunout,*) 'IN STRATAER : flag_verbose_strataer = ',flag_verbose_strataer
+       IF (flag_emit == 1 .OR. flag_emit == 4) THEN
+          WRITE(lunout,*) 'IN STRATAER : flag_H2O2d_nucleation = ',flag_H2O2d_nucleation
+          WRITE(lunout,*) 'IN STRATAER : flag_OH_reduced = ',flag_OH_reduced
+          WRITE(lunout,*) 'IN STRATAER : flag_H2SO4_photolysis = ',flag_H2SO4_photolysis
+          WRITE(lunout,*) 'IN STRATAER : flag_min_rreduce = ',flag_min_rreduce
+          WRITE(lunout,*) 'IN STRATAER : ok_qemiss = ',ok_qemiss
+       ENDIF
+    ENDIF ! if master
+    
+    WRITE(lunout,*) 'IN STRATAER INIT END'
+    
+  END SUBROUTINE strataer_init
+  
+END MODULE strataer_local_var_mod
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_mod.F90	(revision 4584)
+++ 	(revision )
@@ -1,240 +1,0 @@
-! $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 ! Background (default)
-    flag_sulf_emit_distrib = 0 ! Gaussian (default)
-    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)
-       CALL getin_p('flag_sulf_emit_distrib',flag_sulf_emit_distrib)
-
-       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.
-
-       IF (nErupt.GT.0) THEN
-          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))
-       ELSE
-          WRITE(lunout,*) 'ERROR : Using flag_sulf_emit=1 (ie Volcanic eruption) but nErupt (',nErupt,') <=0 !'
-          CALL abort_physic('strataer_mod','No eruption define in physiq_def (nErupt=0). Add one or use background condition.',1)
-       ENDIF
-
-       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)
-    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)
-    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)
-    ENDIF
-
-    !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.
-    nuclat_min=0  ; nuclat_max=0
-    nucpres_min=0 ; nucpres_max=0
-
-    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)
-
-    !IF (is_master) THEN
-    WRITE(lunout,*) 'flag_sulf_emit = ',flag_sulf_emit
-    IF (flag_sulf_emit == 1) THEN
-       WRITE(lunout,*) 'IN STRATAER nErupt: ',nErupt
-       WRITE(lunout,*) 'IN STRATAER injdur: ',injdur
-       WRITE(lunout,*) 'IN STRATAER : year_emit_vol',year_emit_vol
-       WRITE(lunout,*) 'IN STRATAER : mth_emit_vol',mth_emit_vol
-       WRITE(lunout,*) 'IN STRATAER : day_emit_vol',day_emit_vol
-       WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_vol',m_aer_emiss_vol
-       WRITE(lunout,*) 'IN STRATAER : altemiss_vol',altemiss_vol
-       WRITE(lunout,*) 'IN STRATAER : sigma_alt_vol',sigma_alt_vol
-       WRITE(lunout,*) 'IN STRATAER : xlat_min_vol',xlat_min_vol
-       WRITE(lunout,*) 'IN STRATAER : xlat_max_vol',xlat_max_vol
-       WRITE(lunout,*) 'IN STRATAER : xlon_min_vol',xlon_min_vol
-       WRITE(lunout,*) 'IN STRATAER : xlon_max_vol',xlon_max_vol
-    ELSEIF (flag_sulf_emit == 2) THEN
-       WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
-       WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
-       WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
-       WRITE(lunout,*) 'IN STRATAER : xlat_sai',xlat_sai
-       WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
-    ELSEIF (flag_sulf_emit == 3) THEN
-       WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
-       WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
-       WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
-       WRITE(lunout,*) 'IN STRATAER : xlat_min_sai',xlat_min_sai
-       WRITE(lunout,*) 'IN STRATAER : xlat_max_sai',xlat_max_sai
-       WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
-    ENDIF
-    WRITE(lunout,*) 'flag_sulf_emit_distrib = ',flag_sulf_emit_distrib
-    WRITE(lunout,*) 'IN STRATAER : flag_nuc_rate_box = ',flag_nuc_rate_box
-    IF (flag_nuc_rate_box) THEN
-       WRITE(lunout,*) 'IN STRATAER : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max
-       WRITE(lunout,*) 'IN STRATAER : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max
-    ENDIF
-    !ENDIF
-
-    IF (flag_sulf_emit > 0) THEN
-       CALL strataer_ponde_init
-       WRITE(lunout,*) 'IN STRATAER INT2 END'
-    END IF
-    WRITE(lunout,*) 'IN STRATAER 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
-
-    INCLUDE "YOMCST.h"  !--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/Portage_acc/libf/phylmd/StratAer/strataer_nuc_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_nuc_mod.F90	(revision 4743)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/strataer_nuc_mod.F90	(revision 4743)
@@ -0,0 +1,73 @@
+! $Id: strataer_nuc_mod.F90 3930 2021-06-11 19:39:09Z oboucher $
+MODULE strataer_nuc_mod
+! This module contains information about strato microphysic model nucleation parameters
+  
+  IMPLICIT NONE
+  
+CONTAINS
+  
+  ! Init all nucleation params
+  SUBROUTINE strataer_nuc_init()
+    
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE print_control_mod, ONLY : lunout
+    USE mod_phys_lmdz_para, ONLY : is_master
+    USE strataer_local_var_mod, ONLY: flag_nuc_rate_box,nuclat_min,nuclat_max,nucpres_min,nucpres_max
+    
+    !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
+    !
+    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)
+    
+    !============= Print params =============
+    IF (is_master) THEN
+       WRITE(lunout,*) 'IN STRATAER_NUC : flag_nuc_rate_box = ',flag_nuc_rate_box
+       IF (flag_nuc_rate_box) THEN
+          WRITE(lunout,*) 'IN STRATAER_NUC : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max
+          WRITE(lunout,*) 'IN STRATAER_NUC : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max
+       ENDIF
+    ENDIF ! if master
+    
+  END SUBROUTINE strataer_nuc_init
+  
+  ! Init aerosol tracers and large scale scavinging
+  SUBROUTINE tracstrataer_init(aerosol,lessivage)
+    
+    USE infotrac_phy, ONLY: nbtr, nbtr_sulgas, id_H2SO4_strat
+    USE ioipsl, ONLY : getin
+    USE print_control_mod, ONLY : lunout
+    
+    ! Output variables
+    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: aerosol
+    LOGICAL,INTENT(INOUT)                  :: lessivage
+    INTEGER :: it
+    
+    ! Initialization
+    lessivage  =.TRUE.
+    aerosol(:) = .FALSE.
+    
+    DO it= 1, nbtr_sulgas
+       aerosol(it)=.FALSE.
+       IF (it==id_H2SO4_strat) aerosol(it)=.TRUE.
+    ENDDO
+    DO it= nbtr_sulgas+1, nbtr
+       aerosol(it)=.TRUE.
+    ENDDO
+    
+    !!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
+!    CALL getin('lessivage',lessivage)
+    WRITE(lunout,*) 'IN TRACSTRATAER_INIT lessivage: ',lessivage
+    
+  END SUBROUTINE tracstrataer_init
+  
+END MODULE strataer_nuc_mod
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratdistrib.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratdistrib.F90	(revision 4743)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratdistrib.F90	(revision 4743)
@@ -0,0 +1,68 @@
+!
+! $Id: stratemit.F90  2022-07-04 mmarchand $
+!
+SUBROUTINE STRATDISTRIB(altLMDz,altemiss,sigma_alt,f_lay_emiss)
+!
+  USE dimphy, ONLY : klon,klev
+  USE strataer_local_var_mod
+  USE YOMCST, only : RPI
+
+  IMPLICIT NONE
+
+! Input argument
+!---------------
+    REAL,DIMENSION(klev+1)  :: altLMDz    ! altitude of layer interfaces in m
+    REAL                    :: sigma_alt
+    REAL                    :: altemiss
+
+! Local variables
+!----------------
+    REAL                 :: f_lay_sum     ! sum of layer emission fractions
+    REAL                 :: alt
+    INTEGER              :: k,i_int
+    INTEGER,PARAMETER    :: n_int_alt=10  ! subintervals nb for integration over Gaussian emission
+
+! Output argument
+!----------------
+    REAL,DIMENSION(klev) :: f_lay_emiss   ! emiss fraction for every emiss for every vertical layer
+    
+    
+! Inject distrib
+!----------------
+    SELECT CASE(flag_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)* &
+                   & exp(-0.5*((alt-altemiss)/sigma_alt)**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 is considered
+             ! to be half the
+             ! height of the injection, centered around altemiss)
+       f_lay_sum=0.0
+       DO k=1, klev
+          f_lay_emiss(k)=max(min(altemiss+sigma_alt,altLMDz(k+1))- &
+                      & max(altemiss-sigma_alt, &
+                      & altLMDz(k)),0.)/(2.*sigma_alt)
+          f_lay_sum=f_lay_sum+f_lay_emiss(k)
+       ENDDO
+
+       END SELECT        ! End CASE over flag_emit_distrib)
+       
+       WRITE(*,*) "IN stratdistrib f_lay_emiss=",f_lay_emiss
+
+       !correct for step integration error
+       f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
+
+END SUBROUTINE STRATDISTRIB
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratemit.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratemit.F90	(revision 4743)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/stratemit.F90	(revision 4743)
@@ -0,0 +1,162 @@
+!
+! $Id: stratemit.F90  2022-07-04 mmarchand $
+!
+SUBROUTINE STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
+                    m_emiss_vol_daily,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_spec,&
+                    stretchlong,ispeci,id_species_total) 
+
+  USE dimphy, ONLY : klon,klev
+  USE strataer_local_var_mod
+  USE phys_cal_mod
+  USE phys_local_var_mod, ONLY: d_q_emiss
+  USE YOMCST, only : RD, RPI, RG
+  USE geometry_mod, ONLY : cell_area, boundslat
+  USE aerophys
+  USE infotrac_phy
+  USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon
+ 
+  IMPLICIT NONE
+
+! Input argument
+!---------------
+    REAL,INTENT(IN)    :: pdtphys, pdt               ! Pas d'integration pour la physique(seconde)
+    REAL,DIMENSION(klon),INTENT(IN)    :: xlat       ! latitudes pour chaque point
+    REAL,DIMENSION(klon),INTENT(IN)    :: xlon       ! longitudes pour chaque point
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay
+    ! pression pour le mileu de chaque couche (en Pa)
+    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs
+    ! pression pour chaque inter-couche (en Pa)
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri ! Temperature
+    REAL,INTENT(IN)    :: latmin,latmax,lonmin,lonmax
+    REAL,INTENT(IN)    :: m_emiss_vol_daily          ! daily injection mass emission
+    REAL,INTENT(IN)    :: stretchlong                ! length to stretch emission (day)
+    REAL,DIMENSION(klon,klev) :: m_air_gridbox       ! mass ofair in every grid box[kg
+    REAL,INTENT(IN)     :: sigma_alt, altemiss
+    INTEGER,INTENT(IN)  :: id_spec, ispeci, id_species_total
+    
+! Output argument
+!----------------
+    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
+! Local variables
+!----------------
+    INTEGER              :: ieru, i, k, i_int, ilon, ilev
+    REAL                 :: dlat_loc
+    REAL,DIMENSION(klev+1)  :: altLMDz    ! altitude of layer interfaces in m
+    REAL                 :: zrho          ! Density of air [kg/m3]
+    REAL,DIMENSION(klev) :: zdm           ! mass of atm. model layer in kg  
+    REAL                 :: zdz           ! thickness of atm. model layer in m
+    REAL                 :: f_lay_sum     ! sum of layer emission fractions
+    REAL,DIMENSION(klev) :: f_lay_emiss   ! emiss fraction for every emiss for every vertical layer
+    REAL                 :: emission      ! emission
+    REAL                 :: theta_min, theta_max, theta ! for SAI computation between two latitudes
+    
+    IF (is_mpi_root) THEN
+       WRITE(*,*) 'IN STRATEMIT: date from phys_cal_mod=',year_cur,'-',&
+                      &  mth_cur,'-',day_cur,'-',hour,' flh2o=',flh2o
+    ENDIF
+
+!--calculate mass of air in every grid box !NL: DONE ABOVE DIRECTLY IN MAIN LOOP I,K
+!    DO ilon=1, klon
+!       DO ilev=1, klev
+!          m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
+!       ENDDO
+!    ENDDO
+    
+    DO i=1,klon
+       
+       dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3))
+       dlon=360./2./FLOAT(nbp_lon) 
+       ! dlat = half difference of boundary latitudes
+       IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT, dlat=',dlat_loc,xlat(i),xlon(i)
+       
+       theta_min = max(xlat(i)-dlat_loc,latmin)
+       theta_max = min(xlat(i)+dlat_loc,latmax)
+                
+       IF (  xlat(i).GE.latmin-dlat_loc .AND. &
+         &   xlat(i).LT.latmax+dlat_loc .AND. &
+         &   xlon(i).GE.lonmin-dlon .AND. &
+         &   xlon(i).LT.lonmax+dlon ) THEN
+                   !
+          WRITE(*,*) 'coordinates of volcanic injection point=',&
+            &    xlat(i),xlon(i),day_cur,mth_cur,year_cur
+          WRITE(*,*) 'DD m_emiss_vol_daily=', &
+            &    m_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
+          
+          CALL STRATDISTRIB(altLMDz,altemiss,sigma_alt,f_lay_emiss)
+          
+          IF (flag_emit.EQ.3) then
+             theta=(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
+               & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
+          ELSE
+             theta=1.
+          ENDIF
+
+          !vertically distributed emission
+          DO k=1, klev
+            !--calculate mass of air in every grid box
+            m_air_gridbox(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i)
+            ! stretch emission over stretchlong period
+            emission=m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k)/stretchlong/ &
+                 &     (86400.-pdt)*theta
+            
+             IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: emission avant/apres', &
+              &  'alt= ', altLMDz(k), &
+              &  'flh2o= ',flh2o, &
+              &  'id_speac= ',id_spec,id_species_total, &
+              &  'emission= ',emission, &
+              &  'pdtphys= ',pdtphys, &
+              &  'rapport m_emiss/m_air*f_lay= ', m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k), &
+              &   'stretchlong= ', stretchlong, &
+              &   'theta= ', theta
+             
+             IF(emission < 1.E-34) emission = 0.0
+             
+             IF (flh2o.EQ.0) THEN
+                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: tr_ser avant/apres',&
+                  &  'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
+                  &  tr_seri(i,k,id_spec), &
+                  &  tr_seri(i,k,id_spec)+emission*pdtphys
+              
+                tr_seri(i,k,id_spec)=tr_seri(i,k,id_spec)+emission*pdtphys
+                IF (id_species_total.NE.0) THEN
+                   tr_seri(i,k,id_species_total)=tr_seri(i,k,id_species_total)+emission*pdtphys
+                ENDIF
+             ELSE IF(flh2o.EQ.1) THEN
+                d_q_emiss(i,k)=emission*pdtphys
+                IF(d_q_emiss(i,k) > 1.E34) THEN
+                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 = ',d_q_emiss(i,k)
+                ELSE IF(d_q_emiss(i,k) < -1.E34) THEN
+                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 = ',d_q_emiss(i,k)
+                ENDIF
+                
+                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: ',&
+                     &  'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
+                     &  'emission= ',emission, &
+                     &  'd_q_emiss(i,k)= ',d_q_emiss(i,k)
+                
+                IF(d_q_emiss(i,k) > 1.E34) THEN
+                   d_q_emiss(i,k) = 1.E34
+                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 : ',d_q_emiss(i,k)
+                ELSE IF(d_q_emiss(i,k) < -1.E34) THEN
+                   d_q_emiss(i,k) = -1.E34
+                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 : ',d_q_emiss(i,k)
+                ENDIF
+             ENDIF
+             budg_emi(i,ispeci)=budg_emi(i,ispeci)+emission*zdm(k)
+          ENDDO
+          
+       ENDIF ! emission grid cell
+
+    ENDDO ! klon loop
+
+END SUBROUTINE STRATEMIT
Index: LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/traccoag_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/traccoag_mod.F90	(revision 4584)
+++ LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/traccoag_mod.F90	(revision 4743)
@@ -23,11 +23,11 @@
     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     USE mod_phys_lmdz_para, only: gather, scatter
-    USE phys_cal_mod, ONLY : year_len, mth_len, year_cur, mth_cur, day_cur, hour
+    USE phys_cal_mod, ONLY : year_len, year_cur, mth_cur, day_cur, hour
     USE sulfate_aer_mod
     USE phys_local_var_mod, ONLY: stratomask
     USE YOMCST
     USE print_control_mod, ONLY: lunout
-    USE strataer_mod
-
+    USE strataer_local_var_mod
+    
     IMPLICIT NONE
 
@@ -58,4 +58,5 @@
 !----------------
     REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
+    REAL                                   :: m_aer               ! aerosol mass
     INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
     LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
@@ -78,10 +79,14 @@
     REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
     REAL                                   :: dlat_loc
+    REAL                                   :: latmin,latmax,lonmin,lonmax ! lat/lon min/max for injection
+    REAL                                   :: sigma_alt, altemiss ! injection altitude + sigma for distrib
+    REAL                                   :: pdt,stretchlong     ! physic timestep, stretch emission over one day
+    
     INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
     INTEGER                                :: yr, is_bissext
 
-    IF (is_mpi_root) THEN
+    IF (is_mpi_root .AND. flag_verbose_strataer) THEN
        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
+       WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit
     ENDIF
     
@@ -126,20 +131,11 @@
 !--calculate mass of air in every grid box
     DO ilon=1, klon
-    DO ilev=1, klev
-      m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
+       DO ilev=1, klev
+          m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
+       ENDDO
     ENDDO
-    ENDDO
-
-!    IF (debutphy) THEN
-!      CALL gather(tr_seri, tr_seri_glo)
-!      IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
-!--initialising tracer concentrations to zero
-!        DO it=1, nbtr
-!        tr_seri(:,:,it)=0.0
-!        ENDDO
-!      ENDIF
-!    ENDIF
-
+    
 !--initialise emission diagnostics
+    budg_emi(:,1)=0.0
     budg_emi_ocs(:)=0.0
     budg_emi_so2(:)=0.0
@@ -147,6 +143,6 @@
     budg_emi_part(:)=0.0
 
-!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
-    SELECT CASE(flag_sulf_emit)
+!--sulfur emission, depending on chosen scenario (flag_emit)
+    SELECT CASE(flag_emit)
 
     CASE(0) ! background aerosol
@@ -159,71 +155,29 @@
           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), &
+
+             ! daily injection mass emission
+             m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
+             !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
+             m_aer=m_aer*(mSO2mol/mSatom)
+             
+             WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), &
                   '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
+                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'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
-                ENDIF ! emission grid cell
-             ENDDO ! klon loop
-             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
+             
+             latmin=xlat_min_vol(ieru)
+             latmax=xlat_max_vol(ieru)
+             lonmin=xlon_min_vol(ieru)
+             lonmax=xlon_max_vol(ieru)
+             altemiss = altemiss_vol(ieru)
+             sigma_alt = sigma_alt_vol(ieru)
+             pdt=pdtphys
+             ! stretch emission over one day of eruption
+             stretchlong = 1.
+             
+             CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
+                  m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
+                  stretchlong,1,0)
+             
           ENDIF ! emission period
        ENDDO ! eruption number
@@ -305,69 +259,24 @@
         .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
         .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
-
-     DO i=1,klon
-        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
-        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
-!
-!         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)/FLOAT(injdur_sai)/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)/FLOAT(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
-
+        
+       m_aer=m_aer_emiss_sai
+       !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
+       m_aer=m_aer*(mSO2mol/mSatom)
+       
+       latmin=xlat_sai
+       latmax=xlat_sai
+       lonmin=xlon_sai
+       lonmax=xlon_sai
+       altemiss = altemiss_sai
+       sigma_alt = sigma_alt_sai
+       pdt=0.
+       ! stretch emission over whole year (360d)
+       stretchlong=FLOAT(year_len)
+       
+       CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
+            m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
+            stretchlong,1,0)
+       
+       budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
      ENDIF ! Condition over injection dates
 
@@ -375,75 +284,25 @@
             !     lat_min and lat_max
 
-    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
-        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)/ &
-                      & FLOAT(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
-!          ENDDO
-        ENDIF ! emission grid cell
-      ENDDO ! klon loop
-
-    END SELECT ! emission scenario (flag_sulf_emit)
+       m_aer=m_aer_emiss_sai
+       !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
+       m_aer=m_aer*(mSO2mol/mSatom)
+
+       latmin=xlat_min_sai
+       latmax=xlat_max_sai
+       lonmin=xlon_sai
+       lonmax=xlon_sai
+       altemiss = altemiss_sai
+       sigma_alt = sigma_alt_sai
+       pdt=0.
+       ! stretch emission over whole year (360d)
+       stretchlong=FLOAT(year_len)
+
+       CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
+            m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
+            stretchlong,1,0)
+
+       budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
+       
+    END SELECT ! emission scenario (flag_emit)
 
 !--read background concentrations of OCS and SO2 and lifetimes from input file
@@ -463,5 +322,5 @@
     CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
 
-!--call sedimentation routine 
+!--call sedimentation routine
     CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
 
@@ -480,10 +339,5 @@
       ENDDO
     ENDDO
-
-!    CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
-!    DO it=1, nbtr_bin
-!      CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
-!    ENDDO
-
+    
   END SUBROUTINE traccoag
 
