Index: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_emiss_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_emiss_mod.F90	(revision 5110)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_emiss_mod.F90	(revision 5111)
@@ -10,202 +10,203 @@
     USE print_control_mod, ONLY: lunout
     USE lmdz_phys_para, ONLY: is_master
+    USE lmdz_abort_physic, ONLY: abort_physic
 
     ! Local variables
-    INTEGER                  :: ispec
-
-    WRITE(lunout,*) 'IN STRATAER_EMISS INIT WELCOME!'
+    INTEGER :: ispec
+
+    WRITE(lunout, *) 'IN STRATAER_EMISS INIT WELCOME!'
 
     IF (flag_emit==1 .OR. flag_emit==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>0) THEN
-          ALLOCATE(year_emit_vol(nErupt),mth_emit_vol(nErupt),day_emit_vol(nErupt))
-          year_emit_vol=0 ; mth_emit_vol=0 ; day_emit_vol=0
-          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))
-          altemiss_vol=0. ; sigma_alt_vol=0.
-          xlon_min_vol=0. ; xlon_max_vol=0.
-          xlat_min_vol=0. ; xlat_max_vol=0.
-          ! 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
-             ALLOCATE(m_sulf_emiss_vol(nErupt))
-             ALLOCATE(m_aer_emiss_vol(nErupt,nAerErupt))
-             m_aer_emiss_vol=0. ; m_sulf_emiss_vol=0.
-             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))
-             id_species = 0
-             ALLOCATE(m_species_emiss_vol(nErupt,nSpeciesErupt))
-             m_species_emiss_vol=0.
-             ALLOCATE(m_NOx_emiss_vol(nErupt))
-             ALLOCATE(m_H2O_emiss_vol(nErupt))
-             m_Chlore_emiss_vol=0. ; m_Brome_emiss_vol=0.
-             m_NOx_emiss_vol=0. ; m_H2O_emiss_vol=0.
-             ALLOCATE(id_NOx)
-             ALLOCATE(id_H2O)
-             ALLOCATE(budg_emi(klon,nSpeciesErupt))
-             budg_emi(:,:)=0.
+      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>0) THEN
+        ALLOCATE(year_emit_vol(nErupt), mth_emit_vol(nErupt), day_emit_vol(nErupt))
+        year_emit_vol = 0 ; mth_emit_vol = 0 ; day_emit_vol = 0
+        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))
+        altemiss_vol = 0. ; sigma_alt_vol = 0.
+        xlon_min_vol = 0. ; xlon_max_vol = 0.
+        xlat_min_vol = 0. ; xlat_max_vol = 0.
+        ! 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
+          ALLOCATE(m_sulf_emiss_vol(nErupt))
+          ALLOCATE(m_aer_emiss_vol(nErupt, nAerErupt))
+          m_aer_emiss_vol = 0. ; m_sulf_emiss_vol = 0.
+          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
-       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
-
-       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
+          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))
+          id_species = 0
+          ALLOCATE(m_species_emiss_vol(nErupt, nSpeciesErupt))
+          m_species_emiss_vol = 0.
+          ALLOCATE(m_NOx_emiss_vol(nErupt))
+          ALLOCATE(m_H2O_emiss_vol(nErupt))
+          m_Chlore_emiss_vol = 0. ; m_Brome_emiss_vol = 0.
+          m_NOx_emiss_vol = 0. ; m_H2O_emiss_vol = 0.
+          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
+
+      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)
+      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)
+      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
+      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
+      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
+      CALL strataer_ponde_init
+      WRITE(lunout, *) 'IN STRATAER INIT : ponde_lonlat_vol', ponde_lonlat_vol
     ENDIF
 
     !============= Print params =============
     IF (is_master) THEN
-       IF (nErupt > 0) 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
+      IF (nErupt > 0) 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
-          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
+        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
     ENDIF ! if master
 
