Context Navigation

traccoag_mod.F90 @ 3526

Last change on this file since 3526 was 3526, checked in by Laurent Fairhead, 5 years ago
Olivier's modifications + addition of Id keyword
Property svn:keywords set to `Id`
File size: 22.1 KB

Line
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)+emissionzdm(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)+emissionzdm(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)+emissionzdm(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)+emissionzdm(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)+emissionzdm(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_dry4./3.RPI(mdw(it)/2.)*3 &
447	& pplay(i,1)/t_seri(i,1)/RD1.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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format