Context Navigation

phytrac_chimie.F @ 3094

Last change on this file since 3094 was 3035, checked in by amartinez, 16 months ago

======= VENUS PCM ===============

COMMIT BY ANTOINE MARTINEZ

Implementation of vertical ambipolar diffusion in physics

=================================

NEW KEYWORD OF PHYSIQ.DEF

=================================

NEW VERSION OF "physiq.def"

deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def

ok_iondiff: keyword supposed to activate ion ambipolar diffusion

nbapp_chem: replaces nbapp_chim, in order to characterize the Number of calls to the chemistry routines (per Venusian day)

================

phyvenus

================

Iondiff_red.F

Calculation of the Ion Ambipolar Diffusion for 13 ions on 14:

O2+, O+, C+, N+, CO2+,
NO+, CO+, H+, H2O+, H3O+,

HCO+, N2+, OH+

The ion temperature is fixed as the half of the electron temperature
calculated in ion_h.F for the stability of the program and because the
ion temperature is lower than electron temperature.

plasmaphys_venus_mod.F90

parameters of the ambipolar diffusion scheme:

parameter (Pdiffion = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
parameter (naltvert = 168) ! number of level on the altimetric subgrid
parameter (nsubvert = 24000) ! ptimephysiq/nsubvert - minimum sub-timestep allowed
parameter ( nsubmin = 40) ! ptimephysiq/nsubmin - maximum sub-timestep allowed

physic_mod.F

nbapp_chem is not fixed anymore here but deftank/in physiq.def
Ambipolar diffusion activated if (ok_iondiff) is true

conf_phys.f90

add ok_iondiff as parameters to read from physiq.def file set to .false. by default
add nbapp_chem as parameters to characterize the Number of calls to the chemistry routines (per Venusian day) instead of be fixed in physic_mod.F

to read from physiq.def file set to 24000 by default

cleshphys.h

added ok_iondiff & nbapp_chem in COMMON/clesphys_l/
removed nbapp_chim

phytrac_chemistry.F

Added security in the calculation of the sza_local in order to avoid the rare case where the |range| is above 1

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	print*, "SO2 is re-initialised"
180	if (i_so2 /= 0) then
181	trac(:,1:22,i_so2) = 10.e-6
182
183	! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184	! print*, "Tracers are re-initialised"
185	! trac(:,:,:) = 1.0e-30
186
187	! if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
188	! $ .and. (i_so2 /= 0) .and. (i_h2o /= 0)
189	! $ .and. (i_n2 /= 0) .and. (i_co2 /= 0)) then
190
191	! trac(:,1:22,i_ocs) = 3.e-6
192	! trac(:,1:22,i_co) = 25.e-6
193	! trac(:,:,i_hcl) = 0.4e-6
194	! trac(:,1:22,i_so2) = 7.e-6
195	! trac(:,1:22,i_h2o) = 30.e-6
196	! trac(:,:,i_n2) = 0.35e-1
197	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
198
199	! ensure that sum of mixing ratios = 1
200
201	trac_sum(:,:) = 0.
202
203	do iq = 1,nqmax - nmicro
204	if (iq /= i_co2) then
205	trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq)
206	end if
207	end do
208
209	! initialise co2
210
211	trac(:,:,i_co2) = 1. - trac_sum(:,:)
212
213	else
214	write(,) "at least one tracer is missing : stop"
215	stop
216	end if
217
218	! update mmean
219
220	mmean(:,:) = 0.
221	do iq = 1,nqmax - nmicro
222	mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
223	enddo
224	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
225
226	! convert volume to mass mixing ratio
227
228	do iq = 1,nqmax - nmicro
229	trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
230	end do
231
232	end if ! reinit_trac
233
234	end if ! ok_chem
235
236	!-------------------------------------------------------------------
237	! case of detailed microphysics without chemistry
238	!-------------------------------------------------------------------
239
240	if (.not. ok_chem .and. ok_cloud .and. cl_scheme == 2) then
241
242	! convert mass to volume mixing ratio
243
244	do iq = 1,nqmax - nmicro
245	ztrac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
246	end do
247
248	! initialise microphysics
249
250	call vapors4muphy_ini(nlon,nlev,ztrac)
251
252	! convert volume to mass mixing ratio
253
254	do iq = 1,nqmax - nmicro
255	trac(:,:,iq) = ztrac(:,:,iq)*m_tr(iq)/mmean(:,:)
256	end do
257
258	end if
259
260	end if ! debutphy
261
262	!===================================================================
263	! convert mass to volume mixing ratio : gas phase
264	!===================================================================
265
266	do iq = 1,nqmax - nmicro
267	ztrac(:,:,iq) = max(trac(:,:,iq)*mmean(:,:)/m_tr(iq), 1.e-30)
268	end do
269
270	!===================================================================
271	! microphysics: simplified scheme (phd aurelien stolzenbach)
272	!===================================================================
273
274	if (ok_cloud .and. cl_scheme == 1) then
275
276	! convert mass to volume mixing ratio : liquid phase
277
278	ztrac(:,:,i_h2so4liq) = max(trac(:,:,i_h2so4liq)
279	$ *mmean(:,:)/m_tr(i_h2so4liq), 1.e-30)
280	ztrac(:,:,i_h2oliq) = max(trac(:,:,i_h2oliq)
281	$ *mmean(:,:)/m_tr(i_h2oliq), 1.e-30)
282
283	! total water and sulfuric acid (gas + liquid)
284
285	mrtwv(:,:) = ztrac(:,:,i_h2o) + ztrac(:,:,i_h2oliq)
286	mrtsa(:,:) = ztrac(:,:,i_h2so4) + ztrac(:,:,i_h2so4liq)
287
288	! all water and sulfuric acid is put in the gas-phase
289
290	mrwv(:,:) = mrtwv(:,:)
291	mrsa(:,:) = mrtsa(:,:)
292
293	! call microphysics
294
295	call new_cloud_venus(nlev, nlon, temp, pplay,
296	$ mrtwv, mrtsa, mrwv, mrsa)
297
298	! update water vapour and sulfuric acid
299
300	ztrac(:,:,i_h2o) = mrwv(:,:)
301	ztrac(:,:,i_h2oliq) = mrtwv(:,:) - ztrac(:,:,i_h2o)
302
303	ztrac(:,:,i_h2so4) = mrsa(:,:)
304	ztrac(:,:,i_h2so4liq) = mrtsa(:,:) - ztrac(:,:,i_h2so4)
305
306	end if ! simplified scheme
307
308	!===================================================================
309	! microphysics: detailed scheme (phd sabrina guilbon)
310	! !!! to be confirmed whether mad_muphy expects mmr or vmr for h2o and h2so4
311	!===================================================================
312
313	if (ok_cloud .and. cl_scheme == 2) then
314
315	do iq = nqmax-nmicro+1,nqmax
316	ztrac(:,:,iq) = trac(:,:,iq)
317	end do
318
319	do ilon = 1,nlon
320	do ilev = 1, nlev
321	if (temp(ilon,ilev) < 500.) then
322	call mad_muphy(pdtphys, ! timestep
323	$ temp(ilon,ilev),pplay(ilon,ilev), ! temperature and pressure
324	$ ztrac(ilon,ilev,i_h2o),
325	$ ztrac(ilon,ilev,i_h2so4),
326	$ ztrac(ilon,ilev,i_m0_aer),
327	$ ztrac(ilon,ilev,i_m3_aer),
328	$ ztrac(ilon,ilev,i_m0_mode1drop),
329	$ ztrac(ilon,ilev,i_m0_mode1ccn),
330	$ ztrac(ilon,ilev,i_m3_mode1sa),
331	$ ztrac(ilon,ilev,i_m3_mode1w),
332	$ ztrac(ilon,ilev,i_m3_mode1ccn),
333	$ ztrac(ilon,ilev,i_m0_mode2drop),
334	$ ztrac(ilon,ilev,i_m0_mode2ccn),
335	$ ztrac(ilon,ilev,i_m3_mode2sa),
336	$ ztrac(ilon,ilev,i_m3_mode2w),
337	$ ztrac(ilon,ilev,i_m3_mode2ccn))
338	else
339	ztrac(ilon,ilev,i_m0_aer) = 0.
340	ztrac(ilon,ilev,i_m3_aer) = 0.
341	ztrac(ilon,ilev,i_m0_mode1drop) = 0.
342	ztrac(ilon,ilev,i_m0_mode1ccn) = 0.
343	ztrac(ilon,ilev,i_m3_mode1sa) = 0.
344	ztrac(ilon,ilev,i_m3_mode1w) = 0.
345	ztrac(ilon,ilev,i_m3_mode1ccn) = 0.
346	ztrac(ilon,ilev,i_m0_mode2drop) = 0.
347	ztrac(ilon,ilev,i_m0_mode2ccn) = 0.
348	ztrac(ilon,ilev,i_m3_mode2sa) = 0.
349	ztrac(ilon,ilev,i_m3_mode2w) = 0.
350	ztrac(ilon,ilev,i_m3_mode2ccn) = 0.
351	end if
352	end do
353	end do
354
355	end if ! detailed scheme
356
357	!===================================================================
358	! photochemistry
359	!===================================================================
360
361	if (ok_chem) then
362
363	lon_sun = (0.5 - gmtime)2.rpi
364	lon_local(:) = lon(:)*rpi/180.
365	lat_local(:) = lat(:)*rpi/180.
366
367	if (ok_ionchem) then
368
369	vmr_dens_euv(:,:,ix_co2) = ztrac(:,:,i_co2) ! CO2
370	vmr_dens_euv(:,:,ix_co) = ztrac(:,:,i_co) ! CO
371	vmr_dens_euv(:,:,ix_o) = ztrac(:,:,i_o) ! O
372	vmr_dens_euv(:,:,ix_o1d) = ztrac(:,:,i_o1d) ! O(1D)
373	vmr_dens_euv(:,:,ix_o2) = ztrac(:,:,i_o2) ! O2
374	vmr_dens_euv(:,:,ix_o3) = ztrac(:,:,i_o3) ! O3
375	vmr_dens_euv(:,:,ix_h) = ztrac(:,:,i_h) ! H
376	vmr_dens_euv(:,:,ix_h2) = ztrac(:,:,i_h2) ! H2
377	vmr_dens_euv(:,:,ix_oh) = ztrac(:,:,i_oh) ! OH
378	vmr_dens_euv(:,:,ix_ho2) = ztrac(:,:,i_ho2) ! HO2
379	vmr_dens_euv(:,:,ix_h2o2)= ztrac(:,:,i_h2o2) ! H2O2
380	vmr_dens_euv(:,:,ix_h2o) = ztrac(:,:,i_h2o) ! H2O
381	vmr_dens_euv(:,:,ix_n) = ztrac(:,:,i_n) ! N
382	vmr_dens_euv(:,:,ix_n2d) = ztrac(:,:,i_n2d) ! N(2D)
383	vmr_dens_euv(:,:,ix_no) = ztrac(:,:,i_no) ! NO
384	vmr_dens_euv(:,:,ix_no2) = ztrac(:,:,i_no2) ! NO2
385	vmr_dens_euv(:,:,ix_n2) = ztrac(:,:,i_n2) ! N2
386
387	end if
388
389	do ilon = 1,nlon
390	zlocal(:)=zzlay(ilon,:)/1000.
391
392	! solar zenith angle
393	! sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))
394	! $ *cos(lon_sun) + cos(lat_local(ilon))
395	! $ sin(lon_local(ilon))sin(lon_sun))*180./rpi
396
397	sza_local = cos(lat_local(ilon))*cos(lon_local(ilon))
398	$ *cos(lon_sun) + cos(lat_local(ilon))
399	$ sin(lon_local(ilon))sin(lon_sun)
400
401	! Security - Handle rare cases where \|sza_local\| > 1
402	sza_local = min(sza_local,1.)
403	sza_local = max(-1.,sza_local)
404	sza_local = acos(sza_local)*180./rpi
405
406	! electron temperature
407	do ilev = 1, nlev
408	t_elec(ilev) = temp_elect(zlocal(ilev),temp(ilon,ilev),
409	$ sza_local,t_elec_origin)
410	end do
411
412	call photochemistry_venus(nlev, nlon, zlocal, pdtphys,
413	$ ok_jonline,ok_ionchem,tuneupperatm,
414	$ nb_reaction_3_max,nb_reaction_4_max,
415	$ nb_phot_max,nphotion,ilon,
416	$ pplay(ilon,:)/100.,
417	$ temp(ilon,:), t_elec(:),
418	$ ztrac(ilon,:,:),vmr_dens_euv(ilon,:,:),
419	$ mmean(ilon,:),
420	$ sza_local,
421	$ lon(ilon), lat(ilon),
422	$ nqmax, nespeuv, iter(ilon,:),
423	$ prod_tr(ilon,:,:),
424	$ loss_tr(ilon,:,:),
425	$ no_em, o2_em)
426
427	no_emi(ilon,:) = no_em(:)
428	o2_emi(ilon,:) = o2_em(:)
429
430	end do
431	end if ! ok_chem
432
433	!===================================================================
434	! compute tendencies
435	!===================================================================
436
437	! update mmean
438
439	mmean(:,:) = 0.
440	do iq = 1,nqmax - nmicro
441	mmean(:,:) = mmean(:,:)+ztrac(:,:,iq)*m_tr(iq)
442	end do
443	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
444
445	!===================================================================
446	! convert volume to mass mixing ratio / then tendencies in mmr
447	!===================================================================
448
449	! gas phase
450
451	do iq = 1,nqmax - nmicro
452	ztrac(:,:,iq) = max(ztrac(:,:,iq)*m_tr(iq)/mmean(:,:),
453	$ 1.e-30)
454	d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
455	end do
456
457	! liquid phase or moments
458
459	if (ok_cloud .and. cl_scheme == 1) then
460	ztrac(:,:,i_h2so4liq) = max(ztrac(:,:,i_h2so4liq)
461	$ *m_tr(i_h2so4liq)/mmean(:,:),
462	$ 1.e-30)
463	ztrac(:,:,i_h2oliq) = max(ztrac(:,:,i_h2oliq)
464	$ *m_tr(i_h2oliq)/mmean(:,:),
465	$ 1.e-30)
466	d_tr_chem(:,:,i_h2so4liq) = (ztrac(:,:,i_h2so4liq)
467	$ - trac(:,:,i_h2so4liq))/pdtphys
468	d_tr_chem(:,:,i_h2oliq) = (ztrac(:,:,i_h2oliq)
469	$ - trac(:,:,i_h2oliq))/pdtphys
470	end if
471
472
473	if (ok_cloud .and. cl_scheme == 2) then
474	do iq = nqmax-nmicro+1,nqmax
475	d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
476	end do
477	end if
478
479	end subroutine phytrac_chimie

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

Download in other formats:

Original Format