-    WRITE(lunout,*) 'IN STRATAER_EMISS END'
+    WRITE(lunout, *) 'IN STRATAER_EMISS END'
   END SUBROUTINE strataer_emiss_init
 
@@ -222,35 +223,35 @@
 
     ! local var
-    REAL                :: lat_reg_deg,lon_reg_deg ! lat and lon of grid points in degree
-    INTEGER             :: ieru, i, j
+    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>=xlat_min_vol(ieru)-dlat .AND. lat_reg_deg<xlat_max_vol(ieru)+dlat .AND. &
-                  lon_reg_deg>=xlon_min_vol(ieru)-dlon .AND. lon_reg_deg<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
+    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>=xlat_min_vol(ieru) - dlat .AND. lat_reg_deg<xlat_max_vol(ieru) + dlat .AND. &
+                  lon_reg_deg>=xlon_min_vol(ieru) - dlon .AND. lon_reg_deg<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
+    WRITE(lunout, *) 'IN STRATAER_PONDE_INIT ponde_lonlat: ', ponde_lonlat_vol
 
   END SUBROUTINE strataer_ponde_init
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_local_var_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 5110)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 5111)
@@ -1,5 +1,5 @@
 MODULE strataer_local_var_mod
   ! This module contains strato microphysic model parameters & variables
-
+  USE lmdz_abort_physic, ONLY: abort_physic
   IMPLICIT NONE
 
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/traccoag_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/traccoag_mod.F90	(revision 5110)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/traccoag_mod.F90	(revision 5111)
@@ -1,17 +1,17 @@
 MODULE traccoag_mod
 
-! This module calculates the concentration of aerosol particles in certain size bins
-! considering coagulation and sedimentation.
-
+  ! This module calculates the concentration of aerosol particles in certain size bins
+  ! considering coagulation and sedimentation.
+  USE lmdz_abort_physic, ONLY: abort_physic
 CONTAINS
 
   SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
-       presnivs, xlat, xlon, pphis, pphi, &
-       t_seri, pplay, paprs, sh, rh, tr_seri)
-    
-    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & 
-   budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
-   R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode
-    
+          presnivs, xlat, xlon, pphis, pphi, &
+          t_seri, pplay, paprs, sh, rh, tr_seri)
+
+    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
+            budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
+            R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode
+
     USE dimphy
     USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
@@ -19,5 +19,5 @@
     USE geometry_mod, ONLY: cell_area, boundslat
     USE lmdz_grid_phy
-    USE lmdz_phys_mpi_data, ONLY:  is_mpi_root
+    USE lmdz_phys_mpi_data, ONLY: is_mpi_root
     USE lmdz_phys_para, only: gather, scatter
     USE phys_cal_mod, ONLY: year_len, year_cur, mth_cur, day_cur, hour
@@ -27,127 +27,127 @@
     USE print_control_mod, ONLY: lunout
     USE strataer_local_var_mod
-    
+
     IMPLICIT NONE
 
