Context Navigation

calchim_mod.F90 @ 3289

Last change on this file since 3289 was 3289, checked in by jliu, 10 months ago

21/03/2024 == Jliu

Fix several bugs in photochemistry related to the submit-crash-resubmit
problem: Some of the integers in the related routines are not set in
privatethreads. Consequently, the additives such as n=n+1 are accumulated with
threads, causing systermatic problems in the routine. This update fixed the
submit-crash-resubmit problem. The simulated results are same with previous
simulations as tested.

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

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

Download in other formats:

Original Format