Context Navigation

traccoag_mod.F90 @ 5456

Last change on this file since 5456 was 4513, checked in by tlurton, 21 months ago
Strataer: a bug correction in strataer_mod.F90, and an implementation of injection duration for SAI case.
Property svn:keywords set to `Id`
File size: 22.7 KB

Line
1	!
2	! $Id: traccoag_mod.F90 4513 2023-04-20 07:57:22Z aborella $
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_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
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, ONLY : year_len, mth_len, year_cur, mth_cur, day_cur, hour
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
32	IMPLICIT NONE
33
34	! Input argument
35	!---------------
36	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
37	REAL,INTENT(IN) :: gmtime ! Heure courante
38	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
39	INTEGER,INTENT(IN) :: julien ! Jour julien
40
41	REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
42	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
43	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
44	REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol
45	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche
46
47	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
48	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
49	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
50	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
51	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
52
53	! Output argument
54	!----------------
55	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
56
57	! Local variables
58	!----------------
59	REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
60	INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru
61	LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below
62	REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg]
63	REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA]
64	REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m
65	REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer
66	REAL :: f_lay_sum ! sum of layer emission fractions
67	REAL :: alt ! altitude for integral calculation
68	INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile
69	REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m]
70	REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m]
71	REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m]
72	REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg]
73	REAL :: zrho ! Density of air [kg/m3]
74	REAL :: zdz ! thickness of atm. model layer in m
75	REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg
76	REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
77	REAL :: emission ! emission
78	REAL :: theta_min, theta_max ! for SAI computation between two latitudes
79	REAL :: dlat_loc
80	INTEGER :: injdur_sai ! injection duration for SAI case [days]
81	INTEGER :: yr, is_bissext
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	ENDIF
87
88	DO it=1, nbtr_bin
89	r_bin(it)=mdw(it)/2.
90	ENDDO
91
92	!--set boundaries of size bins
93	DO it=1, nbtr_bin
94	IF (it.EQ.1) THEN
95	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
96	r_lower(it)=r_bin(it)**2./r_upper(it)
97	ELSEIF (it.EQ.nbtr_bin) THEN
98	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
99	r_upper(it)=r_bin(it)**2./r_lower(it)
100	ELSE
101	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
102	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
103	ENDIF
104	ENDDO
105
106	IF (debutphy .and. is_mpi_root) THEN
107	DO it=1, nbtr_bin
108	WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
109	ENDDO
110	ENDIF
111
112	!--initialising logical is_strato from stratomask
113	is_strato(:,:)=.FALSE.
114	WHERE (stratomask.GT.0.5) is_strato=.TRUE.
115
116	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
117	! H2SO4 mass fraction in aerosol (%)
118	CALL stracomp(sh,t_seri,pplay)
119
120	! aerosol density (gr/cm3)
121	CALL denh2sa(t_seri)
122
123	! compute factor for converting dry to wet radius (for every grid box)
124	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
125
126	!--calculate mass of air in every grid box
127	DO ilon=1, klon
128	DO ilev=1, klev
129	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
130	ENDDO
131	ENDDO
132
133	! IF (debutphy) THEN
134	! CALL gather(tr_seri, tr_seri_glo)
135	! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
136	!--initialising tracer concentrations to zero
137	! DO it=1, nbtr
138	! tr_seri(:,:,it)=0.0
139	! ENDDO
140	! ENDIF
141	! ENDIF
142
143	!--initialise emission diagnostics
144	budg_emi_ocs(:)=0.0
145	budg_emi_so2(:)=0.0
146	budg_emi_h2so4(:)=0.0
147	budg_emi_part(:)=0.0
148
149	!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
150	SELECT CASE(flag_sulf_emit)
151
152	CASE(0) ! background aerosol
153	! do nothing (no emission)
154
155	CASE(1) ! volcanic eruption
156	!--only emit on day of eruption
157	! stretch emission over one day of Pinatubo eruption
158	DO ieru=1, nErupt
159	IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
160	day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
161	!
162	! daily injection mass emission - NL
163	m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
164	WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
165	'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
166	(injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
167	WRITE(lunout,*) 'IN traccoag, dlon=',dlon
168	DO i=1,klon
169	!Pinatubo eruption at 15.14N, 120.35E
170	dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
171	WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
172	IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
173	xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
174	!
175	WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
176	WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
177	! compute altLMDz
178	altLMDz(:)=0.0
179	DO k=1, klev
180	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
181	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
182	zdz=zdm(k)/zrho !thickness of layer in m
183	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
184	ENDDO
185
186	SELECT CASE(flag_sulf_emit_distrib)
187
188	CASE(0) ! Gaussian distribution
189	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
190	f_lay_sum=0.0
191	DO k=1, klev
192	f_lay_emiss(k)=0.0
193	DO i_int=1, n_int_alt
194	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
195	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_vol(ieru))* &
196	& exp(-0.5((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))2.) &
197	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
198	ENDDO
199	f_lay_sum=f_lay_sum+f_lay_emiss(k)
200	ENDDO
201
202	CASE(1) ! Uniform distribution
203	! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
204	! height of the injection, centered around altemiss_vol(ieru)
205	DO k=1, klev
206	f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
207	& max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
208	f_lay_sum=f_lay_sum+f_lay_emiss(k)
209	ENDDO
210
211	END SELECT ! End CASE over flag_sulf_emit_distrib)
212
213	WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
214	WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
215	!correct for step integration error
216	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
217	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
218	!vertically distributed emission
219	DO k=1, klev
220	! stretch emission over one day of Pinatubo eruption
221	emission=m_aer_emiss_vol_daily(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/1./(86400.-pdtphys)
222	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
223	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
224	ENDDO
225	ENDIF ! emission grid cell
226	ENDDO ! klon loop
227	WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
228	ENDIF ! emission period
229	ENDDO ! eruption number
230
231	CASE(2) ! stratospheric aerosol injections (SAI)
232	!
233	! Computing duration of SAI in days...
234	! ... starting from 0...
235	injdur_sai = 0
236	! ... then adding whole years from first to (n-1)th...
237	DO yr = year_emit_sai_start, year_emit_sai_end-1
238	! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
239	is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
240	injdur_sai = injdur_sai+365+is_bissext
241	ENDDO
242	! ... then subtracting part of the first year where no injection yet...
243	is_bissext = (MOD(year_emit_sai_start,4)==0) .AND. (MOD(year_emit_sai_start,100) /= 0 .OR. MOD(year_emit_sai_start,400) == 0)
244	SELECT CASE(mth_emit_sai_start)
245	CASE(2)
246	injdur_sai = injdur_sai-31
247	CASE(3)
248	injdur_sai = injdur_sai-31-28-is_bissext
249	CASE(4)
250	injdur_sai = injdur_sai-31-28-is_bissext-31
251	CASE(5)
252	injdur_sai = injdur_sai-31-28-is_bissext-31-30
253	CASE(6)
254	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
255	CASE(7)
256	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
257	CASE(8)
258	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
259	CASE(9)
260	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
261	CASE(10)
262	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
263	CASE(11)
264	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
265	CASE(12)
266	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
267	END SELECT
268	injdur_sai = injdur_sai-day_emit_sai_start+1
269	! ... then adding part of the n-th year
270	is_bissext = (MOD(year_emit_sai_end,4)==0) .AND. (MOD(year_emit_sai_end,100) /= 0 .OR. MOD(year_emit_sai_end,400) == 0)
271	SELECT CASE(mth_emit_sai_end)
272	CASE(2)
273	injdur_sai = injdur_sai+31
274	CASE(3)
275	injdur_sai = injdur_sai+31+28+is_bissext
276	CASE(4)
277	injdur_sai = injdur_sai+31+28+is_bissext+31
278	CASE(5)
279	injdur_sai = injdur_sai+31+28+is_bissext+31+30
280	CASE(6)
281	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
282	CASE(7)
283	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
284	CASE(8)
285	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
286	CASE(9)
287	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
288	CASE(10)
289	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
290	CASE(11)
291	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
292	CASE(12)
293	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
294	END SELECT
295	injdur_sai = injdur_sai+day_emit_sai_end
296	! A security: are SAI dates of injection consistent?
297	IF (injdur_sai <= 0) THEN
298	CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
299	ENDIF
300	! Injection in itself
301	IF (( year_emit_sai_start <= year_cur ) &
302	.AND. ( year_cur <= year_emit_sai_end ) &
303	.AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
304	.AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
305	.AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
306	.AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
307
308	DO i=1,klon
309	dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
310	IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
311	& xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
312	!
313	! compute altLMDz
314	altLMDz(:)=0.0
315	DO k=1, klev
316	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
317	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
318	zdz=zdm(k)/zrho !thickness of layer in m
319	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
320	ENDDO
321
322	SELECT CASE(flag_sulf_emit_distrib)
323
324	CASE(0) ! Gaussian distribution
325	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
326	f_lay_sum=0.0
327	DO k=1, klev
328	f_lay_emiss(k)=0.0
329	DO i_int=1, n_int_alt
330	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
331	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_sai)* &
332	& exp(-0.5((alt-altemiss_sai)/sigma_alt_sai)2.) &
333	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
334	ENDDO
335	f_lay_sum=f_lay_sum+f_lay_emiss(k)
336	ENDDO
337
338	CASE(1) ! Uniform distribution
339	f_lay_sum=0.0
340	! In this case, parameter sigma_alt_vol(ieru) is considered to be half
341	! the height of the injection, centered around altemiss_sai
342	DO k=1, klev
343	f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
344	& max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
345	f_lay_sum=f_lay_sum+f_lay_emiss(k)
346	ENDDO
347
348	END SELECT ! Gaussian or uniform distribution
349
350	!correct for step integration error
351	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
352	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
353	!vertically distributed emission
354	DO k=1, klev
355	! stretch emission over whole year (360d)
356	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/FLOAT(injdur_sai)/86400.
357	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
358	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
359	ENDDO
360
361	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
362	! !vertically distributed emission
363	! DO k=1, klev
364	! ! stretch emission over whole year (360d)
365	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/FLOAT(year_len)/86400.
366	! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
367	! budg_emi_part(i)=budg_emi_part(i)+emissionzdm(k)mSatom/mH2SO4mol
368	! ENDDO
369	ENDIF ! emission grid cell
370	ENDDO ! klon loop
371
372	ENDIF ! Condition over injection dates
373
374	CASE(3) ! --- SAI injection over a single band of longitude and between
375	! lat_min and lat_max
376
377	DO i=1,klon
378	! SAI scenario with continuous emission
379	dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
380	theta_min = max(xlat(i)-dlat_loc,xlat_min_sai)
381	theta_max = min(xlat(i)+dlat_loc,xlat_max_sai)
382	IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. &
383	& xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
384	!
385	! compute altLMDz
386	altLMDz(:)=0.0
387	DO k=1, klev
388	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
389	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
390	zdz=zdm(k)/zrho !thickness of layer in m
391	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
392	ENDDO
393
394	SELECT CASE(flag_sulf_emit_distrib)
395
396	CASE(0) ! Gaussian distribution
397	!compute distribution of emission to vertical model layers (based on
398	!Gaussian peak in altitude)
399	f_lay_sum=0.0
400	DO k=1, klev
401	f_lay_emiss(k)=0.0
402	DO i_int=1, n_int_alt
403	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
404	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_sai)* &
405	& exp(-0.5((alt-altemiss_sai)/sigma_alt_sai)2.) &
406	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
407	ENDDO
408	f_lay_sum=f_lay_sum+f_lay_emiss(k)
409	ENDDO
410
411	CASE(1) ! Uniform distribution
412	f_lay_sum=0.0
413	! In this case, parameter sigma_alt_vol(ieru) is considered to be half
414	! the height of the injection, centered around altemiss_sai
415	DO k=1, klev
416	f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
417	& max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
418	f_lay_sum=f_lay_sum+f_lay_emiss(k)
419	ENDDO
420
421	END SELECT ! Gaussian or uniform distribution
422
423	!correct for step integration error
424	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
425	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
426	!vertically distributed emission
427	DO k=1, klev
428	! stretch emission over whole year (360d)
429	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/ &
430	& FLOAT(year_len)/86400.(sin(theta_max/180.RPI)-sin(theta_min/180.*RPI))/ &
431	& (sin(xlat_max_sai/180.RPI)-sin(xlat_min_sai/180.RPI))
432	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
433	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
434	ENDDO
435
436	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
437	! !vertically distributed emission
438	! DO k=1, klev
439	! ! stretch emission over whole year (360d)
440	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/year_len/86400
441	! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
442	! budg_emi_part(i)=budg_emi_part(i)+emissionzdm(k)mSatom/mH2SO4mol
443	! ENDDO
444	ENDIF ! emission grid cell
445	ENDDO ! klon loop
446
447	END SELECT ! emission scenario (flag_sulf_emit)
448
449	!--read background concentrations of OCS and SO2 and lifetimes from input file
450	!--update the variables defined in phys_local_var_mod
451	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
452
453	!--convert OCS to SO2 in the stratosphere
454	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
455
456	!--convert SO2 to H2SO4
457	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
458
459	!--common routine for nucleation and condensation/evaporation with adaptive timestep
460	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
461
462	!--call coagulation routine
463	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
464
465	!--call sedimentation routine
466	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
467
468	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
469	surf_PM25_sulf(:)=0.0
470	DO i=1,klon
471	DO it=1, nbtr_bin
472	IF (mdw(it) .LT. 2.5e-6) THEN
473	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
474	!assume that particles consist of ammonium sulfate at the surface (132g/mol)
475	!and are dry at T = 20 deg. C and 50 perc. humidity
476	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
477	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
478	& pplay(i,1)/t_seri(i,1)/RD1.e9
479	ENDIF
480	ENDDO
481	ENDDO
482
483	! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
484	! DO it=1, nbtr_bin
485	! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
486	! ENDDO
487
488	END SUBROUTINE traccoag
489
490	END MODULE traccoag_mod

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

Download in other formats:

Original Format