-! Input argument
-!---------------
-    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
-    REAL,INTENT(IN)    :: gmtime     ! Heure courante
-    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
-    INTEGER,INTENT(IN) :: julien     ! Jour julien
-
-    REAL,DIMENSION(klev),INTENT(IN)        :: presnivs! pressions approximat. des milieux couches (en PA)
-    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point 
-    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point 
-    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
-    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche
-
-    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
-    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)   :: sh      ! humidite specifique
-    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative   
-
-! Output argument
-!----------------
-    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]  
-
-! Local variables
-!----------------
-    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
-    REAL                                   :: m_aer               ! aerosol mass
-    INTEGER                                :: it, k, i, j, ilon, ilev, itime, i_int, ieru
-    LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
-    REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
-    REAL,DIMENSION(klon_glo,klev,nbtr)     :: tr_seri_glo         ! Concentration Traceur [U/KgA]  
-    REAL,DIMENSION(klev+1)                 :: altLMDz             ! altitude of layer interfaces in m
-    REAL,DIMENSION(klev)                   :: f_lay_emiss         ! fraction of emission for every vertical layer
-    REAL                                   :: f_lay_sum           ! sum of layer emission fractions
-    REAL                                   :: alt                 ! altitude for integral calculation
-    INTEGER,PARAMETER                      :: n_int_alt=10        ! number of subintervals for integration over Gaussian emission profile
-    REAL,DIMENSION(nbtr_bin)               :: r_bin               ! particle radius in size bin [m]
-    REAL,DIMENSION(nbtr_bin)               :: r_lower             ! particle radius at lower bin boundary [m]
-    REAL,DIMENSION(nbtr_bin)               :: r_upper             ! particle radius at upper bin boundary [m]
-    REAL,DIMENSION(nbtr_bin)               :: m_part_dry          ! mass of one dry particle in size bin [kg]
-    REAL                                   :: zrho                ! Density of air [kg/m3]
-    REAL                                   :: zdz                 ! thickness of atm. model layer in m
-    REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
-    REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
-    REAL                                   :: emission            ! emission
-    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
+    ! Input argument
+    !---------------
+    REAL, INTENT(IN) :: pdtphys    ! Pas d'integration pour la physique (seconde)
+    REAL, INTENT(IN) :: gmtime     ! Heure courante
+    LOGICAL, INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
+    INTEGER, INTENT(IN) :: julien     ! Jour julien
+
+    REAL, DIMENSION(klev), INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
+    REAL, DIMENSION(klon), INTENT(IN) :: xlat    ! latitudes pour chaque point
+    REAL, DIMENSION(klon), INTENT(IN) :: xlon    ! longitudes pour chaque point
+    REAL, DIMENSION(klon), INTENT(IN) :: pphis   ! geopotentiel du sol
+    REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi    ! geopotentiel de chaque couche
+
+    REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri  ! Temperature
+    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) :: sh      ! humidite specifique
+    REAL, DIMENSION(klon, klev), INTENT(IN) :: rh      ! humidite relative
+
+    ! Output argument
+    !----------------
+    REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
+
+    ! Local variables
+    !----------------
+    REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
+    REAL :: m_aer               ! aerosol mass
+    INTEGER :: it, k, i, j, ilon, ilev, itime, i_int, ieru
+    LOGICAL, DIMENSION(klon, klev) :: is_strato           ! true = above tropopause, false = below
+    REAL, DIMENSION(klon, klev) :: m_air_gridbox       ! mass of air in every grid box [kg]
+    REAL, DIMENSION(klon_glo, klev, nbtr) :: tr_seri_glo         ! Concentration Traceur [U/KgA]
+    REAL, DIMENSION(klev + 1) :: altLMDz             ! altitude of layer interfaces in m
+    REAL, DIMENSION(klev) :: f_lay_emiss         ! fraction of emission for every vertical layer
+    REAL :: f_lay_sum           ! sum of layer emission fractions
+    REAL :: alt                 ! altitude for integral calculation
+    INTEGER, PARAMETER :: n_int_alt = 10        ! number of subintervals for integration over Gaussian emission profile
+    REAL, DIMENSION(nbtr_bin) :: r_bin               ! particle radius in size bin [m]
+    REAL, DIMENSION(nbtr_bin) :: r_lower             ! particle radius at lower bin boundary [m]
+    REAL, DIMENSION(nbtr_bin) :: r_upper             ! particle radius at upper bin boundary [m]
+    REAL, DIMENSION(nbtr_bin) :: m_part_dry          ! mass of one dry particle in size bin [kg]
+    REAL :: zrho                ! Density of air [kg/m3]
+    REAL :: zdz                 ! thickness of atm. model layer in m
+    REAL, DIMENSION(klev) :: zdm                 ! mass of atm. model layer in kg
+    REAL, DIMENSION(klon, klev) :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
+    REAL :: emission            ! emission
+    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 .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_emit: ',flag_emit
+      WRITE(lunout, *) 'in traccoag: date from phys_cal_mod =', year_cur, '-', mth_cur, '-', day_cur, '-', hour
+      WRITE(lunout, *) 'IN traccoag flag_emit: ', flag_emit
     ENDIF
