Context Navigation

traccoag_mod.F90 @ 4794

Last change on this file since 4794 was 4769, checked in by lguez, 11 months ago

Bug fix: only assign if allocated

Bug fix: only assign to budg_emi if it is allocated.

Property svn:keywords set to Id

File size: 15.3 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, 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_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
17	USE aerophys
18	USE geometry_mod, ONLY : cell_area, boundslat
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, ONLY : year_len, year_cur, mth_cur, day_cur, hour
23	USE sulfate_aer_mod
24	USE phys_local_var_mod, ONLY: stratomask
25	USE YOMCST
26	USE print_control_mod, ONLY: lunout
27	USE strataer_local_var_mod
28
29	IMPLICIT NONE
30
31	! Input argument
32	!---------------
33	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
34	REAL,INTENT(IN) :: gmtime ! Heure courante
35	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
36	INTEGER,INTENT(IN) :: julien ! Jour julien
37
38	REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
39	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
40	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
41	REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol
42	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche
43
44	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
45	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
46	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
47	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
48	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
49
50	! Output argument
51	!----------------
52	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
53
54	! Local variables
55	!----------------
56	REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
57	REAL :: m_aer ! aerosol mass
58	INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru
59	LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below
60	REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg]
61	REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA]
62	REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m
63	REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer
64	REAL :: f_lay_sum ! sum of layer emission fractions
65	REAL :: alt ! altitude for integral calculation
66	INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile
67	REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m]
68	REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m]
69	REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m]
70	REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg]
71	REAL :: zrho ! Density of air [kg/m3]
72	REAL :: zdz ! thickness of atm. model layer in m
73	REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg
74	REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
75	REAL :: emission ! emission
76	REAL :: theta_min, theta_max ! for SAI computation between two latitudes
77	REAL :: dlat_loc
78	REAL :: latmin,latmax,lonmin,lonmax ! lat/lon min/max for injection
79	REAL :: sigma_alt, altemiss ! injection altitude + sigma for distrib
80	REAL :: pdt,stretchlong ! physic timestep, stretch emission over one day
81
82	INTEGER :: injdur_sai ! injection duration for SAI case [days]
83	INTEGER :: yr, is_bissext
84
85	IF (is_mpi_root .AND. flag_verbose_strataer) THEN
86	WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
87	WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit
88	ENDIF
89
90	DO it=1, nbtr_bin
91	r_bin(it)=mdw(it)/2.
92	ENDDO
93
94	!--set boundaries of size bins
95	DO it=1, nbtr_bin
96	IF (it.EQ.1) THEN
97	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
98	r_lower(it)=r_bin(it)**2./r_upper(it)
99	ELSEIF (it.EQ.nbtr_bin) THEN
100	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
101	r_upper(it)=r_bin(it)**2./r_lower(it)
102	ELSE
103	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
104	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
105	ENDIF
106	ENDDO
107
108	IF (debutphy .and. is_mpi_root) THEN
109	DO it=1, nbtr_bin
110	WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
111	ENDDO
112	ENDIF
113
114	!--initialising logical is_strato from stratomask
115	is_strato(:,:)=.FALSE.
116	WHERE (stratomask.GT.0.5) is_strato=.TRUE.
117
118	IF(flag_new_strat_compo) THEN
119	IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Tabazadeh 1994', flag_new_strat_compo
120	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) : binary routine (from reprobus)
121	! H2SO4 mass fraction in aerosol (%) from Tabazadeh et al. (1994).
122	CALL stracomp_bin(sh,t_seri,pplay)
123
124	! aerosol density (gr/cm3) - from Tabazadeh
125	CALL denh2sa_taba(t_seri)
126	ELSE
127	IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Bekki 2D model', flag_new_strat_compo
128	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
129	! H2SO4 mass fraction in aerosol (%)
130	CALL stracomp(sh,t_seri,pplay)
131
132	! aerosol density (gr/cm3)
133	CALL denh2sa(t_seri)
134	ENDIF
135
136	! compute factor for converting dry to wet radius (for every grid box)
137	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
138
139	!--calculate mass of air in every grid box
140	DO ilon=1, klon
141	DO ilev=1, klev
142	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
143	ENDDO
144	ENDDO
145
146	!--initialise emission diagnostics
147	if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:,1)=0.0
148	budg_emi_ocs(:)=0.0
149	budg_emi_so2(:)=0.0
150	budg_emi_h2so4(:)=0.0
151	budg_emi_part(:)=0.0
152
153	!--sulfur emission, depending on chosen scenario (flag_emit)
154	SELECT CASE(flag_emit)
155
156	CASE(0) ! background aerosol
157	! do nothing (no emission)
158
159	CASE(1) ! volcanic eruption
160	!--only emit on day of eruption
161	! stretch emission over one day of Pinatubo eruption
162	DO ieru=1, nErupt
163	IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
164	day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
165
166	! daily injection mass emission
167	m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
168	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
169	m_aer=m_aer*(mSO2mol/mSatom)
170
171	WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), &
172	'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
173	(injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'ieru=',ieru
174	WRITE(lunout,*) 'IN traccoag, dlon=',dlon
175
176	latmin=xlat_min_vol(ieru)
177	latmax=xlat_max_vol(ieru)
178	lonmin=xlon_min_vol(ieru)
179	lonmax=xlon_max_vol(ieru)
180	altemiss = altemiss_vol(ieru)
181	sigma_alt = sigma_alt_vol(ieru)
182	pdt=pdtphys
183	! stretch emission over one day of eruption
184	stretchlong = 1.
185
186	CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
187	m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
188	stretchlong,1,0)
189
190	ENDIF ! emission period
191	ENDDO ! eruption number
192
193	CASE(2) ! stratospheric aerosol injections (SAI)
194	!
195	! Computing duration of SAI in days...
196	! ... starting from 0...
197	injdur_sai = 0
198	! ... then adding whole years from first to (n-1)th...
199	DO yr = year_emit_sai_start, year_emit_sai_end-1
200	! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
201	is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
202	injdur_sai = injdur_sai+365+is_bissext
203	ENDDO
204	! ... then subtracting part of the first year where no injection yet...
205	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)
206	SELECT CASE(mth_emit_sai_start)
207	CASE(2)
208	injdur_sai = injdur_sai-31
209	CASE(3)
210	injdur_sai = injdur_sai-31-28-is_bissext
211	CASE(4)
212	injdur_sai = injdur_sai-31-28-is_bissext-31
213	CASE(5)
214	injdur_sai = injdur_sai-31-28-is_bissext-31-30
215	CASE(6)
216	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
217	CASE(7)
218	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
219	CASE(8)
220	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
221	CASE(9)
222	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
223	CASE(10)
224	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
225	CASE(11)
226	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
227	CASE(12)
228	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
229	END SELECT
230	injdur_sai = injdur_sai-day_emit_sai_start+1
231	! ... then adding part of the n-th year
232	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)
233	SELECT CASE(mth_emit_sai_end)
234	CASE(2)
235	injdur_sai = injdur_sai+31
236	CASE(3)
237	injdur_sai = injdur_sai+31+28+is_bissext
238	CASE(4)
239	injdur_sai = injdur_sai+31+28+is_bissext+31
240	CASE(5)
241	injdur_sai = injdur_sai+31+28+is_bissext+31+30
242	CASE(6)
243	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
244	CASE(7)
245	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
246	CASE(8)
247	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
248	CASE(9)
249	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
250	CASE(10)
251	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
252	CASE(11)
253	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
254	CASE(12)
255	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
256	END SELECT
257	injdur_sai = injdur_sai+day_emit_sai_end
258	! A security: are SAI dates of injection consistent?
259	IF (injdur_sai <= 0) THEN
260	CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
261	ENDIF
262	! Injection in itself
263	IF (( year_emit_sai_start <= year_cur ) &
264	.AND. ( year_cur <= year_emit_sai_end ) &
265	.AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
266	.AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
267	.AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
268	.AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
269
270	m_aer=m_aer_emiss_sai
271	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
272	m_aer=m_aer*(mSO2mol/mSatom)
273
274	latmin=xlat_sai
275	latmax=xlat_sai
276	lonmin=xlon_sai
277	lonmax=xlon_sai
278	altemiss = altemiss_sai
279	sigma_alt = sigma_alt_sai
280	pdt=0.
281	! stretch emission over whole year (360d)
282	stretchlong=FLOAT(year_len)
283
284	CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
285	m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
286	stretchlong,1,0)
287
288	budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
289	ENDIF ! Condition over injection dates
290
291	CASE(3) ! --- SAI injection over a single band of longitude and between
292	! lat_min and lat_max
293
294	m_aer=m_aer_emiss_sai
295	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
296	m_aer=m_aer*(mSO2mol/mSatom)
297
298	latmin=xlat_min_sai
299	latmax=xlat_max_sai
300	lonmin=xlon_sai
301	lonmax=xlon_sai
302	altemiss = altemiss_sai
303	sigma_alt = sigma_alt_sai
304	pdt=0.
305	! stretch emission over whole year (360d)
306	stretchlong=FLOAT(year_len)
307
308	CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
309	m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
310	stretchlong,1,0)
311
312	budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
313
314	END SELECT ! emission scenario (flag_emit)
315
316	!--read background concentrations of OCS and SO2 and lifetimes from input file
317	!--update the variables defined in phys_local_var_mod
318	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
319
320	!--convert OCS to SO2 in the stratosphere
321	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
322
323	!--convert SO2 to H2SO4
324	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
325
326	!--common routine for nucleation and condensation/evaporation with adaptive timestep
327	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
328
329	!--call coagulation routine
330	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
331
332	!--call sedimentation routine
333	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
334
335	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
336	surf_PM25_sulf(:)=0.0
337	DO i=1,klon
338	DO it=1, nbtr_bin
339	IF (mdw(it) .LT. 2.5e-6) THEN
340	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
341	!assume that particles consist of ammonium sulfate at the surface (132g/mol)
342	!and are dry at T = 20 deg. C and 50 perc. humidity
343	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
344	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
345	& pplay(i,1)/t_seri(i,1)/RD1.e9
346	ENDIF
347	ENDDO
348	ENDDO
349
350	END SUBROUTINE traccoag
351
352	END MODULE traccoag_mod

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

Download in other formats:

Original Format