source: LMDZ5/trunk/libf/phylmd/StratAer/traccoag_mod.F90 @ 3564

Last change on this file since 3564 was 2752, checked in by oboucher, 8 years ago

Changes to StratAer? module, including setting up a new output file
hence dimension of output metadata is increased to 10 items.

File size: 13.9 KB
Line 
1MODULE traccoag_mod
2!
3! This module calculates the concentration of aerosol particles in certain size bins
4! considering coagulation and sedimentation.
5!
6CONTAINS
7
8  SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
9       presnivs, xlat, xlon, pphis, pphi, &
10       t_seri, pplay, paprs, sh, rh, tr_seri)
11
12    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
13        & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part
14
15    USE dimphy
16    USE infotrac
17    USE aerophys
18    USE geometry_mod, ONLY : cell_area
19    USE mod_grid_phy_lmdz
20    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
21    USE mod_phys_lmdz_para, only: gather, scatter
22    USE phys_cal_mod
23    USE sulfate_aer_mod
24    USE phys_local_var_mod, ONLY: stratomask
25    USE YOMCST
26
27    IMPLICIT NONE
28
29! Input argument
30!---------------
31    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
32    REAL,INTENT(IN)    :: gmtime     ! Heure courante
33    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
34    INTEGER,INTENT(IN) :: julien     ! Jour julien
35
36    REAL,DIMENSION(klev),INTENT(IN)        :: presnivs! pressions approximat. des milieux couches (en PA)
37    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
38    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point
39    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
40    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche
41
42    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
43    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
44    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
45    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
46    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative   
47
48! Output argument
49!----------------
50    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
51
52! Local variables
53!----------------
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
76    LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
77    REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
78    REAL,DIMENSION(klon_glo,klev,nbtr)     :: tr_seri_glo         ! Concentration Traceur [U/KgA] 
79    REAL,DIMENSION(klev+1)                 :: altLMDz             ! altitude of layer interfaces in m
80    REAL,DIMENSION(klev)                   :: f_lay_emiss         ! fraction of emission for every vertical layer
81    REAL                                   :: f_lay_sum           ! sum of layer emission fractions
82    REAL                                   :: alt                 ! altitude for integral calculation
83    INTEGER,PARAMETER                      :: n_int_alt=10        ! number of subintervals for integration over Gaussian emission profile
84    REAL,DIMENSION(nbtr_bin)               :: r_bin               ! particle radius in size bin [m]
85    REAL,DIMENSION(nbtr_bin)               :: r_lower             ! particle radius at lower bin boundary [m]
86    REAL,DIMENSION(nbtr_bin)               :: r_upper             ! particle radius at upper bin boundary [m]
87    REAL,DIMENSION(nbtr_bin)               :: m_part_dry          ! mass of one dry particle in size bin [kg]
88    REAL                                   :: zrho                ! Density of air [kg/m3]
89    REAL                                   :: zdz                 ! thickness of atm. model layer in m
90    REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
91    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
93    REAL                                   :: emission            ! emission
94
95    IF (is_mpi_root) THEN
96      PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
97    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
102    DO it=1, nbtr_bin
103      r_bin(it)=mdw(it)/2.
104    ENDDO
105
106!--set boundaries of size bins
107    DO it=1, nbtr_bin
108    IF (it.EQ.1) THEN
109      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
110      r_lower(it)=r_bin(it)**2./r_upper(it)
111    ELSEIF (it.EQ.nbtr_bin) THEN
112      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
113      r_upper(it)=r_bin(it)**2./r_lower(it)
114    ELSE
115      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
116      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
117    ENDIF
118    ENDDO
119
120    IF (debutphy .and. is_mpi_root) THEN
121      DO it=1, nbtr_bin
122        PRINT *,'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
123      ENDDO
124    ENDIF
125
126!--initialising logical is_strato from stratomask
127    is_strato(:,:)=.FALSE.
128    WHERE (stratomask.GT.0.5) is_strato=.TRUE.
129
130! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
131! H2SO4 mass fraction in aerosol (%)
132    CALL stracomp(sh,t_seri,pplay)
133
134! aerosol density (gr/cm3)
135    CALL denh2sa(t_seri)
136
137! compute factor for converting dry to wet radius (for every grid box)
138    f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
139
140!--calculate mass of air in every grid box
141    DO ilon=1, klon
142    DO ilev=1, klev
143      m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
144    ENDDO
145    ENDDO
146
147!    IF (debutphy) THEN
148!      CALL gather(tr_seri, tr_seri_glo)
149!      IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
150!--initialising tracer concentrations to zero
151!        DO it=1, nbtr
152!        tr_seri(:,:,it)=0.0
153!        ENDDO
154!      ENDIF
155!    ENDIF
156
157!--initialise emission diagnostics
158    budg_emi_ocs(:)=0.0
159    budg_emi_so2(:)=0.0
160    budg_emi_h2so4(:)=0.0
161    budg_emi_part(:)=0.0
162
163!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
164    SELECT CASE(flag_sulf_emit)
165
166    CASE(0) ! background aerosol
167      ! do nothing (no emission)
168
169    CASE(1) ! volcanic eruption
170      !--only emit on day of eruption
171      ! 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 of Pinatubo eruption
206              emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./86400.
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
214    CASE(2) ! stratospheric aerosol injections (SAI)
215!
216      DO i=1,klon
217!       SAI standard scenario with continuous emission from 1 grid point at the equator
218!       SAI emission on single month
219!       IF  ((mth_cur==4 .AND. &
220!       SAI continuous emission o
221        IF  ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
222          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
223!
224          PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
225!         compute altLMDz
226          altLMDz(:)=0.0
227          DO k=1, klev
228            zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
229            zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
230            zdz=zdm(k)/zrho                      !thickness of layer in m
231            altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
232          ENDDO
233          !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
234          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
245          !correct for step integration error
246          f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
247          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
248          !vertically distributed emission
249          DO k=1, klev
250            ! 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. 
252            tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
253            budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
254          ENDDO
255!          !emission as monodisperse particles with 0.1um dry radius (BIN21)
256!          !vertically distributed emission
257!          DO k=1, klev
258!            ! 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
260!            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
261!            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
262!          ENDDO
263        ENDIF ! emission grid cell
264      ENDDO ! klon loop
265
266    END SELECT ! emission scenario (flag_sulf_emit)
267
268!--read background concentrations of OCS and SO2 and lifetimes from input file
269!--update the variables defined in phys_local_var_mod
270    CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
271
272!--convert OCS to SO2 in the stratosphere
273    CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
274
275!--convert SO2 to H2SO4
276    CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
277
278!--common routine for nucleation and condensation/evaporation with adaptive timestep
279    CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
280
281!--call coagulation routine
282    CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
283
284!--call sedimentation routine
285    CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
286
287!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
288    surf_PM25_sulf(:)=0.0
289    DO i=1,klon
290      DO it=1, nbtr_bin
291        IF (mdw(it) .LT. 2.5e-6) THEN
292          !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
294          surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
295                           & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
296                           & *pplay(i,1)/t_seri(i,1)/RD*1.e9
297        ENDIF
298      ENDDO
299    ENDDO
300
301!    CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
302!    DO it=1, nbtr_bin
303!      CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
304!    ENDDO
305
306  END SUBROUTINE traccoag
307
308END MODULE traccoag_mod
Note: See TracBrowser for help on using the repository browser.