-    
+
     !   radius [m]
-    DO it=1, nbtr_bin
-      r_bin(it)=mdw(it)/2.
-    ENDDO
-
-!--set boundaries of size bins
-    DO it=1, nbtr_bin
-    IF (it==1) THEN
-      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
-      r_lower(it)=r_bin(it)**2./r_upper(it)
-    ELSEIF (it==nbtr_bin) THEN
-      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
-      r_upper(it)=r_bin(it)**2./r_lower(it)
+    DO it = 1, nbtr_bin
+      r_bin(it) = mdw(it) / 2.
+    ENDDO
+
+    !--set boundaries of size bins
+    DO it = 1, nbtr_bin
+      IF (it==1) THEN
+        r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it))
+        r_lower(it) = r_bin(it)**2. / r_upper(it)
+      ELSEIF (it==nbtr_bin) THEN
+        r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1))
+        r_upper(it) = r_bin(it)**2. / r_lower(it)
+      ELSE
+        r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1))
+        r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it))
+      ENDIF
+    ENDDO
+
+    IF (debutphy .and. is_mpi_root) THEN
+      DO it = 1, nbtr_bin
+        WRITE(lunout, *) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
+      ENDDO
+    ENDIF
+
+    !--initialising logical is_strato from stratomask
+    is_strato(:, :) = .FALSE.
+    WHERE (stratomask>0.5) is_strato = .TRUE.
+
+    IF(flag_new_strat_compo) THEN
+      IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
+      ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
+      CALL stracomp_kelvin(sh, t_seri, pplay)
     ELSE
-      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
-      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
+      IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
+      ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
+      ! H2SO4 mass fraction in aerosol (%)
+      CALL stracomp(sh, t_seri, pplay)
+
+      ! aerosol density (gr/cm3)
+      CALL denh2sa(t_seri)
+
+      ! compute factor for converting dry to wet radius (for every grid box)
+      f_r_wet(:, :) = (dens_aer_dry / (DENSO4(:, :) * 1000.) / (R2SO4(:, :) / 100.))**(1. / 3.)
     ENDIF
-    ENDDO
-
-    IF (debutphy .and. is_mpi_root) THEN
-      DO it=1, nbtr_bin
-        WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
-      ENDDO
-    ENDIF
-
-!--initialising logical is_strato from stratomask
-    is_strato(:,:)=.FALSE. 
-    WHERE (stratomask>0.5) is_strato=.TRUE.
-
-    IF(flag_new_strat_compo) THEN
-       IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
-       ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
-       CALL stracomp_kelvin(sh,t_seri,pplay)
-    ELSE
-       IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
-       ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) 
-       ! H2SO4 mass fraction in aerosol (%)
-       CALL stracomp(sh,t_seri,pplay)
-       
-       ! aerosol density (gr/cm3)
-       CALL denh2sa(t_seri)
-       
-       ! compute factor for converting dry to wet radius (for every grid box)
-       f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
-    ENDIF
-    
-!--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)
-       ENDDO
-    ENDDO
-    
-!--initialise emission diagnostics
-    if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:,1)=0.0
-    budg_emi_ocs(:)=0.0
-    budg_emi_so2(:)=0.0
-    budg_emi_h2so4(:)=0.0
-    budg_emi_part(:)=0.0
-
-!--sulfur emission, depending on chosen scenario (flag_emit)
+
+    !--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)
+      ENDDO
+    ENDDO
+
+    !--initialise emission diagnostics
+    if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:, 1) = 0.0
+    budg_emi_ocs(:) = 0.0
+    budg_emi_so2(:) = 0.0
+    budg_emi_h2so4(:) = 0.0
+    budg_emi_part(:) = 0.0
+
+    !--sulfur emission, depending on chosen scenario (flag_emit)
     SELECT CASE(flag_emit)
 
