1 | ! |
---|
2 | ! $Id: traccoag_mod.F90 3526 2019-05-28 13:00:44Z fairhead $ |
---|
3 | ! |
---|
4 | MODULE traccoag_mod |
---|
5 | ! |
---|
6 | ! This module calculates the concentration of aerosol particles in certain size bins |
---|
7 | ! considering coagulation and sedimentation. |
---|
8 | ! |
---|
9 | CONTAINS |
---|
10 | |
---|
11 | SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, & |
---|
12 | presnivs, xlat, xlon, pphis, pphi, & |
---|
13 | t_seri, pplay, paprs, sh, rh, tr_seri) |
---|
14 | |
---|
15 | USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & |
---|
16 | & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part |
---|
17 | |
---|
18 | USE dimphy |
---|
19 | USE infotrac |
---|
20 | USE aerophys |
---|
21 | USE geometry_mod, ONLY : cell_area, boundslat |
---|
22 | USE mod_grid_phy_lmdz |
---|
23 | USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root |
---|
24 | USE mod_phys_lmdz_para, only: gather, scatter |
---|
25 | USE phys_cal_mod |
---|
26 | USE sulfate_aer_mod |
---|
27 | USE phys_local_var_mod, ONLY: stratomask |
---|
28 | USE YOMCST |
---|
29 | USE print_control_mod, ONLY: lunout |
---|
30 | USE strataer_mod |
---|
31 | USE phys_cal_mod, ONLY : year_len |
---|
32 | |
---|
33 | IMPLICIT NONE |
---|
34 | |
---|
35 | ! Input argument |
---|
36 | !--------------- |
---|
37 | REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) |
---|
38 | REAL,INTENT(IN) :: gmtime ! Heure courante |
---|
39 | LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique |
---|
40 | INTEGER,INTENT(IN) :: julien ! Jour julien |
---|
41 | |
---|
42 | REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA) |
---|
43 | REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point |
---|
44 | REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point |
---|
45 | REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol |
---|
46 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche |
---|
47 | |
---|
48 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
---|
49 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
---|
50 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) |
---|
51 | REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique |
---|
52 | REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative |
---|
53 | |
---|
54 | ! Output argument |
---|
55 | !---------------- |
---|
56 | REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
---|
57 | |
---|
58 | ! Local variables |
---|
59 | !---------------- |
---|
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 |
---|
63 | LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below |
---|
64 | REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] |
---|
65 | REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] |
---|
66 | REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m |
---|
67 | REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer |
---|
68 | REAL :: f_lay_sum ! sum of layer emission fractions |
---|
69 | REAL :: alt ! altitude for integral calculation |
---|
70 | INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile |
---|
71 | REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m] |
---|
72 | REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m] |
---|
73 | REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m] |
---|
74 | REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg] |
---|
75 | REAL :: zrho ! Density of air [kg/m3] |
---|
76 | REAL :: zdz ! thickness of atm. model layer in m |
---|
77 | REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg |
---|
78 | REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction |
---|
79 | REAL :: emission ! emission |
---|
80 | REAL :: theta_min, theta_max ! for SAI computation between two latitudes |
---|
81 | REAL :: dlat_loc |
---|
82 | |
---|
83 | IF (is_mpi_root) THEN |
---|
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 |
---|
119 | ENDIF |
---|
120 | |
---|
121 | DO it=1, nbtr_bin |
---|
122 | r_bin(it)=mdw(it)/2. |
---|
123 | ENDDO |
---|
124 | |
---|
125 | !--set boundaries of size bins |
---|
126 | DO it=1, nbtr_bin |
---|
127 | IF (it.EQ.1) THEN |
---|
128 | r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) |
---|
129 | r_lower(it)=r_bin(it)**2./r_upper(it) |
---|
130 | ELSEIF (it.EQ.nbtr_bin) THEN |
---|
131 | r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) |
---|
132 | r_upper(it)=r_bin(it)**2./r_lower(it) |
---|
133 | ELSE |
---|
134 | r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) |
---|
135 | r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) |
---|
136 | ENDIF |
---|
137 | ENDDO |
---|
138 | |
---|
139 | IF (debutphy .and. is_mpi_root) THEN |
---|
140 | DO it=1, nbtr_bin |
---|
141 | WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' |
---|
142 | ENDDO |
---|
143 | ENDIF |
---|
144 | |
---|
145 | !--initialising logical is_strato from stratomask |
---|
146 | is_strato(:,:)=.FALSE. |
---|
147 | WHERE (stratomask.GT.0.5) is_strato=.TRUE. |
---|
148 | |
---|
149 | ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) |
---|
150 | ! H2SO4 mass fraction in aerosol (%) |
---|
151 | CALL stracomp(sh,t_seri,pplay) |
---|
152 | |
---|
153 | ! aerosol density (gr/cm3) |
---|
154 | CALL denh2sa(t_seri) |
---|
155 | |
---|
156 | ! compute factor for converting dry to wet radius (for every grid box) |
---|
157 | f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.) |
---|
158 | |
---|
159 | !--calculate mass of air in every grid box |
---|
160 | DO ilon=1, klon |
---|
161 | DO ilev=1, klev |
---|
162 | m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon) |
---|
163 | ENDDO |
---|
164 | ENDDO |
---|
165 | |
---|
166 | ! IF (debutphy) THEN |
---|
167 | ! CALL gather(tr_seri, tr_seri_glo) |
---|
168 | ! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN |
---|
169 | !--initialising tracer concentrations to zero |
---|
170 | ! DO it=1, nbtr |
---|
171 | ! tr_seri(:,:,it)=0.0 |
---|
172 | ! ENDDO |
---|
173 | ! ENDIF |
---|
174 | ! ENDIF |
---|
175 | |
---|
176 | !--initialise emission diagnostics |
---|
177 | budg_emi_ocs(:)=0.0 |
---|
178 | budg_emi_so2(:)=0.0 |
---|
179 | budg_emi_h2so4(:)=0.0 |
---|
180 | budg_emi_part(:)=0.0 |
---|
181 | |
---|
182 | !--sulfur emission, depending on chosen scenario (flag_sulf_emit) |
---|
183 | SELECT CASE(flag_sulf_emit) |
---|
184 | |
---|
185 | CASE(0) ! background aerosol |
---|
186 | ! do nothing (no emission) |
---|
187 | |
---|
188 | CASE(1) ! volcanic eruption |
---|
189 | !--only emit on day of eruption |
---|
190 | ! stretch emission over one day of Pinatubo eruption |
---|
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 | |
---|
269 | CASE(2) ! stratospheric aerosol injections (SAI) |
---|
270 | ! |
---|
271 | DO i=1,klon |
---|
272 | ! SAI standard scenario with continuous emission from 1 grid point at the equator |
---|
273 | ! SAI emission on single month |
---|
274 | ! SAI continuous emission o |
---|
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. & |
---|
279 | & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN |
---|
280 | ! |
---|
281 | WRITE(lunout,*) 'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur |
---|
282 | ! compute altLMDz |
---|
283 | altLMDz(:)=0.0 |
---|
284 | DO k=1, klev |
---|
285 | zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 |
---|
286 | zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg |
---|
287 | zdz=zdm(k)/zrho !thickness of layer in m |
---|
288 | altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface |
---|
289 | ENDDO |
---|
290 | |
---|
291 | SELECT CASE(flag_sulf_emit_distrib) |
---|
292 | |
---|
293 | CASE(0) ! Gaussian distribution |
---|
294 | !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) |
---|
295 | f_lay_sum=0.0 |
---|
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 | |
---|
319 | !correct for step integration error |
---|
320 | f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum |
---|
321 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
---|
322 | !vertically distributed emission |
---|
323 | DO k=1, klev |
---|
324 | ! stretch emission over whole year (360d) |
---|
325 | emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400. |
---|
326 | tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys |
---|
327 | budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol |
---|
328 | ENDDO |
---|
329 | |
---|
330 | ! !emission as monodisperse particles with 0.1um dry radius (BIN21) |
---|
331 | ! !vertically distributed emission |
---|
332 | ! DO k=1, klev |
---|
333 | ! ! stretch emission over whole year (360d) |
---|
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 |
---|
410 | ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys |
---|
411 | ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol |
---|
412 | ! ENDDO |
---|
413 | ENDIF ! emission grid cell |
---|
414 | ENDDO ! klon loop |
---|
415 | |
---|
416 | END SELECT ! emission scenario (flag_sulf_emit) |
---|
417 | |
---|
418 | !--read background concentrations of OCS and SO2 and lifetimes from input file |
---|
419 | !--update the variables defined in phys_local_var_mod |
---|
420 | CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri) |
---|
421 | |
---|
422 | !--convert OCS to SO2 in the stratosphere |
---|
423 | CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
---|
424 | |
---|
425 | !--convert SO2 to H2SO4 |
---|
426 | CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
---|
427 | |
---|
428 | !--common routine for nucleation and condensation/evaporation with adaptive timestep |
---|
429 | CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) |
---|
430 | |
---|
431 | !--call coagulation routine |
---|
432 | CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) |
---|
433 | |
---|
434 | !--call sedimentation routine |
---|
435 | CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) |
---|
436 | |
---|
437 | !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies |
---|
438 | surf_PM25_sulf(:)=0.0 |
---|
439 | DO i=1,klon |
---|
440 | DO it=1, nbtr_bin |
---|
441 | IF (mdw(it) .LT. 2.5e-6) THEN |
---|
442 | !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & |
---|
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 |
---|
445 | surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & |
---|
446 | & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 & |
---|
447 | & *pplay(i,1)/t_seri(i,1)/RD*1.e9 |
---|
448 | ENDIF |
---|
449 | ENDDO |
---|
450 | ENDDO |
---|
451 | |
---|
452 | ! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2') |
---|
453 | ! DO it=1, nbtr_bin |
---|
454 | ! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ') |
---|
455 | ! ENDDO |
---|
456 | |
---|
457 | END SUBROUTINE traccoag |
---|
458 | |
---|
459 | END MODULE traccoag_mod |
---|