Context Navigation

traccoag_mod.F90 @ 5411

Last change on this file since 5411 was 4727, checked in by idelkadi, 14 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
Property svn:keywords set to `Id`
File size: 14.6 KB

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

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

Download in other formats:

Original Format