Context Navigation

calchim_mod.F90 @ 2599

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

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

Download in other formats:

Original Format