@@ -158,238 +158,238 @@
       !--only emit on day of eruption
       ! stretch emission over one day of Pinatubo eruption
-       DO ieru=1, nErupt
-          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
-             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,'ieru=',ieru
-             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
-             
-             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
-       
+      DO ieru = 1, nErupt
+        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
+          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, 'ieru=', ieru
+          WRITE(lunout, *) 'IN traccoag, dlon=', dlon
+
+          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
+
     CASE(2) ! stratospheric aerosol injections (SAI)
 
-     ! Computing duration of SAI in days...
-     ! ... starting from 0...
-     injdur_sai = 0
-     ! ... then adding whole years from first to (n-1)th...
-     DO yr = year_emit_sai_start, year_emit_sai_end-1
-       ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
-       is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
-       injdur_sai = injdur_sai+365+is_bissext
-     ENDDO
-     ! ... then subtracting part of the first year where no injection yet...
-     is_bissext = (MOD(year_emit_sai_start,4)==0) .AND. (MOD(year_emit_sai_start,100) /= 0 .OR. MOD(year_emit_sai_start,400) == 0)
-     SELECT CASE(mth_emit_sai_start)
-     CASE(2)
-        injdur_sai = injdur_sai-31
-     CASE(3)
-        injdur_sai = injdur_sai-31-28-is_bissext
-     CASE(4)
-        injdur_sai = injdur_sai-31-28-is_bissext-31
-     CASE(5)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30
-     CASE(6)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
-     CASE(7)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
-     CASE(8)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
-     CASE(9)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
-     CASE(10)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
-     CASE(11)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
-     CASE(12)
-        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
-     END SELECT
-     injdur_sai = injdur_sai-day_emit_sai_start+1
-     ! ... then adding part of the n-th year 
-     is_bissext = (MOD(year_emit_sai_end,4)==0) .AND. (MOD(year_emit_sai_end,100) /= 0 .OR. MOD(year_emit_sai_end,400) == 0)
-     SELECT CASE(mth_emit_sai_end)
-     CASE(2)
-        injdur_sai = injdur_sai+31
-     CASE(3)
-        injdur_sai = injdur_sai+31+28+is_bissext
-     CASE(4)
-        injdur_sai = injdur_sai+31+28+is_bissext+31
-     CASE(5)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30
-     CASE(6)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
-     CASE(7)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
-     CASE(8)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
-     CASE(9)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
-     CASE(10)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
-     CASE(11)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
-     CASE(12)
-        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
-     END SELECT
-     injdur_sai = injdur_sai+day_emit_sai_end
-     ! A security: are SAI dates of injection consistent?
-     IF (injdur_sai <= 0) THEN
-        CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
-     ENDIF
-     ! Injection in itself
-     IF (( year_emit_sai_start <= year_cur ) &
-        .AND. ( year_cur <= year_emit_sai_end ) &
-        .AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
-        .AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
-        .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
-        
-       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
+      ! Computing duration of SAI in days...
+      ! ... starting from 0...
+      injdur_sai = 0
+      ! ... then adding whole years from first to (n-1)th...
+      DO yr = year_emit_sai_start, year_emit_sai_end - 1
+        ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
+        is_bissext = (MOD(yr, 4)==0) .AND. (MOD(yr, 100) /= 0 .OR. MOD(yr, 400) == 0)
+        injdur_sai = injdur_sai + 365 + is_bissext
+      ENDDO
+      ! ... then subtracting part of the first year where no injection yet...
+      is_bissext = (MOD(year_emit_sai_start, 4)==0) .AND. (MOD(year_emit_sai_start, 100) /= 0 .OR. MOD(year_emit_sai_start, 400) == 0)
+      SELECT CASE(mth_emit_sai_start)
+      CASE(2)
+        injdur_sai = injdur_sai - 31
+      CASE(3)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext
+      CASE(4)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31
+      CASE(5)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30
+      CASE(6)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31
+      CASE(7)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30
+      CASE(8)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31
+      CASE(9)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31
+      CASE(10)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30
+      CASE(11)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31
+      CASE(12)
+        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 - 30
+      END SELECT
+      injdur_sai = injdur_sai - day_emit_sai_start + 1
+      ! ... then adding part of the n-th year
+      is_bissext = (MOD(year_emit_sai_end, 4)==0) .AND. (MOD(year_emit_sai_end, 100) /= 0 .OR. MOD(year_emit_sai_end, 400) == 0)
+      SELECT CASE(mth_emit_sai_end)
+      CASE(2)
+        injdur_sai = injdur_sai + 31
+      CASE(3)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext
+      CASE(4)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31
+      CASE(5)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30
+      CASE(6)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31
+      CASE(7)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30
+      CASE(8)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31
+      CASE(9)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31
+      CASE(10)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30
+      CASE(11)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31
+      CASE(12)
+        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30
+      END SELECT
+      injdur_sai = injdur_sai + day_emit_sai_end
+      ! A security: are SAI dates of injection consistent?
+      IF (injdur_sai <= 0) THEN
+        CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.', 1)
+      ENDIF
+      ! Injection in itself
+      IF ((year_emit_sai_start <= year_cur) &
+              .AND. (year_cur <= year_emit_sai_end) &
+              .AND. (mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur) &
+              .AND. (mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end) &
+              .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
+
+        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
 
     CASE(3) ! --- SAI injection over a single band of longitude and between
-            !     lat_min and lat_max
-
-       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
-       
+      !     lat_min and lat_max
+
+      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
-!--update the variables defined in phys_local_var_mod
-    CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
-
-!--convert OCS to SO2 in the stratosphere
-    CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
-
-!--convert SO2 to H2SO4
-    CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
-
-!--common routine for nucleation and condensation/evaporation with adaptive timestep
-    CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
-
-!--CALL coagulation routine
-    CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
-
-!--CALL sedimentation routine
+    !--read background concentrations of OCS and SO2 and lifetimes from input file
+    !--update the variables defined in phys_local_var_mod
+    CALL interp_sulf_input(debutphy, pdtphys, paprs, tr_seri)
+
+    !--convert OCS to SO2 in the stratosphere
+    CALL ocs_to_so2(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato)
+
+    !--convert SO2 to H2SO4
+    CALL so2_to_h2so4(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato)
+
+    !--common routine for nucleation and condensation/evaporation with adaptive timestep
+    CALL micphy_tstep(pdtphys, tr_seri, t_seri, pplay, paprs, rh, is_strato)
+
+    !--CALL coagulation routine
+    CALL coagulate(pdtphys, mdw, tr_seri, t_seri, pplay, dens_aer, is_strato)
+
+    !--CALL sedimentation routine
     CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
 
-!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
-    surf_PM25_sulf(:)=0.0
-    DO i=1,klon
-      DO it=1, nbtr_bin
+    !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
+    surf_PM25_sulf(:) = 0.0
+    DO i = 1, klon
+      DO it = 1, nbtr_bin
         IF (mdw(it) < 2.5e-6) THEN
           !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
           !assume that particles consist of ammonium sulfate at the surface (132g/mol) 
           !and are dry at T = 20 deg. C and 50 perc. humidity
-          surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
-   *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
-   *pplay(i,1)/t_seri(i,1)/RD*1.e9
+          surf_PM25_sulf(i) = surf_PM25_sulf(i) + tr_seri(i, 1, it + nbtr_sulgas) &
+                  * 132. / 98. * dens_aer_dry * 4. / 3. * RPI * (mdw(it) / 2.)**3 &
+                  * pplay(i, 1) / t_seri(i, 1) / RD * 1.e9
         ENDIF
       ENDDO
     ENDDO
-    
-!--compute 
-!     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
-!     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
-!     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
-!     nd_mode: DRY(?) particle concentration in different modes (part/m3)
-     sulfmmr(:,:)=0.0
-     SAD_sulfate(:,:)=0.0
-     sulfmmr_mode(:,:,:)=0.0
-     nd_mode(:,:,:)=0.0
-     
-     DO i=1,klon
-        DO j=1,klev
-           DO it=1, nbtr_bin
-              !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
-              !assume that particles consist of ammonium sulfate at the surface (132g/mol) 
-              !and are dry at T = 20 deg. C and 50 perc. humidity
-              
-              !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
-              !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
-              sulfmmr_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &        ! [DRY part/kgA in bin it] 
-    *(4./3.)*RPI*(mdw(it)/2.)**3.   &                   ! [mdw: dry diameter in m]
-    *dens_aer_dry                                       ! [dry aerosol mass density in kg/m3]
-              
-              !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio) 
-              !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
-              sulfmmr(i,j) = sulfmmr(i,j) + sulfmmr_mode(i,j,it)
-              
-              !     nd_mode: particle concentration in different modes (DRY part/m3)
-              nd_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &             ! [DRY part/kgA in bin it] 
-   *pplay(i,j)/t_seri(i,j)/RD                           ! [air mass concentration in kg air /m3A]
-              
-              IF(flag_new_strat_compo) THEN
-                 !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
-                 SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)] 
-    *4.*RPI*( mdw(it)*f_r_wetB(i,j,it)/2. )**2. &       ! [WET SA of part it in m2]
-    *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
-              ELSE
-                 !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
-                 SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)] 
-    *4.*RPI*( mdw(it)*f_r_wet(i,j)/2. )**2. &           ! [WET SA of part it in m2]
-    *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
-              ENDIF
-           ENDDO
+
+    !--compute
+    !     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
+    !     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
+    !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
+    !     nd_mode: DRY(?) particle concentration in different modes (part/m3)
+    sulfmmr(:, :) = 0.0
+    SAD_sulfate(:, :) = 0.0
+    sulfmmr_mode(:, :, :) = 0.0
+    nd_mode(:, :, :) = 0.0
+
+    DO i = 1, klon
+      DO j = 1, klev
+        DO it = 1, nbtr_bin
+          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
+          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
+          !and are dry at T = 20 deg. C and 50 perc. humidity
+
+          !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
+          !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
+          sulfmmr_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) &        ! [DRY part/kgA in bin it]
+                  * (4. / 3.) * RPI * (mdw(it) / 2.)**3.   &                   ! [mdw: dry diameter in m]
+                  * dens_aer_dry                                       ! [dry aerosol mass density in kg/m3]
+
+          !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio)
+          !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
+          sulfmmr(i, j) = sulfmmr(i, j) + sulfmmr_mode(i, j, it)
+
+          !     nd_mode: particle concentration in different modes (DRY part/m3)
+          nd_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) &             ! [DRY part/kgA in bin it]
+                  * pplay(i, j) / t_seri(i, j) / RD                           ! [air mass concentration in kg air /m3A]
+
+          IF(flag_new_strat_compo) THEN
+            !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
+            SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
+                    * 4. * RPI * (mdw(it) * f_r_wetB(i, j, it) / 2.)**2. &       ! [WET SA of part it in m2]
+                    * 1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
+          ELSE
+            !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
+            SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
+                    * 4. * RPI * (mdw(it) * f_r_wet(i, j) / 2.)**2. &           ! [WET SA of part it in m2]
+                    * 1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
+          ENDIF
         ENDDO
-     ENDDO
-     
+      ENDDO
+    ENDDO
+
   END SUBROUTINE traccoag
 
