Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

phytrac_chimie.F

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

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

Download in other formats:

Original Format