Index: LMDZ6/branches/Ocean_skin/libf/phylmd/StratAer/traccoag_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/StratAer/traccoag_mod.F90	(revision 3428)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/StratAer/traccoag_mod.F90	(revision 3605)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 MODULE traccoag_mod
 !
@@ -16,5 +19,5 @@
     USE infotrac
     USE aerophys
-    USE geometry_mod, ONLY : cell_area
+    USE geometry_mod, ONLY : cell_area, boundslat
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
@@ -24,4 +27,7 @@
     USE phys_local_var_mod, ONLY: stratomask
     USE YOMCST
+    USE print_control_mod, ONLY: lunout
+    USE strataer_mod
+    USE phys_cal_mod, ONLY : year_len
 
     IMPLICIT NONE
@@ -52,26 +58,6 @@
 ! Local variables
 !----------------
-! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI)
-    INTEGER,PARAMETER  :: flag_sulf_emit=2 
-!
-!--flag_sulf_emit=1 --example Pinatubo
-    INTEGER,PARAMETER :: year_emit_vol=1991          ! year of emission date
-    INTEGER,PARAMETER :: mth_emit_vol=6              ! month of emission date
-    INTEGER,PARAMETER :: day_emit_vol=15             ! day of emission date 
-    REAL,PARAMETER    :: m_aer_emiss_vol=7.e9   ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
-    REAL,PARAMETER    :: altemiss_vol=17.e3     ! emission altitude in m
-    REAL,PARAMETER    :: sigma_alt_vol=1.e3     ! standard deviation of emission altitude in m
-    REAL,PARAMETER    :: xlat_vol=15.14         ! latitude of volcano in degree
-    REAL,PARAMETER    :: xlon_vol=120.35        ! longitude of volcano in degree
-
-!--flag_sulf_emit=2 --SAI
-    REAL,PARAMETER    :: m_aer_emiss_sai=1.e10  ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
-    REAL,PARAMETER    :: altemiss_sai=17.e3     ! emission altitude in m
-    REAL,PARAMETER    :: sigma_alt_sai=1.e3     ! standard deviation of emission altitude in m
-    REAL,PARAMETER    :: xlat_sai=0.01          ! latitude of SAI in degree
-    REAL,PARAMETER    :: xlon_sai=120.35        ! longitude of SAI in degree
-
-!--other local variables
-    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int
+    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
+    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
     LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
     REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
@@ -90,14 +76,13 @@
     REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
     REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
-    REAL                                   :: dlat, dlon          ! d latitude and d longitude of grid in degree
     REAL                                   :: emission            ! emission
+    REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
+    REAL                                   :: dlat_loc
 
     IF (is_mpi_root) THEN
-      PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
+       WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
+       WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit
     ENDIF
-
-    dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
-    dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
-
+    
     DO it=1, nbtr_bin
       r_bin(it)=mdw(it)/2.
@@ -120,5 +105,5 @@
     IF (debutphy .and. is_mpi_root) THEN
       DO it=1, nbtr_bin
-        PRINT *,'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
+        WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
       ENDDO
     ENDIF
@@ -170,46 +155,77 @@
       !--only emit on day of eruption
       ! stretch emission over one day of Pinatubo eruption
