Context Navigation

calchim_mod.F90 @ 2325

Last change on this file since 2325 was 2321, checked in by emillour, 5 years ago

Mars GCM:

Finalized water ion chemisty (added H3O+ and OH+ ions). Added an example of relevent traceur.def file in deftank.
Added the computation of NO nightglow.

FGG

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

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

Download in other formats:

Original Format