Context Navigation

calchim_mod.F90 @ 2573

Last change on this file since 2573 was 2563, checked in by emillour, 4 years ago

Mars GCM:
Improve wstats: the callstats flag is now embedded within wstats and checked
internaly so no need to have some "if (callstats)" conditions around calls
to wstats anymore.
In addition: wstats now looks for an (optional) stats.def file in the
directory where the GCM is run to know which variable should be included
in the stats.nc file. The stats.def ASCII file should simply contain
one variable name per line, in the same way as the diagfi.def file for
diagfi outputs. If there is no stats.def file then all variables sent to
wstats will be in the stats.nc file (which matches the behaviour prior to
this improvement).
EM

File size: 37.9 KB

Line
1	MODULE calchim_mod
2
3	IMPLICIT NONE
4
5	INTEGER, SAVE :: ichemistry ! compute chemistry every ichemistry physics step
6	REAL,SAVE,ALLOCATABLE :: zdqchim(:,:,:) ! Tendancy on pq due to photochemistry
7	REAL,SAVE,ALLOCATABLE :: zdqschim(:,:) ! Tendancy on qsurf due to photochemistry
8
9	CONTAINS
10
11	subroutine calchim(ngrid,nlayer,nq, &
12	ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, &
13	zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud, &
14	dqscloud,tau,co2ice, &
15	pu,pdu,pv,pdv,surfdust,surfice)
16
17	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
18	igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2, &
19	igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap, &
20	igcm_no, igcm_n, igcm_no2, igcm_n2d, &
21	igcm_o2plus, igcm_co2plus, igcm_oplus, &
22	igcm_coplus, igcm_cplus, igcm_nplus, &
23	igcm_noplus, igcm_n2plus, igcm_hplus, &
24	igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, &
25	igcm_h3oplus, igcm_ohplus, igcm_elec, &
26	igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, &
27	igcm_do2, igcm_hdo2, mmol
28
29	use conc_mod, only: mmean ! mean molecular mass of the atmosphere
30	use comcstfi_h, only: pi
31	use photolysis_mod, only: init_photolysis, nphot
32	use iono_h, only: temp_elect
33	use wstats_mod, only: wstats
34
35	implicit none
36
37	!=======================================================================
38	!
39	! subject:
40	! --------
41	!
42	! Prepare the call for the photochemical module, and send back the
43	! tendencies from photochemistry in the chemical species mass mixing ratios
44	!
45	! Arguments:
46	! ----------
47	!
48	! Input:
49	!
50	! ptimestep timestep (s)
51	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
52	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa)
53	! pt(ngrid,nlayer) Temperature (K)
54	! pdt(ngrid,nlayer) Temperature tendency (K)
55	! pu(ngrid,nlayer) u component of the wind (ms-1)
56	! pdu(ngrid,nlayer) u component tendency (K)
57	! pv(ngrid,nlayer) v component of the wind (ms-1)
58	! pdv(ngrid,nlayer) v component tendency (K)
59	! dist_sol distance of the sun (AU)
60	! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith)
61	! pq(ngrid,nlayer,nq) advected fields, ie chemical species here
62	! pdq(ngrid,nlayer,nq) previous tendencies on pq
63	! tau(ngrid) dust optical depth
64	! co2ice(ngrid) co2 ice surface layer (kg.m-2)
65	! surfdust(ngrid,nlayer) dust surface area (m2/m3)
66	! surfice(ngrid,nlayer) ice surface area (m2/m3)
67	!
68	! Output:
69	!
70	! dqchim(ngrid,nlayer,nq) tendencies on pq due to chemistry
71	! dqschim(ngrid,nq) tendencies on qsurf
72	!
73	!=======================================================================
74
75	include "callkeys.h"
76
77	! input:
78
79	integer,intent(in) :: ngrid ! number of atmospheric columns
80	integer,intent(in) :: nlayer ! number of atmospheric layers
81	integer,intent(in) :: nq ! number of tracers
82	real :: ptimestep
83	real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers
84	real :: zzlay(ngrid,nlayer) ! pressure at the middle of the layers
85	real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels
86	real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
87	real :: pt(ngrid,nlayer) ! temperature
88	real :: pdt(ngrid,nlayer) ! temperature tendency
89	real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1)
90	real :: pdu(ngrid,nlayer) ! u component tendency
91	real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1)
92	real :: pdv(ngrid,nlayer) ! v component tendency
93	real :: dist_sol ! distance of the sun (AU)
94	real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith)
95	real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio
96	real :: pdq(ngrid,nlayer,nq) ! previous tendencies
97	real :: zday ! date (time since Ls=0, in martian days)
98	real :: tau(ngrid) ! dust optical depth
99	real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
100	real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
101	real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3)
102
103	! output:
104
105	real :: dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
106	real :: dqschim(ngrid,nq) ! tendencies on qsurf
107	real :: dqcloud(ngrid,nlayer,nq) ! tendencies on pq due to condensation
108	real :: dqscloud(ngrid,nq) ! tendencies on qsurf
109
110	! local variables:
111
112	integer,save :: nbq ! number of tracers used in the chemistry
113	integer,allocatable,save :: niq(:) ! array storing the indexes of the tracers
114	integer :: iloc(1) ! index of major species
115	integer :: ig,l,i,iq,iqmax
116	integer :: foundswitch, lswitch
117	integer,save :: chemthermod
118
119	integer,save :: i_co2 = 0
120	integer,save :: i_co = 0
121	integer,save :: i_o = 0
122	integer,save :: i_o1d = 0
123	integer,save :: i_o2 = 0
124	integer,save :: i_o3 = 0
125	integer,save :: i_h = 0
126	integer,save :: i_h2 = 0
127	integer,save :: i_oh = 0
128	integer,save :: i_ho2 = 0
129	integer,save :: i_h2o2 = 0
130	integer,save :: i_ch4 = 0
131	integer,save :: i_n2 = 0
132	integer,save :: i_h2o = 0
133	integer,save :: i_n = 0
134	integer,save :: i_no = 0
135	integer,save :: i_no2 = 0
136	integer,save :: i_n2d = 0
137	integer,save :: i_co2plus=0
138	integer,save :: i_oplus=0
139	integer,save :: i_o2plus=0
140	integer,save :: i_coplus=0
141	integer,save :: i_cplus=0
142	integer,save :: i_nplus=0
143	integer,save :: i_noplus=0
144	integer,save :: i_n2plus=0
145	integer,save :: i_hplus=0
146	integer,save :: i_hco2plus=0
147	integer,save :: i_hcoplus=0
148	integer,save :: i_h2oplus=0
149	integer,save :: i_h3oplus=0
150	integer,save :: i_ohplus=0
151	integer,save :: i_elec=0
152	integer,save :: i_hdo=0
153	integer,save :: i_od=0
154	integer,save :: i_d=0
155	integer,save :: i_hd=0
156	integer,save :: i_do2=0
157	integer,save :: i_hdo2=0
158
159	integer :: ig_vl1
160
161	integer :: nb_reaction_3_max ! number of quadratic reactions
162	integer :: nb_reaction_4_max ! number of bimolecular reactions
163	integer :: nquench ! number of quenching + heterogeneous reactions
164	integer :: nphotion ! number of photoionizations
165	integer :: nb_reaction_4_ion ! quadratic reactions for ionosphere
166	integer :: nb_reaction_4_deut ! quadratic reactions for deuterium chem
167	integer :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions
168
169
170	real :: latvl1, lonvl1
171	real :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry
172	! new mole fraction after
173	real :: zt(ngrid,nlayer) ! temperature
174	real :: zu(ngrid,nlayer) ! u component of the wind
175	real :: zv(ngrid,nlayer) ! v component of the wind
176	real :: taucol ! dust optical depth at the surface
177	real :: kb ! boltzmann constant
178
179	logical, save :: firstcall = .true.
180	logical, save :: depos ! switch for dry deposition
181	logical, save :: ionchem ! switch for ionospheric chemistry
182	logical, save :: deutchem ! switch for deuterium chemistry
183	logical, save :: jonline ! switch for online photodissociation rates or lookup table
184	logical, save :: unichim ! only one unified chemical scheme at all
185	! layers (default), or upper atmospheric
186	! scheme in the thermosphere
187
188	! for each column of atmosphere:
189
190	real :: zpress(nlayer) ! Pressure (mbar)
191	real :: zdens(nlayer) ! Density (cm-3)
192	real :: ztemp(nlayer) ! Temperature (K)
193	real :: ztelec(nlayer) ! Electronic temperature (K)
194	real :: zlocal(nlayer) ! Altitude (km)
195	real :: zycol(nlayer,nq) ! Composition (mole fractions)
196	real :: zmmean(nlayer) ! Mean molecular mass (g.mole-1)
197	real :: szacol ! Solar zenith angle
198	real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3)
199	real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3)
200	real :: jo3(nlayer) ! Photodissociation rate O3->O1D (s-1)
201	real :: jh2o(nlayer) ! Photodissociation rate H2O->H+OH (s-1)
202	real :: em_no(nlayer) ! NO nightglow emission rate
203	real :: em_o2(nlayer) ! O2 nightglow emission rate
204
205	integer :: iter(nlayer) ! Number of chemical iterations
206	! within one physical timestep
207
208	! for output:
209
210	logical,save :: output ! to issue calls to writediagfi and stats
211	real :: jo3_3d(ngrid,nlayer) ! Photodissociation rate O3->O1D (s-1)
212	real :: jh2o_3d(ngrid,nlayer) ! Photodissociation rate H2O->H+OH (s-1)
213	real :: emission_no(ngrid,nlayer) !NO emission rate
214	real :: emission_o2(ngrid,nlayer) !O2 emission rate
215	real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
216	! within one physical timestep
217
218	!=======================================================================
219	! initialization of the chemistry (first call only)
220	!=======================================================================
221
222	if (firstcall) then
223
224	!=======================================================================
225	! main dashboard for the chemistry
226	!=======================================================================
227
228	unichim = .true. ! true : unified chemistry ! false : separate models in lower and upper atmosphere
229	jonline = .true. ! true : on-line calculation of photodissociation rates ! false : lookup table
230	ionchem = .false. ! switch for ionospheric chemistry
231	deutchem= .false. ! switch for deuterium chemistry
232	depos = .false. ! switch for dry deposition
233	output = .true. ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc)
234
235	! if (photochem) then
236	! if (jonline) then
237	! print*,'calchim: Read UV absorption cross-sections'
238	! call init_photolysis ! on-line photolysis
239	! else
240	! print*,'calchim: Read photolysis lookup table'
241	! call read_phototable ! off-line photolysis
242	! end if
243	! end if
244
245	! if(.not.unichim) then
246	!Read reaction rates from external file if the upper atmospheric
247	!chemistry is called
248	! call chemthermos_readini
249	! endif
250
251	! find index of chemical tracers to use
252	allocate(niq(nq))
253	! Listed here are all tracers that can go into photochemistry
254	nbq = 0 ! to count number of tracers
255	! Species ALWAYS present if photochem=.T.
256	i_co2 = igcm_co2
257	if (i_co2 == 0) then
258	write(,) "calchim: Error; no CO2 tracer !!!"
259	call abort_physic("calchim","missing co2 tracer",1)
260	else
261	nbq = nbq + 1
262	niq(nbq) = i_co2
263	end if
264	i_co = igcm_co
265	if (i_co == 0) then
266	write(,) "calchim: Error; no CO tracer !!!"
267	call abort_physic("calchim","missing co tracer",1)
268	else
269	nbq = nbq + 1
270	niq(nbq) = i_co
271	end if
272	i_o = igcm_o
273	if (i_o == 0) then
274	write(,) "calchim: Error; no O tracer !!!"
275	call abort_physic("calchim","missing o tracer",1)
276	else
277	nbq = nbq + 1
278	niq(nbq) = i_o
279	end if
280	i_o1d = igcm_o1d
281	if (i_o1d == 0) then
282	write(,) "calchim: Error; no O1D tracer !!!"
283	call abort_physic("calchim","missing o1d tracer",1)
284	else
285	nbq = nbq + 1
286	niq(nbq) = i_o1d
287	end if
288	i_o2 = igcm_o2
289	if (i_o2 == 0) then
290	write(,) "calchim: Error; no O2 tracer !!!"
291	call abort_physic("calchim","missing o2 tracer",1)
292	else
293	nbq = nbq + 1
294	niq(nbq) = i_o2
295	end if
296	i_o3 = igcm_o3
297	if (i_o3 == 0) then
298	write(,) "calchim: Error; no O3 tracer !!!"
299	call abort_physic("calchim","missing o3 tracer",1)
300	else
301	nbq = nbq + 1
302	niq(nbq) = i_o3
303	end if
304	i_h = igcm_h
305	if (i_h == 0) then
306	write(,) "calchim: Error; no H tracer !!!"
307	call abort_physic("calchim","missing h tracer",1)
308	else
309	nbq = nbq + 1
310	niq(nbq) = i_h
311	end if
312	i_h2 = igcm_h2
313	if (i_h2 == 0) then
314	write(,) "calchim: Error; no H2 tracer !!!"
315	call abort_physic("calchim","missing h2 tracer",1)
316	else
317	nbq = nbq + 1
318	niq(nbq) = i_h2
319	end if
320	i_oh = igcm_oh
321	if (i_oh == 0) then
322	write(,) "calchim: Error; no OH tracer !!!"
323	call abort_physic("calchim","missing oh tracer",1)
324	else
325	nbq = nbq + 1
326	niq(nbq) = i_oh
327	end if
328	i_ho2 = igcm_ho2
329	if (i_ho2 == 0) then
330	write(,) "calchim: Error; no HO2 tracer !!!"
331	call abort_physic("calchim","missing ho2 tracer",1)
332	else
333	nbq = nbq + 1
334	niq(nbq) = i_ho2
335	end if
336	i_h2o2 = igcm_h2o2
337	if (i_h2o2 == 0) then
338	write(,) "calchim: Error; no H2O2 tracer !!!"
339	call abort_physic("calchim","missing h2o2 tracer",1)
340	else
341	nbq = nbq + 1
342	niq(nbq) = i_h2o2
343	end if
344	i_ch4 = igcm_ch4
345	if (i_ch4 == 0) then
346	write(,) "calchim: Error; no CH4 tracer !!!"
347	write(,) "CH4 will be ignored in the chemistry"
348	else
349	nbq = nbq + 1
350	niq(nbq) = i_ch4
351	end if
352	i_n2 = igcm_n2
353	if (i_n2 == 0) then
354	write(,) "calchim: Error; no N2 tracer !!!"
355	call abort_physic("calchim","missing n2 tracer",1)
356	else
357	nbq = nbq + 1
358	niq(nbq) = i_n2
359	end if
360	i_n = igcm_n
361	if (i_n == 0) then
362	if (photochem) then
363	write(,) "calchim: Error; no N tracer !!!"
364	call abort_physic("calchim","missing n tracer",1)
365	end if
366	else
367	nbq = nbq + 1
368	niq(nbq) = i_n
369	end if
370	i_n2d = igcm_n2d
371	if (i_n2d == 0) then
372	if (photochem) then
373	write(,) "calchim: Error; no N2D tracer !!!"
374	call abort_physic("calchim","missing n2d tracer",1)
375	end if
376	else
377	nbq = nbq + 1
378	niq(nbq) = i_n2d
379	end if
380	i_no = igcm_no
381	if (i_no == 0) then
382	if (photochem) then
383	write(,) "calchim: Error; no NO tracer !!!"
384	call abort_physic("calchim","missing no tracer",1)
385	end if
386	else
387	nbq = nbq + 1
388	niq(nbq) = i_no
389	end if
390	i_no2 = igcm_no2
391	if (i_no2 == 0) then
392	if (photochem) then
393	write(,) "calchim: Error; no NO2 tracer !!!"
394	call abort_physic("calchim","missing no2 tracer",1)
395	end if
396	else
397	nbq = nbq + 1
398	niq(nbq) = i_no2
399	end if
400	i_h2o = igcm_h2o_vap
401	if (i_h2o == 0) then
402	write(,) "calchim: Error; no water vapor tracer !!!"
403	call abort_physic("calchim","missing h2o_vap tracer",1)
404	else
405	nbq = nbq + 1
406	niq(nbq) = i_h2o
407	end if
408	i_hdo=igcm_hdo_vap
409	if (i_hdo == 0) then
410	write(,) "calchim: no HDO tracer !!!"
411	! call abort_physic("calchim","missing hdo_vap tracer",1)
412	write(,) "No deuterium chemistry considered"
413	else
414	nbq = nbq + 1
415	niq(nbq) = i_hdo
416	deutchem = .true.
417	write(,) "calchim: HDO tracer found in traceur.def"
418	write(,) "Deuterium chemistry included"
419	end if
420	i_od=igcm_od
421	if(deutchem) then
422	if (i_od == 0) then
423	write(,) "calchim: Error, no OD tracer !!!"
424	write(,) "OD is needed if HDO is in traceur.def"
425	call abort_physic("calchim","missing od tracer",1)
426	else
427	nbq = nbq + 1
428	niq(nbq) = i_od
429	end if
430	else
431	if (i_oplus /= 0) then
432	write(,) "calchim: Error: OD is present, but HDO is not!!!"
433	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
434	call abort_physic("calchim","missing hdo tracer",1)
435	endif
436	endif
437	i_d=igcm_d
438	if(deutchem) then
439	if (i_d == 0) then
440	write(,) "calchim: Error, no D tracer !!!"
441	write(,) "D is needed if HDO is in traceur.def"
442	call abort_physic("calchim","missing d tracer",1)
443	else
444	nbq = nbq + 1
445	niq(nbq) = i_d
446	end if
447	else
448	if (i_d /= 0) then
449	write(,) "calchim: Error: D is present, but HDO is not!!!"
450	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
451	call abort_physic("calchim","missing hdo tracer",1)
452	endif
453	endif
454	i_hd=igcm_hd
455	if(deutchem) then
456	if (i_hd == 0) then
457	write(,) "calchim: Error, no HD tracer !!!"
458	write(,) "HD is needed if HDO is in traceur.def"
459	call abort_physic("calchim","missing hd tracer",1)
460	else
461	nbq = nbq + 1
462	niq(nbq) = i_hd
463	end if
464	else
465	if (i_hd /= 0) then
466	write(,) "calchim: Error: HD is present, but HDO is not!!!"
467	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
468	call abort_physic("calchim","missing hd tracer",1)
469	endif
470	endif
471	i_do2=igcm_do2
472	if(deutchem) then
473	if (i_do2 == 0) then
474	write(,) "calchim: Error, no DO2 tracer !!!"
475	write(,) "DO2 is needed if HDO is in traceur.def"
476	call abort_physic("calchim","missing do2 tracer",1)
477	else
478	nbq = nbq + 1
479	niq(nbq) = i_do2
480	end if
481	else
482	if (i_do2 /= 0) then
483	write(,) "calchim: Error: DO2 is present, but HDO is not!!!"
484	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
485	call abort_physic("calchim","missing do2 tracer",1)
486	endif
487	endif
488	i_hdo2=igcm_hdo2
489	if(deutchem) then
490	if (i_hdo2 == 0) then
491	write(,) "calchim: Error, no HDO2 tracer !!!"
492	write(,) "HDO2 is needed if HDO is in traceur.def"
493	call abort_physic("calchim","missing hdo2 tracer",1)
494	else
495	nbq = nbq + 1
496	niq(nbq) = i_hdo2
497	end if
498	else
499	if (i_hdo2 /= 0) then
500	write(,) "calchim: Error: HDO2 is present, but HDO is not!!!"
501	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
502	call abort_physic("calchim","missing hdo2 tracer",1)
503	endif
504	endif
505	i_o2plus = igcm_o2plus
506	if (i_o2plus == 0) then
507	write(,) "calchim: no O2+ tracer "
508	write(,) "Only neutral chemistry"
509	else
510	nbq = nbq + 1
511	niq(nbq) = i_o2plus
512	ionchem = .true.
513	write(,) "calchim: O2+ tracer found in traceur.def"
514	write(,) "Ion chemistry included"
515	end if
516	i_co2plus = igcm_co2plus
517	if(ionchem) then
518	if (i_co2plus == 0) then
519	write(,) "calchim: Error, no CO2+ tracer !!!"
520	write(,) "CO2+ is needed if O2+ is in traceur.def"
521	call abort_physic("calchim","missing co2plus tracer",1)
522	else
523	nbq = nbq + 1
524	niq(nbq) = i_co2plus
525	end if
526	else
527	if (i_co2plus /= 0) then
528	write(,) "calchim: Error: CO2+ is present, but O2+ is not!!!"
529	write(,) "Both must be in traceur.def if ion chemistry wanted"
530	call abort_physic("calchim","missing o2plus tracer",1)
531	endif
532	endif
533	i_oplus=igcm_oplus
534	if(ionchem) then
535	if (i_oplus == 0) then
536	write(,) "calchim: Error, no O+ tracer !!!"
537	write(,) "O+ is needed if O2+ is in traceur.def"
538	call abort_physic("calchim","missing oplus tracer",1)
539	else
540	nbq = nbq + 1
541	niq(nbq) = i_oplus
542	end if
543	else
544	if (i_oplus /= 0) then
545	write(,) "calchim: Error: O+ is present, but O2+ is not!!!"
546	write(,) "Both must be in traceur.def if ion chemistry wanted"
547	call abort_physic("calchim","missing o2plus tracer",1)
548	endif
549	endif
550	i_noplus=igcm_noplus
551	if(ionchem) then
552	if (i_noplus == 0) then
553	write(,) "calchim: Error, no NO+ tracer !!!"
554	write(,) "NO+ is needed if O2+ is in traceur.def"
555	call abort_physic("calchim","missing noplus tracer",1)
556	else
557	nbq = nbq + 1
558	niq(nbq) = i_noplus
559	end if
560	else
561	if (i_noplus /= 0) then
562	write(,) "calchim: Error: NO+ is present, but O2+ is not!!!"
563	write(,) "Both must be in traceur.def if ion chemistry wanted"
564	endif
565	endif
566	i_coplus=igcm_coplus
567	if(ionchem) then
568	if (i_coplus == 0) then
569	write(,) "calchim: Error, no CO+ tracer !!!"
570	write(,) "CO+ is needed if O2+ is in traceur.def"
571	call abort_physic("calchim","missing coplus tracer",1)
572	else
573	nbq = nbq + 1
574	niq(nbq) = i_coplus
575	end if
576	else
577	if (i_coplus /= 0) then
578	write(,) "calchim: Error: CO+ is present, but O2+ is not!!!"
579	write(,) "Both must be in traceur.def if ion chemistry wanted"
580	endif
581	endif
582	i_cplus=igcm_cplus
583	if(ionchem) then
584	if (i_cplus == 0) then
585	write(,) "calchim: Error, no C+ tracer !!!"
586	write(,) "C+ is needed if O2+ is in traceur.def"
587	call abort_physic("calchim","missing cplus tracer",1)
588	else
589	nbq = nbq + 1
590	niq(nbq) = i_cplus
591	end if
592	else
593	if (i_cplus /= 0) then
594	write(,) "calchim: Error: C+ is present, but O2+ is not!!!"
595	write(,) "Both must be in traceur.def if ion chemistry wanted"
596	endif
597	endif
598	i_n2plus=igcm_n2plus
599	if(ionchem) then
600	if (i_n2plus == 0) then
601	write(,) "calchim: Error, no N2+ tracer !!!"
602	write(,) "N2+ is needed if O2+ is in traceur.def"
603	call abort_physic("calchim","missing n2plus tracer",1)
604	else
605	nbq = nbq + 1
606	niq(nbq) = i_n2plus
607	end if
608	else
609	if (i_n2plus /= 0) then
610	write(,) "calchim: Error: N2+ is present, but O2+ is not!!!"
611	write(,) "Both must be in traceur.def if ion chemistry wanted"
612	endif
613	endif
614	i_nplus=igcm_nplus
615	if(ionchem) then
616	if (i_nplus == 0) then
617	write(,) "calchim: Error, no N+ tracer !!!"
618	write(,) "N+ is needed if O2+ is in traceur.def"
619	call abort_physic("calchim","missing nplus tracer",1)
620	else
621	nbq = nbq + 1
622	niq(nbq) = i_nplus
623	end if
624	else
625	if (i_nplus /= 0) then
626	write(,) "calchim: Error: N+ is present, but O2+ is not!!!"
627	write(,) "Both must be in traceur.def if ion chemistry wanted"
628	endif
629	endif
630	i_hplus=igcm_hplus
631	if(ionchem) then
632	if (i_hplus == 0) then
633	write(,) "calchim: Error, no H+ tracer !!!"
634	write(,) "H+ is needed if O2+ is in traceur.def"
635	call abort_physic("calchim","missing hplus tracer",1)
636	else
637	nbq = nbq + 1
638	niq(nbq) = i_hplus
639	end if
640	else
641	if (i_hplus /= 0) then
642	write(,) "calchim: Error: H+ is present, but O2+ is not!!!"
643	write(,) "Both must be in traceur.def if ion chemistry wanted"
644	endif
645	endif
646	i_hco2plus=igcm_hco2plus
647	if(ionchem) then
648	if (i_hco2plus == 0) then
649	write(,) "calchim: Error, no HCO2+ tracer !!!"
650	write(,) "HCO2+ is needed if O2+ is in traceur.def"
651	call abort_physic("calchim","missing hco2plus tracer",1)
652	else
653	nbq = nbq + 1
654	niq(nbq) = i_hco2plus
655	end if
656	else
657	if (i_hco2plus /= 0) then
658	write(,) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
659	write(,) "Both must be in traceur.def if ion chemistry wanted"
660	endif
661	endif
662	i_hcoplus=igcm_hcoplus
663	if(ionchem) then
664	if (i_hcoplus == 0) then
665	write(,) "calchim: Error, no HCO+ tracer !!!"
666	write(,) "HCO+ is needed if O2+ is in traceur.def"
667	call abort_physic("calchim","missing hcoplus tracer",1)
668	else
669	nbq = nbq + 1
670	niq(nbq) = i_hcoplus
671	end if
672	else
673	if (i_hcoplus /= 0) then
674	write(,) "calchim: Error: HCO+ is present, but O2+ is not!!!"
675	write(,) "Both must be in traceur.def if ion chemistry wanted"
676	endif
677	endif
678	i_h2oplus=igcm_h2oplus
679	if(ionchem) then
680	if (i_h2oplus == 0) then
681	write(,) "calchim: Error, no H2O+ tracer !!!"
682	write(,) "H2O+ is needed if O2+ is in traceur.def"
683	call abort_physic("calchim","missing h2oplus tracer",1)
684	else
685	nbq = nbq + 1
686	niq(nbq) = i_h2oplus
687	end if
688	else
689	if (i_h2oplus /= 0) then
690	write(,) "calchim: Error: H2O+ is present, but O2+ is not!!!"
691	write(,) "Both must be in traceur.def if ion chemistry wanted"
692	endif
693	endif
694	i_h3oplus=igcm_h3oplus
695	if(ionchem) then
696	if (i_h3oplus == 0) then
697	write(,) "calchim: Error, no H3O+ tracer !!!"
698	write(,) "H3O+ is needed if O2+ is in traceur.def"
699	call abort_physic("calchim","missing h3oplus tracer",1)
700	else
701	nbq = nbq + 1
702	niq(nbq) = i_h3oplus
703	end if
704	else
705	if (i_h3oplus /= 0) then
706	write(,) "calchim: Error: H3O+ is present, but O2+ is not!!!"
707	write(,) "Both must be in traceur.def if ion chemistry wanted"
708	endif
709	endif
710	i_ohplus=igcm_ohplus
711	if(ionchem) then
712	if (i_ohplus == 0) then
713	write(,) "calchim: Error, no OH+ tracer !!!"
714	write(,) "OH+ is needed if O2+ is in traceur.def"
715	call abort_physic("calchim","missing ohplus tracer",1)
716	else
717	nbq = nbq + 1
718	niq(nbq) = i_ohplus
719	end if
720	else
721	if (i_ohplus /= 0) then
722	write(,) "calchim: Error: OH+ is present, but O2+ is not!!!"
723	write(,) "Both must be in traceur.def if ion chemistry wanted"
724	endif
725	endif
726	i_elec = igcm_elec
727	if(ionchem) then
728	if (i_elec == 0) then
729	write(,) "calchim: Error, no e- tracer !!!"
730	write(,) "e- is needed if O2+ is in traceur.def"
731	call abort_physic("calchim","missing elec tracer",1)
732	else
733	nbq = nbq + 1
734	niq(nbq) = i_elec
735	end if
736	else
737	if (i_elec /= 0) then
738	write(,) "calchim: Error: e- is present, but O2+ is not!!!"
739	write(,) "Both must be in traceur.def if ion chemistry wanted"
740	endif
741	endif
742	write(,) 'calchim: found nbq = ',nbq,' tracers'
743
744	write(,) 'calchim: tracer indices=',niq(:)
745
746
747	if (photochem) then
748	if (jonline) then
749	print*,'calchim: Read UV absorption cross-sections'
750	!Add two photodissociations in deuterium chemistry included
751	if(deutchem) nphot = nphot + 2
752	call init_photolysis ! on-line photolysis
753	else
754	print*,'calchim: Read photolysis lookup table'
755	call read_phototable ! off-line photolysis
756	end if
757	end if
758
759	if(.not.unichim) then
760	!Read reaction rates from external file if the upper atmospheric
761	!chemistry is called
762	call chemthermos_readini
763	endif
764
765	! stop
766	firstcall = .false.
767	end if ! if (firstcall)
768
769	! Initializations
770
771	zycol(:,:) = 0.
772	dqchim(:,:,:) = 0.
773	dqschim(:,:) = 0.
774
775	kb = 1.3806e-23
776
777	! latvl1= 22.27
778	! lonvl1= -47.94
779	! ig_vl1= 1+ int( (1.5-(latvl1-90.)48./180.) -2 )64. + &
780	! int(1.5+(lonvl1+180)*64./360.)
781
782	!=======================================================================
783	! loop over grid
784	!=======================================================================
785
786	do ig = 1,ngrid
787
788	foundswitch = 0
789	do l = 1,nlayer
790	do i = 1,nbq
791	iq = niq(i) ! get tracer index
792	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
793	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
794	end do
795	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
796	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
797	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
798	zpress(l) = pplay(ig,l)/100.
799	ztemp(l) = zt(ig,l)
800	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
801	zlocal(l) = zzlay(ig,l)/1000.
802	zmmean(l) = mmean(ig,l)
803	ztelec(l) = temp_elect(zlocal(l),ztemp(l),1)
804	!Electronic temperature. Index 1 -> Viking; Index 2-> MAVEN
805
806	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
807
808	surfdust1d(l) = surfdust(ig,l)*1.e-2
809	surfice1d(l) = surfice(ig,l)*1.e-2
810
811	! search for switch index between regions
812
813	if (unichim) then
814	lswitch = nlayer + 1
815	else
816	if (foundswitch == 0 .and. pplay(ig,l) < 10.) then
817	lswitch = l
818	foundswitch = 1
819	end if
820	endif
821	end do ! of do l=1,nlayer
822
823	szacol = acos(mu0(ig))*180./pi
824	taucol = tau(ig)
825
826	!=======================================================================
827	! call chemical subroutines
828	!=======================================================================
829
830	if (photochem) then
831	! set number of reactions, depending on ion chemistry or not
832	nb_reaction_4_ion = 64
833	nb_reaction_4_deut = 35
834
835	!Default numbers if no ion and no deuterium chemistry included
836
837	nb_reaction_4_max = 31 ! set number of bimolecular reactions
838	nb_reaction_3_max = 6 ! set number of quadratic reactions
839	nquench = 9 ! set number of quenching + heterogeneous
840	nphotion = 0 ! set number of photoionizations
841
842	if (ionchem) then
843	nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_ion
844	nphotion = 18 ! set number of photoionizations
845	endif
846	if(deutchem) then
847	nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut
848	end if
849
850	! nb_phot_max is the total number of processes that are treated
851	! numerically as a photolysis:
852
853	nb_phot_max = nphot + nphotion + nquench
854
855	! call main photochemistry routine
856
857	call photochemistry(nlayer,nq,nbq,ionchem,deutchem, &
858	nb_reaction_3_max,nb_reaction_4_max,nphot, &
859	nb_phot_max,nphotion, &
860	jonline,ig,lswitch,zycol,szacol,ptimestep, &
861	zpress,zlocal,ztemp,ztelec,zdens,zmmean, &
862	dist_sol,zday,surfdust1d,surfice1d, &
863	jo3,jh2o,em_no,em_o2,taucol,iter)
864
865	! ozone photolysis, for output
866
867	do l = 1,nlayer
868	jo3_3d(ig,l) = jo3(l)
869	jh2o_3d(ig,l) = jh2o(l)
870	iter_3d(ig,l) = iter(l)
871	end do
872
873	! condensation of h2o2
874
875	call perosat(ngrid, nlayer, nq, &
876	ig,ptimestep,pplev,pplay, &
877	ztemp,zycol,dqcloud,dqscloud)
878
879	! case of separate photochemical model in the thermosphere
880
881	if (.not.unichim) then
882	chemthermod = 3 !C/O/H/N/ions
883	call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
884	zdens,zpress,zlocal,szacol,ptimestep,zday,&
885	em_no,em_o2)
886	end if
887
888	do l = 1,nlayer
889	emission_no(ig,l) = em_no(l)
890	emission_o2(ig,l) = em_o2(l)
891	end do
892	end if ! photochem
893
894	! dry deposition
895
896	if (depos) then
897	call deposition(ngrid, nlayer, nq, &
898	ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
899	zu, zv, zt, zycol, ptimestep, co2ice)
900	end if
901
902	!=======================================================================
903	! tendencies
904	!=======================================================================
905
906	! index of the most abundant species at each level
907
908	! major(:) = maxloc(zycol, dim = 2)
909
910	! tendency for the most abundant species = - sum of others
911
912	do l = 1,nlayer
913	iloc=maxloc(zycol(l,:))
914	iqmax=iloc(1)
915	do i = 1,nbq
916	iq = niq(i) ! get tracer index
917	if (iq /= iqmax) then
918	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
919	- zq(ig,l,iq))/ptimestep
920	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
921	- dqchim(ig,l,iq)
922	end if
923	end do
924	end do ! of do l = 1,nlayer
925
926	!=======================================================================
927	! end of loop over grid
928	!=======================================================================
929
930	end do ! of do ig=1,ngrid
931
932	!=======================================================================
933	! write outputs
934	!=======================================================================
935
936	! value of parameter 'output' to trigger writting of outputs
937	! is set above at the declaration of the variable.
938
939	if (photochem .and. output) then
940	if (ngrid > 1) then
941	call writediagfi(ngrid,'jo3','j o3->o1d', &
942	's-1',3,jo3_3d(1,1))
943	call writediagfi(ngrid,'jh2o','jh2o', &
944	's-1',3,jh2o_3d(1,1))
945	call writediagfi(ngrid,'iter','iterations', &
946	' ',3,iter_3d(1,1))
947
948	! if (.not. unichim) then
949	call writediagfi(ngrid,'emission_no', &
950	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
951	call writediagfi(ngrid,'emission_o2', &
952	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
953	! endif
954
955	call wstats(ngrid,'jo3','j o3->o1d', &
956	's-1',3,jo3_3d(1,1))
957	call wstats(ngrid,'emission_no', &
958	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
959	call wstats(ngrid,'emission_o2', &
960	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
961	call wstats(ngrid,'mmean','mean molecular mass', &
962	'g.mole-1',3,mmean(1,1))
963	end if ! of if (ngrid.gt.1)
964	end if ! of if (output)
965
966	end subroutine calchim
967
968
969	subroutine ini_calchim_mod(ngrid,nlayer,nq)
970
971	implicit none
972
973	integer,intent(in) :: ngrid ! number of atmospheric columns
974	integer,intent(in) :: nlayer ! number of atmospheric layers
975	integer,intent(in) :: nq ! number of tracers
976
977	allocate(zdqchim(ngrid,nlayer,nq))
978	zdqchim(:,:,:)=0
979	allocate(zdqschim(ngrid,nq))
980	zdqschim(:,:)=0
981
982	end subroutine ini_calchim_mod
983
984
985	subroutine end_calchim_mod
986
987	implicit none
988
989	if (allocated(zdqchim)) deallocate(zdqchim)
990	if (allocated(zdqschim)) deallocate(zdqschim)
991
992	end subroutine end_calchim_mod
993
994	END MODULE calchim_mod
995

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

Download in other formats:

Original Format