Context Navigation

phytrac_chimie.F @ 3567

Last change on this file since 3567 was 3323, checked in by flefevre, 8 months ago
1) revised cloud microphysical parameters (this changes the results) 2) the number of modes is passed as an argument to cloud routines (this does not change the results) 3) some cosmetics to chemparam_mod.F90 (this does not change the results)
Property svn:executable set to ``*
File size: 18.5 KB

Line
1	subroutine phytrac_chimie (
2	$ debutphy,
3	$ gmtime,
4	$ nqmax,
5	$ nlon,
6	$ lat,
7	$ lon,
8	$ zzlay,
9	$ nlev,
10	$ pdtphys,
11	$ temp,
12	$ pplay,
13	$ trac,
14	$ d_tr_chem,
15	$ iter,
16	$ prod_tr,
17	$ loss_tr,
18	$ no_emi,
19	$ o2_emi)
20
21	use chemparam_mod
22	use conc, only: mmean,rnew
23	use photolysis_mod, only: init_photolysis, nphot
24	use iono_h, only: temp_elect
25	#ifdef CPP_XIOS
26	use xios_output_mod, only: send_xios_field
27	#endif
28
29	implicit none
30
31	#include "clesphys.h"
32	#include "YOMCST.h"
33
34	!===================================================================
35	! input
36	!===================================================================
37
38	integer :: nlon, nlev ! number of gridpoints and levels
39	integer :: nqmax ! number of tracers
40
41	real :: gmtime ! day fraction
42	real :: pdtphys ! phytrac_chimie timestep (s)
43	real, dimension(nlon,nlev) :: temp ! temperature (k)
44	real, dimension(nlon,nlev) :: pplay ! pressure (pa)
45	real, dimension(nlon,nlev) :: zzlay ! altitude (m)
46	real, dimension(nlon,nlev,nqmax) :: trac ! tracer mass mixing ratio
47
48	logical :: debutphy ! first call flag
49
50	!===================================================================
51	! output
52	!===================================================================
53
54	integer, dimension(nlon,nlev) :: iter ! chemical iterations
55	real, dimension(nlon,nlev,nqmax) :: d_tr_chem ! chemical tendency for each tracer
56	real, dimension(nlon,nlev,nqmax) :: prod_tr, loss_tr ! production and loss terms (for info)
57	real, dimension(nlon,nlev) :: no_emi ! no emission
58	real, dimension(nlon,nlev) :: o2_emi ! o2 emission
59
60	!===================================================================
61	! local
62	!===================================================================
63
64	real :: sza_local ! solar zenith angle (deg)
65	real :: lon_sun
66	real :: zlocal(nlev) ! altitude for photochem (km)
67	real :: t_elec(nlev) ! electron temperature [K]
68
69	integer, parameter :: t_elec_origin=2
70	!Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data
71
72	integer :: i, iq
73	integer :: ilon, ilev
74
75	real lat(nlon), lat_local(nlon)
76	real lon(nlon), lon_local(nlon)
77
78	real, dimension(nlon,nlev) :: mrtwv, mrtsa ! total water and total sulfuric acid
79	real, dimension(nlon,nlev) :: mrwv, mrsa ! gas-phase water and gas-phase sulfuric acid
80	real, dimension(nlon,nlev) :: trac_sum
81	real, dimension(nlon,nlev,nqmax) :: ztrac ! local tracer mixing ratio
82	real, dimension(nlev) :: no_em
83	real, dimension(nlev) :: o2_em
84
85	integer, save :: nb_reaction_3_max ! number of quadratic reactions
86	integer, save :: nb_reaction_4_max ! number of bimolecular reactions
87	integer, save :: nquench ! number of quenching + heterogeneous reactions
88	integer, save :: nphotion ! number of photoionizations
89	integer, save :: nb_reaction_4_ion ! quadratic reactions for ionosphere
90	! integer, save :: nb_reaction_4_deut ! quadratic reactions for deuterium chem
91	integer, save :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions
92
93	! tracer indexes for the EUV heating:
94	!!! ATTENTION. These values have to be identical to those in euvheat.F90
95	!!! If the values are changed there, the same has to be done here !!!
96
97	! integer,parameter :: i_co2=1
98	! integer,parameter :: i_n2=13
99	! integer,parameter :: i_n=14
100	! integer,parameter :: i_o=3
101	! integer,parameter :: i_co=4
102
103	integer,parameter :: ix_co2 = 1
104	integer,parameter :: ix_co = 2
105	integer,parameter :: ix_o = 3
106	integer,parameter :: ix_o1d = 4
107	integer,parameter :: ix_o2 = 5
108	integer,parameter :: ix_o3 = 6
109	integer,parameter :: ix_h = 7
110	integer,parameter :: ix_h2 = 8
111	integer,parameter :: ix_oh = 9
112	integer,parameter :: ix_ho2 = 10
113	integer,parameter :: ix_h2o2 = 11
114	integer,parameter :: ix_h2o = 12
115	integer,parameter :: ix_n = 13
116	integer,parameter :: ix_n2d = 14
117	integer,parameter :: ix_no = 15
118	integer,parameter :: ix_no2 = 16
119	integer,parameter :: ix_n2 = 17
120
121	! NEED TO BE THE SAME THAN IN EUVHEAT.F90
122	integer,parameter :: nespeuv=17 ! Number of species considered (11, 12 or 17 (with nitrogen))
123
124	real :: vmr_dens_euv(nlon,nlev,nespeuv) ! local species density for EUV heating
125
126	!===================================================================
127	! first call : initialisations
128	!===================================================================
129
130	if (debutphy) then
131
132	!--- Adjustment of Helium amount
133	! if (i_he/=0) then
134	! trac(:,:,i_he)=trac(:,:,i_he)/2.
135	! endif
136	!---
137
138	!-------------------------------------------------------------------
139	! Determination of chemistry reaction number
140	!-------------------------------------------------------------------
141
142	if (ok_chem) then
143	! set number of reactions, depending on ion chemistry or not
144	nb_reaction_4_ion = 64
145	!nb_reaction_4_deut = 35
146
147	!Default numbers if no ion and no deuterium chemistry included
148
149	nb_reaction_4_max = 98 ! set number of bimolecular reactions
150	nb_reaction_3_max = 12 ! set number of quadratic reactions
151	nquench = 13 ! set number of quenching + heterogeneous
152	nphotion = 0 ! set number of photoionizations
153
154	if (ok_ionchem) then
155	nb_reaction_4_max = nb_reaction_4_max+nb_reaction_4_ion
156	nphotion = 18 ! set number of photoionizations
157	endif
158	!if(deutchem) then
159	! nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut
160	!end if
161
162	!nb_phot_max is the total number of processes that are treated
163	!numerically as a photolysis:
164
165	nb_phot_max = nphot + nphotion + nquench
166
167	!-------------------------------------------------------------------
168	! case of tracers re-initialisation with chemistry
169	!-------------------------------------------------------------------
170
171	if (reinit_trac .and. ok_chem) then
172
173	!!! in this reinitialisation, trac is VOLUME mixing ratio
174	! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
175	! convert mass to volume mixing ratio
176	do iq = 1,nqmax - nmicro
177	trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
178	end do
179
180	print*, "SO2 is re-initialised"
181
182	if (i_so2 /= 0) then
183	trac(:,:,i_so2) = 0.
184	trac(:,1:22,i_so2) = 10.e-6
185
186	! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187	! print*, "Tracers are re-initialised"
188	! trac(:,:,:) = 1.0e-30
189
190	! if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
191	! $ .and. (i_so2 /= 0) .and. (i_h2o /= 0)
192	! $ .and. (i_n2 /= 0) .and. (i_co2 /= 0)) then
193
194	! trac(:,1:22,i_ocs) = 3.e-6
195	! trac(:,1:22,i_co) = 25.e-6
196	! trac(:,:,i_hcl) = 0.4e-6
197	! trac(:,1:22,i_so2) = 7.e-6
198	! trac(:,1:22,i_h2o) = 30.e-6
199	! trac(:,:,i_n2) = 0.35e-1
200	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201
202	! ensure that sum of mixing ratios = 1
203
204	trac_sum(:,:) = 0.
205
206	do iq = 1,nqmax - nmicro
207	if (iq /= i_co2) then
208	trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq)
209	end if
210	end do
211
212	! initialise co2
213
214	trac(:,:,i_co2) = 1. - trac_sum(:,:)
215
216	else
217	write(,) "at least one tracer is missing : stop"
218	stop
219	end if
220
221	! update mmean
222
223	mmean(:,:) = 0.
224	do iq = 1,nqmax - nmicro
225	mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
226	enddo
227	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
228
229	! convert volume to mass mixing ratio
230
231	do iq = 1,nqmax - nmicro
232	trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
233	end do
234
235	end if ! reinit_trac
236
237	end if ! ok_chem
238
239	!-------------------------------------------------------------------
240	! case of detailed microphysics without chemistry
241	!-------------------------------------------------------------------
242
243	if (.not. ok_chem .and. ok_cloud .and. cl_scheme == 2) then
244
245	! convert mass to volume mixing ratio
246
247	do iq = 1,nqmax - nmicro
248	ztrac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
249	end do
250
251	! initialise microphysics
252
253	call vapors4muphy_ini(nlon,nlev,ztrac)
254
255	! convert volume to mass mixing ratio
256
257	do iq = 1,nqmax - nmicro
258	trac(:,:,iq) = ztrac(:,:,iq)*m_tr(iq)/mmean(:,:)
259	end do
260
261	end if
262
263	end if ! debutphy
264
265	!===================================================================
266	! convert mass to volume mixing ratio : gas phase
267	!===================================================================
268
269	do iq = 1,nqmax - nmicro
270	ztrac(:,:,iq) = max(trac(:,:,iq)*mmean(:,:)/m_tr(iq), 1.e-30)
271	end do
272
273	!===================================================================
274	! microphysics: simplified scheme (phd aurelien stolzenbach)
275	!===================================================================
276
277	if (ok_cloud .and. cl_scheme == 1) then
278
279	! convert mass to volume mixing ratio : liquid phase
280
281	ztrac(:,:,i_h2so4liq) = max(trac(:,:,i_h2so4liq)
282	$ *mmean(:,:)/m_tr(i_h2so4liq), 1.e-30)
283	ztrac(:,:,i_h2oliq) = max(trac(:,:,i_h2oliq)
284	$ *mmean(:,:)/m_tr(i_h2oliq), 1.e-30)
285
286	! total water and sulfuric acid (gas + liquid)
287
288	mrtwv(:,:) = ztrac(:,:,i_h2o) + ztrac(:,:,i_h2oliq)
289	mrtsa(:,:) = ztrac(:,:,i_h2so4) + ztrac(:,:,i_h2so4liq)
290
291	! all water and sulfuric acid is put in the gas-phase
292
293	mrwv(:,:) = mrtwv(:,:)
294	mrsa(:,:) = mrtsa(:,:)
295
296	! call microphysics
297
298	call new_cloud_venus(nb_mode, nlev, nlon, temp, pplay,
299	$ mrtwv, mrtsa, mrwv, mrsa)
300
301	! update water vapour and sulfuric acid
302
303	ztrac(:,:,i_h2o) = mrwv(:,:)
304	ztrac(:,:,i_h2oliq) = mrtwv(:,:) - ztrac(:,:,i_h2o)
305
306	ztrac(:,:,i_h2so4) = mrsa(:,:)
307	ztrac(:,:,i_h2so4liq) = mrtsa(:,:) - ztrac(:,:,i_h2so4)
308
309	end if ! simplified scheme
310
311	!===================================================================
312	! microphysics: detailed scheme (phd sabrina guilbon)
313	! !!! to be confirmed whether mad_muphy expects mmr or vmr for h2o and h2so4
314	!===================================================================
315
316	if (ok_cloud .and. cl_scheme == 2) then
317
318	do iq = nqmax-nmicro+1,nqmax
319	ztrac(:,:,iq) = trac(:,:,iq)
320	end do
321
322	do ilon = 1,nlon
323	do ilev = 1, nlev
324	if (temp(ilon,ilev) < 500.) then
325	call mad_muphy(pdtphys, ! timestep
326	$ temp(ilon,ilev),pplay(ilon,ilev), ! temperature and pressure
327	$ ztrac(ilon,ilev,i_h2o),
328	$ ztrac(ilon,ilev,i_h2so4),
329	$ ztrac(ilon,ilev,i_m0_aer),
330	$ ztrac(ilon,ilev,i_m3_aer),
331	$ ztrac(ilon,ilev,i_m0_mode1drop),
332	$ ztrac(ilon,ilev,i_m0_mode1ccn),
333	$ ztrac(ilon,ilev,i_m3_mode1sa),
334	$ ztrac(ilon,ilev,i_m3_mode1w),
335	$ ztrac(ilon,ilev,i_m3_mode1ccn),
336	$ ztrac(ilon,ilev,i_m0_mode2drop),
337	$ ztrac(ilon,ilev,i_m0_mode2ccn),
338	$ ztrac(ilon,ilev,i_m3_mode2sa),
339	$ ztrac(ilon,ilev,i_m3_mode2w),
340	$ ztrac(ilon,ilev,i_m3_mode2ccn))
341	else
342	ztrac(ilon,ilev,i_m0_aer) = 0.
343	ztrac(ilon,ilev,i_m3_aer) = 0.
344	ztrac(ilon,ilev,i_m0_mode1drop) = 0.
345	ztrac(ilon,ilev,i_m0_mode1ccn) = 0.
346	ztrac(ilon,ilev,i_m3_mode1sa) = 0.
347	ztrac(ilon,ilev,i_m3_mode1w) = 0.
348	ztrac(ilon,ilev,i_m3_mode1ccn) = 0.
349	ztrac(ilon,ilev,i_m0_mode2drop) = 0.
350	ztrac(ilon,ilev,i_m0_mode2ccn) = 0.
351	ztrac(ilon,ilev,i_m3_mode2sa) = 0.
352	ztrac(ilon,ilev,i_m3_mode2w) = 0.
353	ztrac(ilon,ilev,i_m3_mode2ccn) = 0.
354	end if
355	end do
356	end do
357
358	end if ! detailed scheme
359
360	!===================================================================
361	! photochemistry
362	!===================================================================
363
364	if (ok_chem) then
365
366	lon_sun = (0.5 - gmtime)2.rpi
367	lon_local(:) = lon(:)*rpi/180.
368	lat_local(:) = lat(:)*rpi/180.
369
370	if (ok_ionchem) then
371
372	vmr_dens_euv(:,:,ix_co2) = ztrac(:,:,i_co2) ! CO2
373	vmr_dens_euv(:,:,ix_co) = ztrac(:,:,i_co) ! CO
374	vmr_dens_euv(:,:,ix_o) = ztrac(:,:,i_o) ! O
375	vmr_dens_euv(:,:,ix_o1d) = ztrac(:,:,i_o1d) ! O(1D)
376	vmr_dens_euv(:,:,ix_o2) = ztrac(:,:,i_o2) ! O2
377	vmr_dens_euv(:,:,ix_o3) = ztrac(:,:,i_o3) ! O3
378	vmr_dens_euv(:,:,ix_h) = ztrac(:,:,i_h) ! H
379	vmr_dens_euv(:,:,ix_h2) = ztrac(:,:,i_h2) ! H2
380	vmr_dens_euv(:,:,ix_oh) = ztrac(:,:,i_oh) ! OH
381	vmr_dens_euv(:,:,ix_ho2) = ztrac(:,:,i_ho2) ! HO2
382	vmr_dens_euv(:,:,ix_h2o2)= ztrac(:,:,i_h2o2) ! H2O2
383	vmr_dens_euv(:,:,ix_h2o) = ztrac(:,:,i_h2o) ! H2O
384	vmr_dens_euv(:,:,ix_n) = ztrac(:,:,i_n) ! N
385	vmr_dens_euv(:,:,ix_n2d) = ztrac(:,:,i_n2d) ! N(2D)
386	vmr_dens_euv(:,:,ix_no) = ztrac(:,:,i_no) ! NO
387	vmr_dens_euv(:,:,ix_no2) = ztrac(:,:,i_no2) ! NO2
388	vmr_dens_euv(:,:,ix_n2) = ztrac(:,:,i_n2) ! N2
389
390	end if
391
392	do ilon = 1,nlon
393	zlocal(:)=zzlay(ilon,:)/1000.
394
395	! solar zenith angle
396	! sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))
397	! $ *cos(lon_sun) + cos(lat_local(ilon))
398	! $ sin(lon_local(ilon))sin(lon_sun))*180./rpi
399
400	sza_local = cos(lat_local(ilon))*cos(lon_local(ilon))
401	$ *cos(lon_sun) + cos(lat_local(ilon))
402	$ sin(lon_local(ilon))sin(lon_sun)
403
404	! Security - Handle rare cases where \|sza_local\| > 1
405	sza_local = min(sza_local,1.)
406	sza_local = max(-1.,sza_local)
407	sza_local = acos(sza_local)*180./rpi
408
409	! electron temperature
410	do ilev = 1, nlev
411	t_elec(ilev) = temp_elect(zlocal(ilev),temp(ilon,ilev),
412	$ sza_local,t_elec_origin)
413	end do
414
415	call photochemistry_venus(nlev, nlon, zlocal, pdtphys,
416	$ ok_jonline,ok_ionchem,tuneupperatm,
417	$ nb_reaction_3_max,nb_reaction_4_max,
418	$ nb_phot_max,nphotion,ilon,
419	$ pplay(ilon,:)/100.,
420	$ temp(ilon,:), t_elec(:),
421	$ ztrac(ilon,:,:),vmr_dens_euv(ilon,:,:),
422	$ mmean(ilon,:),
423	$ sza_local,
424	$ lon(ilon), lat(ilon),
425	$ nqmax, nespeuv, iter(ilon,:),
426	$ prod_tr(ilon,:,:),
427	$ loss_tr(ilon,:,:),
428	$ no_em, o2_em)
429
430	no_emi(ilon,:) = no_em(:)
431	o2_emi(ilon,:) = o2_em(:)
432
433	end do
434	end if ! ok_chem
435
436	!===================================================================
437	! compute tendencies
438	!===================================================================
439
440	! update mmean
441
442	mmean(:,:) = 0.
443	do iq = 1,nqmax - nmicro
444	mmean(:,:) = mmean(:,:)+ztrac(:,:,iq)*m_tr(iq)
445	end do
446	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
447
448	!===================================================================
449	! convert volume to mass mixing ratio / then tendencies in mmr
450	!===================================================================
451
452	! gas phase
453
454	do iq = 1,nqmax - nmicro
455	ztrac(:,:,iq) = max(ztrac(:,:,iq)*m_tr(iq)/mmean(:,:),
456	$ 1.e-30)
457	d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
458	end do
459
460	! liquid phase or moments
461
462	if (ok_cloud .and. cl_scheme == 1) then
463	ztrac(:,:,i_h2so4liq) = max(ztrac(:,:,i_h2so4liq)
464	$ *m_tr(i_h2so4liq)/mmean(:,:),
465	$ 1.e-30)
466	ztrac(:,:,i_h2oliq) = max(ztrac(:,:,i_h2oliq)
467	$ *m_tr(i_h2oliq)/mmean(:,:),
468	$ 1.e-30)
469	d_tr_chem(:,:,i_h2so4liq) = (ztrac(:,:,i_h2so4liq)
470	$ - trac(:,:,i_h2so4liq))/pdtphys
471	d_tr_chem(:,:,i_h2oliq) = (ztrac(:,:,i_h2oliq)
472	$ - trac(:,:,i_h2oliq))/pdtphys
473	end if
474
475
476	if (ok_cloud .and. cl_scheme == 2) then
477	do iq = nqmax-nmicro+1,nqmax
478	d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
479	end do
480	end if
481
482	end subroutine phytrac_chimie

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

Download in other formats:

Original Format