Ignore:
Timestamp:
May 28, 2019, 2:52:20 PM (6 years ago)
Author:
Laurent Fairhead
Message:

Modifs necessaires à la version 6.0.10
OB

Location:
LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/micphy_tstep.F90

    • Property svn:keywords set to Id
    r3098 r3525  
     1!
     2! $Id$
     3!
    14SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
    25
     6  USE geometry_mod, ONLY : latitude_deg !NL- latitude corr. to local domain
    37  USE dimphy, ONLY : klon,klev
    48  USE aerophys
    … …  
    913  USE sulfate_aer_mod, ONLY : STRAACT
    1014  USE YOMCST, ONLY : RPI, RD, RG
    11 
     15  USE print_control_mod, ONLY: lunout
     16  USE strataer_mod
     17 
    1218  IMPLICIT NONE
    1319
    … …  
    8995      ! compute nucleation rate in kg(H2SO4)/kgA/s
    9096      CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), &
    91              & a_xm,b_xm,c_xm,nucl_rate,ntot,x)
     97           & a_xm,b_xm,c_xm,nucl_rate,ntot,x)
     98      !NL - add nucleation box (if flag on)
     99      IF (flag_nuc_rate_box) THEN
     100         IF (latitude_deg(ilon).LE.(nuclat_min) .OR. latitude_deg(ilon).GE.(nuclat_max) &
     101              .OR. pplay(ilon,ilev).GE.nucpres_max .AND. pplay(ilon,ilev) .LE. nucpres_min ) THEN
     102            nucl_rate=0.0
     103         ENDIF
     104      ENDIF
    92105      ! compute cond/evap rate in kg(H2SO4)/kgA/s
    93106      CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
    … …  
    160173    DO it=1, nbtr
    161174      IF (tr_seri(ilon,ilev,it).LT.0.0) THEN
    162         PRINT *, 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it
     175         WRITE(lunout,*) 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it
    163176      ENDIF
    164177    ENDDO

  • LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/StratAer/traccoag_mod.F90

    • Property svn:keywords set to Id
    r3114 r3525  
     1!
     2! $Id$
     3!
    14MODULE traccoag_mod
    25!
    … …  
    1619    USE infotrac
    1720    USE aerophys
    18     USE geometry_mod, ONLY : cell_area
     21    USE geometry_mod, ONLY : cell_area, boundslat
    1922    USE mod_grid_phy_lmdz
    2023    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    … …  
    2427    USE phys_local_var_mod, ONLY: stratomask
    2528    USE YOMCST
     29    USE print_control_mod, ONLY: lunout
     30    USE strataer_mod
     31    USE phys_cal_mod, ONLY : year_len
    2632
    2733    IMPLICIT NONE
    … …  
    5258! Local variables
    5359!----------------
    54 ! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI)
    55     INTEGER,PARAMETER  :: flag_sulf_emit=2
    56 !
    57 !--flag_sulf_emit=1 --example Pinatubo
    58     INTEGER,PARAMETER :: year_emit_vol=1991          ! year of emission date
    59     INTEGER,PARAMETER :: mth_emit_vol=6              ! month of emission date
    60     INTEGER,PARAMETER :: day_emit_vol=15             ! day of emission date
    61     REAL,PARAMETER    :: m_aer_emiss_vol=7.e9   ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
    62     REAL,PARAMETER    :: altemiss_vol=17.e3     ! emission altitude in m
    63     REAL,PARAMETER    :: sigma_alt_vol=1.e3     ! standard deviation of emission altitude in m
    64     REAL,PARAMETER    :: xlat_vol=15.14         ! latitude of volcano in degree
    65     REAL,PARAMETER    :: xlon_vol=120.35        ! longitude of volcano in degree
    66 
    67 !--flag_sulf_emit=2 --SAI
    68     REAL,PARAMETER    :: m_aer_emiss_sai=1.e10  ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
    69     REAL,PARAMETER    :: altemiss_sai=17.e3     ! emission altitude in m
    70     REAL,PARAMETER    :: sigma_alt_sai=1.e3     ! standard deviation of emission altitude in m
    71     REAL,PARAMETER    :: xlat_sai=0.01          ! latitude of SAI in degree
    72     REAL,PARAMETER    :: xlon_sai=120.35        ! longitude of SAI in degree
    73 
    74 !--other local variables
    75     INTEGER                                :: it, k, i, ilon, ilev, itime, i_int
     60    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
     61    REAL                                   :: sum_emi_so2         ! Test sum of all LON for budg_emi_so2
     62    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
    7663    LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
    7764    REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
    … …  
    9077    REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
    9178    REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
    92     REAL                                   :: dlat, dlon          ! d latitude and d longitude of grid in degree
    9379    REAL                                   :: emission            ! emission
     80    REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
     81    REAL                                   :: dlat_loc
    9482
    9583    IF (is_mpi_root) THEN
    96       PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
     84       WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
     85       WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit
     86       IF (flag_sulf_emit == 1) THEN
     87          WRITE(lunout,*) 'IN traccoag nErupt: ',nErupt
     88          WRITE(lunout,*) 'IN traccoag injdur: ',injdur
     89          WRITE(lunout,*) 'IN traccoag : year_emit_vol',year_emit_vol
     90          WRITE(lunout,*) 'IN traccoag : mth_emit_vol',mth_emit_vol
     91          WRITE(lunout,*) 'IN traccoag : day_emit_vol',day_emit_vol
     92          WRITE(lunout,*) 'IN traccoag : m_aer_emiss_vol',m_aer_emiss_vol
     93          WRITE(lunout,*) 'IN traccoag : altemiss_vol',altemiss_vol
     94          WRITE(lunout,*) 'IN traccoag : sigma_alt_vol',sigma_alt_vol
     95          WRITE(lunout,*) 'IN traccoag : ponde_lonlat_vol',ponde_lonlat_vol
     96          WRITE(lunout,*) 'IN traccoag : xlat_min_vol',xlat_min_vol
     97          WRITE(lunout,*) 'IN traccoag : xlat_max_vol',xlat_max_vol
     98          WRITE(lunout,*) 'IN traccoag : xlon_min_vol',xlon_min_vol
     99          WRITE(lunout,*) 'IN traccoag : xlon_max_vol',xlon_max_vol
     100       ELSEIF (flag_sulf_emit == 2) THEN
     101          WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai
     102          WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai
     103          WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai
     104          WRITE(lunout,*) 'IN traccoag : xlat_sai',xlat_sai
     105          WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai
     106       ELSEIF (flag_sulf_emit == 3) THEN
     107          WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai
     108          WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai
     109          WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai
     110          WRITE(lunout,*) 'IN traccoag : xlat_min_sai',xlat_min_sai
     111          WRITE(lunout,*) 'IN traccoag : xlat_max_sai',xlat_max_sai
     112          WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai
     113       ENDIF
     114       WRITE(lunout,*) 'IN traccoag : flag_nuc_rate_box = ',flag_nuc_rate_box
     115       IF (flag_nuc_rate_box) THEN
     116          WRITE(lunout,*) 'IN traccoag : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max
     117          WRITE(lunout,*) 'IN traccoag : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max
     118       ENDIF
    97119    ENDIF
    98 
    99     dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
    100     dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
    101 
     120   
    102121    DO it=1, nbtr_bin
    103122      r_bin(it)=mdw(it)/2.
    … …  
    120139    IF (debutphy .and. is_mpi_root) THEN
    121140      DO it=1, nbtr_bin
    122         PRINT *,'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
     141        WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
    123142      ENDDO
    124143    ENDIF
    … …  
    170189      !--only emit on day of eruption
    171190      ! stretch emission over one day of Pinatubo eruption
    172       IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN
    173 !
    174         DO i=1,klon
    175           !Pinatubo eruption at 15.14N, 120.35E
    176           IF  ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
    177                 xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
    178 !
    179           PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
    180 !         compute altLMDz
    181             altLMDz(:)=0.0
    182             DO k=1, klev
    183               zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
    184               zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
    185               zdz=zdm(k)/zrho                      !thickness of layer in m
    186               altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
    187             ENDDO
    188             !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
    189             f_lay_sum=0.0
    190             DO k=1, klev
    191               f_lay_emiss(k)=0.0
    192               DO i_int=1, n_int_alt
    193                 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    194                 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* &
    195                 &              exp(-0.5*((alt-altemiss_vol)/sigma_alt_vol)**2.)*   &
    196                 &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    197               ENDDO
    198               f_lay_sum=f_lay_sum+f_lay_emiss(k)
    199             ENDDO
    200             !correct for step integration error
    201             f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
    202             !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    203             !vertically distributed emission
    204             DO k=1, klev
    205               ! stretch emission over one day (minus one timestep) of Pinatubo eruption
    206               emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
    207               tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
    208               budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
    209             ENDDO
    210           ENDIF ! emission grid cell
    211         ENDDO ! klon loop
    212       ENDIF ! emission period
    213 
     191       DO ieru=1, nErupt
     192          IF (is_mpi_root) THEN
     193             sum_emi_so2 = 0.0 ! Init sum
     194          ENDIF
     195          IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
     196               day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
     197             !
     198             ! daily injection mass emission - NL
     199             m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
     200             WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
     201                  'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
     202                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
     203             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
     204             DO i=1,klon
     205                !Pinatubo eruption at 15.14N, 120.35E
     206                dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
     207                WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
     208                IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
     209                     xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
     210                   !
     211                   WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
     212                   WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
     213                   !         compute altLMDz
     214                   altLMDz(:)=0.0
     215                   DO k=1, klev
     216                      zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
     217                      zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
     218                      zdz=zdm(k)/zrho                      !thickness of layer in m
     219                      altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
     220                   ENDDO
     221
     222                   SELECT CASE(flag_sulf_emit_distrib)
     223                   
     224                   CASE(0) ! Gaussian distribution
     225                   !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
     226                   f_lay_sum=0.0
     227                   DO k=1, klev
     228                      f_lay_emiss(k)=0.0
     229                      DO i_int=1, n_int_alt
     230                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
     231                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* &
     232                              &              exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)*   &
     233                              &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
     234                      ENDDO
     235                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
     236                   ENDDO
     237                   
     238                   CASE(1) ! Uniform distribution
     239                   ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
     240                   ! height of the injection, centered around altemiss_vol(ieru)
     241                   DO k=1, klev
     242                      f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
     243                      & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
     244                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
     245                   ENDDO
     246
     247                   END SELECT        ! End CASE over flag_sulf_emit_distrib)
     248
     249                   WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
     250                   WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
     251                   !correct for step integration error
     252                   f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
     253                   !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     254                   !vertically distributed emission
     255                   DO k=1, klev
     256                      ! stretch emission over one day of Pinatubo eruption
     257                      emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
     258                      tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
     259                      budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
     260                   ENDDO
     261                   sum_emi_so2 = sum_emi_so2 + budg_emi_so2(i) ! Sum all LON
     262                ENDIF ! emission grid cell
     263             ENDDO ! klon loop
     264             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") global sum_emi_so2=",sum_emi_so2
     265             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
     266          ENDIF ! emission period
     267       ENDDO ! eruption number
     268       
    214269    CASE(2) ! stratospheric aerosol injections (SAI)
    215270!
    … …  
    217272!       SAI standard scenario with continuous emission from 1 grid point at the equator
    218273!       SAI emission on single month
    219 !       IF  ((mth_cur==4 .AND. &
    220274!       SAI continuous emission o
    221         IF  ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
     275        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
     276        WRITE(lunout,*) "IN traccoag, dlon=",dlon
     277        WRITE(lunout,*) "IN traccoag, dlat=",dlat_loc
     278        IF  ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
    222279          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
    223280!
    224           PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
     281          WRITE(lunout,*) 'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
    225282!         compute altLMDz
    226283          altLMDz(:)=0.0
    … …  
    231288            altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
    232289          ENDDO
     290
     291          SELECT CASE(flag_sulf_emit_distrib)
     292
     293          CASE(0) ! Gaussian distribution
    233294          !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
    234295          f_lay_sum=0.0
    235           DO k=1, klev
    236             f_lay_emiss(k)=0.0
    237             DO i_int=1, n_int_alt
    238               alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    239               f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
    240               &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
    241               &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    242             ENDDO
    243             f_lay_sum=f_lay_sum+f_lay_emiss(k)
    244           ENDDO
     296               DO k=1, klev
     297                     f_lay_emiss(k)=0.0
     298                     DO i_int=1, n_int_alt
     299                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
     300                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
     301                         &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
     302                         &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
     303                     ENDDO
     304                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
     305               ENDDO
     306
     307          CASE(1) ! Uniform distribution
     308          f_lay_sum=0.0
     309          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
     310          ! the height of the injection, centered around altemiss_sai
     311               DO k=1, klev
     312                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
     313                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
     314                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
     315               ENDDO
     316
     317          END SELECT ! Gaussian or uniform distribution
     318
    245319          !correct for step integration error
    246320          f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
    … …  
    249323          DO k=1, klev
    250324            ! stretch emission over whole year (360d)
    251             emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400. 
     325            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400. 
    252326            tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
    253327            budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
    254328          ENDDO
     329
    255330!          !emission as monodisperse particles with 0.1um dry radius (BIN21)
    256331!          !vertically distributed emission
    257332!          DO k=1, klev
    258333!            ! stretch emission over whole year (360d)
    259 !            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400
     334!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400
     335!            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
     336!            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
     337!          ENDDO
     338        ENDIF ! emission grid cell
     339      ENDDO ! klon loop
     340
     341    CASE(3) ! --- SAI injection over a single band of longitude and between
     342            !     lat_min and lat_max
     343
     344    WRITE(lunout,*) 'IN traccoag, dlon=',dlon
     345    DO i=1,klon
     346!       SAI scenario with continuous emission
     347        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
     348        WRITE(lunout,*) 'IN traccoag, dlat = ',dlat_loc
     349        theta_min = max(xlat(i)-dlat_loc,xlat_min_sai)
     350        theta_max = min(xlat(i)+dlat_loc,xlat_max_sai)
     351        IF  ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. &
     352          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
     353!
     354!         compute altLMDz
     355          altLMDz(:)=0.0
     356          DO k=1, klev
     357            zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
     358            zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
     359            zdz=zdm(k)/zrho                      !thickness of layer in m
     360            altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
     361          ENDDO
     362
     363          SELECT CASE(flag_sulf_emit_distrib)
     364
     365          CASE(0) ! Gaussian distribution
     366          !compute distribution of emission to vertical model layers (based on
     367          !Gaussian peak in altitude)
     368          f_lay_sum=0.0
     369               DO k=1, klev
     370                     f_lay_emiss(k)=0.0
     371                     DO i_int=1, n_int_alt
     372                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
     373                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
     374                         & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
     375                         & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
     376                     ENDDO
     377                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
     378               ENDDO
     379
     380          CASE(1) ! Uniform distribution
     381          f_lay_sum=0.0
     382          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
     383          ! the height of the injection, centered around altemiss_sai
     384               DO k=1, klev
     385                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
     386                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
     387                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
     388               ENDDO
     389
     390          END SELECT ! Gaussian or uniform distribution
     391
     392          !correct for step integration error
     393          f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
     394          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     395          !vertically distributed emission
     396          DO k=1, klev
     397            ! stretch emission over whole year (360d)
     398            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ &
     399                      & year_len/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
     400                      & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
     401            tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
     402            budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
     403          ENDDO
     404
     405!          !emission as monodisperse particles with 0.1um dry radius (BIN21)
     406!          !vertically distributed emission
     407!          DO k=1, klev
     408!            ! stretch emission over whole year (360d)
     409!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400
    260410!            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
    261411!            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
    … …  
    291441        IF (mdw(it) .LT. 2.5e-6) THEN
    292442          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
    293           !assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity
     443          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
     444          !and are dry at T = 20 deg. C and 50 perc. humidity
    294445          surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
    295446                           & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
Note: See TracChangeset for help on using the changeset viewer.