-      IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN 
-!
-        DO i=1,klon
-          !Pinatubo eruption at 15.14N, 120.35E
-          IF  ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
-                xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
-! 
-          PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
-!         compute altLMDz
-            altLMDz(:)=0.0
-            DO k=1, klev
-              zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
-              zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
-              zdz=zdm(k)/zrho                      !thickness of layer in m
-              altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
-            ENDDO
-            !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
-            f_lay_sum=0.0
-            DO k=1, klev
-              f_lay_emiss(k)=0.0
-              DO i_int=1, n_int_alt
-                alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
-                f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* &
-                &              exp(-0.5*((alt-altemiss_vol)/sigma_alt_vol)**2.)*   & 
-                &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 
-              ENDDO
-              f_lay_sum=f_lay_sum+f_lay_emiss(k)
-            ENDDO
-            !correct for step integration error
-            f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
-            !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
-            !vertically distributed emission
-            DO k=1, klev
-              ! stretch emission over one day (minus one timestep) of Pinatubo eruption
-              emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
-              tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 
-              budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
-            ENDDO
-          ENDIF ! emission grid cell
-        ENDDO ! klon loop
-      ENDIF ! emission period
-
+       DO ieru=1, nErupt
+          IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
+               day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
+             !
+             ! daily injection mass emission - NL
+             m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
+             WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
+                  'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
+                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
+             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
+             DO i=1,klon
+                !Pinatubo eruption at 15.14N, 120.35E
+                dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
+                WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
+                IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
+                     xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
+                   !
+                   WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
+                   WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
+                   !         compute altLMDz
+                   altLMDz(:)=0.0
+                   DO k=1, klev
+                      zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
+                      zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
+                      zdz=zdm(k)/zrho                      !thickness of layer in m
+                      altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
+                   ENDDO
+
+                   SELECT CASE(flag_sulf_emit_distrib)
+                   
+                   CASE(0) ! Gaussian distribution
+                   !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
+                   f_lay_sum=0.0
+                   DO k=1, klev
+                      f_lay_emiss(k)=0.0
+                      DO i_int=1, n_int_alt
+                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* &
+                              &              exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)*   &
+                              &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                      ENDDO
+                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
+                   ENDDO
+                   
+                   CASE(1) ! Uniform distribution
+                   ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
+                   ! height of the injection, centered around altemiss_vol(ieru)
+                   DO k=1, klev
+                      f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
+                      & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
+                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
+                   ENDDO
+
+                   END SELECT        ! End CASE over flag_sulf_emit_distrib)
+
+                   WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
+                   WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
+                   !correct for step integration error
+                   f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
+                   !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
+                   !vertically distributed emission
+                   DO k=1, klev
+                      ! stretch emission over one day of Pinatubo eruption
+                      emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
+                      tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
+                      budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
+                   ENDDO
+                ENDIF ! emission grid cell
+             ENDDO ! klon loop
+             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
+          ENDIF ! emission period
+       ENDDO ! eruption number
+       
     CASE(2) ! stratospheric aerosol injections (SAI)
 !
@@ -217,10 +233,9 @@
 !       SAI standard scenario with continuous emission from 1 grid point at the equator
 !       SAI emission on single month
-!       IF  ((mth_cur==4 .AND. &
 !       SAI continuous emission o
-        IF  ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
+        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
+        IF  ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
           &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
 !
-          PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
 !         compute altLMDz
           altLMDz(:)=0.0
@@ -231,16 +246,33 @@
             altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
           ENDDO
+
+          SELECT CASE(flag_sulf_emit_distrib)
+
+          CASE(0) ! Gaussian distribution
           !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
           f_lay_sum=0.0
-          DO k=1, klev
-            f_lay_emiss(k)=0.0
-            DO i_int=1, n_int_alt
-              alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
-              f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
-              &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   & 
-              &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 
-            ENDDO
-            f_lay_sum=f_lay_sum+f_lay_emiss(k)
-          ENDDO
+               DO k=1, klev
+                     f_lay_emiss(k)=0.0
+                     DO i_int=1, n_int_alt
+                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
+                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
+                         &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   & 
+                         &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 
+                     ENDDO
+                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
+               ENDDO
+
+          CASE(1) ! Uniform distribution
+          f_lay_sum=0.0
+          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
+          ! the height of the injection, centered around altemiss_sai
+               DO k=1, klev
+                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
+                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
+                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
+               ENDDO
+
+          END SELECT ! Gaussian or uniform distribution
+
           !correct for step integration error
           f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
@@ -249,13 +281,87 @@
           DO k=1, klev
             ! stretch emission over whole year (360d)
-            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400.  
+            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400.  
             tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
             budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
           ENDDO
+
 !          !emission as monodisperse particles with 0.1um dry radius (BIN21)
 !          !vertically distributed emission
 !          DO k=1, klev
 !            ! stretch emission over whole year (360d)
-!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400 
+!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/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
+
+    CASE(3) ! --- SAI injection over a single band of longitude and between
+            !     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
@@ -291,5 +397,6 @@
         IF (mdw(it) .LT. 2.5e-6) THEN
           !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
-          !assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity
+          !assume that particles consist of ammonium sulfate at the surface (132g/mol) 
+          !and are dry at T = 20 deg. C and 50 perc. humidity
           surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
                            & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
