Context Navigation

traccoag_mod.F90 @ 2694

Last change on this file since 2694 was 2691, checked in by oboucher, 8 years ago
Cleaning up the routine.
File size: 14.4 KB

Line
1	MODULE traccoag_mod
2	!
3	! This module calculates the concentration of aerosol particles in certain size bins
4	! considering coagulation and sedimentation.
5	!
6	CONTAINS
7
8	SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
9	presnivs, xlat, xlon, pphis, pphi, &
10	t_seri, pplay, paprs, sh, rh , &
11	tr_seri)
12
13	USE phys_local_var_mod, ONLY: mdw, sulf_convert, sulf_nucl, sulf_cond_evap, &
14	& sfluxaer, ocs_convert, R2SO4, DENSO4, f_r_wet, SO2_backgr_tend, OCS_backgr_tend, &
15	& OCS_lifetime, SO2_lifetime, surf_PM25_sulf
16
17	USE dimphy
18	USE infotrac
19	USE aerophys
20	USE geometry_mod, ONLY : cell_area
21	USE mod_grid_phy_lmdz
22	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
23	USE mod_phys_lmdz_para, only: gather, scatter
24	USE phys_cal_mod
25	USE sulfate_aer_mod
26	USE nucleate_mod
27	USE phys_local_var_mod, ONLY: stratomask
28	USE sulfate_aer_mod, ONLY : stracomp, denh2sa
29	USE YOMCST
30
31	IMPLICIT NONE
32
33	! Input argument
34	!---------------
35	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
36	REAL,INTENT(IN) :: gmtime ! Heure courante
37	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
38	INTEGER,INTENT(IN) :: julien ! Jour julien
39
40	REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
41	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
42	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
43	REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol
44	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche
45
46	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
47	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
48	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
49	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
50	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
51
52	! Output argument
53	!----------------
54	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
55
56	! Local variables
57	!----------------
58	! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI)
59	INTEGER,PARAMETER :: flag_sulf_emit=2
60	!
61	!--flag_sulf_emit=1 --example Pinatubo
62	INTEGER,PARAMETER :: year_emit_vol=1991 ! year of emission date
63	INTEGER,PARAMETER :: mth_emit_vol=6 ! month of emission date
64	INTEGER,PARAMETER :: day_emit_vol=15 ! day of emission date
65	REAL,PARAMETER :: m_aer_emiss_vol=7.e9 ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
66	REAL,PARAMETER :: altemiss_vol=17.e3 ! emission altitude in m
67	REAL,PARAMETER :: sigma_alt_vol=1.e3 ! standard deviation of emission altitude in m
68	REAL,PARAMETER :: xlat_vol=15.14 ! latitude of volcano in degree
69	REAL,PARAMETER :: xlon_vol=120.35 ! longitude of volcano in degree
70
71	!--flag_sulf_emit=2 --SAI
72	REAL,PARAMETER :: m_aer_emiss_sai=1.e10 ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
73	REAL,PARAMETER :: altemiss_sai=17.e3 ! emission altitude in m
74	REAL,PARAMETER :: sigma_alt_sai=1.e3 ! standard deviation of emission altitude in m
75	REAL,PARAMETER :: xlat_sai=0.0 ! latitude of SAI in degree
76	REAL,PARAMETER :: xlon_sai=120.35 ! longitude of SAI in degree
77
78	!--be careful - this needs to be changed with resolution - here for 96x95
79	REAL,PARAMETER :: dlat=0.9474 ! d latitude in degree
80	REAL,PARAMETER :: dlon=1.875 ! d longitude in degree
81
82	!--other local variables
83	INTEGER :: it, k, i, ilon, ilev, itime, i_int
84	LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below
85	REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg]
86	REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA]
87	REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m
88	REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer
89	REAL :: f_lay_sum ! sum of layer emission fractions
90	REAL :: alt_lower ! integral boundary for Gaussian emission profile in altitude
91	REAL :: alt_upper ! integral boundary for Gaussian emission profile in altitude
92	INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile
93	REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m]
94	REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m]
95	REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m]
96	REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg]
97	REAL :: zrho ! Density of air [kg/m3]
98	REAL :: zdz ! thickness of atm. model layer in m
99	REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
100	REAL,DIMENSION(klon,klev) :: sulf_convert_sub ! SO2 converted to H2SO4 [kg(S)/m2/layer/s]
101	REAL,DIMENSION(klon,klev) :: sulf_nucl_sub ! H2SO4 converted from gas to particles [kg(S)/m2/layer/s]
102	REAL,DIMENSION(klon,klev) :: ocs_convert_sub ! OCS converted to SO2 [kg(S)/m2/layer/s]
103	REAL,DIMENSION(klon,klev) :: SO2_backgr_tend_sub ! SO2 from background [kg(S)/m2/layer/s]
104	REAL,DIMENSION(klon,klev) :: OCS_backgr_tend_sub ! OCS from background [kg(S)/m2/layer/s]
105
106	IF (is_mpi_root) THEN
107	PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
108	ENDIF
109
110	DO it=1, nbtr_bin
111	r_bin(it)=mdw(it)/2.
112	ENDDO
113
114	!--set boundaries of size bins
115	DO it=1, nbtr_bin
116	IF (it.EQ.1) THEN
117	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
118	r_lower(it)=r_bin(it)**2./r_upper(it)
119	ELSEIF (it.EQ.nbtr_bin) THEN
120	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
121	r_upper(it)=r_bin(it)**2./r_lower(it)
122	ELSE
123	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
124	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
125	ENDIF
126	ENDDO
127
128	IF (debutphy .and. is_mpi_root) THEN
129	DO it=1, nbtr_bin
130	PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
131	ENDDO
132	ENDIF
133
134	!--initialising logical is_strato from stratomask
135	is_strato(:,:)=.FALSE.
136	WHERE (stratomask(:,:).GT.0.5) is_strato(:,:)=.TRUE.
137
138	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
139	! H2SO4 mass fraction in aerosol (%)
140	CALL stracomp(sh,t_seri,pplay)
141
142	! aerosol density (gr/cm3)
143	CALL denh2sa(t_seri)
144
145	! compute factor for converting dry to wet radius (for every grid box)
146	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
147
148	!--calculate mass of air in every grid box
149	DO ilon=1, klon
150	DO ilev=1, klev
151	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1)) / RG * cell_area(ilon)
152	ENDDO
153	ENDDO
154
155	IF (debutphy) THEN
156	CALL gather(tr_seri, tr_seri_glo)
157	IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
158	!--initialising tracer concentrations to zero
159	DO it=1, nbtr
160	tr_seri(:,:,it)=0.0
161	ENDDO
162	ENDIF
163	ENDIF
164
165	!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
166	SELECT CASE(flag_sulf_emit)
167
168	CASE(0) ! background aerosol
169	! do nothing (no emission)
170
171	CASE(1) ! volcanic eruption
172	!--only emit on day of eruption
173	! stretch emission over one day of Pinatubo eruption
174	IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN
175	!
176	DO i=1,klon
177	!Pinatubo eruption at 15.14N, 120.35E
178	IF ( xlat(i).GT.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
179	xlon(i).GT.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
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	zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG !thickness of layer in m
185	altLMDz(k+1)=altLMDz(k)+zdz
186	ENDDO
187	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
188	f_lay_sum=0.0
189	DO k=1, klev
190	f_lay_emiss(k)=0.0
191	alt_lower=altLMDz(k)
192	DO i_int=1, n_int_alt
193	alt_upper=alt_lower+(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((altemiss_vol-alt_lower)/sigma_alt_vol)*2.)+ &
196	& exp(-0.5((alt_upper-altemiss_vol)/sigma_alt_vol)2.)) &
197	& (alt_upper-alt_lower)/2.
198	ENDDO
199	f_lay_sum=f_lay_sum+f_lay_emiss(k)
200	ENDDO
201	!correct for step integration error
202	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
203	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
204	!vertically distributed emission
205	DO k=1, klev
206	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+ &
207	& m_aer_emiss_vol(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k) &
208	& /(1.*86400./pdtphys) ! stretch emission over one day of Pinatubo eruption
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).GT.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
222	& xlon(i).GT.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
223	! compute altLMDz
224	altLMDz(:)=0.0
225	DO k=1, klev
226	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
227	zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG !thickness of layer in m
228	altLMDz(k+1)=altLMDz(k)+zdz
229	ENDDO
230	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
231	f_lay_sum=0.0
232	DO k=1, klev
233	f_lay_emiss(k)=0.0
234	alt_lower=altLMDz(k)
235	DO i_int=1, n_int_alt
236	alt_upper=alt_lower+(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
237	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_sai)* &
238	& (exp(-0.5((altemiss_sai-alt_lower)/sigma_alt_sai)*2.)+ &
239	& exp(-0.5((alt_upper-altemiss_sai)/sigma_alt_sai)2.)) &
240	& (alt_upper-alt_lower)/2.
241	ENDDO
242	f_lay_sum=f_lay_sum+f_lay_emiss(k)
243	ENDDO
244	!correct for step integration error
245	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
246	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
247	!vertically distributed emission
248	DO k=1, klev
249	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+ &
250	& m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k) &
251	& /(360.*86400./pdtphys) ! stretch emission over whole year (360d)
252	! & /(60.*86400./pdtphys) ! stretch emission over 2 months (seasonal emission)
253	! & /7. ! distribute equally over 7 emission grid points
254	ENDDO
255	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
256	! !vertically distributed emission
257	! DO k=1, klev
258	! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+ &
259	! & m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k) &
260	! & /(360.*86400./pdtphys) & ! stretch emission over whole year (360d)
261	! & /7. ! distribute equally over 7 emission grid points
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	SO2_backgr_tend_sub(:,:)=0.0
270	OCS_backgr_tend_sub(:,:)=0.0
271	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri,SO2_backgr_tend_sub, &
272	& OCS_backgr_tend_sub,SO2_lifetime,OCS_lifetime)
273	SO2_backgr_tend(:,:)=SO2_backgr_tend_sub(:,:)
274	OCS_backgr_tend(:,:)=OCS_backgr_tend_sub(:,:)
275
276	!--convert OCS to SO2 in the stratosphere
277	ocs_convert_sub(:,:)=0.0
278	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,OCS_lifetime,is_strato,ocs_convert_sub)
279	ocs_convert(:,:)=ocs_convert_sub(:,:)
280
281	!--convert SO2 to H2SO4
282	sulf_convert_sub(:,:)=0.0
283	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,SO2_lifetime,is_strato,sulf_convert_sub)
284	sulf_convert(:,:)=sulf_convert_sub(:,:)
285
286	!--common routine for nucleation and condensation/evaporation with adaptive timestep
287	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
288
289	!--call coagulation routine
290	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
291
292	!--call sedimentation routine
293	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
294
295	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
296	surf_PM25_sulf(:)=0.0
297	DO i=1,klon
298	DO it=1, nbtr_bin
299	IF (mdw(it) .LT. 2.5e-6) THEN
300	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
301	!assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity
302	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
303	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
304	& pplay(i,1)/t_seri(i,1)/RD1e9
305	ENDIF
306	ENDDO
307	ENDDO
308
309	END SUBROUTINE traccoag
310
311	END MODULE traccoag_mod

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

Download in other formats:

Original Format