Context Navigation

calchim_mod.F90 @ 2808

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

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

Download in other formats:

Original Format