Context Navigation

paramfoto_compact.F @ 3026

Last change on this file since 3026 was 2615, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars
File size: 175.4 KB

Line
1	c**********************************************************************
2
3	subroutine paramfoto_compact
4	$(ig,nlayer,chemthermod,lswitch,tx,timestep,zenit,zx,rm,nesptherm)
5
6	c Main thermospheric photochemistry routine.
7
8	c may 2008 FGG+MALV,GG
9	c**********************************************************************
10
11	use iono_h
12	use param_v4_h
13	implicit none
14
15	c arguments
16
17	integer lswitch,ig,nesptherm,chemthermod,nlayer
18	real zdens(nlayer)
19	real tx(nlayer)
20	real zenit
21	real zx(nlayer)
22	real rm(nlayer,nesptherm)
23	real timestep
24
25
26	c local variables
27
28	real*8 deltat,timefrac_sec
29	real*8 co2xnew,o2xnew,o3pxnew,coxnew,hxnew,ohxnew
30	real*8 ho2xnew,h2xnew
31	real*8 h2o2xnew,o1dxnew,o3xnew,h2oxnew
32	real*8 noxnew, nxnew, n2xnew, n2dxnew,no2xnew
33	real*8 oplusxnew, o2plusxnew,co2plusxnew
34	real*8 nplusxnew, n2plusxnew,noplusxnew, hplusxnew
35	real*8 electxnew,coplusxnew, cplusxnew, hco2plusxnew
36
37	real*8 co2xoutput,o2xoutput,o3pxoutput,coxoutput
38	real*8 ho2xoutput,h2xoutput,hxoutput,ohxoutput
39	real*8 h2o2xoutput,o1dxoutput,o3xoutput,h2oxoutput
40	real*8 nxoutput,noxoutput,n2xoutput,n2dxoutput,no2xoutput
41	real*8 co2plusxoutput,coplusxoutput,oplusxoutput,o2plusxoutput
42	real*8 cplusxoutput,noplusxoutput,n2plusxoutput,hplusxoutput
43	real*8 electxoutput,nplusxoutput,hco2plusxoutput
44	real*8 electxoutput_timemarching, electxoutput_neutrality
45
46	real*8 co2xinput,o2xinput,o3pxinput,coxinput
47	real*8 ho2xinput,h2xinput,hxinput,ohxinput
48	real*8 h2o2xinput,o1dxinput,o3xinput,h2oxinput
49	real*8 nxinput,noxinput,n2xinput,n2dxinput,no2xinput
50	real*8 co2plusxinput,coplusxinput,oplusxinput,o2plusxinput
51	real*8 cplusxinput,noplusxinput,n2plusxinput,hplusxinput
52	real*8 electxinput,nplusxinput,hco2plusxinput
53
54	real*8 co2xini,o2xini,o3pxini,coxini
55	real*8 ho2xini,h2xini,hxini,ohxini
56	real*8 h2o2xini,o1dxini,o3xini,h2oxini
57	real*8 nxini,noxini,n2xini,n2dxini,no2xini
58	real*8 co2plusxini,coplusxini,oplusxini,o2plusxini
59	real*8 cplusxini,noplusxini,n2plusxini,hplusxini
60	real*8 electxini,nplusxini,hco2plusxini
61
62	real*8 dco2x,do2x,do3px,dcox,dhx,dohx,dho2x,dh2x
63	real*8 dh2ox,dh2o2x,do1dx,do3x,dnx,dnox,dn2x,dn2dx,dno2x
64	real*8 dco2plusx,dcoplusx,doplusx,do2plusx
65	real*8 dcplusx,dnoplusx,dn2plusx,dhplusx,dhco2plusx
66	real*8 delectx,dnplusx
67
68	real*8 jdistot8(nabs,nlayer)
69	real*8 jdistot8_b(nabs,nlayer)
70	real*8 jion8(nabs,nlayer,4)
71	real*8 tx8
72
73
74
75	real*8 alfa_laststep, IonMostAbundant
76
77	real*8 tmin(nlayer)
78	real*8 fmargin1,critere
79
80	integer compmin(nlayer)
81	integer i,j,k
82	integer numpasos
83	integer n_comp_en_EQ(nlayer), paso
84
85	! Tracer indexes in the thermospheric chemistry:
86	!!! ATTENTION. These values have to be identical to those in chemthermos.F90
87	!!! If the values are changed there, the same has to be done here !!!
88	! integer,parameter :: i_co2=1
89	! integer,parameter :: i_o2=2
90	! integer,parameter :: i_o=3
91	! integer,parameter :: i_co=4
92	! integer,parameter :: i_h=5
93	! integer,parameter :: i_oh=6
94	! integer,parameter :: i_ho2=7
95	! integer,parameter :: i_h2=8
96	! integer,parameter :: i_h2o=9
97	! integer,parameter :: i_h2o2=10
98	! integer,parameter :: i_o1d=11
99	! integer,parameter :: i_o3=12
100	! integer,parameter :: i_n2=13
101	! integer,parameter :: i_n=14
102	! integer,parameter :: i_no=15
103	! integer,parameter :: i_n2d=16
104	! integer,parameter :: i_no2=17
105	integer,parameter :: i_co2 = 1
106	integer,parameter :: i_co = 2
107	integer,parameter :: i_o = 3
108	integer,parameter :: i_o1d = 4
109	integer,parameter :: i_o2 = 5
110	integer,parameter :: i_o3 = 6
111	integer,parameter :: i_h = 7
112	integer,parameter :: i_h2 = 8
113	integer,parameter :: i_oh = 9
114	integer,parameter :: i_ho2 = 10
115	integer,parameter :: i_h2o2 = 11
116	integer,parameter :: i_h2o = 12
117	integer,parameter :: i_n = 13
118	integer,parameter :: i_n2d = 14
119	integer,parameter :: i_no = 15
120	integer,parameter :: i_no2 = 16
121	integer,parameter :: i_n2 = 17
122	integer,parameter :: i_co2plus=18
123	integer,parameter :: i_oplus=19
124	integer,parameter :: i_o2plus=20
125	integer,parameter :: i_coplus=21
126	integer,parameter :: i_cplus=22
127	integer,parameter :: i_nplus=23
128	integer,parameter :: i_noplus=24
129	integer,parameter :: i_n2plus=25
130	integer,parameter :: i_hplus=26
131	integer,parameter :: i_hco2plus=27
132	integer,parameter :: i_elec=28
133
134	c formats
135
136	c**********************************************************************
137
138
139	c external timestep
140	timefrac_sec=dble(timestep)
141
142	C Start: altitude loop
143	do i=nlayer,lswitch,-1
144	c Temperature and concentrations to real*8
145	tx8=dble(tx(i))
146	co2xini=dble(rm(i,i_co2))
147	o2xini=dble(rm(i,i_o2))
148	o3pxini=dble(rm(i,i_o))
149	coxini=dble(rm(i,i_co))
150	hxini=dble(rm(i,i_h))
151	ohxini=dble(rm(i,i_oh))
152	ho2xini=dble(rm(i,i_ho2))
153	h2xini=dble(rm(i,i_h2))
154	h2oxini=dble(rm(i,i_h2o))
155	h2o2xini=dble(rm(i,i_h2o2))
156	o1dxini=dble(rm(i,i_o1d))
157	!Only if O3, N or ion chemistry requested
158	if(chemthermod.ge.1) o3xini=dble(rm(i,i_o3))
159	!Only if N or ion chemistry requested
160	if(chemthermod.ge.2) then
161	n2xini=dble(rm(i,i_n2))
162	nxini=dble(rm(i,i_n))
163	noxini=dble(rm(i,i_no))
164	n2dxini=dble(rm(i,i_n2d))
165	no2xini=dble(rm(i,i_no2))
166	endif
167	!Only if ion chemistry requested
168	if(chemthermod.eq.3) then
169	co2plusxini=dble(rm(i,i_co2plus))
170	oplusxini=dble(rm(i,i_oplus))
171	o2plusxini=dble(rm(i,i_o2plus))
172	coplusxini=dble(rm(i,i_coplus))
173	cplusxini=dble(rm(i,i_cplus))
174	nplusxini=dble(rm(i,i_nplus))
175	n2plusxini=dble(rm(i,i_n2plus))
176	noplusxini=dble(rm(i,i_noplus))
177	hplusxini=dble(rm(i,i_hplus))
178	hco2plusxini=dble(rm(i,i_hco2plus))
179	electxini=dble(rm(i,i_elec))
180	endif
181
182	!Calculation of photodissociation and photoionization rates
183	!from photoabsorption rates and ionization-to-dissociation
184	!branching ratios
185	call phdisrate(ig,nlayer,chemthermod,zenit,i)
186	! Conversion to double precision
187	do j=1,nabs
188	jdistot8(j,i) = dble(jdistot(j,i))
189	jdistot8_b(j,i) = dble(jdistot_b(j,i))
190	do k=1,4
191	jion8(j,i,k)=dble(jion(j,i,k))
192	enddo
193	end do
194
195	!Reaction rates
196	call getch( ig, chemthermod,tx8, zx(i))
197
198	!Lifetimes and temporal integration
199	call lifetimes(ig,i,nlayer,chemthermod,zenit,zx,
200	$ jdistot8,jdistot8_b,jion8,
201	$ tmin(i),compmin(i),
202	$ n_comp_en_EQ(i),co2xini,o2xini,o3pxini,coxini,hxini,
203	$ ohxini,ho2xini,h2xini,h2oxini,h2o2xini,o1dxini,o3xini,
204	$ n2xini,nxini,noxini,no2xini,n2dxini,co2plusxini,oplusxini,
205	$ o2plusxini,coplusxini,cplusxini,nplusxini,noplusxini,
206	$ n2plusxini,hplusxini,hco2plusxini,electxini )
207
208	!Calculation of the internal timestep and revision of the
209	!validity of the photochemical equilibrium approximation
210	!for each species
211
212	! JYC criteria added to avoid instabilities in (H) + (O+) <-> (H+) + (O) reactions when H+ is important
213	fmargin1=5
214	!Only if ion chemistry requested
215	if(chemthermod.eq.3) then
216	critere=hplusxini/(o3pxini+hxini+h2xini)
217	if (critere .gt. 5d-4) then
218	fmargin1=2000.*critere
219	if (fmargin1 .gt. 50.) fmargin1=50
220	endif
221	endif !Of chemthermod.eq.3
222
223	call timemarching ( ig,i,nlayer,chemthermod,n_comp_en_EQ,
224	. compmin,tmin,timefrac_sec, deltat,fmargin1)
225
226	!Number of timesteps
227	numpasos = int( timefrac_sec / deltat )
228	alfa_laststep = 1.d0 + timefrac_sec/deltat - dble(numpasos)
229	do paso=1,numpasos
230
231	!Concentrations at the first step
232	if(paso.eq.1) then
233	co2xinput=co2xini
234	o2xinput=o2xini
235	o3pxinput=o3pxini
236	coxinput=coxini
237	hxinput=hxini
238	ohxinput=ohxini
239	ho2xinput=ho2xini
240	h2xinput=h2xini
241	h2oxinput=h2oxini
242	h2o2xinput=h2o2xini
243	o1dxinput=o1dxini
244	o3xinput=o3xini
245	nxinput=nxini
246	noxinput=noxini
247	n2xinput=n2xini
248	n2dxinput=n2dxini
249	no2xinput=no2xini
250	!
251	co2plusxinput = co2plusxini
252	oplusxinput = oplusxini
253	o2plusxinput = o2plusxini
254	coplusxinput = coplusxini
255	cplusxinput = cplusxini
256	nplusxinput = nplusxini
257	n2plusxinput = n2plusxini
258	noplusxinput = noplusxini
259	hplusxinput = hplusxini
260	hco2plusxinput= hco2plusxini
261	electxinput = electxini
262	else
263	!Concentrations for the new step
264	co2xinput=co2xinput+dco2x
265	o2xinput=o2xinput+do2x
266	o3pxinput=o3pxinput+do3px
267	coxinput=coxinput+dcox
268	hxinput=hxinput+dhx
269	ohxinput=ohxinput+dohx
270	ho2xinput=ho2xinput+dho2x
271	h2xinput=h2xinput+dh2x
272	h2oxinput=h2oxinput+dh2ox
273	h2o2xinput=h2o2xinput+dh2o2x
274	o1dxinput=o1dxinput+do1dx
275	!Only if O3, N or ion chemistry requested
276	if(chemthermod.ge.1) o3xinput=o3xinput+do3x
277	!Only if N or ion chemistry requested
278	if(chemthermod.ge.2) then
279	nxinput=nxinput+dnx
280	noxinput=noxinput+dnox
281	n2xinput=n2xinput+dn2x
282	n2dxinput=n2dxinput+dn2dx
283	no2xinput=no2xinput+dno2x
284	endif
285	!Only if ion chemistry requested
286	if(chemthermod.eq.3) then
287	co2plusxinput = co2plusxinput + dco2plusx
288	oplusxinput = oplusxinput + doplusx
289	o2plusxinput = o2plusxinput + do2plusx
290	coplusxinput = coplusxinput + dcoplusx
291	cplusxinput = cplusxinput + dcplusx
292	nplusxinput = nplusxinput + dnplusx
293	n2plusxinput = n2plusxinput + dn2plusx
294	noplusxinput = noplusxinput + dnoplusx
295	hplusxinput = hplusxinput + dhplusx
296	hco2plusxinput= hco2plusxinput+ dhco2plusx
297	electxinput = electxinput + delectx
298	endif
299
300	end if
301	!Calculation of productions and losses
302	call prodsandlosses (ig,i,nlayer,chemthermod,zenit, zx,
303	& jdistot8, jdistot8_b, jion8,
304	& co2xinput, o2xinput, o3pxinput,
305	& coxinput, h2xinput, o3xinput,
306	& h2oxinput, nxinput, noxinput,
307	& h2o2xinput, n2xinput,
308	& o1dxinput, ohxinput, ho2xinput,
309	& hxinput, n2dxinput, no2xinput,
310	& co2plusxinput, o2plusxinput, coplusxinput,
311	& oplusxinput, cplusxinput, noplusxinput,
312	& n2plusxinput, hplusxinput, nplusxinput,
313	& hco2plusxinput,electxinput )
314
315
316	!New abundances, implicit scheme for the timemarching
317
318	!First, for the 11 species that can not be in PE
319	!( CO2, O2, O3P, CO, H2, H2O, H2O2, O3, N, NO, N2 )
320
321	call implicito ( ig, co2xoutput, ' CO2',
322	& co2xinput, Pco2tot(i), Lco2tot(i), deltat )
323	call implicito ( ig, o2xoutput, ' O2',
324	& o2xinput, Po2tot(i), Lo2tot(i), deltat )
325	call implicito ( ig, o3pxoutput, ' O3P',
326	& o3pxinput, Po3ptot(i), Lo3ptot(i), deltat )
327	call implicito ( ig, coxoutput, ' CO',
328	& coxinput, Pcotot(i), Lcotot(i), deltat )
329	call implicito ( ig, h2xoutput, ' H2',
330	& h2xinput, Ph2tot(i), Lh2tot(i), deltat )
331	call implicito ( ig, h2oxoutput, ' H2O',
332	& h2oxinput, Ph2otot(i), Lh2otot(i), deltat )
333	call implicito ( ig, h2o2xoutput, 'H2O2',
334	& h2o2xinput, Ph2o2tot(i), Lh2o2tot(i), deltat )
335	!only if O3, N or ion chemistry requested
336	if(chemthermod.ge.1)
337	$ call implicito ( ig, o3xoutput, ' O3',
338	& o3xinput, Po3tot(i), Lo3tot(i), deltat )
339	!Only if N or ion chemistry requested
340	if(chemthermod.ge.2) then
341	call implicito ( ig, nxoutput, ' N',
342	& nxinput, Pntot(i), Lntot(i), deltat )
343	call implicito ( ig, noxoutput, ' NO',
344	& noxinput, Pnotot(i), Lnotot(i), deltat )
345	call implicito ( ig, n2xoutput, ' N2',
346	& n2xinput, Pn2tot(i), Ln2tot(i), deltat )
347	endif
348
349
350	!Second, 6+10 species that can be in PE, but are not
351	! 6 neutral , O1D, OH, HO2, H, N2D, NO2
352	if(o1d_eq(i).eq.'N') then
353	call implicito ( ig, o1dxoutput, ' O1D',
354	& o1dxinput, Po1dtot(i), Lo1dtot(i), deltat )
355	end if
356	if(oh_eq(i).eq.'N') then
357	call implicito ( ig, ohxoutput, ' OH',
358	& ohxinput, Pohtot(i), Lohtot(i), deltat )
359	end if
360	if(ho2_eq(i).eq.'N') then
361	call implicito ( ig, ho2xoutput, ' HO2',
362	& ho2xinput, Pho2tot(i), Lho2tot(i), deltat )
363	end if
364	if(h_eq(i).eq.'N') then
365	call implicito ( ig, hxoutput, ' H',
366	& hxinput, Phtot(i), Lhtot(i), deltat )
367	end if
368	!Only if N or ion chemistry requested
369	if(chemthermod.ge.2) then
370	if(n2d_eq(i).eq.'N') then
371	call implicito ( ig, n2dxoutput, ' N2D',
372	& n2dxinput, Pn2dtot(i), Ln2dtot(i), deltat )
373	end if
374	if(no2_eq(i).eq.'N') then
375	call implicito ( ig, no2xoutput, ' NO2',
376	& no2xinput, Pno2tot(i), Lno2tot(i), deltat )
377	end if
378	endif
379
380	! 9 ions (all of them) and electrons
381	!Only if ion chemistry requested
382	if(chemthermod.ge.3) then
383	if(n2plus_eq(i).eq.'N') then
384	call implicito ( ig, n2plusxoutput, ' N2+',
385	& n2plusxinput,Pn2plustot(i),Ln2plustot(i),deltat)
386	end if
387	if(cplus_eq(i).eq.'N') then
388	call implicito ( ig, cplusxoutput, ' C+',
389	& cplusxinput,Pcplustot(i),Lcplustot(i),deltat)
390	end if
391	if(coplus_eq(i).eq.'N') then
392	call implicito ( ig, coplusxoutput, ' CO+',
393	& coplusxinput,Pcoplustot(i),Lcoplustot(i),deltat)
394	end if
395	if(co2plus_eq(i).eq.'N') then
396	call implicito ( ig, co2plusxoutput, 'CO2+',
397	& co2plusxinput,Pco2plustot(i),Lco2plustot(i),deltat)
398	end if
399	if(oplus_eq(i).eq.'N') then
400	call implicito ( ig, oplusxoutput, ' O+',
401	& oplusxinput,Poplustot(i),Loplustot(i),deltat)
402	end if
403	if(hplus_eq(i).eq.'N') then
404	call implicito ( ig, hplusxoutput, ' H+',
405	& hplusxinput,Phplustot(i),Lhplustot(i),deltat)
406	end if
407	if(o2plus_eq(i).eq.'N') then
408	call implicito ( ig, o2plusxoutput, ' O2+',
409	& o2plusxinput,Po2plustot(i),Lo2plustot(i),deltat)
410	end if
411	if(noplus_eq(i).eq.'N') then
412	call implicito ( ig, noplusxoutput, ' NO+',
413	& noplusxinput,Pnoplustot(i),Lnoplustot(i),deltat)
414	end if
415	if(nplus_eq(i).eq.'N') then
416	call implicito ( ig, nplusxoutput, ' N+',
417	& nplusxinput,Pnplustot(i),Lnplustot(i),deltat)
418	end if
419	if(hco2plus_eq(i).eq.'N') then
420	call implicito ( ig, hco2plusxoutput, 'CO2+',
421	& hco2plusxinput,Phco2plustot(i),Lhco2plustot(i),
422	$ deltat)
423	end if
424	! elect
425	call implicito ( ig, electxoutput_timemarching, 'elec',
426	& electxinput,Pelecttot(i),Lelecttot(i),deltat)
427	endif !Of chemthermod.eq.3
428
429
430	!Third, those species (among the 16 that can be in PE) that are in PE
431	call EF_oscilacion
432	& ( ig,i,nlayer, paso,chemthermod,zenit, zx,
433	& jdistot8, jdistot8_b,jion8,
434	& deltat,
435	$ co2xoutput, co2xinput,
436	$ o2xoutput, o2xinput,
437	$ o3pxoutput, o3pxinput,
438	$ coxoutput, coxinput,
439	$ h2xoutput, h2xinput,
440	$ h2oxoutput, h2oxinput,
441	$ h2o2xoutput, h2o2xinput,
442	$ o3xoutput, o3xinput,
443	$ nxoutput, nxinput,
444	$ noxoutput, noxinput,
445	$ n2xoutput, n2xinput,
446	& o1dxoutput, o1dxinput,
447	& ohxoutput, ohxinput,
448	& ho2xoutput, ho2xinput,
449	& hxoutput, hxinput,
450	& n2dxoutput, n2dxinput,
451	& no2xoutput, no2xinput,
452	& co2plusxoutput, co2plusxinput,
453	& o2plusxoutput, o2plusxinput,
454	& coplusxoutput, coplusxinput,
455	& oplusxoutput, oplusxinput,
456	& cplusxoutput, cplusxinput,
457	& noplusxoutput, noplusxinput,
458	& n2plusxoutput, n2plusxinput,
459	& hplusxoutput, hplusxinput,
460	& nplusxoutput, nplusxinput,
461	$ hco2plusxoutput,hco2plusxinput,
462	& electxoutput, electxinput,
463	& electxoutput_timemarching )
464
465	!Electrons given by the condition of global neutrality
466	!Only if ion chemistry requested
467	if(chemthermod.eq.3) then
468	electxoutput = o2plusxoutput +
469	@ co2plusxoutput +
470	@ coplusxoutput +
471	@ oplusxoutput +
472	@ cplusxoutput +
473	@ n2plusxoutput +
474	@ nplusxoutput +
475	@ noplusxoutput +
476	@ hplusxoutput +
477	$ hco2plusxoutput
478	electxoutput_neutrality = electxoutput
479	!
480	IonMostAbundant = o2plusxoutput
481	IonMostAbundant = max( co2plusxoutput, IonMostAbundant)
482	IonMostAbundant = max( coplusxoutput, IonMostAbundant)
483	IonMostAbundant = max( oplusxoutput, IonMostAbundant)
484	IonMostAbundant = max( cplusxoutput, IonMostAbundant)
485	IonMostAbundant = max( n2plusxoutput, IonMostAbundant)
486	IonMostAbundant = max( noplusxoutput, IonMostAbundant)
487	IonMostAbundant = max( nplusxoutput, IonMostAbundant)
488	IonMostAbundant = max( hplusxoutput, IonMostAbundant)
489	IonMostAbundant = max( hco2plusxoutput, IonMostAbundant)
490	IonMostAbundant = IonMostAbundant / electxoutput
491	endif !Of chemthermod.eq.3
492
493	!Concentration changes for this time step
494	dco2x=co2xoutput-co2xinput
495	do2x=o2xoutput-o2xinput
496	do3px=o3pxoutput-o3pxinput
497	dcox=coxoutput-coxinput
498	dhx=hxoutput-hxinput
499	dohx=ohxoutput-ohxinput
500	dho2x=ho2xoutput-ho2xinput
501	dh2x=h2xoutput-h2xinput
502	dh2ox=h2oxoutput-h2oxinput
503	dh2o2x=h2o2xoutput-h2o2xinput
504	do1dx=o1dxoutput-o1dxinput
505	!Only if O3, N or ion chemistry requested
506	if(chemthermod.ge.1) do3x=o3xoutput-o3xinput
507	!Only if N or ion chemistry requested
508	if(chemthermod.ge.2) then
509	dnx=nxoutput-nxinput
510	dnox=noxoutput-noxinput
511	dn2x=n2xoutput-n2xinput
512	dn2dx=n2dxoutput-n2dxinput
513	dno2x=no2xoutput-no2xinput
514	endif
515	!Only if ion chemistry requested
516	if(chemthermod.eq.3) then
517	dco2plusx=co2plusxoutput-co2plusxinput
518	do2plusx=o2plusxoutput-o2plusxinput
519	doplusx=oplusxoutput-oplusxinput
520	dcoplusx=coplusxoutput-coplusxinput
521	dcplusx=cplusxoutput-cplusxinput
522	dnplusx=nplusxoutput-nplusxinput
523	dn2plusx=n2plusxoutput-n2plusxinput
524	dnoplusx=noplusxoutput-noplusxinput
525	dhplusx=hplusxoutput-hplusxinput
526	dhco2plusx=hco2plusxoutput-hco2plusxinput
527	delectx=electxoutput- electxinput
528	endif
529	if(paso.eq.numpasos) then
530	!Final concentrations after last time step
531	co2xnew = co2xinput + dco2x * alfa_laststep
532	if(co2xnew.lt.0)co2xnew=1.e-30
533	o2xnew = o2xinput + do2x * alfa_laststep
534	if(o2xnew.lt.0)o2xnew=1.e-30
535	o3pxnew = o3pxinput + do3px * alfa_laststep
536	if(o3pxnew.lt.0)o3pxnew=1.e-30
537	coxnew = coxinput + dcox * alfa_laststep
538	if(coxnew.lt.0)coxnew=1.e-30
539	hxnew = hxinput + dhx * alfa_laststep
540	if(hxnew.lt.0)hxnew=1.e-30
541	ohxnew = ohxinput + dohx * alfa_laststep
542	if(ohxnew.lt.0)ohxnew=1.e-30
543	ho2xnew = ho2xinput + dho2x * alfa_laststep
544	if(ho2xnew.lt.0)ho2xnew=1.e-30
545	h2xnew = h2xinput + dh2x * alfa_laststep
546	if(h2xnew.lt.0)h2xnew=1.e-30
547	h2oxnew = h2oxinput + dh2ox * alfa_laststep
548	if(h2oxnew.lt.0)h2oxnew=1.e-30
549	h2o2xnew = h2o2xinput + dh2o2x * alfa_laststep
550	if(h2o2xnew.lt.0)h2o2xnew=1.e-30
551	o1dxnew = o1dxinput + do1dx * alfa_laststep
552	if(o1dxnew.lt.0)o1dxnew=1.e-30
553	!Only if O3, N or ion chemistry requested
554	if(chemthermod.ge.1) then
555	o3xnew = o3xinput + do3x * alfa_laststep
556	if(o3xnew.lt.0)o3xnew=1.e-30
557	endif
558	!Only if N or ion chemistry requested
559	if(chemthermod.ge.2) then
560	nxnew = nxinput + dnx * alfa_laststep
561	if(nxnew.lt.0)nxnew=1.e-30
562	noxnew = noxinput + dnox * alfa_laststep
563	if(noxnew.lt.0)noxnew=1.e-30
564	n2xnew = n2xinput + dn2x * alfa_laststep
565	if(n2xnew.lt.0)n2xnew=1.e-30
566	n2dxnew = n2dxinput + dn2dx * alfa_laststep
567	if(n2dxnew.lt.0)n2dxnew=1.e-30
568	no2xnew = no2xinput + dno2x * alfa_laststep
569	if(no2xnew.lt.0)no2xnew=1.e-30
570	endif
571	!Only if ion chemistry requested
572	if(chemthermod.ge.3) then
573	co2plusxnew = co2plusxinput+dco2plusx*alfa_laststep
574	if(co2plusxnew.lt.0)co2plusxnew=1.e-30
575	o2plusxnew = o2plusxinput+do2plusx*alfa_laststep
576	if(o2plusxnew.lt.0)o2plusxnew=1.e-30
577	oplusxnew = oplusxinput+doplusx*alfa_laststep
578	if(oplusxnew.lt.0)oplusxnew=1.e-30
579	coplusxnew = coplusxinput+dcoplusx*alfa_laststep
580	if(coplusxnew.lt.0)coplusxnew=1.e-30
581	nplusxnew = nplusxinput +dnplusx*alfa_laststep
582	if(nplusxnew.lt.0)nplusxnew=1.e-30
583	n2plusxnew = n2plusxinput+dn2plusx*alfa_laststep
584	if(n2plusxnew.lt.0)n2plusxnew=1.e-30
585	noplusxnew = noplusxinput+dnoplusx*alfa_laststep
586	if(noplusxnew.lt.0)noplusxnew=1.e-30
587	hplusxnew = hplusxinput+dhplusx*alfa_laststep
588	if(hplusxnew.lt.0)hplusxnew=1.e-30
589	cplusxnew = cplusxinput+dcplusx*alfa_laststep
590	if(cplusxnew.lt.0)cplusxnew=1.e-30
591	hco2plusxnew = hco2plusxinput+dhco2plusx*alfa_laststep
592	if(hco2plusxnew.lt.0)hco2plusxnew=1.e-30
593	electxnew = electxinput+delectx*alfa_laststep
594	if(electxnew.lt.0)electxnew=1.e-30
595	endif !Of chemthermod.ge.3
596	endif !Of paso.eq.numpasos
597
598
599	end do
600
601	!New concentrations to be returned
602	rm(i,i_co2) = real(co2xnew)
603	rm(i,i_o2) = real(o2xnew)
604	rm(i,i_o) = real(o3pxnew)
605	rm(i,i_co) = real(coxnew)
606	rm(i,i_h) = real(hxnew)
607	rm(i,i_oh) = real(ohxnew)
608	rm(i,i_ho2) = real(ho2xnew)
609	rm(i,i_h2) = real(h2xnew)
610	rm(i,i_h2o) = real(h2oxnew)
611	rm(i,i_h2o2) = real(h2o2xnew)
612	rm(i,i_o1d) = real(o1dxnew)
613	!Only if O3, N or ion chemistry requested
614	if(chemthermod.ge.1)
615	$ rm(i,i_o3) = real(o3xnew)
616	!Only if N or ion chemistry requested
617	if(chemthermod.ge.2) then
618	rm(i,i_n) = real(nxnew)
619	rm(i,i_n2) = real(n2xnew)
620	rm(i,i_no) = real(noxnew)
621	rm(i,i_n2d) = real(n2dxnew)
622	rm(i,i_no2) = real(no2xnew)
623	endif
624	!Only if ion chemistry requested
625	if(chemthermod.eq.3) then
626	rm(i,i_co2plus) = real(co2plusxnew)
627	rm(i,i_oplus) = real(oplusxnew)
628	rm(i,i_o2plus) = real(o2plusxnew)
629	rm(i,i_coplus) = real(coplusxnew)
630	rm(i,i_cplus) = real(cplusxnew)
631	rm(i,i_nplus) = real(nplusxnew)
632	rm(i,i_n2plus) = real(n2plusxnew)
633	rm(i,i_noplus) = real(noplusxnew)
634	rm(i,i_hplus) = real(hplusxnew)
635	rm(i,i_hco2plus)= real(hco2plusxnew)
636	rm(i,i_elec) = real(electxnew)
637	endif
638	end do
639	cccccc End altitude loop
640
641	return
642
643
644	end
645
646
647	c**********************************************************************
648	c**********************************************************************
649
650	subroutine implicito ( ig,c_output, text,
651	& c_input, Prod, Loss, tstep )
652
653
654	c Given the productions and losses, calculates new concentrations using
655	c an implicit integration scheme. Checks if there are negative values and
656	c avoids underflows.
657
658	c jul 2008 MA Version en subrutina
659	c**********************************************************************
660
661	implicit none
662
663	c arguments
664	c
665	integer ig
666	real*8 c_output ! O.
667	real*8 c_input ! I.
668	real*8 tstep ! I.
669	real*8 Prod ! I.
670	real*8 Loss ! I.
671	character*4 text ! I.
672
673	ccccccccccccccc CODE STARTS
674
675	c_output = (c_input + Prod * tstep) / (1.d0 + Loss * tstep)
676
677	! Stop is negative prods, losses, concentrations or times
678	!
679	if ( c_output.lt.0.d0 ) then
680	write(,) text//' < 0 !!!'
681	write (,) ' Terms of the implicit equation: '
682	write (,) ' c_input =', c_input
683	write (,) ' Prod = ', Prod
684	write (,) ' Loss = ', Loss
685	write (,) ' tstep = ', tstep
686	write (,) ' c_output =', c_output
687	write (,) ' ig = ', ig
688	stop ' Stop at IMPLICIT , PHCHEM '
689	endif
690
691	! Avoid underflow
692	!
693	if ( c_output.lt.1.d-30) c_output=1.d-30
694
695
696	return
697	c END
698	end
699
700
701
702	c***********************************************************************
703	function ionsec_nplus (zenit, alt)
704
705	c Calculates the N+ production by photoelectrons, following
706	c Nicholson et al. 2009
707
708	c FGG sep 2010 first version
709	c***********************************************************************
710
711	implicit none
712
713	! Arguments
714	real*8 ionsec_nplus
715	real zenit
716	real alt
717
718	! Local variables
719	real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14
720	real*8 b0,b1,b2,b3,b4
721	real*8 altaux
722	real*8 zenit_rad
723
724	!!!!!!! Program starts
725
726	zenit_rad=dble(zenit*3.141592/180.)
727
728	if(zenit.le.90.) then
729	altaux=dble(alt)+
730	$ 15.cos(zenit_rad)-40.sqrt(cos(zenit_rad))+25.
731	else
732	altaux=dble(alt)
733	endif
734
735	if(altaux.gt.108.4) then
736	a0 = 1.139925703d3
737	a1 = -4.742398258d1
738	a2 = 8.404232989d-1
739	a3 = -8.108229906d-3
740	a4 = 4.420892285d-5
741	a5 = -1.154901432d-7
742	a6 = -3.514073816d-11
743	a7 = 8.790819159d-13
744	a8 = -1.320788149d-16
745	a9 = -8.202233732d-18
746	a10 = -1.256480521d-22
747	a11 = 1.329336168e-22
748	a12 = -4.403185142d-25
749	a13 = 6.098474897d-28
750	a14 = -3.256951018d-31
751	ionsec_nplus = a0 + a1altaux + a2altaux*2 + a3altaux**3 +
752	$ a4altaux4 + a5altaux*5 + a6altaux*6 + a7altaux**7
753	$ + a8altaux8 + a9altaux*9 + a10altaux**10 +
754	$ a11altaux11 + a12altaux*12 + a13altaux**13 +
755	$ a14altaux*14
756	ionsec_nplus = 10**(ionsec_nplus-2.)
757	elseif(altaux.gt.80..and.altaux.le.108.4) then
758	b0 = 6.346190854d4
759	b1 = -2.623253212d3
760	b2 = 4.050319629d1
761	b3 = -2.767987276d-1
762	b4 = 7.064439029d-4
763	ionsec_nplus = b0 + b1altaux + b2altaux*2 + b3altaux**3 +
764	$ b4altaux*4
765	else
766	ionsec_nplus=0.d0
767	endif
768	if(ionsec_nplus.gt.100.d0.or.ionsec_nplus.lt.0.d0)
769	$ ionsec_nplus=0.d0
770
771
772	return
773
774	end
775
776
777	c***********************************************************************
778	function ionsec_n2plus (zenit, alt)
779
780	c N2+ production by photoelectrons, following Nicholson et al. 2009
781
782	c FGG sep 2010 first version
783	c***********************************************************************
784
785	implicit none
786
787	! Arguments
788	real*8 ionsec_n2plus
789	real zenit
790	real alt
791
792	! Local variables
793
794	real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9
795	real*8 b0,b1,b2,b3,b4,b5
796	real*8 altaux
797	real*8 zenit_rad
798
799	!!!!!!! Program starts
800
801	zenit_rad=dble(zenit*3.141592/180.)
802	if(zenit.le.90.) then
803	altaux=dble(alt)+
804	$ 15.cos(zenit_rad)-40.sqrt(cos(zenit_rad))+25.
805	else
806	altaux=dble(alt)
807	endif
808
809
810	if(altaux.gt.108.4) then
811	a0 = 9.843804026d2
812	a1 = -3.978696855d1
813	a2 = 7.028448262d-1
814	a3 = -7.11195117d-3
815	a4 = 4.545683986d-5
816	a5 = -1.905046447d-7
817	a6 = 5.240068127d-10
818	a7 = -9.130399894d-13
819	a8 = 9.151792207d-16
820	a9 = -4.023230948d-19
821	ionsec_n2plus = a0 + a1altaux + a2altaux*2 + a3altaux**3 +
822	$ a4altaux4 + a5altaux*5 + a6altaux**6 +
823	$ a7altaux7 + a8altaux*8 + a9altaux**9
824	ionsec_n2plus = 10**(ionsec_n2plus-2.)
825	elseif(altaux.gt.80..and.altaux.le.108.4) then
826	b0 = 5.146111566d4
827	b1 = -1.771736158d3
828	b2 = 1.811156914d1
829	b3 = 3.583722498d-3
830	b4 = -9.891151731d-4
831	b5 = 3.994474097d-6
832	ionsec_n2plus = b0 + b1altaux + b2altaux*2 + b3altaux**3 +
833	$ b4altaux4 + b5altaux**5
834	else
835	ionsec_n2plus = 0.d0
836	endif
837	if(ionsec_n2plus.gt.100.d0.or.ionsec_n2plus.lt.0.d0)
838	$ ionsec_n2plus=0.d0
839
840	return
841
842	end
843
844
845
846	c***********************************************************************
847	function ionsec_oplus (zenit, alt)
848
849	c O+ production by photoelectrons, from Nicholson et al. 2009
850
851	c FGG sep 2010 first version
852	c***********************************************************************
853
854	implicit none
855
856	! Arguments
857	real*8 ionsec_oplus
858	real zenit
859	real alt
860
861	! Local variables
862
863	real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12
864	real*8 b0,b1,b2,b3,b4,b5,b6
865	real*8 altaux
866	real*8 zenit_rad
867
868	!!!!!!! Program starts
869
870	zenit_rad=dble(zenit*3.141592/180.)
871
872	if(zenit.le.90.) then
873	altaux=dble(alt)+
874	$ 15.cos(zenit_rad)-40.sqrt(cos(zenit_rad))+25.
875	else
876	altaux=dble(alt)
877	endif
878
879	if(altaux.gt.112.9) then
880	a0 = 6.453740731d2
881	a1 = -2.547798991d1
882	a2 = 4.384613636d-1
883	a3 = -4.288387072d-3
884	a4 = 2.596437447d-5
885	a5 = -9.750300694d-8
886	a6 = 1.986722344d-10
887	a7 = -2.293667699d-14
888	a8 = -1.080547019d-15
889	a9 = 3.171787989d-18
890	a10 = -4.566493384d-21
891	a11 = 3.472393897d-24
892	a12 = -1.115699979d-27
893	ionsec_oplus = a0 + a1altaux + a2altaux*2 + a3altaux**3 +
894	$ a4altaux4 + a5altaux*5 + a6altaux*6 + a7altaux**7
895	$ + a8altaux8 + a9altaux*9 + a10altaux**10 +
896	$ a11altaux11 +a12altaux**12
897	ionsec_oplus = 10**(ionsec_oplus-2.)
898	elseif(altaux.gt.80..and.altaux.le.112.9) then
899	b0 = -5.934881676d5
900	b1 = 3.546095705d4
901	b2 = -8.806801303d2
902	b3 = 1.163735173d1
903	b4 = -8.62992438d-2
904	b5 = 3.40547333d-4
905	b6 = -5.587037506d-7
906	ionsec_oplus = b0 + b1altaux + b2altaux*2 + b3altaux**3 +
907	$ b4altaux4 + b5altaux*5 + b6altaux**6
908	else
909	ionsec_oplus=0.d0
910	endif
911
912	if(ionsec_oplus.gt.100.d0.or.ionsec_oplus.lt.0.d0)
913	$ ionsec_oplus=0.d0
914
915	return
916
917	end
918
919
920
921	c***********************************************************************
922	function ionsec_coplus (zenit, alt)
923
924	c CO+ production by photoelectrons from Nicholson et al. 2009
925
926	c FGG sep 2010 first version
927	c***********************************************************************
928
929	implicit none
930
931	! Arguments
932	real*8 ionsec_coplus
933	real zenit
934	real alt
935
936	! Local variables
937
938	real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9
939	real*8 b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
940	real*8 altaux
941	real*8 zenit_rad
942
943	!!!!!!! Program starts
944
945	zenit_rad=dble(zenit*3.141592/180.)
946
947	if(zenit.le.90.) then
948	altaux=dble(alt)+
949	$ 15.cos(zenit_rad)-40.sqrt(cos(zenit_rad))+25.
950	else
951	altaux=dble(alt)
952	endif
953
954	if(altaux.gt.110.6) then
955	a0 = 7.33258229d2
956	a1 = -2.919984139d1
957	a2 = 5.079651482d-1
958	a3 = -5.057170037d-3
959	a4 = 3.178156709d-5
960	a5 = -1.309076957d-7
961	a6 = 3.53858799d-10
962	a7 = -6.060315762d-13
963	a8 = 5.973573923d-16
964	a9 = -2.584454326d-19
965	ionsec_coplus = a0 + a1altaux + a2altaux*2 + a3altaux**3 +
966	$ a4altaux4 + a5altaux*5 + a6altaux*6 + a7altaux**7
967	$ + a8altaux8 + a9altaux**9
968	ionsec_coplus = 10**(ionsec_coplus-2.)
969	elseif(altaux.gt.80..and.altaux.le.110.6) then
970	b0 = -1.165107657d6
971	b1 = 4.315606169d4
972	b2 = -3.480483017d2
973	b3 = -3.831253024d0
974	b4 = 4.33316742d-2
975	b5 = 2.43075291d-4
976	b6 = -7.835732322d-8
977	b7 = -3.80859833d-8
978	b8 = -1.947628467d-10
979	b9 = 3.440753726d-12
980	b10 = 2.336227916d-14
981	b11 = -3.575877198d-16
982	b12 = 1.030801684d-18
983	ionsec_coplus = b0 + b1altaux + b2altaux*2 + b3altaux**3 +
984	$ b4altaux4 + b5altaux*5 + b6altaux**6 +
985	$ b7altaux7 + b8altaux*8 + b9altaux**9 +
986	$ b10altaux10 + b11altaux*11 + b12altaux**12
987	else
988	ionsec_coplus=0.d0
989	endif
990	if(ionsec_coplus.gt.100..or.ionsec_coplus.lt.0.d0)
991	$ ionsec_coplus=0.d0
992
993	return
994
995	end
996
997
998
999	c***********************************************************************
1000	function ionsec_co2plus (zenit, alt)
1001
1002	c CO2+ production by photoelectrons, from Nicholson et al. 2009
1003
1004	c FGG sep 2010 first version
1005	c***********************************************************************
1006
1007	implicit none
1008
1009	! Arguments
1010	real*8 ionsec_co2plus
1011	real zenit
1012	real alt
1013
1014	! Local variables
1015
1016	real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9
1017	real*8 b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10
1018	real*8 altaux
1019	real*8 zenit_rad
1020
1021	!!!!!!!Program starts
1022
1023	zenit_rad=dble(zenit*3.141592/180.)
1024
1025	if(zenit.le.90.) then
1026	altaux=dble(alt)+
1027	$ 15.cos(zenit_rad)-40.sqrt(cos(zenit_rad))+25.
1028	else
1029	altaux=dble(alt)
1030	endif
1031
1032	if(altaux.gt.112.9) then
1033	a0 = 8.64096192d2
1034	a1 = -3.471713960d1
1035	a2 = 6.072614479d-1
1036	a3 = -6.050002721d-3
1037	a4 = 3.779639483d-5
1038	a5 = -1.533626303d-7
1039	a6 = 4.032987841d-10
1040	a7 = -6.602964344d-13
1041	a8 = 6.067681784d-16
1042	a9 = -2.356829271d-19
1043	ionsec_co2plus = a0 + a1altaux + a2altaux*2 + a3altaux**3 +
1044	$ a4altaux4 + a5altaux*5 + a6altaux**6 +
1045	$ a7altaux7 + a8altaux*8 + a9altaux**9
1046	ionsec_co2plus = 10**(ionsec_co2plus-2.)
1047	elseif(altaux.ge.80..and.altaux.le.112.9) then
1048	b0 = 1.159404818d6
1049	b1 = -5.617523193d4
1050	b2 = 8.535640078d2
1051	b3 = -5.956128524d-1
1052	b4 = -8.532255532d-2
1053	b5 = 1.217829692d-4
1054	b6 = 9.641535217d-6
1055	b7 = -4.502450788d-8
1056	b8 = -4.9684920146d-10
1057	b9 = 4.830572346d-12
1058	b10 = -1.168635127d-14
1059	ionsec_co2plus = b0 + b1altaux + b2altaux*2 + b3altaux**3 +
1060	$ b4altaux4 + b5altaux*5 + b6altaux**6 +
1061	$ b7altaux7 + b8altaux*8 + b9altaux**9 +
1062	$ b10altaux*10
1063	else
1064	ionsec_co2plus = 0.d0
1065	endif
1066	if(ionsec_co2plus.gt.100.d0.or.ionsec_co2plus.lt.0.d0)
1067	$ ionsec_co2plus=0.d0
1068
1069	return
1070
1071	end
1072
1073
1074	c***********************************************************************
1075	function ionsec_o2plus (zenit, alt)
1076
1077	c O2+ production by photoelectrons, from Nicholson et al. 2009
1078
1079	c FGG sep 2010 first version
1080	c***********************************************************************
1081
1082	implicit none
1083
1084	! Arguments
1085	real*8 ionsec_o2plus
1086	real zenit
1087	real alt
1088
1089	! Local variables
1090
1091	real*8 a0,a1,a2,a3,a4,a5,a6,a7
1092	real*8 b0,b1,b2,b3,b4,b5,b6,b7,b8
1093	real*8 altaux
1094	real*8 zenit_rad
1095
1096	!!!!!!! Program starts
1097
1098	zenit_rad=dble(zenit*3.141592/180.)
1099
1100	if(zenit.le.90.) then
1101	altaux=dble(alt)+
1102	$ 15.cos(zenit_rad)-40.sqrt(cos(zenit_rad))+25.
1103	else
1104	altaux=dble(alt)
1105	endif
1106
1107	if(altaux.gt.112.9) then
1108	a0 = 7.265142765d2
1109	a1 = -2.714716832d1
1110	a2 = 4.315022800d-1
1111	a3 = -3.774025573d-3
1112	a4 = 1.962771814d-5
1113	a5 = -6.076128732d-8
1114	a6 = 1.037835637d-10
1115	a7 = -7.552930040d-14
1116	ionsec_o2plus = a0 + a1altaux + a2altaux*2 + a3altaux**3 +
1117	$ a4altaux4 + a5altaux*5 + a6altaux*6 + a7altaux**7
1118	ionsec_o2plus = 10**(ionsec_o2plus-2.)
1119	elseif(altaux.gt.80..and.altaux.le.112.9) then
1120	b0 = 3.622091694d6
1121	b1 = -1.45914419d5
1122	b2 = 1.764604914d3
1123	b3 = 1.186141031d0
1124	b4 = -1.331821089d-1
1125	b5 = -3.661686584d-4
1126	b6 = 1.934372959d-5
1127	b7 = -1.307294421d-7
1128	b8 = 2.846288872d-10
1129	ionsec_o2plus = b0 + b1altaux + b2altaux*2 + b3altaux**3 +
1130	$ b4altaux4 + b5altaux*5 + b6altaux*6 + b7altaux**7
1131	$ + b8altaux*8
1132	else
1133	ionsec_o2plus = 0.d0
1134	endif
1135	if(ionsec_o2plus.gt.100.d0.or.ionsec_o2plus.lt.0.d0)
1136	$ ionsec_o2plus=0.d0
1137
1138
1139	return
1140
1141	end
1142
1143
1144
1145
1146	c**********************************************************************
1147	c**********************************************************************
1148
1149	subroutine phdisrate(ig,nlayer,chemthermod,zenit,i)
1150
1151	c Calculates photoionization and photodissociation rates from the
1152	c photoabsorption rates calculated in jthermcalc_e107 and the
1153	c ionization/dissociation branching ratios in param_read_e107
1154
1155	c apr 2002 fgg first version
1156	c**********************************************************************
1157
1158	use param_v4_h, only: ninter,nabs,
1159	. jfotsout,fluxtop,
1160	. jion,jdistot,jdistot_b,
1161	. efdisco2,efdiso2,efdish2o,
1162	. efdish2o2,efdish2,efdiso3,
1163	. efdiso,efdisn,efdish,
1164	. efdisno,efdisn2,efdisno2,
1165	. efdisco,efionco2,efiono2,efionn2,
1166	. efionco,efiono3p,efionn,
1167	. efionno,efionh
1168
1169	implicit none
1170
1171	c arguments
1172
1173	integer i !altitude
1174	integer ig,chemthermod,nlayer
1175	real zenit
1176
1177	c local variables
1178
1179	integer inter,iz,j
1180	real lambda
1181	real jdis(nabs,ninter,nlayer)
1182	character*1 dn
1183
1184
1185	c**********************************************************************
1186
1187	c photodissociation and photoionization rates initialization
1188
1189	jion(:,i,:) = 0.d0
1190	jdistot(:,i) = 0.d0
1191	jdistot_b(:,i) = 0.d0
1192
1193	! jion(1,i,1) = 0.d0 ! CO2 channel 1 ( --> CO2+ + e- )
1194	! jion(1,i,2) = 0.d0 ! CO2 channel 2 ( --> O+ + CO + e- )
1195	! jion(1,i,3) = 0.d0 ! CO2 channel 3 ( --> CO+ + O + e- )
1196	! jion(1,i,4) = 0.d0 ! CO2 channel 4 ( --> C+ + O2 + e- )
1197	! jion(2,i,1) = 0.d0 ! O2 (only one ionization channel)
1198	! jion(3,i,1) = 0.d0 ! O3P (only one ionization channel)
1199	! jion(4,i,1) = 0.d0 ! H2O (no ionization)
1200	! jion(5,i,1) = 0.d0 ! H2 (no ionization)
1201	! jion(6,i,1) = 0.d0 ! H2O2 (no ionization)
1202	! jion(7,i,1) = 0.d0 ! O3 (no ionization)
1203	! jion(8,i,1) = 0.d0 ! N2 channel 1 ( --> n2+ + e- )
1204	! jion(8,i,2) = 0.d0 ! N2 channel 2 ( --> n+ + n + e- )
1205	! jion(9,i,1) = 0.d0 ! N ( --> N+ + e- )
1206	! jion(10,i,1)= 0.d0 ! We do not mind its ionization
1207	! jion(11,i,1) = 0.d0 ! CO channel 1 ( --> co+ + e- )
1208	! jion(11,i,2) = 0.d0 ! CO channel 2 ( --> c+ + o + e- )
1209	! jion(11,i,3) = 0.d0 ! CO channel 3 ( --> o+ + c + e- )
1210	! jion(12,i,1) = 0.d0 ! H ( --> H+ + e- )
1211	! jion(13,i,1) = 0.d0
1212	! do j=1,nabs
1213	! jdistot(j,i) = 0.
1214	! jdistot_b(j,i) = 0.
1215	! end do
1216
1217	if(zenit.gt.140.) then
1218	dn='n'
1219	else
1220	dn='d'
1221	end if
1222
1223	if(dn.eq.'n') return
1224
1225	c photodissociation and photoionization rates for each species
1226
1227	do inter=1,ninter
1228	c CO2
1229	jdis(1,inter,i) = jfotsout(inter,1,i) * fluxtop(inter)
1230	$ * efdisco2(inter)
1231	if(inter.gt.29.and.inter.le.32) then
1232	jdistot(1,i) = jdistot(1,i) + jdis(1,inter,i)
1233	else if(inter.le.29) then
1234	jdistot_b(1,i) = jdistot_b(1,i) + jdis(1,inter,i)
1235	end if
1236	jion(1,i,1)=jion(1,i,1) +
1237	$ jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,1)
1238	jion(1,i,2)=jion(1,i,2) +
1239	$ jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,2)
1240	jion(1,i,3)=jion(1,i,3) +
1241	$ jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,3)
1242	jion(1,i,4)=jion(1,i,4) +
1243	$ jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,4)
1244
1245
1246	c O2
1247	jdis(2,inter,i) = jfotsout(inter,2,i) * fluxtop(inter)
1248	$ * efdiso2(inter)
1249	if(inter.ge.31) then
1250	jdistot(2,i) = jdistot(2,i) + jdis(2,inter,i)
1251	else if(inter.eq.30) then
1252	jdistot(2,i)=jdistot(2,i)+0.02*jdis(2,inter,i)
1253	jdistot_b(2,i)=jdistot_b(2,i)+0.98*jdis(2,inter,i)
1254	else if(inter.lt.31) then
1255	jdistot_b(2,i) = jdistot_b(2,i) + jdis(2,inter,i)
1256	end if
1257	jion(2,i,1)=jion(2,i,1) +
1258	$ jfotsout(inter,2,i) * fluxtop(inter) * efiono2(inter,1)
1259	jion(2,1,2)=jion(2,1,2) +
1260	$ jfotsout(inter,2,i) * fluxtop(inter) * efiono2(inter,2)
1261	!(1.-efdiso2(inter))
1262
1263	c O3P
1264	jion(3,i,1)=jion(3,i,1) +
1265	$ jfotsout(inter,3,i) * fluxtop(inter) * efiono3p(inter)
1266
1267	c H2O
1268	jdis(4,inter,i) = jfotsout(inter,4,i) * fluxtop(inter)
1269	$ * efdish2o(inter)
1270	jdistot(4,i) = jdistot(4,i) + jdis(4,inter,i)
1271
1272	c H2
1273	jdis(5,inter,i) = jfotsout(inter,5,i) * fluxtop(inter)
1274	$ * efdish2(inter)
1275	jdistot(5,i) = jdistot(5,i) + jdis(5,inter,i)
1276
1277	c H2O2
1278	jdis(6,inter,i) = jfotsout(inter,6,i) * fluxtop(inter)
1279	$ * efdish2o2(inter)
1280	jdistot(6,i) = jdistot(6,i) + jdis(6,inter,i)
1281
1282	!Only if O3, N or ion chemistry requested
1283	if(chemthermod.ge.1) then
1284	c O3
1285	jdis(7,inter,i) = jfotsout(inter,7,i) * fluxtop(inter)
1286	$ * efdiso3(inter)
1287	if(inter.eq.34) then
1288	jdistot(7,i) = jdistot(7,i) + jdis(7,inter,i)
1289	else if (inter.eq.35) then
1290	jdistot(7,i) = jdistot(7,i) + 0.997 * jdis(7,inter,i)
1291	jdistot_b(7,i) = jdistot_b(7,i) + 0.003 * jdis(7,inter,i)
1292	else if (inter.eq.36) then
1293	jdistot_b(7,i) = jdistot_b(7,i) + jdis(7,inter,i)
1294	endif
1295	endif !Of chemthermod.ge.1
1296
1297	!Only if N or ion chemistry requested
1298	if(chemthermod.ge.2) then
1299	c N2
1300	jdis(8,inter,i) = jfotsout(inter,8,i) * fluxtop(inter)
1301	$ * efdisn2(inter)
1302	jdistot(8,i) = jdistot(8,i) + jdis(8,inter,i)
1303	jion(8,i,1) = jion(8,i,1) + jfotsout(inter,8,i) *
1304	$ fluxtop(inter) * efionn2(inter,1)
1305	jion(8,i,2) = jion(8,i,2) + jfotsout(inter,8,i) *
1306	$ fluxtop(inter) * efionn2(inter,2)
1307
1308	c N
1309	jion(9,i,1) = jion(9,i,1) + jfotsout(inter,9,i) *
1310	$ fluxtop(inter) * efionn(inter)
1311
1312	c NO
1313	jdis(10,inter,i) = jfotsout(inter,10,i) * fluxtop(inter)
1314	$ * efdisno(inter)
1315	jdistot(10,i) = jdistot(10,i) + jdis(10,inter,i)
1316	jion(10,i,1) = jion(10,i,1) + jfotsout(inter,10,i) *
1317	$ fluxtop(inter) * efionno(inter)
1318
1319	c NO2
1320	jdis(13,inter,i) = jfotsout(inter,13,i) * fluxtop(inter)
1321	$ * efdisno2(inter)
1322	jdistot(13,i) = jdistot(13,i) + jdis(13,inter,i)
1323
1324	endif !Of chemthermod.ge.2
1325
1326	!Only if ion chemistry requested
1327	if(chemthermod.eq.3) then
1328	c CO
1329	jdis(11,inter,i) = jfotsout(inter,11,i) * fluxtop(inter)
1330	$ * efdisco(inter)
1331	jdistot(11,i) = jdistot(11,i) + jdis(11,inter,i)
1332	jion(11,i,1) = jion(11,i,1) + jfotsout(inter,11,i) *
1333	$ fluxtop(inter) * efionco(inter,1)
1334	jion(11,i,2) = jion(11,i,2) + jfotsout(inter,11,i) *
1335	$ fluxtop(inter) * efionco(inter,2)
1336	jion(11,i,3) = jion(11,i,3) + jfotsout(inter,11,i) *
1337	$ fluxtop(inter) * efionco(inter,3)
1338
1339	c H
1340	jion(12,i,1) = jion(12,i,1) + jfotsout(inter,12,i) *
1341	$ fluxtop(inter) * efionh(inter)
1342	endif !Of chemthermod.eq.3
1343
1344
1345	end do
1346
1347
1348	return
1349
1350
1351	end
1352
1353
1354
1355
1356	c**********************************************************************
1357	c***************************************************************************
1358	subroutine getch (ig,chemthermod,tt,zkm)
1359
1360
1361	c Reaction rates. The parameters rcoef are read from the
1362	c chemthermos_reactionrates.def file
1363
1364
1365	c***************************************************************************
1366
1367	use param_v4_h, only: rcoef,
1368	. ch2, ch3, ch4, ch5, ch7,ch9,ch10,ch11,ch13,ch14,ch15,ch18,
1369	. ch19,ch20,ch21,ch22,ch23,ch24,ch30,ch31,ch32,ch33,ch34,
1370	. ch35,ch36,ch37,ch38,ch39,ch40,ch41,ch42,ch43,ch45,
1371	. ch46,ch47,ch48,ch49,ch50,ch55,ch56,ch57,ch58,ch59,ch62,
1372	. ch63,ch64,ch65,ch66,ch67,ch68,ch69,ch70,ch71,
1373	. ch72,ch73,ch74,ch75,ch76,ch85,ch86,ch87
1374
1375
1376
1377	implicit none
1378
1379	c Arguments
1380	integer ig,chemthermod
1381	real*8 tt ! Temperature
1382	real zkm ! Altitude in km
1383
1384	c local variables:
1385	real*8 tcte
1386	real*8 t_elect ! electronic temperatures
1387	real*8 val ! valores de alturas corresp a t_elect
1388	real*8 zhanson(9),tehanson(9)
1389	real*8 incremento
1390	integer ii, i1, i2
1391
1392	c**************************************************************************
1393
1394	tcte = tt
1395
1396	! goto 151
1397	!Electronic temperatures
1398	! (Hanson et al. 1977) approx. from Mars p. 107
1399	zhanson(1) = 120.
1400	zhanson(2) = 130.
1401	zhanson(3) = 150.
1402	zhanson(4) = 175.
1403	zhanson(5) = 200.
1404	zhanson(6) = 225.
1405	zhanson(7) = 250.
1406	zhanson(8) = 275.
1407	zhanson(9) = 300.
1408	tehanson(1) = tt
1409	tehanson(2) = 200.
1410	tehanson(3) = 300.
1411	tehanson(4) = 500.
1412	tehanson(5) = 1250.
1413	tehanson(6) = 2000.
1414	tehanson(7) = 2200.
1415	tehanson(8) = 2400.
1416	tehanson(9) = 2500.
1417	if ( zkm .le. 120. ) then
1418	t_elect = tt
1419	else if(zkm .ge.300.) then
1420	t_elect=tehanson(9)
1421	else
1422	do ii=9,2,-1
1423	if ( zkm .lt. zhanson(ii) ) then
1424	i1 = ii - 1
1425	i2 = ii
1426	endif
1427	enddo
1428	incremento = ( tehanson(i2)-tehanson(i1) ) /
1429	$ ( zhanson(i2)-zhanson(i1) )
1430	t_elect = tehanson(i1) + (zkm-zhanson(i1)) * incremento
1431
1432	! t_elect = t_elect * 2.
1433	endif
1434	! 151 continue
1435	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
1436	! t_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
1437
1438	!Initializations
1439	ch2 = 0.d0
1440	ch3 = 0.0
1441	ch4 = 0.d0
1442	ch5 = 0.d0
1443	ch7 = 0.d0
1444	ch9 = 0.d0
1445	ch10 = 0.d0
1446	ch11 = 0.d0
1447	ch13 = 0.d0
1448	ch14 = 0.d0
1449	ch15 = 0.d0
1450	ch18 = 0.d0
1451	ch19 = 0.d0
1452	ch20 = 0.d0
1453	ch21 = 0.d0
1454	ch22 = 0.d0
1455	ch23 = 0.d0
1456	ch24 = 0.d0
1457	ch30 = 0.d0
1458	ch31 = 0.d0
1459	ch32 = 0.d0
1460	ch33 = 0.d0
1461	ch34 = 0.d0
1462	ch35 = 0.d0
1463	ch36 = 0.d0
1464	ch37 = 0.d0
1465	ch38 = 0.d0
1466	ch39 = 0.d0
1467	ch40 = 0.d0
1468	ch41 = 0.d0
1469	ch42 = 0.d0
1470	ch43 = 0.d0
1471	ch45 = 0.d0
1472	ch46 = 0.d0
1473	ch47 = 0.d0
1474	ch48 = 0.d0
1475	ch49 = 0.d0
1476	ch50 = 0.d0
1477	ch55 = 0.d0
1478	ch56 = 0.d0
1479	ch57 = 0.d0
1480	ch58 = 0.d0
1481	ch59 = 0.d0
1482	ch62 = 0.d0
1483	ch63 = 0.d0
1484	ch64 = 0.d0
1485	ch65 = 0.d0
1486	ch66 = 0.d0
1487	ch67 = 0.d0
1488	ch68 = 0.d0
1489	ch69 = 0.d0
1490	ch70 = 0.d0
1491	ch71 = 0.d0
1492	ch72 = 0.d0
1493	ch73 = 0.d0
1494	ch74 = 0.d0
1495	ch75 = 0.d0
1496	ch76 = 0.d0
1497	ch85 = 0.d0
1498	ch86 = 0.d0
1499
1500
1501	!Reaction rates
1502	!ch2: h + o2 + co2 --> ho2 + co2
1503	! JPL 2003 (low pressure limit)*2.5
1504	! ch2 = 1.425d-31 * (tcte / 300.)**(-1.6d0)
1505	! JPL 2011 (low pressure limit)*2.5
1506	! ch2 = 1.1e-31 * (tcte / 300.)**(-1.3)
1507	ch2=rcoef(1,1)((tcte/300.)rcoef(1,2))exp(rcoef(1,3)/tcte)
1508
1509	!ch3: o + ho2 --> oh + o2
1510	! JPL 2011:
1511	! ch3 = 3.0d-11 * exp(200.d0 / tcte)
1512	ch3=rcoef(2,1)((tcte/300.)rcoef(2,2))exp(rcoef(2,3)/tcte)
1513
1514	ch4: co + oh --> co2 + h
1515	!Nair et al, 1994:
1516	!ch4 = 3.2d-13 * exp(-300.d0 / tcte)
1517	!mccabe et al., grl, 28, 3135, 2001
1518	!ch4 = 1.57d-13 + 3.54d-33*concco2
1519	!JPL 2011 (low pressure limit):
1520	! ch4 = 1.5d-13 * (tcte/300.)**0.6
1521	ch4=rcoef(3,1)((tcte/300.)rcoef(3,2))exp(rcoef(3,3)/tcte)
1522
1523	ch5: ho2 + ho2 --> h2o2 + o2
1524	!JPL 2003:
1525	!ch5 = 2.3d-13 * exp(600.d0 / tcte)
1526	!JPL 2011:
1527	! ch5 = 3.0d-13 * exp(460.d0 / tcte)
1528	ch5=rcoef(4,1)((tcte/300.)rcoef(4,2))exp(rcoef(4,3)/tcte)
1529
1530	ch7: oh + ho2 --> h2o + o2
1531	!JPL 2011:
1532	! ch7 = 4.8d-11 * exp(250.d0 / tcte)
1533	ch7=rcoef(5,1)((tcte/300.)rcoef(5,2))exp(rcoef(5,3)/tcte)
1534
1535	ch9: o(1d) + h2o --> 2oh
1536	!JPL 2003:
1537	!ch9 = 2.2d-10
1538	!JPL 2011:
1539	! ch9 = 1.63d-10 * exp(60.d0 / tcte)
1540	ch9=rcoef(6,1)((tcte/300.)rcoef(6,2))exp(rcoef(6,3)/tcte)
1541
1542	ch10: o + o + co2 --> o2 + co2
1543	!JPL 1990:
1544	! ch10 = 1.1d-27 * (tcte **(-2.0d0)) !Estandard en el 1-D
1545	!Tsang and Hampson, 1986:
1546	! ch10 = 1.3d-34 * exp(900.d0/tcte)
1547	ch10=rcoef(7,1)((tcte/300.)rcoef(7,2))exp(rcoef(7,3)/tcte)
1548
1549	ch11: o + oh --> o2 + h
1550	!JPL 2003:
1551	!ch11 = 2.2d-11 * exp(120.d0 / tcte)
1552	!JPL 2011:
1553	! ch11 = 1.8d-11 * exp(180.d0 /tcte)
1554	ch11=rcoef(8,1)((tcte/300.)rcoef(8,2))exp(rcoef(8,3)/tcte)
1555
1556	ch13: h + ho2 --> h2 + o2
1557	!JPL 2003:
1558	!ch13 = 6.5d-12
1559	!JPL 2011:
1560	! ch13 = 6.9d-12
1561	ch13=rcoef(9,1)((tcte/300.)rcoef(9,2))exp(rcoef(9,3)/tcte)
1562
1563	ch14: o(1d) + h2 --> h + oh
1564	!JPL 2003:
1565	!ch14 = 1.1d-10
1566	!JPL 2011:
1567	! ch14 = 1.2d-10
1568	ch14=rcoef(10,1)((tcte/300.)rcoef(10,2))
1569	$ exp(rcoef(10,3)/tcte)
1570
1571	ch15: oh + h2 --> h + h2o
1572	!JPL 2003:
1573	!ch15 = 5.5d-12 * exp (-2000.d0 / tcte)
1574	!JPL 2011:
1575	! ch15 = 2.8d-12 * exp (-1800.d0 / tcte)
1576	ch15=rcoef(11,1)((tcte/300.)rcoef(11,2))
1577	$ exp(rcoef(11,3)/tcte)
1578
1579	ch18: oh + h2o2 --> h2o + ho2
1580	!JPL 2003:
1581	!ch18 = 2.9d-12 * exp (-160.d0 / tcte)
1582	!JPL 2011:
1583	! ch18 = 1.8d-12
1584	ch18=rcoef(12,1)((tcte/300.)rcoef(12,2))
1585	$ exp(rcoef(12,3)/tcte)
1586
1587	ch19: o(1d) + co2 --> o + co2
1588	!JPL 2003:
1589	!ch19 = 7.4d-11 * exp(120.d0 / tcte)
1590	!JPL 2011:
1591	! ch19 = 7.5d-11 * exp(115.d0 / tcte)
1592	ch19=rcoef(13,1)((tcte/300.)rcoef(13,2))
1593	$ exp(rcoef(13,3)/tcte)
1594
1595	ch20: o(1d) + o2 --> o + o2
1596	!JPL 2003:
1597	!ch20 = 3.2d-11 * exp (70.d0 / tcte)
1598	!JPL 2011:
1599	! ch20 = 3.3d-11 * exp(55.d0 / tcte)
1600	ch20=rcoef(14,1)((tcte/300.)rcoef(14,2))
1601	$ exp(rcoef(14,3)/tcte)
1602
1603	ch21: o + o2 + co2 --> o3 + co2
1604	!JPL 2011 * 2.5:
1605	! ch21 = 1.5d-33 * ((tcte / 300.d0) ** (-2.4d0))
1606	ch21=rcoef(15,1)((tcte/300.)rcoef(15,2))
1607	$ exp(rcoef(15,3)/tcte)
1608
1609
1610	!Only if O3, N or ion chemistry requested
1611	if(chemthermod.ge.1) then
1612
1613	ch22: o3 + h --> o2 + oh
1614	!JPL 2011:
1615	! ch22 = 1.4d-10 * exp (-470.d0 / tcte)
1616	ch22=rcoef(16,1)((tcte/300.)rcoef(16,2))
1617	$ exp(rcoef(16,3)/tcte)
1618
1619	ch23: o3 + oh --> ho2 + o2
1620	!JPL 2011:
1621	! ch23 = 1.7d-12 * exp (-940.d0 / tcte)
1622	ch23=rcoef(17,1)((tcte/300.)rcoef(17,2))
1623	$ exp(rcoef(17,3)/tcte)
1624
1625	ch24: o3 + ho2 --> oh + 2o2
1626	!JPL 2011:
1627	! ch24 = 1.0d-14 * exp (-490.d0 / tcte)
1628	ch24=rcoef(18,1)((tcte/300.)rcoef(18,2))
1629	$ exp(rcoef(18,3)/tcte)
1630
1631	endif
1632
1633	!Only if N or ion chemistry requested
1634	if(chemthermod.ge.2) then
1635	c<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1636	c N chemistry
1637	ch30: n + no --> n2 + o
1638	!JPL 2011:
1639	! ch30 = 2.1d-11 * exp (100.d0 / tcte)
1640	ch30=rcoef(19,1)((tcte/300.)rcoef(19,2))
1641	$ exp(rcoef(19,3)/tcte)
1642
1643	ch31: n2 + o(1d) --> n2 + o
1644	!JPL 2011:
1645	! ch31 = 2.15d-11 * exp (110.d0 / tcte)
1646	ch31=rcoef(20,1)((tcte/300.)rcoef(20,2))
1647	$ exp(rcoef(20,3)/tcte)
1648
1649	ch32: n + o2 --> no + o
1650	!JPL 2011:
1651	! ch32 = 1.5d-11 * exp (-3600.d0 / tcte)
1652	ch32=rcoef(21,1)((tcte/300.)rcoef(21,2))
1653	$ exp(rcoef(21,3)/tcte)
1654
1655	ch33: n + oh --> no + h
1656	!Atkinson et al., 1989 (usado en Nair et al., 1994)
1657	! ch33 = 3.8d-11 * exp (85.d0 / tcte)
1658	ch33=rcoef(22,1)((tcte/300.)rcoef(22,2))
1659	$ exp(rcoef(22,3)/tcte)
1660
1661	ch34: n + o3 --> no + o2
1662	!JPL 2011 (it is an upper limit):
1663	! ch34 = 1.0d-16
1664	ch34=rcoef(23,1)((tcte/300.)rcoef(23,2))
1665	$ exp(rcoef(23,3)/tcte)
1666
1667	ch35: n + ho2 --> no + oh
1668	!Brune et al., 1983 (from Nair et al., 1994)
1669	! ch35 = 2.2d-11
1670	ch35=rcoef(24,1)((tcte/300.)rcoef(24,2))
1671	$ exp(rcoef(24,3)/tcte)
1672
1673	ch36: n(2d) + o --> n + o
1674	!Fell et al., 1990 (from Nair et al., 1994)
1675	!ch36 = 6.9d-13
1676	!Herron, 1999:
1677	! ch36 = 3.3d-12 * exp(-260.d0 / tcte)
1678	ch36=rcoef(25,1)((tcte/300.)rcoef(25,2))
1679	$ exp(rcoef(25,3)/tcte)
1680
1681	ch37: n(2d) + n2 --> n + n2
1682	!Herron, 1999:
1683	!Coincides with Nair et al., 1994:
1684	! ch37 = 1.7d-14
1685	ch37=rcoef(26,1)((tcte/300.)rcoef(26,2))
1686	$ exp(rcoef(26,3)/tcte)
1687
1688	ch38: n(2d) + co2 --> no + co
1689	!Pipper et al., 1989 (from Nair et al., 1994):
1690	!ch38 = 3.5d-13
1691	!Herron, 1999:
1692	! ch38 = 3.6d-13
1693	ch38=rcoef(27,1)((tcte/300.)rcoef(27,2))
1694	$ exp(rcoef(27,3)/tcte)
1695
1696	ch39: no + ho2 --> no2+oh
1697	!JPL 2006:
1698	!ch39 = 3.5d-12 * exp (250.d0 / tcte)
1699	!JPL 2011:
1700	! ch39 = 3.3d-12 * exp(270.d0 / tcte)
1701	ch39=rcoef(28,1)((tcte/300.)rcoef(28,2))
1702	$ exp(rcoef(28,3)/tcte)
1703
1704	ch40: o + no + co2 --> no2 + co2
1705	!JPL 2011 * 2.5 (low pressure limit)
1706	! ch40 = 2.5d0 * 9.0d-32 * ((tcte / 300.d0) ** (-1.5d0))
1707	ch40=rcoef(29,1)((tcte/300.)rcoef(29,2))
1708	$ exp(rcoef(29,3)/tcte)
1709
1710	ch41: o + no2 --> no + o2
1711	!JPL 2011:
1712	! ch41 = 5.1d-12 * exp (210.d0 / tcte)
1713	ch41=rcoef(30,1)((tcte/300.)rcoef(30,2))
1714	$ exp(rcoef(30,3)/tcte)
1715
1716	ch42: no + o3 --> no2 + o2
1717	!JPL 2011:
1718	! ch42 = 3.0d-12 * exp (-1500.d0 / tcte)
1719	ch42=rcoef(31,1)((tcte/300.)rcoef(31,2))
1720	$ exp(rcoef(31,3)/tcte)
1721
1722	ch43: h + no2 --> no + oh
1723	!JPL 2011:
1724	! ch43 = 4.0d-10 * exp (-340.d0 / tcte)
1725	ch43=rcoef(32,1)((tcte/300.)rcoef(32,2))
1726	$ exp(rcoef(32,3)/tcte)
1727
1728	ch45: n + o --> no
1729	!Du and Dalgarno, 1990
1730	! ch45 = 2.8d-17 * ((300.d0 / tcte) ** 0.5)
1731	ch45=rcoef(33,1)((tcte/300.)rcoef(33,2))
1732	$ exp(rcoef(33,3)/tcte)
1733
1734	endif !of if(chemthermod.ge.2)
1735
1736
1737	!Only if ion chemistry requested
1738	if(chemthermod.eq.3) then
1739	c
1740	c Ionosphere
1741	c
1742	ch46: co2+ + O2 --> O2+ + CO2
1743	!Moffat et al., 2005 (given by GG):
1744	!ch46 = 5.0d-11
1745	!Copp et al., 1982:
1746	!ch46 = 5.3d-11
1747	!Aninich 1993 (from Fox and Sung 2001):
1748	! ch46 = 5.5d-11 * (300.d0/t_elect)**0.82
1749	ch46=rcoef(34,1)((t_elect/300.)rcoef(34,2))
1750	$ exp(rcoef(34,3)/t_elect)
1751
1752
1753	ch47: CO2+ + O --> O+ + CO2
1754	!Original (incorrect) value (corresponds to ch48):
1755	!ch47 = 1.64d-10
1756	!Fehsenfeld et al., 1970 (from UMIST,
1757	!Fox and Sung 2001, Krasnopolsky 2002):
1758	! ch47 = 9.6d-11
1759	ch47=rcoef(35,1)((t_elect/300.)rcoef(35,2))
1760	$ exp(rcoef(35,3)/t_elect)
1761
1762	ch48: CO2+ + O --> O2+ + CO
1763	!Original (incorrect) value (corresponds to ch47):
1764	!ch48 = 9.6d-11
1765	!Fehsenfeld et al., 1970 (from UMIST,
1766	!Fox and Sung 2001, Krasnopolsky 2002):
1767	! ch48 = 1.64d-10
1768	ch48=rcoef(36,1)((t_elect/300.)rcoef(36,2))
1769	$ exp(rcoef(36,3)/t_elect)
1770
1771	ch49: O2+ + elect --> O + O
1772	!Alge et al., 1983:
1773	!Here we do not divide into reaction producing different
1774	!O atomic states. O + O(1d) seems to be the dominant products
1775	!(see Fox and Sung 2002). We should consider dividing
1776	!into two different reactions
1777	! ch49 = 2.0d-7(300.d0/t_elect)*(0.7d0)
1778	ch49=rcoef(37,1)((t_elect/300.)rcoef(37,2))
1779	$ exp(rcoef(37,3)/t_elect)
1780
1781	ch50: O+ + CO2 --> O2+ + CO
1782	!Adams et al., 1980 (from UMIST):
1783	! ch50 = 9.4d-10
1784	!Anicich 1993 (from Fox and Sung 2001):
1785	!ch50 = 1.1d-9
1786	ch50=rcoef(38,1)((t_elect/300.)rcoef(38,2))
1787	$ exp(rcoef(38,3)/t_elect)
1788
1789	ch55: CO2+ + e ----> CO + O
1790	!Mitchell, 1990 (from UMIST):
1791	! ch55 = 3.8d-7(300.d0/t_elect)*(0.5d0)
1792	!Gougousi et al., 1997 (from Fox and Sung 2001):
1793	!ch55 = 3.5d-7(300.d0/t_elect)*0.5d0
1794	ch55=rcoef(39,1)((t_elect/300.)rcoef(39,2))
1795	$ exp(rcoef(39,3)/t_elect)
1796
1797	ch56: O+ + CO2 ---> O2 + CO+
1798	!Original, Kim et al., 1989:
1799	!ch56 = 9.4d-10
1800	!It does not appear in any other paper. Its presence in
1801	!Kim et al., 1989 is probably a confusion with ch50.
1802	! ch56 = 0.d0
1803	ch56=rcoef(40,1)((t_elect/300.)rcoef(40,2))
1804	$ exp(rcoef(40,3)/t_elect)
1805
1806	ch57: CO+ + CO2 ---> CO2+ + CO
1807	!Adams et al., 1978 (from UMIST):
1808	! ch57 = 1.0d-9
1809	!Anicich 1993 (from Fox and Sung 2001):
1810	!ch57 = 1.1d-9
1811	ch57=rcoef(41,1)((t_elect/300.)rcoef(41,2))
1812	$ exp(rcoef(41,3)/t_elect)
1813
1814
1815	ch58: CO+ + O ---> O+ + CO
1816	!Fenhsenfeld et al. 1970 (from UMIST, F&S2001, K2002):
1817	! ch58 = 1.4d-10
1818	ch58=rcoef(42,1)((t_elect/300.)rcoef(42,2))
1819	$ exp(rcoef(42,3)/t_elect)
1820
1821	ch59: C+ + CO2 ---> CO+ + CO !!!! NEW !!!
1822	!Fahey et al., 1981 (from UMIST, F&S2001, K2002):
1823	! ch59 = 1.1d-9
1824	ch59=rcoef(43,1)((t_elect/300.)rcoef(43,2))
1825	$ exp(rcoef(43,3)/t_elect)
1826
1827	ch62: CO2+ + NO --> NO+ + CO2
1828	!Copp et al., 1982 (from UMIST):
1829	! ch62 = 1.2d-10
1830	!Anicich 1993 (from Fox and Sung 2001):
1831	!ch62 = 1.23d-10
1832	ch62=rcoef(44,1)((t_elect/300.)rcoef(44,2))
1833	$ exp(rcoef(44,3)/t_elect)
1834
1835	ch63: CO2+ + N --> NO + CO+
1836	!Kim et al., 1989:
1837	!ch63 = 1.0d-11
1838	!Scott et al., 1998 (from Fox and Sung 2001):
1839	! ch63 = 3.4d-10
1840	ch63=rcoef(45,1)((t_elect/300.)rcoef(45,2))
1841	$ exp(rcoef(45,3)/t_elect)
1842
1843	ch64: O2+ + NO --> NO+ + O2
1844	!Middey and Vigiano 1999 (from Fox and Sung 2001):
1845	! ch64 = 4.5d-10
1846	!Aninich 1993 (from UMIST):
1847	!ch64 = 4.6d-10
1848	ch64=rcoef(46,1)((t_elect/300.)rcoef(46,2))
1849	$ exp(rcoef(46,3)/t_elect)
1850
1851	ch65: O2+ + N2 --> NO+ + NO
1852	!Original from GG, Moffat 2005:
1853	!ch65 = 1.0d-16
1854	!Ferguson 1973 (from Fox and Sung 2001):
1855	! ch65 = 1.0d-15
1856	ch65=rcoef(47,1)((t_elect/300.)rcoef(47,2))
1857	$ exp(rcoef(47,3)/t_elect)
1858
1859	ch66: O2+ + N --> NO+ + O
1860	!Kim et al., 1989:
1861	!ch66 = 1.2d-10
1862	!Scott et al., 1998 (from Fox and Sung 2001):
1863	! ch66 = 1.0d-10
1864	!Goldan et al., 1966 (from UMIST):
1865	!ch66 = 1.8d-10
1866	ch66=rcoef(48,1)((t_elect/300.)rcoef(48,2))
1867	$ exp(rcoef(48,3)/t_elect)
1868
1869
1870	ch67: O+ + N2 --> NO+ + N
1871	!Moffat 2005:
1872	!ch67 = 1.2d-12 * (300.d0/t_elect)**(0.41d0)
1873	!Hierl et al. 1997 (from Fox and Sung 2001):
1874	! ch67 = 1.2d-12 * (300.d0/t_elect)**0.45d0
1875	!Adams et al., 1980 (from UMIST):
1876	!ch67=2.42d-12 * (300.d0/t_elec)*(-0.21)exp(44./t_elec)
1877	ch67=rcoef(49,1)((t_elect/300.)rcoef(49,2))
1878	$ exp(rcoef(49,3)/t_elect)
1879
1880	ch68: N2+ + CO2 --> CO2+ + N2
1881	!Adams et al. 1980 (from UMIST):
1882	!ch68 = 7.7d-10
1883	!Dotan et al. 2000 (from F&S2001):
1884	! ch68 = 9.0d-10 * (300./t_elect)**0.23
1885	ch68=rcoef(50,1)((t_elect/300.)rcoef(50,2))
1886	$ exp(rcoef(50,3)/t_elect)
1887
1888	ch69: N2+ + O3p --> NO+ + N
1889	!McFarland et al., 1974 (from UMIST):
1890	!ch69 = 1.3d-10
1891	!Scott et al. 1999 (from F&S2001):
1892	! ch69 = 1.33d-10 * (300./t_elect)**0.44
1893	ch69=rcoef(51,1)((t_elect/300.)rcoef(51,2))
1894	$ exp(rcoef(51,3)/t_elect)
1895
1896	ch70: N2+ + CO --> N2 + CO+
1897	!Adams et al., 1980 (from UMIST):
1898	! ch70 = 7.4d-11
1899	!Frost et al., 1998 (from F&S2001):
1900	!ch70 = 7.6d-11
1901	ch70=rcoef(52,1)((t_elect/300.)rcoef(52,2))
1902	$ exp(rcoef(52,3)/t_elect)
1903
1904	ch71: N2+ + e- --> N + N
1905	!Moffat 2005
1906	!ch71 = 3.5d-7 * (300.d0/t_elect)**(0.5d0)
1907	!Peterson et al. 1998 (from UMIST):
1908	! ch71 = 1.7d-7 * (300.d0/t_elect)**0.3
1909	!Zipf 1980+Kella et al 1996 (from F&S2001):
1910	!ch71 = 2.2d-7 * (300.d0/t_elect)**0.39
1911	ch71=rcoef(53,1)((t_elect/300.)rcoef(53,2))
1912	$ exp(rcoef(53,3)/t_elect)
1913
1914
1915	ch72: N2+ + O3p --> O+ + N2
1916	!Moffat 2005:
1917	!ch72 = 4.1d-10
1918	!McFarland et al. 1974 (From UMIST):
1919	!ch72 = 1.1d-11
1920	!Scott et al., 1999 (from F&S2001):
1921	! ch72 = 7.0d-12 * (300.d0/t_elect)**0.23
1922	ch72=rcoef(54,1)((t_elect/300.)rcoef(54,2))
1923	$ exp(rcoef(54,3)/t_elect)
1924
1925
1926	ch73 CO+ + H --> H+ + CO
1927	!Scott et al., 1997 (from F&S2001):
1928	! ch73 = 4.0d-10
1929	!Federer et al. 1984 (from UMIST):
1930	!ch73 = 7.5d-10
1931	ch73=rcoef(55,1)((t_elect/300.)rcoef(55,2))
1932	$ exp(rcoef(55,3)/t_elect)
1933
1934
1935	ch74: O+ + H --> H+ + O
1936	!Krasnopolsky 2002:
1937	!ch74 = 5.7d-10 * (tt/300.d0)**(0.36d0)
1938	!Stancil et al. 1999 (from UMIST):
1939	! ch74 = 5.66d-10(tt/300.)0.36exp(8.6/tt)
1940	!Aninich 1993 (from F&S2001):
1941	!ch74 = 6.4e-10
1942	ch74=rcoef(56,1)((tcte/300.)rcoef(56,2))
1943	$ exp(rcoef(56,3)/tcte)
1944
1945	ch75: NO+ + e- --> N + O
1946	!Mitchel 1990 (from UMIST):
1947	! ch75 = 4.3d-7 * (300.d0/t_elect)**(0.37d0)
1948	!Vejby-Christensen et al. 1996 (from F&S2001):
1949	!ch75=4.0d-7 * (300.d0/t_elect)**0.5d0
1950	ch75=rcoef(57,1)((t_elect/300.)rcoef(57,2))
1951	$ exp(rcoef(57,3)/t_elect)
1952
1953	ch76: H+ + O3p --> O+ + H
1954	!Krasnopolsky et al. 2002:
1955	!ch76 = 7.3d-10 * (tt/300.d0)*(0.23d0) exp(-226./tt)
1956	!Stancil et al. 1999 (from UMIST):
1957	! ch76 = 6.86e-10* (tt/300.)*0.26exp(-224.3/tt)
1958	ch76=rcoef(58,1)((tcte/300.)rcoef(58,2))
1959	$ exp(rcoef(58,3)/tcte)
1960
1961	ch85: N+ + CO2 --> CO2+ + N
1962	!Krasnopolsky 2002:
1963	!ch85 = 1.d-9
1964	!Adams et al. 1980 (from UMIST):
1965	! ch85 = 7.5d-10
1966	!Aninich et al. 1993 (from F&S2001):
1967	!ch85 = 9.2d-10
1968	ch85=rcoef(59,1)((t_elect/300.)rcoef(59,2))
1969	$ exp(rcoef(59,3)/t_elect)
1970
1971	ch86: H2 + CO2+ --> H + HCO2+
1972	!Scott et al. 1998 (from F&S2001 and K2002):
1973	!ch86 = 8.7d-10
1974	!Copp et al. 1982 (from UMIST):
1975	! ch86 = 9.5d-10
1976	ch86=rcoef(60,1)((t_elect/300.)rcoef(60,2))
1977	$ exp(rcoef(60,3)/t_elect)
1978
1979
1980	c h87: HCO2+ + e -> CO2 + H
1981	!Krasnopolsky 2002:
1982	! ch87 = 3.4d-7 * (300.d0/t_elect)**(0.5d0)
1983	!UMIST 2012: the reactions has 3 different sets of products: CO2+H,
1984	!CO+O+H (dominante) y CO+OH. Habria que tener esto en cuenta
1985	ch87=rcoef(61,1)((t_elect/300.)rcoef(61,2))
1986	$ exp(rcoef(61,3)/t_elect)
1987
1988
1989	endif !Of if(chemthermod.eq.3)
1990
1991	return
1992	end
1993
1994
1995
1996	c**********************************************************************
1997	c**********************************************************************
1998
1999	subroutine lifetimes
2000	& ( ig,i,nlayer,chemthermod,zenit,zx,jdistot8, jdistot8_b, jion8,
2001	$ xtmin, xcompmin, xn_comp_en_EQ,
2002	$ co2xini,o2xini,o3pxini,coxini,hxini,ohxini,ho2xini,h2xini,
2003	$ h2oxini,h2o2xini,o1dxini,o3xini,n2xini,nxini,noxini,no2xini,
2004	$ n2dxini,co2plusxini,oplusxini,o2plusxini,coplusxini,
2005	$ cplusxini,nplusxini,noplusxini,n2plusxini,hplusxini,
2006	$ hco2plusxini,electxini)
2007
2008
2009	c Calculates the lifetime of each species at each time step (itime)
2010	c and each altitude (i), and the minimum of them (tmin)
2011	c It also computes the number of species in PE
2012	c
2013	c
2014	c jul 2008 MA Version en subrutina
2015	c 2009 FGG Adaptation to GCM
2016	c**********************************************************************
2017
2018	use iono_h
2019	use param_v4_h
2020
2021	implicit none
2022
2023	include 'callkeys.h'
2024
2025	c arguments
2026	c
2027	integer i,ig,nlayer ! I. Layer
2028	integer chemthermod
2029	real zenit
2030	real zx(nlayer)
2031	real*8 jdistot8(nabs,nlayer) ! I.
2032	real*8 jdistot8_b(nabs,nlayer) ! I.
2033	real*8 jion8(nabs,nlayer,4) ! I.
2034
2035	real*8 xtmin ! O.
2036	integer xcompmin ! O.
2037	integer xn_comp_en_EQ ! O.
2038
2039	real*8 co2xini,o2xini,o3pxini,coxini
2040	real*8 ho2xini,h2xini,hxini,ohxini
2041	real*8 h2o2xini,o1dxini,o3xini,h2oxini
2042	real*8 nxini,noxini,n2xini,n2dxini,no2xini
2043	real*8 co2plusxini,coplusxini,oplusxini,o2plusxini
2044	real*8 cplusxini,noplusxini,n2plusxini,hplusxini
2045	real*8 electxini,nplusxini,hco2plusxini
2046
2047	c local variables
2048	c
2049	integer j
2050
2051
2052	external ionsec_nplus
2053	real*8 ionsec_nplus
2054
2055	external ionsec_n2plus
2056	real*8 ionsec_n2plus
2057
2058	external ionsec_oplus
2059	real*8 ionsec_oplus
2060
2061	external ionsec_coplus
2062	real*8 ionsec_coplus
2063
2064	external ionsec_co2plus
2065	real*8 ionsec_co2plus
2066
2067	external ionsec_o2plus
2068	real*8 ionsec_o2plus
2069
2070
2071
2072	ccccccccccccccc CODE STARTS
2073
2074	!Initialization of lifetimes
2075	do j = 1, nreact
2076	tauco2(j,i) = 1.d30
2077	tauo2(j,i) = 1.d30
2078	tauo3p(j,i) = 1.d30
2079	tauco(j,i) = 1.d30
2080	tauh(j,i) = 1.d30
2081	tauoh(j,i) = 1.d30
2082	tauho2(j,i) = 1.d30
2083	tauh2(j,i) = 1.d30
2084	tauh2o(j,i) = 1.d30
2085	tauo1d(j,i) = 1.d30
2086	tauh2o2(j,i) = 1.d30
2087	tauo3(j,i) = 1.d30
2088	taun2(j,i) = 1.d30
2089	taun(j,i) = 1.d30
2090	tauno(j,i) = 1.d30
2091	taun2d(j,i) = 1.d30
2092	tauno2(j,i) = 1.d30
2093	tauco2plus(j,i) = 1.d30
2094	tauoplus(j,i) = 1.d30
2095	tauo2plus(j,i) = 1.d30
2096	taucoplus(j,i) = 1.d30
2097	taucplus(j,i) = 1.d30
2098	taunplus(j,i) = 1.d30
2099	taun2plus(j,i) = 1.d30
2100	taunoplus(j,i) = 1.d30
2101	tauhplus(j,i) = 1.d30
2102	tauhco2plus(j,i)= 1.d30
2103	end do
2104
2105	!Lifetime of each species in each reaction
2106	if(jdistot8(1,i).gt.1.d-30) tauco2(1,i) = 1.d0 / jdistot8(1,i)
2107
2108	if(ch2o2xinico2xini.gt.1.d-30)
2109	$ tauh(2,i) = 1.d0 / (ch2 * o2xini * co2xini)
2110	if(ch2hxinico2xini.gt.1.d-30)
2111	$ tauo2(2,i) = 1.d0 / (ch2 * hxini * co2xini)
2112
2113	if(ch3*o3pxini.gt.1.d-30) tauho2(3,i) = 1.d0 /
2114	$ (ch3 * o3pxini)
2115	if(ch3*ho2xini.gt.1.d-30) tauo3p(3,i) = 1.d0 /
2116	$ (ch3 * ho2xini)
2117
2118	if(ch4*coxini.gt.1.d-30) tauoh(4,i) = 1.d0 /
2119	$ (ch4 * coxini)
2120	if(ch4*ohxini.gt.1.d-30) tauco(4,i) = 1.d0 /
2121	$ (ch4 * ohxini)
2122
2123	if(ch5*ho2xini.gt.1.d-30)tauho2(5,i)=1.d0 /
2124	$ (2.d0ch5ho2xini)
2125
2126
2127	if(jdistot8(6,i).gt.1.d-30) tauh2o2(6,i) = 1.d0 / jdistot8(6,i)
2128
2129	if(ch7*ohxini.gt.1.d-30) tauho2(7,i) = 1.d0 /
2130	$ (ch7 * ohxini)
2131	if(ch7*ho2xini.gt.1.d-30) tauoh(7,i) = 1.d0 /
2132	$ (ch7 * ho2xini)
2133
2134	if(jdistot8(4,i).gt.1.d-30) tauh2o(8,i) = 1.d0 / jdistot8(4,i)
2135
2136	if(ch9*o1dxini.gt.1.d-30) tauh2o(9,i) = 1.d0 /
2137	$ (ch9 * o1dxini)
2138	if(ch9*h2oxini.gt.1.d-30) tauo1d(9,i) = 1.d0 /
2139	$ (ch9 * h2oxini)
2140
2141	if(ch10o3pxinico2xini.gt.1.d-30)
2142	$ tauo3p(10,i) = 1.d0 /
2143	$ (2.d0 * ch10 * o3pxini * co2xini)
2144
2145	if(ch11*o3pxini.gt.1.d-30) tauoh(11,i)=1.d0 /
2146	$ (ch11 * o3pxini)
2147	if(ch11*ohxini.gt.1.d-30) tauo3p(11,i) = 1.d0 /
2148	$ (ch11 * ohxini)
2149
2150	if(jdistot8(2,i).gt.1.d-30) tauo2(12,i) = 1.d0 / jdistot8(2,i)
2151
2152	if(ch13*hxini.gt.1.d-30) tauho2(13,i) = 1.d0 /
2153	$ (ch13 * hxini)
2154	if(ch13*ho2xini.gt.1.d-30) tauh(13,i) = 1.d0 /
2155	$ (ch13 * ho2xini)
2156
2157	if(ch14*o1dxini.gt.1.d-30) tauh2(14,i) = 1.d0 /
2158	$ (ch14 * o1dxini)
2159	if(ch14*h2xini.gt.1.d-30) tauo1d(14,i) = 1.d0 /
2160	$ (ch14 * h2xini)
2161
2162	if(ch15*ohxini.gt.1.d-30) tauh2(15,i) = 1.d0 /
2163	$ (ch15 * ohxini)
2164	if(ch15*h2xini.gt.1.d-30) tauoh(15,i) = 1.d0 /
2165	$ (ch15 * h2xini)
2166
2167	if(jdistot8_b(1,i).gt.1.d-30) tauco2(16,i)=1.d0/jdistot8_b(1,i)
2168
2169	if(jdistot8_b(2,i).gt.1.d-30) tauo2(17,i)=1.d0/jdistot8_b(2,i)
2170
2171	if(ch18*ohxini.gt.1.d-30) tauh2o2(18,i)=1.d0 /
2172	$ (ch18 * ohxini)
2173	if(ch18*h2o2xini.gt.1.d-30) tauoh(18,i)=1.d0 /
2174	$ (ch18 * h2o2xini)
2175
2176	if(ch19*co2xini.gt.1.d-30)tauo1d(19,i)=1.d0 /
2177	$ (ch19 * co2xini)
2178
2179	if(ch20*o2xini.gt.1.d-30)tauo1d(20,i)= 1.d0 /
2180	$ (ch20 * o2xini)
2181
2182	if(ch21o2xinico2xini.gt.1.d-30)tauo3p(21,i)= 1.d0 /
2183	$ (ch21 * o2xini * co2xini)
2184	if(ch21o3pxinico2xini.gt.1.d-30) tauo2(21,i) = 1.d0 /
2185	$ (ch21 * o3pxini * co2xini)
2186
2187	!Only if O3, N or ion chemistry requested
2188	if(chemthermod.ge.1) then
2189	if(ch22*hxini.gt.1.d-30) tauo3(22,i) = 1.d0 /
2190	$ (ch22 * hxini)
2191	if(ch22*o3xini.gt.1.d-30) tauh(22,i) = 1.d0 /
2192	$ (ch22 * o3xini)
2193
2194	if(ch23*ohxini.gt.1.d-30) tauo3(23,i) = 1.d0 /
2195	$ (ch23 * ohxini)
2196	if(ch23*o3xini.gt.1.d-30) tauoh(23,i) = 1.d0 /
2197	$ (ch23 * o3xini)
2198
2199	if(ch24 * ho2xini.gt.1.d-30)tauo3(24,i)= 1.d0 /
2200	$ (ch24 * ho2xini)
2201	if(ch24 * o3xini.gt.1.d-30)tauho2(24,i)= 1.d0 /
2202	$ (ch24 * o3xini)
2203
2204	if(jdistot8(7,i).gt.1.d-30) tauo3(25,i)=1.d0 /
2205	$ jdistot8(7,i)
2206
2207	if(jdistot8_b(7,i).gt.1.d-30) tauo3(26,i)= 1.d0 /
2208	$ jdistot8_b(7,i)
2209
2210	endif !Of chemthermod.ge.1
2211
2212	if(jdistot8(5,i).gt.1.d-30) tauh2(27,i)= 1.d0/jdistot8(5,i)
2213
2214	!Only if N or ion chemistry requested
2215	if(chemthermod.ge.2) then
2216	if(jdistot8(8,i).gt.1.d-30) taun2(28,i) = 1.d0 /
2217	$ jdistot8(8,i)
2218
2219	if(jdistot8(10,i).gt.1.d-30) tauno(29,i) = 1.d0 /
2220	$ jdistot8(10,i)
2221
2222	if(ch30 * noxini.gt.1.d-30) taun(30,i) = 1.d0 /
2223	$ (ch30 * noxini)
2224	if(ch30 * nxini.gt.1.d-30) tauno(30,i) = 1.d0 /
2225	$ (ch30 * nxini)
2226
2227	if(ch31 * o1dxini.gt.1.d-30) taun2(31,i) = 1.d0 /
2228	$ (ch31 * o1dxini)
2229	if(ch31 * n2xini.gt.1.d-30) tauo1d(31,i) = 1.d0 /
2230	$ (ch31 * n2xini)
2231
2232	if(ch32 * o2xini.gt.1.d-30) taun(32,i) = 1.d0 /
2233	$ (ch32 * o2xini)
2234	if(ch32 * nxini.gt.1.d-30) tauo2(32,i) = 1.d0 /
2235	$ (ch32 * nxini)
2236
2237	if(ch33 * ohxini.gt.1.d-30) taun(33,i) = 1.d0 /
2238	$ (ch33 * ohxini)
2239	if(ch33 * nxini.gt.1.d-30) tauoh(33,i) = 1.d0 /
2240	$ (ch33 * nxini)
2241
2242	if(ch34 * o3xini.gt.1.d-30) taun(34,i) = 1.d0 /
2243	$ (ch34 * o3xini)
2244	if(ch34 * nxini.gt.1.d-30) tauo3(34,i) = 1.d0 /
2245	$ (ch34 * nxini)
2246
2247	if(ch35 * ho2xini.gt.1.d-30) taun(35,i) = 1.d0 /
2248	$ (ch35 * ho2xini)
2249	if(ch35 * nxini.gt.1.d-30) tauho2(35,i) = 1.d0 /
2250	$ (ch35 * nxini)
2251
2252	if(ch36 * o3pxini.gt.1.d-30) taun2d(36,i) = 1.d0 /
2253	$ (ch36 * o3pxini)
2254	if(ch36 * n2dxini.gt.1.d-30) tauo3p(36,i) = 1.d0 /
2255	$ (ch36 * n2dxini)
2256
2257	if(ch37 * n2xini.gt.1.d-30) taun2d(37,i) = 1.d0 /
2258	$ (ch37 * n2xini)
2259	if(ch37 * n2dxini.gt.1.d-30) taun2(37,i) = 1.d0 /
2260	$ (ch37 * n2dxini)
2261
2262	if(ch38 * co2xini.gt.1.d-30) taun2d(38,i) = 1.d0 /
2263	$ (ch38 * co2xini)
2264	if(ch38 * n2dxini.gt.1.d-30) tauco2(38,i) = 1.d0 /
2265	$ (ch38 * n2dxini)
2266
2267	if(ch39 * ho2xini.gt.1.d-30) tauno(39,i) = 1.d0 /
2268	$ (ch39 * ho2xini)
2269	if(ch39 * noxini.gt.1.d-30) tauho2(39,i) = 1.d0 /
2270	$ (ch39 * noxini)
2271
2272	if(ch40 * noxini * co2xini.gt.1.d-30) tauo3p(40,i) = 1.d0 /
2273	$ (ch40 * noxini * co2xini)
2274	if(ch40 * o3pxini * co2xini.gt.1.d-30) tauno(40,i) = 1.d0 /
2275	$ (ch40 * o3pxini * co2xini)
2276
2277	if(ch41 * no2xini.gt.1.d-30) tauo3p(41,i) = 1.d0 /
2278	$ (ch41 * no2xini)
2279	if(ch41 * o3pxini.gt.1.d-30) tauno2(41,i) = 1.d0 /
2280	$ (ch41 * o3pxini)
2281
2282	if(ch42 * noxini.gt.1.d-30) tauo3(42,i) = 1.d0 /
2283	$ (ch42 * noxini)
2284	if(ch42 * o3xini.gt.1.d-30) tauno(42,i) = 1.d0 /
2285	$ (ch42 * o3xini)
2286
2287	if(ch43 * no2xini.gt.1.d-30) tauh(43,i) = 1.d0 /
2288	$ (ch43 * no2xini)
2289	if(ch43 * hxini.gt.1.d-30) tauno2(43,i) = 1.d0 /
2290	$ (ch43 * hxini)
2291
2292	if(jdistot8(13,i).gt.1.d-30) tauno2(44,i) = 1.d0 /
2293	$ jdistot8(13,i)
2294
2295	if(ch45 * nxini.gt.1.d-30) tauo3p(45,i) = 1.d0 /
2296	$ (ch45 * nxini)
2297	if(ch45 * o3pxini.gt.1.d-30) taun(45,i) = 1.d0 /
2298	$ (ch45 * o3pxini)
2299
2300	endif !Of chemthermod.ge.2
2301
2302	c>>>>>>>>>>>>>>>>>>>>>>>> IONOSPHERE >>>>>>>>>>>>>>>>
2303
2304	!Only if ion chemistry requested
2305	if(chemthermod.eq.3) then
2306	if(ch46 * co2plusxini .gt.1.d-30) tauo2(46,i) =
2307	@ 1.d0/(ch46*co2plusxini)
2308	if(ch46 * o2xini .gt.1.d-30) tauco2plus(46,i) =
2309	@ 1.d0/(ch46*o2xini)
2310
2311	if ( ch47*o3pxini .gt. 1.d-30 ) tauco2plus(47,i) =
2312	@ 1.d0/( ch47*o3pxini )
2313	if ( ch47*co2plusxini .gt. 1.d-30 ) tauo3p(47,i) =
2314	@ 1.d0/( ch47*co2plusxini )
2315
2316	if ( ch48*o3pxini .gt. 1.d-30 ) tauco2plus(48,i) =
2317	@ 1.d0/(ch48*o3pxini)
2318	if ( ch48*co2plusxini .gt. 1.d-30 ) tauo3p(48,i) =
2319	@ 1.d0/(ch48*co2plusxini)
2320
2321	if ( ch49*electxini .gt. 1.d-30 ) tauo2plus(49,i) =
2322	@ 1.d0/(ch49*electxini)
2323
2324	if ( ch50*co2xini .gt. 1.d-30 ) tauoplus(50,i) =
2325	@ 1.d0/(ch50*co2xini)
2326	if ( ch50*oplusxini .gt. 1.d-30 ) tauco2(50,i) =
2327	@ 1.d0/(ch50*oplusxini)
2328
2329	if ( jion8(1,i,1).gt.1.d-30 ) tauco2(51,i) =
2330	$ 1.d0 / jion8(1,i,1)
2331
2332	if ( jion8(1,i,2).gt.1.d-30 ) tauco2(52,i) =
2333	$ 1.d0 / jion8(1,i,2)
2334
2335	if ( jion8(1,i,3).gt.1.d-30 ) tauco2(53,i) =
2336	$ 1.d0 / jion8(1,i,3)
2337
2338	if ( jion8(1,i,4).gt.1.d-30 ) tauco2(54,i) =
2339	$ 1.d0 / jion8(1,i,4)
2340
2341	if ( ch55*electxini .gt. 1.d-30 ) tauco2plus(55,i) =
2342	@ 1.d0/(ch55*electxini)
2343
2344	if ( ch56*oplusxini .gt. 1.d-30 ) tauco2(56,i) =
2345	@ 1.d0/(ch56*oplusxini)
2346	if ( ch56*co2xini .gt. 1.d-30 ) tauoplus(56,i) =
2347	@ 1.d0/(ch56*co2xini)
2348
2349	if ( ch57*coplusxini .gt. 1.d-30 ) tauco2(57,i) =
2350	@ 1.d0/(ch57*coplusxini)
2351	if ( ch57*co2xini .gt. 1.d-30 ) taucoplus(57,i) =
2352	@ 1.d0/(ch57*co2xini)
2353
2354	if ( ch58*coplusxini .gt. 1.d-30 ) tauo3p(58,i) =
2355	@ 1.d0/(ch58*coplusxini)
2356	if ( ch58*o3pxini .gt. 1.d-30 ) taucoplus(58,i) =
2357	@ 1.d0/(ch58*o3pxini)
2358
2359	if ( ch59*cplusxini .gt. 1.d-30 ) tauco2(59,i) =
2360	@ 1.d0/(ch59*cplusxini)
2361	if ( ch59*co2xini .gt. 1.d-30 ) taucplus(59,i) =
2362	@ 1.d0/(ch59*co2xini)
2363
2364	if ( jion8(2,i,1).gt.1.d-30 ) tauo2(60,i) =
2365	$ 1.d0 / jion8(2,i,1)
2366
2367	if ( jion8(3,i,1).gt.1.d-30 ) tauo3p(61,i) =
2368	$ 1.d0 / jion8(3,i,1)
2369
2370	if ( ch62*co2plusxini .gt. 1.d-30 ) tauno(62,i) =
2371	@ 1.d0/(ch62*co2plusxini)
2372	if ( ch62*noxini .gt. 1.d-30 ) tauco2plus(62,i) =
2373	@ 1.d0/(ch62*noxini)
2374
2375	if ( ch63*co2plusxini .gt. 1.d-30 ) taun(63,i) =
2376	@ 1.d0/(ch63*cplusxini)
2377	if ( ch63*nxini .gt. 1.d-30 ) tauco2plus(63,i) =
2378	@ 1.d0/(ch63*nxini)
2379
2380	if ( ch64*o2plusxini .gt. 1.d-30 ) tauno(64,i) =
2381	@ 1.d0/(ch64*o2plusxini)
2382	if ( ch64*noxini .gt. 1.d-30 ) tauo2plus(64,i) =
2383	@ 1.d0/(ch64*noxini)
2384
2385	if ( ch65*o2plusxini .gt. 1.d-30 ) taun2(65,i) =
2386	@ 1.d0/(ch65*o2plusxini)
2387	if ( ch65*n2xini .gt. 1.d-30 ) tauo2plus(65,i) =
2388	@ 1.d0/(ch65*n2xini)
2389
2390	if ( ch66*o2plusxini .gt. 1.d-30 ) taun(66,i) =
2391	@ 1.d0/(ch66*o2plusxini)
2392	if ( ch66*nxini .gt. 1.d-30 ) tauo2plus(66,i) =
2393	@ 1.d0/(ch66*nxini)
2394
2395	if ( ch67*oplusxini .gt. 1.d-30 ) taun2(67,i) =
2396	@ 1.d0/(ch67*oplusxini)
2397	if ( ch67*n2xini .gt. 1.d-30 ) tauoplus(67,i) =
2398	@ 1.d0/(ch67*n2xini)
2399
2400	if ( ch68*n2plusxini .gt. 1.d-30 ) tauco2(68,i) =
2401	@ 1.d0/(ch68*n2plusxini)
2402	if ( ch68*co2xini .gt. 1.d-30 ) taun2plus(68,i) =
2403	@ 1.d0/(ch68*co2xini)
2404
2405	if ( ch69*cplusxini .gt. 1.d-30 ) tauco2(69,i) =
2406	@ 1.d0/(ch69*cplusxini)
2407	if ( ch69*co2xini .gt. 1.d-30 ) taucplus(69,i) =
2408	@ 1.d0/(ch69*co2xini)
2409
2410	if ( ch70*n2plusxini .gt. 1.d-30 ) tauco(70,i) =
2411	@ 1.d0/(ch70*n2plusxini)
2412	if ( ch70*coxini .gt. 1.d-30 ) taun2plus(70,i) =
2413	@ 1.d0/(ch70*coxini)
2414
2415	if ( ch71*electxini .gt. 1.d-30 ) taun2plus(71,i) =
2416	@ 1.d0/(ch71*electxini)
2417
2418	if ( ch72*n2plusxini .gt. 1.d-30 ) tauo3p(72,i) =
2419	@ 1.d0/(ch72*n2plusxini)
2420	if ( ch72*o3pxini .gt. 1.d-30 ) taun2plus(72,i) =
2421	@ 1.d0/(ch72*o3pxini)
2422
2423	if ( ch73*coplusxini .gt. 1.d-30 ) tauh(73,i) =
2424	@ 1.d0/(ch73*coplusxini)
2425	if ( ch73*hxini .gt. 1.d-30 ) taucoplus(73,i) =
2426	@ 1.d0/(ch73*hxini)
2427
2428	if ( ch74*oplusxini .gt. 1.d-30 ) tauh(74,i) =
2429	@ 1.d0/(ch74*oplusxini)
2430	if ( ch74*hxini .gt. 1.d-30 ) tauoplus(74,i) =
2431	@ 1.d0/(ch74*hxini)
2432
2433	if ( ch75*electxini .gt. 1.d-30 ) taunoplus(75,i) =
2434	@ 1.d0/(ch75*electxini)
2435
2436	if ( ch76*hplusxini .gt. 1.d-30 ) tauo3p(76,i) =
2437	@ 1.d0/(ch76*hplusxini)
2438	if ( ch76*o3pxini .gt. 1.d-30 ) tauhplus(76,i) =
2439	@ 1.d0/(ch76*o3pxini)
2440
2441	if( jion8(11,i,1).gt.1.d-30 ) tauco(77,i) =
2442	$ 1.d0 / jion8(11,i,1)
2443
2444	if( jion8(11,i,2).gt.1.d-30 ) tauco(78,i) =
2445	$ 1.d0 / jion8(11,i,2)
2446
2447	!if( jion8(11,i,3).gt.1.d-30 ) tauco(79,i) =
2448	! $ 1.d0 / jion8(11,i,3)
2449
2450	if( jion8(10,i,1).gt.1.d-30 ) tauno(80,i) =
2451	$ 1.d0 / jion8(10,i,1)
2452
2453	if( jion8(8,i,1).gt.1.d-30 ) taun2(81,i) =
2454	$ 1.d0 / jion8(8,i,1)
2455
2456	if( jion8(8,i,2).gt.1.d-30 ) taun2(82,i) =
2457	$ 1.d0 / jion8(8,i,2)
2458
2459	if( jion8(12,i,1).gt.1.d-30 ) tauh(83,i) =
2460	$ 1.d0 / jion8(12,i,1)
2461
2462	if( jion8(9,i,1).gt.1.d-30 ) taun(84,i) =
2463	$ 1.d0 / jion8(9,i,1)
2464
2465	if ( ch85*nplusxini .gt. 1.d-30 ) tauco2(85,i) =
2466	@ 1.d0/(ch85*nplusxini)
2467	if ( ch85*co2xini .gt. 1.d-30 ) taunplus(85,i) =
2468	@ 1.d0/(ch85*co2xini)
2469
2470	if ( ch86*co2plusxini .gt. 1.d-30) tauh2(86,i) =
2471	$ 1.d0/(ch86*co2plusxini)
2472	if ( ch86*h2xini .gt. 1.d-30) tauco2plus(86,i) =
2473	$ 1.d0/(ch86*h2xini)
2474
2475	if ( ch87*electxini .gt. 1.d-30) tauhco2plus(87,i) =
2476	$ 1.d0/(ch87*electxini)
2477
2478	if ( jion8(9,i,1)*ionsec_nplus(zenit,zx(i)).gt.1.d-30)
2479	$ taun(88,i) = 1.d0 /
2480	$ (jion8(9,i,1)*ionsec_nplus(zenit,zx(i)))
2481
2482	if ( jion8(8,i,1)*ionsec_n2plus(zenit,zx(i)).gt.1.d-30)
2483	$ taun2(89,i) = 1.d0 /
2484	$ (jion8(8,i,1)*ionsec_n2plus(zenit,zx(i)))
2485
2486	if ( jion8(3,i,1)*ionsec_oplus(zenit,zx(i)).gt.1.d-30)
2487	$ tauo3p(90,i) = 1.d0 /
2488	$ (jion8(3,i,1)*ionsec_oplus(zenit,zx(i)))
2489
2490	if (jion8(11,i,1)*ionsec_coplus(zenit,zx(i)).gt.1.d-30)
2491	$ tauco(91,i) = 1.d0 /
2492	$ (jion8(11,i,1)*ionsec_coplus(zenit,zx(i)))
2493
2494	if (jion8(1,i,1)*ionsec_co2plus(zenit,zx(i)).gt.1.d-30)
2495	$ tauco2(92,i) = 1.d0 /
2496	$ (jion8(1,i,1)*ionsec_co2plus(zenit,zx(i)))
2497
2498	if ( jion8(2,i,1)*ionsec_o2plus(zenit,zx(i)).gt.1.d-30)
2499	$ tauo2(93,i) = 1.d0 /
2500	$ (jion8(2,i,1)*ionsec_o2plus(zenit,zx(i)))
2501
2502	endif !Of chemthermod.eq.3
2503	c>>>>>>>>>>>>>>>>>>>>>>>>
2504
2505	!Minimum lifetime for each species
2506	tminco2(i) = 1.d30
2507	tmino2(i) = 1.d30
2508	tmino3p(i) = 1.d30
2509	tminco(i) = 1.d30
2510	tminh(i) = 1.d30
2511	tminoh(i) = 1.d30
2512	tminho2(i) = 1.d30
2513	tminh2(i) = 1.d30
2514	tminh2o(i) = 1.d30
2515	tmino1d(i) = 1.d30
2516	tminh2o2(i) = 1.d30
2517	tmino3(i) = 1.d30
2518	tminn(i) = 1.d30
2519	tminno(i) = 1.d30
2520	tminn2(i) = 1.d30
2521	tminn2d(i) = 1.d30
2522	tminno2(i) = 1.d30
2523	tminco2plus(i) = 1.d30
2524	tminoplus(i) = 1.d30
2525	tmino2plus(i) = 1.d30
2526	tmincoplus(i) = 1.d30
2527	tmincplus(i) = 1.d30
2528	tminnplus(i) = 1.d30
2529	tminn2plus(i) = 1.d30
2530	tminnoplus(i) = 1.d30
2531	tminhplus(i) = 1.d30
2532	tminhco2plus(i)=1.d30
2533
2534	do j=1,nreact
2535	tminco2(i) = min(tminco2(i),tauco2(j,i))
2536	tmino2(i) = min(tmino2(i),tauo2(j,i))
2537	tmino3p(i) = min(tmino3p(i),tauo3p(j,i))
2538	tminco(i) = min(tminco(i),tauco(j,i))
2539	tminh(i) = min(tminh(i),tauh(j,i))
2540	tminoh(i) = min(tminoh(i),tauoh(j,i))
2541	tminho2(i) = min(tminho2(i),tauho2(j,i))
2542	tminh2(i) = min(tminh2(i),tauh2(j,i))
2543	tminh2o(i) = min(tminh2o(i),tauh2o(j,i))
2544	tmino1d(i) = min(tmino1d(i),tauo1d(j,i))
2545	tminh2o2(i) = min(tminh2o2(i),tauh2o2(j,i))
2546	tmino3(i) = min(tmino3(i),tauo3(j,i))
2547	tminn(i) = min(tminn(i),taun(j,i))
2548	tminno(i) = min(tminno(i),tauno(j,i))
2549	tminn2(i) = min(tminn2(i),taun2(j,i))
2550	tminn2d(i) = min(tminn2d(i),taun2d(j,i))
2551	tminno2(i) = min(tminno2(i),tauno2(j,i))
2552	!
2553	tminco2plus(i) = min(tminco2plus(i),tauco2plus(j,i))
2554	tminoplus(i) = min(tminoplus(i),tauoplus(j,i))
2555	tmino2plus(i) = min(tmino2plus(i),tauo2plus(j,i))
2556	tmincoplus(i) = min(tmincoplus(i),taucoplus(j,i))
2557	tmincplus(i) = min(tmincplus(i),taucplus(j,i))
2558	tminnplus(i) = min(tminnplus(i),taunplus(j,i))
2559	tminn2plus(i) = min(tminn2plus(i),taun2plus(j,i))
2560	tminnoplus(i) = min(tminnoplus(i),taunoplus(j,i))
2561	tminhplus(i) = min(tminhplus(i),tauhplus(j,i))
2562	tminhco2plus(i)= min(tminhco2plus(i),tauhco2plus(j,i))
2563	end do
2564
2565	!!! Minimum lifetime for all species
2566
2567	xtmin = 1.d30
2568
2569	! Neutrals that can not be in PE
2570
2571	xtmin = min(xtmin,tminco2(i))
2572	if(xtmin.eq.tminco2(i)) xcompmin=1
2573	xtmin = min(xtmin,tmino2(i))
2574	if(xtmin.eq.tmino2(i)) xcompmin=2
2575	xtmin = min(xtmin,tmino3p(i))
2576	if(xtmin.eq.tmino3p(i)) xcompmin=3
2577	xtmin = min(xtmin,tminco(i))
2578	if(xtmin.eq.tminco(i)) xcompmin=4
2579
2580	xtmin = min(xtmin,tminh2(i))
2581	if(xtmin.eq.tminh2(i)) xcompmin=8
2582	xtmin = min(xtmin,tminh2o(i))
2583	if(xtmin.eq.tminh2o(i)) xcompmin=9
2584
2585	xtmin = min(xtmin,tminh2o2(i))
2586	if(xtmin.eq.tminh2o2(i)) xcompmin=11
2587	xtmin = min(xtmin,tmino3(i))
2588	if(xtmin.eq.tmino3(i)) xcompmin=12
2589	xtmin = min(xtmin,tminn(i))
2590	if(xtmin.eq.tminn(i)) xcompmin=13
2591	xtmin = min(xtmin,tminno(i))
2592	if(xtmin.eq.tminno(i)) xcompmin=14
2593	xtmin = min(xtmin,tminn2(i))
2594	if(xtmin.eq.tminn2(i)) xcompmin=15
2595
2596	! Neutrals that can be in PE
2597	! [ O1D , OH , HO2 , H , N2D , NO2 ]
2598
2599	xn_comp_en_EQ = 0
2600
2601	! H
2602	h_eq(i)='Y'
2603	xn_comp_en_EQ = xn_comp_en_EQ + 1
2604
2605	! OH
2606	oh_eq(i)='Y'
2607	xn_comp_en_EQ = xn_comp_en_EQ + 1
2608
2609	! HO2
2610	ho2_eq(i)='Y'
2611	xn_comp_en_EQ = xn_comp_en_EQ + 1
2612
2613	! O1D
2614	o1d_eq(i)='Y'
2615	xn_comp_en_EQ = xn_comp_en_EQ + 1
2616
2617	! O3
2618	o3_eq(i)='N'
2619	! o3_eq(i)='Y'
2620	! xn_comp_en_EQ = xn_comp_en_EQ + 1
2621
2622
2623	!Only if N or ion chemistry requested
2624	if(chemthermod.ge.2) then
2625	! N2D
2626	n2d_eq(i)='Y'
2627	xn_comp_en_EQ = xn_comp_en_EQ + 1
2628
2629	! NO2
2630	no2_eq(i)='Y'
2631	xn_comp_en_EQ = xn_comp_en_EQ + 1
2632
2633
2634	! NO
2635	no_eq(i)='N'
2636
2637	! no_eq(i)='Y'
2638	! xn_comp_en_EQ = xn_comp_en_EQ + 1
2639
2640	endif !Of chemthermod.ge.2
2641
2642	! Ions
2643	!Only if ion chemistry requested
2644
2645	if(chemthermod.eq.3) then
2646	! C+
2647	cplus_eq(i)='Y'
2648	xn_comp_en_EQ = xn_comp_en_EQ + 1
2649
2650	! CO+
2651	coplus_eq(i)='Y'
2652	xn_comp_en_EQ = xn_comp_en_EQ + 1
2653
2654	! O+
2655	! oplus_eq(i)='N'
2656	oplus_eq(i)='Y'
2657	xn_comp_en_EQ = xn_comp_en_EQ + 1
2658
2659	! N2+
2660	n2plus_eq(i)='Y'
2661	xn_comp_en_EQ = xn_comp_en_EQ + 1
2662
2663	! H+
2664	! hplus_eq(i)='N'
2665
2666	hplus_eq(i)='Y'
2667	xn_comp_en_EQ = xn_comp_en_EQ + 1
2668
2669	! CO2+
2670	co2plus_eq(i)='Y'
2671	xn_comp_en_EQ = xn_comp_en_EQ + 1
2672
2673	! O2+
2674	o2plus_eq(i)='Y'
2675	xn_comp_en_EQ = xn_comp_en_EQ + 1
2676
2677	! NO+
2678	noplus_eq(i)='Y'
2679	xn_comp_en_EQ = xn_comp_en_EQ + 1
2680
2681	! N+
2682	nplus_eq(i)='Y'
2683	xn_comp_en_EQ = xn_comp_en_EQ + 1
2684
2685	! HCO2+
2686	hco2plus_eq(i)='Y'
2687	xn_comp_en_EQ = xn_comp_en_EQ + 1
2688	endif !Of chemthermod.eq.3
2689
2690	return
2691	c END
2692	end
2693
2694
2695
2696	c**********************************************************************
2697	c**********************************************************************
2698
2699	subroutine timemarching( ig,i,nlayer,chemthermod,n_comp_en_EQ,
2700	$ compmin,tmin,timefrac_sec, deltat,fmargin1 )
2701
2702
2703	c Calculates the timestep of the model, deltat, and verifies if the species
2704	c that can be in PE actually verify or not the PE condition
2705	c
2706	c jul 2008 MA Version en subrutina
2707	c 2009 FGG Adaptation to GCM
2708	c**********************************************************************
2709
2710	use iono_h
2711	use param_v4_h, only: tminco2,tmino2,tmino3p,tminco,tminh,tminoh,
2712	. tminho2,tminh2,tminh2o,tmino1d,tminh2o2,tmino3,tminn,tminno,
2713	. tminno2,tminn2,tminn2d,tminco2plus,tminoplus,tmino2plus,
2714	. tmincoplus,tmincplus,tminnplus,tminnoplus,tminn2plus,
2715	. tminhplus,tminhco2plus
2716
2717	implicit none
2718
2719	c arguments
2720	c
2721	integer i,ig,nlayer ! I. Layer
2722	integer chemthermod
2723	integer n_comp_en_EQ(nlayer) ! Number of species in PE
2724	integer compmin(nlayer) ! Species with minimum lifetime
2725	real*8 tmin(nlayer) ! Minimum lifetime
2726	real*8 timefrac_sec ! I.
2727	real*8 deltat ! O. TimeMarching step
2728
2729	c local variables
2730	c
2731	integer j
2732	real*8 tminaux
2733
2734	real*8 fmargin1, fmargin2
2735
2736	ccccccccccccccc CODE STARTS
2737
2738	! fmargin1=1.
2739	fmargin2=5.
2740
2741	!Internal timestep as the minimum of the external (GCM) timestep
2742	!and the minimum lifetime of the species not in PE divided by a
2743	!given factor
2744	tminaux = min( timefrac_sec, tmin(i)/fmargin1 )
2745
2746	!Given the internal timestep, we verify if the species that can be
2747	!in PE verify or not the PE condition (lifetime < deltat/fmargin2)
2748
2749	! 6 neutral species that can be in PE
2750	! [ O1D , OH , HO2 , H , N2D , NO2 ]
2751
2752	! O1D
2753	if ( (o1d_eq(i).eq.'Y') .and.
2754	& (tmino1d(i).gt.tminaux/fmargin2) ) then
2755	o1d_eq(i)='N'
2756	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2757	endif
2758	! OH
2759	if ( (oh_eq(i).eq.'Y') .and.
2760	& (tminoh(i).gt.tminaux/fmargin2) ) then
2761	oh_eq(i)='N'
2762	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2763	endif
2764	! HO2
2765	if ( (ho2_eq(i).eq.'Y') .and.
2766	& (tminho2(i).gt.tminaux/fmargin2) ) then
2767	ho2_eq(i)='N'
2768	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2769	endif
2770	! H
2771	if ( (h_eq(i).eq.'Y') .and.
2772	& (tminh(i).gt.tminaux/fmargin2) ) then
2773	h_eq(i)='N'
2774	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2775	endif
2776
2777	!Only if N or ion chemistry requested
2778	if(chemthermod.ge.2) then
2779	! N2D
2780	if ( (n2d_eq(i).eq.'Y') .and.
2781	& (tminn2d(i).gt.tminaux/fmargin2) ) then
2782	n2d_eq(i)='N'
2783	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2784	endif
2785	! NO2
2786	if ( (no2_eq(i).eq.'Y') .and.
2787	& (tminno2(i).gt.tminaux/fmargin2) ) then
2788	no2_eq(i)='N'
2789	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2790	endif
2791
2792	endif !Of chemthermod.ge.2
2793
2794
2795	!
2796	! 9 ions
2797	!
2798
2799	!Only if ion chemistry requested
2800	if(chemthermod.eq.3) then
2801	! C+
2802	if ( (cplus_eq(i).eq.'Y') .and.
2803	& (tmincplus(i).gt.tminaux/fmargin2) ) then
2804	cplus_eq(i)='N'
2805	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2806	endif
2807	! CO+
2808	if ( (coplus_eq(i).eq.'Y') .and.
2809	& (tmincoplus(i).gt.tminaux/fmargin2) ) then
2810	coplus_eq(i)='N'
2811	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2812	endif
2813	! O+
2814	if ( (oplus_eq(i).eq.'Y') .and.
2815	& (tminoplus(i).gt.tminaux/fmargin2) ) then
2816	oplus_eq(i)='N'
2817	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2818	endif
2819	! N2+
2820	if ( (n2plus_eq(i).eq.'Y') .and.
2821	& (tminn2plus(i).gt.tminaux/fmargin2) ) then
2822	n2plus_eq(i)='N'
2823	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2824	endif
2825	! H+
2826	if ( (hplus_eq(i).eq.'Y') .and.
2827	& (tminhplus(i).gt.tminaux/fmargin2) ) then
2828	hplus_eq(i)='N'
2829	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2830	endif
2831	! CO2+
2832	if ( (co2plus_eq(i).eq.'Y') .and.
2833	& (tminco2plus(i).gt.tminaux/fmargin2) ) then
2834	co2plus_eq(i)='N'
2835	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2836	endif
2837	! O2+
2838	if ( (o2plus_eq(i).eq.'Y') .and.
2839	& (tmino2plus(i).gt.tminaux/fmargin2) ) then
2840	o2plus_eq(i)='N'
2841	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2842	endif
2843	! NO+
2844	if ( (noplus_eq(i).eq.'Y') .and.
2845	& (tminnoplus(i).gt.tminaux/fmargin2) ) then
2846	noplus_eq(i)='N'
2847	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2848	endif
2849	! N+
2850	if ( (nplus_eq(i).eq.'Y') .and.
2851	& (tminnplus(i).gt.tminaux/fmargin2) ) then
2852	nplus_eq(i)='N'
2853	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2854	endif
2855	! HCO2+
2856	if ( (hco2plus_eq(i).eq.'Y') .and.
2857	& (tminhco2plus(i).gt.tminaux/fmargin2) ) then
2858	hco2plus_eq(i)='N'
2859	n_comp_en_EQ(i) = n_comp_en_EQ(i) - 1
2860	endif
2861
2862	endif !Of chemthermod.eq.3
2863
2864
2865
2866	! And we set the internal timestep
2867	!
2868	deltat = tminaux
2869
2870	return
2871	c END
2872	end
2873
2874
2875
2876	c**********************************************************************
2877	c**********************************************************************
2878
2879	subroutine prodsandlosses ( ig,i,nlayer,chemthermod,zenit,zx,
2880	& jdistot8, jdistot8_b, jion8,
2881	& co2xinput, o2xinput, o3pxinput,
2882	& coxinput, h2xinput, o3xinput,
2883	& h2oxinput, nxinput, noxinput,
2884	& h2o2xinput, n2xinput,
2885	& o1dxinput, ohxinput, ho2xinput,
2886	& hxinput, n2dxinput, no2xinput,
2887	& co2plusxinput, o2plusxinput, coplusxinput,
2888	& oplusxinput, cplusxinput, noplusxinput,
2889	& n2plusxinput, hplusxinput, nplusxinput,
2890	& hco2plusxinput,electxinput )
2891
2892
2893
2894	c Calculates the productions and losses of each species in each reaction
2895	c and each altitude
2896	c
2897	c jul 2008 MA Version en subrutina
2898	c apr 2009 FGG Compact version to save CPU time, adaptation
2899	c to GCM
2900	c**********************************************************************
2901
2902	use param_v4_h
2903	implicit none
2904
2905	c arguments
2906	c
2907	integer ig,nlayer
2908	integer i ! I. Layer
2909	integer chemthermod
2910	real zx(nlayer)
2911	real zenit
2912	real*8 jdistot8(nabs,nlayer)
2913	real*8 jdistot8_b(nabs,nlayer)
2914	real*8 jion8(nabs,nlayer,4)
2915	real*8 co2xinput,o2xinput,o3pxinput,coxinput
2916	real*8 ho2xinput,h2xinput,hxinput,ohxinput
2917	real*8 h2o2xinput,o1dxinput,o3xinput,h2oxinput
2918	real*8 nxinput,noxinput,n2xinput,n2dxinput,no2xinput
2919	real*8 co2plusxinput,coplusxinput,oplusxinput,o2plusxinput
2920	real*8 cplusxinput,noplusxinput,n2plusxinput,hplusxinput
2921	real*8 electxinput,nplusxinput,hco2plusxinput
2922
2923	c local variables
2924	c
2925	integer j
2926
2927	external ionsec_nplus
2928	real*8 ionsec_nplus
2929
2930	external ionsec_n2plus
2931	real*8 ionsec_n2plus
2932
2933	external ionsec_oplus
2934	real*8 ionsec_oplus
2935
2936	external ionsec_coplus
2937	real*8 ionsec_coplus
2938
2939	external ionsec_co2plus
2940	real*8 ionsec_co2plus
2941
2942	external ionsec_o2plus
2943	real*8 ionsec_o2plus
2944
2945	logical firstcall
2946	save firstcall
2947
2948	!$OMP THREADPRIVATE(firstcall)
2949
2950	data firstcall /.true./
2951
2952	ccccccccccccccc CODE STARTS
2953
2954	!Initialization'
2955	! if (firstcall) then
2956	! firstcall= .false.
2957	do j=1,nreact
2958	Lco2(i,j) = 0.d0
2959	Pco2(i,j) = 0.d0
2960	Lo2(i,j) = 0.d0
2961	Po2(i,j) = 0.d0
2962	Lo3p(i,j) = 0.d0
2963	Po3p(i,j) = 0.d0
2964	Lco(i,j) = 0.d0
2965	Pco(i,j) = 0.d0
2966	Ph(i,j) = 0.d0
2967	Lh(i,j) = 0.d0
2968	Poh(i,j) = 0.d0
2969	Loh(i,j) = 0.d0
2970	Pho2(i,j) = 0.d0
2971	Lho2(i,j) = 0.d0
2972	Ph2(i,j) = 0.d0
2973	Lh2(i,j) = 0.d0
2974	Ph2o(i,j) = 0.d0
2975	Lh2o(i,j) = 0.d0
2976	Po1d(i,j) = 0.d0
2977	Lo1d(i,j) = 0.d0
2978	Ph2o2(i,j) = 0.d0
2979	Lh2o2(i,j) = 0.d0
2980	Po3(i,j) = 0.d0
2981	Lo3(i,j) = 0.d0
2982	Pn(i,j) = 0.d0
2983	Ln(i,j) = 0.d0
2984	Pno(i,j) = 0.d0
2985	Lno(i,j) = 0.d0
2986	Pn2(i,j) = 0.d0
2987	Ln2(i,j) = 0.d0
2988	Pn2d(i,j) = 0.d0
2989	Ln2d(i,j) = 0.d0
2990	Pno2(i,j) = 0.d0
2991	Lno2(i,j) = 0.d0
2992	Lco2plus(i,j) = 0.d0
2993	Pco2plus(i,j) = 0.d0
2994	Loplus(i,j) = 0.d0
2995	Poplus(i,j) = 0.d0
2996	Lo2plus(i,j) = 0.d0
2997	Po2plus(i,j) = 0.d0
2998	Pcoplus(i,j) = 0.d0
2999	Lcoplus(i,j) = 0.d0
3000	Pcplus(i,j) = 0.d0
3001	Lcplus(i,j) = 0.d0
3002	Pnplus(i,j) = 0.d0
3003	Lnplus(i,j) = 0.d0
3004	Pnoplus(i,j) = 0.d0
3005	Lnoplus(i,j) = 0.d0
3006	Pn2plus(i,j) = 0.d0
3007	Ln2plus(i,j) = 0.d0
3008	Phplus(i,j) = 0.d0
3009	Lhplus(i,j) = 0.d0
3010	Phco2plus(i,j) = 0.d0
3011	Lhco2plus(i,j) = 0.d0
3012	Pelect(i,j) = 0.d0
3013	Lelect(i,j) = 0.d0
3014	end do
3015	Pco2tot(i) = 0.d0
3016	Lco2tot(i) = 0.d0
3017	Po2tot(i) = 0.d0
3018	Lo2tot(i) = 0.d0
3019	Po3ptot(i) = 0.d0
3020	Lo3ptot(i) = 0.d0
3021	Pcotot(i) = 0.d0
3022	Lcotot(i) = 0.d0
3023	Phtot(i) = 0.d0
3024	Lhtot(i) = 0.d0
3025	Pohtot(i) = 0.d0
3026	Lohtot(i) = 0.d0
3027	Pho2tot(i) = 0.d0
3028	Lho2tot(i) = 0.d0
3029	Ph2tot(i) = 0.d0
3030	Lh2tot(i) = 0.d0
3031	Ph2otot(i) = 0.d0
3032	Lh2otot(i) = 0.d0
3033	Po1dtot(i) = 0.d0
3034	Lo1dtot(i) = 0.d0
3035	Ph2o2tot(i) = 0.d0
3036	Lh2o2tot(i) = 0.d0
3037	Po3tot(i) = 0.d0
3038	Lo3tot(i) = 0.d0
3039	Pntot(i) = 0.d0
3040	Lntot(i) = 0.d0
3041	Pnotot(i) = 0.d0
3042	Lnotot(i) = 0.d0
3043	Pn2tot(i) = 0.d0
3044	Ln2tot(i) = 0.d0
3045	Pn2dtot(i) = 0.d0
3046	Ln2dtot(i) = 0.d0
3047	Pno2tot(i) = 0.d0
3048	Lno2tot(i) = 0.d0
3049	!
3050	Pco2plustot(i) = 0.d0
3051	Lco2plustot(i) = 0.d0
3052	Loplustot(i) = 0.d0
3053	Poplustot(i) = 0.d0
3054	Lo2plustot(i) = 0.d0
3055	Po2plustot(i) = 0.d0
3056	Lcoplustot(i) = 0.d0
3057	Pcoplustot(i) = 0.d0
3058	Lcplustot(i) = 0.d0
3059	Pcplustot(i) = 0.d0
3060	Lnplustot(i) = 0.d0
3061	Pnplustot(i) = 0.d0
3062	Lnoplustot(i) = 0.d0
3063	Pnoplustot(i) = 0.d0
3064	Ln2plustot(i) = 0.d0
3065	Pn2plustot(i) = 0.d0
3066	Lhplustot(i) = 0.d0
3067	Phplustot(i) = 0.d0
3068	Lhco2plustot(i) = 0.d0
3069	Phco2plustot(i) = 0.d0
3070	Pelecttot(i) = 0.0d0
3071	Lelecttot(i) = 0.0d0
3072	! endif
3073
3074
3075
3076	!!! Productions and losses reaction by reaction
3077
3078	c R1: CO2 + hv -> CO + O
3079
3080	Lco2(i,1) = jdistot8(1,i)
3081	Pco(i,1) = co2xinput * jdistot8(1,i)
3082	Po3p(i,1) = Pco(i,1)!co2xinput * jdistot8(1,i)
3083
3084
3085	c R16(1b): CO2 + hv -> CO + O1D
3086
3087	Lco2(i,16) = jdistot8_b(1,i)
3088	Pco(i,16) = co2xinput * jdistot8_b(1,i)
3089	Po1d(i,16) = Pco(i,16)!co2xinput * jdistot8_b(1,i)
3090
3091
3092	c R2: H + O2 + CO2 -> HO2 + CO2
3093
3094	Lh(i,2) = ch2 * o2xinput * co2xinput
3095	Lo2(i,2) = ch2 * hxinput * co2xinput
3096	Pho2(i,2) = ch2 * hxinput * o2xinput * co2xinput
3097
3098
3099	c R3: O + HO2 -> OH + O2
3100
3101	Lo3p(i,3) = ch3 * ho2xinput
3102	Lho2(i,3) = ch3 * o3pxinput
3103	Poh(i,3) = ch3 * o3pxinput * ho2xinput
3104	Po2(i,3) = Poh(i,3)!ch3 * o3pxinput * ho2xinput
3105
3106
3107	c R4: CO + OH -> CO2 + H
3108
3109	Lco(i,4) = ch4 * ohxinput
3110	Loh(i,4) = ch4 * coxinput
3111	Pco2(i,4) = ch4 * coxinput * ohxinput
3112	Ph(i,4) = Pco2(i,4)!ch4 * coxinput * ohxinput
3113
3114
3115	c R5: 2HO2 -> H2O2 + O2
3116
3117	Lho2(i,5) = 2.d0 * ch5 * ho2xinput
3118	Po2(i,5) = ch5 * ho2xinput * ho2xinput
3119	Ph2o2(i,5) = Po2(i,5)!ch5 * ho2xinput * ho2xinput
3120
3121
3122	c R6: H2O2 + hv -> 2OH
3123
3124	Lh2o2(i,6) = jdistot8(6,i)
3125	Poh(i,6) = 2.d0 * jdistot8(6,i) * h2o2xinput
3126
3127
3128	c R7: OH + HO2 -> H2O + O2
3129
3130	Loh(i,7) = ch7 * ho2xinput
3131	Lho2(i,7) = ch7 * ohxinput
3132	Po2(i,7) = ch7 * ohxinput * ho2xinput
3133	Ph2o(i,7) = Po2(i,7)!ch7 * ohxinput * ho2xinput
3134
3135
3136	c R8: H20 + hv -> H + OH
3137
3138	Lh2o(i,8) = jdistot8(4,i)
3139	Ph(i,8) = jdistot8(4,i) * h2oxinput
3140	Poh(i,8) = Ph(i,8)!jdistot8(4,i) * h2oxinput
3141
3142
3143	c R9: O(1D) + H2O -> 2OH
3144
3145	Lo1d(i,9) = ch9 * h2oxinput
3146	Lh2o(i,9) = ch9 * o1dxinput
3147	Poh(i,9) = 2.d0 * ch9 * o1dxinput * h2oxinput
3148
3149
3150	c R10: 2O + CO2 -> O2 + CO2
3151
3152	Lo3p(i,10) = 2.d0 * ch10 * o3pxinput * co2xinput
3153	Po2(i,10) = ch10 * o3pxinput * o3pxinput * co2xinput
3154
3155
3156	c R11: O + OH -> O2 + H
3157
3158	Lo3p(i,11) = ch11 * ohxinput
3159	Loh(i,11) = ch11 * o3pxinput
3160	Po2(i,11) = ch11 * o3pxinput * ohxinput
3161	Ph(i,11) = Po2(i,11)!ch11 * o3pxinput * ohxinput
3162
3163
3164	c R12: O2 + hv -> 2O
3165
3166	Lo2(i,12) = jdistot8(2,i)
3167	Po3p(i,12) = 2.d0 * jdistot8(2,i) * o2xinput
3168
3169
3170	c R17(12b): O2 + hv -> O + O1D
3171
3172	Lo2(i,17) = jdistot8_b(2,i)
3173	Po3p(i,17) = jdistot8_b(2,i) * o2xinput
3174	Po1d(i,17) = Po3p(i,17)!jdistot8_b(2,i) * o2xinput
3175
3176
3177	c R13: H + HO2 -> H2 + O2
3178
3179	Lh(i,13) = ch13 * ho2xinput
3180	Lho2(i,13) = ch13 * hxinput
3181	Ph2(i,13) = ch13 * hxinput * ho2xinput
3182	Po2(i,13) = Ph2(i,13)!ch13 * hxinput * ho2xinput
3183
3184
3185	c R14: O(1D) + H2 -> H + OH
3186
3187	Lo1d(i,14) = ch14 * h2xinput
3188	Lh2(i,14) = ch14 * o1dxinput
3189	Ph(i,14) = ch14 * o1dxinput * h2xinput
3190	Poh(i,14) = Ph(i,14)!ch14 * o1dxinput * h2xinput
3191
3192
3193	c R15: OH + H2 -> H + H20
3194
3195	Loh(i,15) = ch15 * h2xinput
3196	Lh2(i,15) = ch15 * ohxinput
3197	Ph(i,15) = ch15 * ohxinput * h2xinput
3198	Ph2o(i,15) = Ph(i,15)!ch15 * ohxinput * h2xinput
3199
3200
3201	c R18: OH + H2O2 -> H2O + HO2
3202
3203	Loh(i,18) = ch18 * h2o2xinput
3204	Lh2o2(i,18) = ch18 * ohxinput
3205	Ph2o(i,18) = ch18 * ohxinput * h2o2xinput
3206	Pho2(i,18) = Ph2o(i,18)!ch18 * ohxinput * h2o2xinput
3207
3208
3209	c R19: O(1D) + CO2 -> O + CO2
3210
3211	Lo1d(i,19) = ch19 * co2xinput
3212	Po3p(i,19) = ch19 * o1dxinput * co2xinput
3213
3214
3215	c R20: O(1D) + O2 -> O + O2
3216
3217	Lo1d(i,20) = ch20 * o2xinput
3218	Po3p(i,20) = ch20 * o1dxinput * o2xinput
3219
3220
3221	c R21: O + O2 + CO2 -> O3 + CO2
3222
3223	Lo3p(i,21) = ch21 * o2xinput * co2xinput
3224	Lo2(i,21) = ch21 * o3pxinput * co2xinput
3225	Po3(i,21) = ch21 * o3pxinput * o2xinput * co2xinput
3226
3227
3228	!Only if O3, N or ion chemistry requested
3229	if(chemthermod.ge.1) then
3230
3231	c R22: O3 + H -> OH + O2
3232
3233	Lo3(i,22) = ch22 * hxinput
3234	Lh(i,22) = ch22 * o3xinput
3235	Poh(i,22) = ch22 * o3xinput * hxinput
3236	Po2(i,22) = Poh(i,22) !ch22 * o3xinput * hxinput
3237
3238
3239	c R23: O3 + OH -> HO2 + O2
3240
3241	Lo3(i,23) = ch23 * ohxinput
3242	Loh(i,23) = ch23 * o3xinput
3243	Pho2(i,23) = ch23 * ohxinput * o3xinput
3244	Po2(i,23) = Pho2(i,23) !ch23 * ohxinput * o3xinput
3245
3246
3247	c R24: O3 + HO2 -> OH + 2O2
3248
3249	Lo3(i,24) = ch24 * ho2xinput
3250	Lho2(i,24) = ch24 * o3xinput
3251	Poh(i,24) = ch24 * o3xinput * ho2xinput
3252	Po2(i,24) = Poh(i,24) !2.d0 * ch24 * o3xinput * ho2xinput
3253
3254
3255	c R25: O3 + hv -> O2 + O3P
3256
3257	Lo3(i,25) = jdistot8(7,i)
3258	Po2(i,25) = jdistot8(7,i) * o3xinput
3259	Po3p(i,25) = Po2(i,25) !jdistot8(7,i) * o3xinput
3260
3261
3262	c R26 (R25_b): O3 + hv -> O2 + O1D
3263
3264	Lo3(i,26) = jdistot8_b(7,i)
3265	Po2(i,26) = jdistot8_b(7,i) * o3xinput
3266	Po1d(i,26) = Po2(i,26) !jdistot8_b(7,i) * o3xinput
3267
3268	endif !Of chemthermod.ge.1
3269
3270
3271	c R27: H2 + hv -> 2H
3272
3273	Lh2(i,27) = jdistot8(5,i)
3274	Ph(i,27) = 2.d0 * jdistot8(5,i) * h2xinput
3275
3276
3277	!Only if N or ion chemistry requested
3278	if(chemthermod.ge.2) then
3279
3280	c R28: N2 + hv -> N + N2D
3281
3282	Ln2(i,28) = jdistot8(8,i)
3283	Pn(i,28) = jdistot8(8,i) * n2xinput
3284	Pn2d(i,28) = Pn(i,28) !jdistot8(8,i) * n2xinput
3285
3286
3287	c R29: NO + hv -> N + O
3288
3289	Lno(i,29) = jdistot8(10,i)
3290	Pn(i,29) = jdistot8(10,i) * noxinput
3291	Po3p(i,29) = Pn(i,29) !jdistot8(10,i) * noxinput
3292
3293
3294	c R30: N + NO -> N2 + O
3295
3296	Ln(i,30) = ch30 * noxinput
3297	Lno(i,30) = ch30 * nxinput
3298	Pn2(i,30) = ch30 * nxinput * noxinput
3299	Po3p(i,30) = Pn2(i,30) !ch30 * nxinput * noxinput
3300
3301
3302	c R31: N2 + O1D -> N2 + O
3303
3304	Lo1d(i,31) = ch31 * n2xinput
3305	Po3p(i,31) = ch31 * n2xinput * o1dxinput
3306
3307
3308	c R32: N + O2 -> NO + O
3309
3310	Ln(i,32) = ch32 * o2xinput
3311	Lo2(i,32) = ch32 * nxinput
3312	Pno(i,32) = ch32 * o2xinput * nxinput
3313	Po3p(i,32) = Pno(i,32) !ch32 * o2xinput * nxinput
3314
3315
3316	c R33: N + OH -> NO + H
3317
3318	Ln(i,33) = ch33 * ohxinput
3319	Loh(i,33) = ch33 * nxinput
3320	Pno(i,33) = ch33 * nxinput * ohxinput
3321	Ph(i,33) = Pno(i,33) !Pch33 * nxinput * ohxinput
3322
3323
3324	c R34: N + O3 -> NO + O2
3325
3326	Ln(i,34) = ch34 * o3xinput
3327	Lo3(i,34) = ch34 * nxinput
3328	Pno(i,34) = ch34 * nxinput * o3xinput
3329	Po2(i,34) = Pno(i,34) !ch34 * nxinput * o3xinput
3330
3331
3332	c R35: N + HO2 -> NO + OH
3333
3334	Ln(i,35) = ch35 * ho2xinput
3335	Lho2(i,35) = ch35 * nxinput
3336	Pno(i,35) = ch35 * nxinput * ho2xinput
3337	Poh(i,35) = Pno(i,35) !ch35 * nxinput * ho2xinput
3338
3339
3340	c R36: N2D + O -> N + O
3341
3342	Ln2d(i,36) = ch36 * o3pxinput
3343	Pn(i,36) = ch36 * n2dxinput * o3pxinput
3344
3345
3346	c R37: N2D + N2 -> N + N2
3347
3348	Ln2d(i,37) = ch37 * n2xinput
3349	Pn(i,37) = ch37 * n2dxinput * n2xinput
3350
3351
3352	c R38: N2D + CO2 -> NO + CO
3353
3354	Ln2d(i,38) = ch38 * co2xinput
3355	Lco2(i,38) = ch38 * n2dxinput
3356	Pno(i,38) = ch38 * n2dxinput * co2xinput
3357	Pco(i,38) = Pno(i,38) !ch38 * n2dxinput * co2xinput
3358
3359
3360	c R39: NO + HO2 -> NO2 + OH
3361
3362	Lno(i,39) = ch39 * ho2xinput
3363	Lho2(i,39) = ch39 * noxinput
3364	Pno2(i,39) = ch39 * noxinput * ho2xinput
3365	Poh(i,39) = Pno2(i,39) !ch39 * noxinput * ho2xinput
3366
3367
3368	c R40: O + NO + CO2 -> NO2 + CO2
3369
3370	Lo3p(i,40) = ch40 * noxinput * co2xinput
3371	Lno(i,40) = ch40 * o3pxinput * co2xinput
3372	Pno2(i,40) = ch40 * o3pxinput * noxinput * co2xinput
3373
3374
3375	c R41: O + NO2 -> NO + O2
3376
3377	Lo3p(i,41) = ch41 * no2xinput
3378	Lno2(i,41) = ch41 * o3pxinput
3379	Pno(i,41) = ch41 * o3pxinput * no2xinput
3380	Po2(i,41) = Pno(i,41) !ch41 * o3pxinput * no2xinput
3381
3382
3383	c R42: NO + O3 -> NO2 + O2
3384
3385	Lno(i,42) = ch42 * o3xinput
3386	Lo3(i,42) = ch42 * noxinput
3387	Pno2(i,42) = ch42 * noxinput * o3xinput
3388	Po2(i,42) = Pno2(i,42) !ch42 * noxinput * o3xinput
3389
3390
3391	c R43: H + NO2 -> NO + OH
3392
3393	Lh(i,43) = ch43 * no2xinput
3394	Lno2(i,43) = ch43 * hxinput
3395	Pno(i,43) = ch43 * no2xinput * hxinput
3396	Poh(i,43) = Pno(i,43) !ch43 * no2xinput * hxinput
3397
3398
3399	c R44: NO2 + hv -> NO + O
3400
3401	Lno2(i,44) = jdistot8(13,i)
3402	Pno(i,44) = jdistot8(13,i) * no2xinput
3403	Po3p(i,44) = Pno(i,44) !jdistot8(13,i) * no2xinput
3404
3405
3406	c R45: N + O -> NO
3407
3408	Ln(i,45) = ch45 * o3pxinput
3409	Lo3p(i,45) = ch45 * nxinput
3410	Pno(i,45) = ch45 * o3pxinput * nxinput
3411
3412	endif !Of chemthermod.ge.2
3413
3414
3415	! >>>> IONOSFERA
3416
3417	!Only if ion chemistry requested
3418	if(chemthermod.eq.3) then
3419
3420	c R46: CO2+ + O2 -> CO2 + O2+
3421
3422	Lco2plus(i,46) = ch46 * o2xinput
3423	Lo2(i,46) = ch46 * co2plusxinput
3424	Pco2(i,46) = ch46 * co2plusxinput * o2xinput
3425	Po2plus(i,46) = Pco2(i,46) !ch46 * co2plusxinput * o2xinput
3426
3427
3428	c R47: CO2+ + O -> O+ + CO2
3429
3430	Lco2plus(i,47) = ch47 * o3pxinput
3431	Lo3p(i,47) = ch47 * co2plusxinput
3432	Pco2(i,47) = ch47 * co2plusxinput * o3pxinput
3433	Poplus(i,47) = Pco2(i,47) !ch47 * co2plusxinput * o3pxinput
3434
3435
3436	c R48: CO2+ O -> O2+ + CO
3437
3438	Lco2plus(i,48) = ch48 * o3pxinput
3439	Lo3p(i,48) = ch48 * co2plusxinput
3440	Po2plus(i,48) = ch48 * co2plusxinput * o3pxinput
3441	Pco(i,48) = Po2plus(i,48) !ch48 * co2plusxinput * o3pxinput
3442
3443
3444	c R49: O2+ + e -> O + O
3445
3446	Lo2plus(i,49) = ch49 * electxinput
3447	Lelect(i,49) = ch49 * o2plusxinput
3448	Po3p(i,49) = ch49 * o2plusxinput * electxinput *2.d0
3449
3450
3451	c R50: O+ + CO2 -> O2+ + CO
3452
3453	Loplus(i,50) = ch50 * co2xinput
3454	Lco2(i,50) = ch50 * oplusxinput
3455	Po2plus(i,50)= ch50 * oplusxinput * co2xinput
3456	Pco(i,50) = Po2plus(i,50) !ch50 * oplusxinput * co2xinput
3457
3458
3459	c R51: CO2 + hv -> CO2+ + e
3460
3461	Lco2(i,51) = jion8(1,i,1)
3462	Pco2plus(i,51) = co2xinput * jion8(1,i,1)
3463	Pelect(i,51) = Pco2plus(i,51) !co2xinput * jion8(1,i,1)
3464
3465
3466	c R52: CO2 + hv --> O+ + CO + e
3467
3468	Lco2(i,52) = jion8(1,i,2)
3469	Pco(i,52) = co2xinput * jion8(1,i,2)
3470	Poplus(i,52) = Pco(i,52) !co2xinput * jion8(1,i,2)
3471	Pelect(i,52) = Pco(i,52) !co2xinput * jion8(1,i,2)
3472
3473
3474	c R53: CO2 + hv --> CO+ + O + e
3475
3476	Lco2(i,53) = jion8(1,i,3)
3477	Pcoplus(i,53) = co2xinput * jion8(1,i,3)
3478	Po3p(i,53) = Pcoplus(i,53) !co2xinput * jion8(1,i,3)
3479	Pelect(i,53) = Pcoplus(i,53) !co2xinput * jion8(1,i,3)
3480
3481
3482	c R54: CO2 + hv --> C+ + O2 + e
3483
3484	Lco2(i,54) = jion8(1,i,4)
3485	Pcplus(i,54) = co2xinput * jion8(1,i,4)
3486	Po2(i,54) = Pcplus(i,54) !co2xinput * jion8(1,i,4)
3487	Pelect(i,54) = Pcplus(i,54) !co2xinput * jion8(1,i,4)
3488
3489
3490	c R55: CO2+ + e --> CO + O
3491
3492	Lco2plus(i,55) = ch55 * electxinput
3493	Lelect(i,55) = ch55 * co2plusxinput
3494	Pco(i,55) = ch55 * co2plusxinput * electxinput
3495	Po3p(i,55) = Pco(i,55) !ch55 * co2plusxinput * electxinput
3496
3497
3498	c R56: O+ + CO2 --> O2 + CO+ Probablemente no se dá
3499
3500	Lco2(i,56) = ch56 * oplusxinput
3501	Loplus(i,56) = ch56 * co2xinput
3502	Po2(i,56) = ch56 * co2xinput * oplusxinput
3503	Pcoplus(i,56)= Po2(i,56) !ch56 * co2xinput * oplusxinput
3504
3505
3506	c R57: CO+ + CO2 --> CO2+ + CO
3507
3508	Lco2(i,57) = ch57 * coplusxinput
3509	Lcoplus(i,57) = ch57 * co2xinput
3510	Pco2plus(i,57)= ch57 * coplusxinput * co2xinput
3511	Pco(i,57) = Pco2plus(i,57) !ch57 * coplusxinput * co2xinput
3512
3513
3514	c R58: CO+ + O --> O+ + CO
3515
3516	Lo3p(i,58) = ch58 * coplusxinput
3517	Lcoplus(i,58) = ch58 * o3pxinput
3518	Poplus(i,58) = ch58 * coplusxinput * o3pxinput
3519	Pco(i,58) = Poplus(i,58) !ch58 * coplusxinput * o3pxinput
3520
3521
3522	c R59: C+ + CO2 --> CO+ + CO
3523
3524	Lco2(i,59) = ch59 * cplusxinput
3525	Lcplus(i,59) = ch59 * co2xinput
3526	Pcoplus(i,59)= ch59 * cplusxinput * co2xinput
3527	Pco(i,59) = Pcoplus(i,59) !ch59 * cplusxinput * co2xinput
3528
3529
3530	c R60: O2 + hv --> O2+ + e
3531
3532	Lo2(i,60) = jion8(2,i,1)
3533	Po2plus(i,60) = o2xinput * jion8(2,i,1)
3534	Pelect(i,60) = Po2plus(i,60) !o2xinput * jion8(2,i,1)
3535
3536
3537	c R61: O + hv --> O+ + e
3538
3539	Lo3p(i,61) = jion8(3,i,1)
3540	Poplus(i,61) = o3pxinput * jion8(3,i,1)
3541	Pelect(i,61) = Poplus(i,61) !o3pxinput * jion8(3,i,1)
3542
3543
3544	c R62 : CO2+ + NO --> NO+ + CO2
3545
3546	Lco2plus(i,62) = ch62 * noxinput
3547	Lno(i,62) = ch62 * co2plusxinput
3548	Pnoplus(i,62) = ch62 * co2plusxinput * noxinput
3549	Pco2(i,62) = Pnoplus(i,62) !ch62 * co2plusxinput * noxinput
3550
3551
3552	c R63 : CO2+ + N --> NO + CO+
3553
3554	Lco2plus(i,63) = ch63 * nxinput
3555	Ln(i,63) = ch63 * co2plusxinput
3556	Pcoplus(i,63) = ch63 * co2plusxinput * nxinput
3557	Pno(i,63) = Pcoplus(i,63) !ch63 * co2plusxinput * nxinput
3558
3559
3560	c R64 : O2+ + NO --> NO+ + O2
3561
3562	Lo2plus(i,64) = ch64 * noxinput
3563	Lno(i,64) = ch64 * o2plusxinput
3564	Pnoplus(i,64) = ch64 * o2plusxinput * noxinput
3565	Po2(i,64) = Pnoplus(i,64) !ch64 * o2plusxinput * noxinput
3566
3567
3568	c R65 : O2+ + N2 --> NO+ + NO
3569
3570	Lo2plus(i,65) = ch65 * n2xinput
3571	Ln2(i,65) = ch65 * o2plusxinput
3572	Pnoplus(i,65) = ch65 * o2plusxinput * n2xinput
3573	Pno(i,65) = Pnoplus(i,65) !ch65 * o2plusxinput * n2xinput
3574
3575
3576	c R66: O2+ + N --> NO+ + O
3577
3578	Lo2plus(i,66) = ch66 * nxinput
3579	Ln(i,66) = ch66 * o2plusxinput
3580	Pnoplus(i,66) = ch66 * o2plusxinput * nxinput
3581	Po3p(i,66) = Pnoplus(i,66) !ch66 * o2plusxinput * nxinput
3582
3583
3584	c R67 : O+ + N2 --> NO+ + N
3585
3586	Loplus(i,67) = ch67 * n2xinput
3587	Ln2(i,67) = ch67 * oplusxinput
3588	Pnoplus(i,67) = ch67 * oplusxinput * n2xinput
3589	Pn(i,67) = Pnoplus(i,67) !ch67 * oplusxinput * n2xinput
3590
3591
3592	c R68 : N2+ + CO2 --> CO2+ + N2
3593
3594	Ln2plus(i,68) = ch68 * co2xinput
3595	Lco2(i,68) = ch68 * n2plusxinput
3596	Pco2plus(i,68) = ch68 * n2plusxinput * co2xinput
3597	Pn2(i,68) = Pco2plus(i,68) !ch68 * n2plusxinput * co2xinput
3598
3599
3600	c R69 : N2+ + O3p --> NO+ + N
3601
3602	Ln2plus(i,69) = ch69 * o3pxinput
3603	Lo3p(i,69) = ch69 * n2plusxinput
3604	Pnoplus(i,69) = ch69 * n2plusxinput * o3pxinput
3605	Pn(i,69) = Pnoplus(i,69) !ch69 * n2plusxinput * o3pxinput
3606
3607
3608	c R70 : N2+ + CO --> N2 + CO+
3609
3610	Ln2plus(i,70) = ch70 * coxinput
3611	Lco(i,70) = ch70 * n2plusxinput
3612	Pcoplus(i,70) = ch70 * n2plusxinput * coxinput
3613	Pn2(i,70) = Pcoplus(i,70) !ch70 * n2plusxinput * coxinput
3614
3615
3616	c R71 : N2+ + e- --> N + N
3617
3618	Ln2plus(i,71) = ch71 * electxinput
3619	Lelect(i,71) = ch71 * n2plusxinput
3620	Pn(i,71) = 2. * ch71 * n2plusxinput * electxinput
3621
3622
3623	c R72 : N2+ + O3p --> O+ + N2
3624
3625	Ln2plus(i,72) = ch72 * o3pxinput
3626	Lo3p(i,72) = ch72 * n2plusxinput
3627	Poplus(i,72) = ch72 * n2plusxinput * o3pxinput
3628	Pn2(i,72) = Poplus(i,72) !ch72 * n2plusxinput * o3pxinput
3629
3630
3631	c R73 : CO+ + H --> H+ + CO
3632
3633	Lcoplus(i,73) = ch73 * hxinput
3634	Lh(i,73) = ch73 * coplusxinput
3635	Phplus(i,73) = ch73 * coplusxinput * hxinput
3636	Pco(i,73) = Phplus(i,73) !ch73 * coplusxinput * hxinput
3637
3638	c R74 : O+ + H --> H+ + O
3639
3640	Loplus(i,74) = ch74 * hxinput
3641	Lh(i,74) = ch74 * oplusxinput
3642	Phplus(i,74) = ch74* oplusxinput * hxinput
3643	Po3p(i,74) = Phplus(i,74) !ch74 * oplusxinput * hxinput
3644
3645
3646	c R75: NO+ + e- --> N + O
3647
3648	Lnoplus(i,75) = ch75 * electxinput
3649	Lelect(i,75) = ch75 * noplusxinput
3650	Pn(i,75)= ch75 * noplusxinput * electxinput
3651	Po3p(i,75)= Pn(i,75) !ch75 * noplusxinput * electxinput
3652
3653
3654	c R76: H+ + O3p --> O+ + H
3655
3656	Lhplus(i,76) = ch76 * o3pxinput
3657	Lo3p(i,76) = ch76 * hplusxinput
3658	Poplus(i,76) = ch76 * hplusxinput * o3pxinput
3659	Ph(i,76) = Poplus(i,76) !ch76 * hplusxinput * o3pxinput
3660
3661
3662	c R77: CO + hv --> CO+ + e-
3663
3664	Lco(i,77) = jion8(11,i,1)
3665	Pcoplus(i,77) = coxinput * jion8(11,i,1)
3666	Pelect(i,77) = Pcoplus(i,77) !coxinput * jion8(11,i,1)
3667
3668
3669	c R78: CO + hv --> C+ + O + e-
3670
3671	Lco(i,78) = jion8(11,i,2)
3672	Pcplus(i,78) = coxinput * jion8(11,i,2)
3673	Po3p(i,78) = Pcplus(i,78) !coxinput * jion8(11,i,2)
3674	Pelect(i,78) = Pcplus(i,78) !coxinput * jion8(11,i,2)
3675
3676
3677	c R79: CO + hv --> C + O+ + e-
3678
3679	! Lco(i,79) = jion8(11,i,3)
3680	! Pc(i,79) = coxinput * jion8(11,i,3)
3681	! Poplus(i,79) = Pc(i,79)!Pc(i,79)coxinput * jion8(11,i,3)
3682	! Pelect(i,79) = Pc(i,79)!coxinput * jion8(11,i,3)
3683
3684
3685	c R80: NO + hv --> NO+ + e-
3686
3687	Lno(i,80) = jion8(10,i,1)
3688	Pnoplus(i,80) = noxinput * jion8(10,i,1)
3689	Pelect(i,80) = Pnoplus(i,80) !noxinput * jion8(10,i,1)
3690
3691
3692	c R81: N2 + hv --> N2+ + e-
3693
3694	Ln2(i,81) = jion8(8,i,1)
3695	Pn2plus(i,81) = n2xinput * jion8(8,i,1)
3696	Pelect(i,81) = Pn2plus(i,81) !n2xinput * jion8(8,i,1)
3697
3698
3699	c R82: N2 + hv --> N+ + N + e-
3700
3701	Ln2(i,82) = jion8(8,i,2)
3702	Pnplus(i,82) = n2xinput * jion8(8,i,2)
3703	Pn(i,82) = Pnplus(i,82) !n2xinput * jion8(8,i,2)
3704	Pelect(i,82) = Pnplus(i,82) !n2xinput * jion8(8,i,2)
3705
3706
3707	c R83: H + hv --> H+ + e
3708
3709	Lh(i,83) = jion8(12,i,1)
3710	Phplus(i,83) = hxinput * jion8(12,i,1)
3711	Pelect(i,83) = Phplus(i,83) !hxinput * jion8(12,i,1)
3712
3713
3714	c R84: N + hv --> N+ + e
3715
3716	Ln(i,84) = jion8(9,i,1)
3717	Pnplus(i,84) = nxinput * jion8(9,i,1)
3718	Pelect(i,84) = Pnplus(i,84) !nxinput * jion8(9,i,1)
3719
3720
3721	c R85: N+ + CO2 --> CO2+ + N
3722
3723	Pn(i,85) = ch85 * co2xinput * nplusxinput
3724	Pco2plus(i,85) = Pn(i,85) !ch85 * co2xinput * nplusxinput
3725	Lnplus(i,85) = ch85 * co2xinput
3726	Lco2(i,85) = ch85 * nplusxinput
3727
3728
3729	c R86: H2 + CO2+ --> H + HCO2+
3730	Ph(i,86) = ch86 * co2plusxinput * h2xinput
3731	Lco2plus(i,86) = ch86 * h2xinput
3732	Lh2(i,86) = ch86 * co2plusxinput
3733	!!!
3734
3735	c R87: HCO2+ + e --> H + CO2
3736	Ph(i,87) = ch87 * hco2plusxinput * electxinput
3737	Pco2(i,87) = Ph(i,87)
3738	Lhco2plus(i,87) = ch87 * electxinput
3739	Lelect(i,87) = ch87 * hco2plusxinput
3740
3741
3742	c R88: N + e --> N+ + e + e
3743
3744	Pnplus(i,88) = ionsec_nplus(zenit,zx(i))*Pnplus(i,84)
3745	Pelect(i,88) = Pnplus(i,88)
3746	if(nxinput.ge.1.d-20) then
3747	Ln(i,88) = Pnplus(i,88)/nxinput
3748	else
3749	Ln(i,88) = 1.d-30
3750	endif
3751
3752
3753	c R89: N2 + e --> N2+ + e + e
3754
3755	Pn2plus(i,89) = ionsec_n2plus(zenit,zx(i))*Pn2plus(i,81)
3756	Pelect(i,89) = Pn2plus(i,89)
3757	if(n2xinput.ge.1.d-20) then
3758	Ln2(i,89) = Pn2plus(i,89)/n2xinput
3759	else
3760	Ln2(i,89) = 1.d-30
3761	endif
3762
3763
3764	c R90: O + e --> O+ + e + e
3765
3766	Poplus(i,90) = ionsec_oplus(zenit,zx(i))*Poplus(i,61)
3767	Pelect(i,90) = Poplus(i,90)
3768	if(o3pxinput.ge.1.d-20) then
3769	Lo3p(i,90) = Poplus(i,90)/o3pxinput
3770	else
3771	Lo3p(i,90) = 1.d-30
3772	endif
3773
3774
3775	c R91: CO + e --> CO+ + e + e
3776
3777	Pcoplus(i,91) = ionsec_coplus(zenit,zx(i))*Pcoplus(i,77)
3778	Pelect(i,91) = Pcoplus(i,91)
3779	if(coxinput.ge.1.d-20) then
3780	Lco(i,91) = Pcoplus(i,91)/coxinput
3781	else
3782	Lco(i,91) = 1.d-30
3783	endif
3784
3785
3786	c R92: CO2 + e --> CO2+ + e + e
3787
3788	Pco2plus(i,92) = ionsec_co2plus(zenit,zx(i))*
3789	$ Pco2plus(i,51)
3790	Pelect(i,92) = Pco2plus(i,92)
3791	if(co2xinput.ge. 1.d-20) then
3792	Lco2(i,92) = Pco2plus(i,92)/co2xinput
3793	else
3794	Lco2(i,92) = 1.d-30
3795	endif
3796
3797
3798	c R93: O2 + e --> O2+ + e
3799	Po2plus(i,93) = ionsec_o2plus(zenit,zx(i))*Po2plus(i,60)
3800	Pelect(i,93) = Po2plus(i,93)
3801	if(o2xinput.ge.1.d-20) then
3802	Lo2(i,93) = Po2plus(i,93)/o2xinput
3803	else
3804	Lo2(i,93) = 1.d-30
3805	endif
3806
3807	endif !Of chemthermod.eq.3
3808
3809
3810
3811	c Total production and loss for each species. To save CPU time, we
3812	c do not sum over all reactions, but only over those where the species
3813	c participates
3814
3815	c CO2:
3816	c Prod: R4, R46, R47, R62, R87
3817	c Loss: R1, R16, R38, R50, R51, R52, R53, R54, R57, R59, R68, R85, R92
3818
3819	Pco2tot(i) = Pco2(i,4) + Pco2(i,46) + Pco2(i,47) + Pco2(i,62)
3820	$ + Pco2(i,87)
3821	Lco2tot(i) = Lco2(i,1) + Lco2(i,16) + Lco2(i,38) + Lco2(i,50) +
3822	$ Lco2(i,51) + Lco2(i,52) + Lco2(i,53) + Lco2(i,54) +
3823	$ Lco2(i,56) + Lco2(i,57) + Lco2(i,59) + Lco2(i,68) +
3824	$ Lco2(i,85) + Lco2(i,92)
3825
3826	c O2
3827	c Prod: R3, R5, R7, R10, R11, R13, R22, R23, R24, R25, R26, R34, R41, R42,
3828	c R54, R64
3829	c Loss: R2, R12, R17, R21, R32, R46, R60, R93
3830
3831	Po2tot(i) = Po2(i,3) + Po2(i,5) + Po2(i,7) + Po2(i,10) +
3832	$ Po2(i,11) + Po2(i,13) + Po2(i,22) + Po2(i,23) + Po2(i,24) +
3833	$ Po2(i,25) + Po2(i,26) + Po2(i,34) + Po2(i,41) + Po2(i,42) +
3834	$ Po2(i,54) + Po2(i,56) + Po2(i,64)
3835	Lo2tot(i) = Lo2(i,2) + Lo2(i,12) + Lo2(i,17) + Lo2(i,21) +
3836	$ Lo2(i,32) + Lo2(i,46) + Lo2(i,60) + Lo2(i,93)
3837
3838
3839	c O(3p)
3840	c Prod: R1, R12, R17, R19, R20, R25, R29, R30, R31, R32, R44, R49, R53,
3841	c R55, R66, R74, R75, R78
3842	c Loss: R3, R10, R11, R21, R40, R41, R45, R47, R48, R58, R61, R69, R72,
3843	c R76, R90
3844
3845	Po3ptot(i) = Po3p(i,1) + Po3p(i,12) + Po3p(i,17) + Po3p(i,19) +
3846	$ Po3p(i,20) + Po3p(i,25) + Po3p(i,29) + Po3p(i,30) +
3847	$ Po3p(i,31) + Po3p(i,32) + Po3p(i,44) + Po3p(i,49) +
3848	$ Po3p(i,53) + Po3p(i,55) + Po3p(i,66) + Po3p(i,74) +
3849	$ Po3p(i,75) + Po3p(i,78)
3850	Lo3ptot(i) = Lo3p(i,3) + Lo3p(i,10) + Lo3p(i,11) + Lo3p(i,21) +
3851	$ Lo3p(i,40) + Lo3p(i,41) + Lo3p(i,45) + Lo3p(i,47) +
3852	$ Lo3p(i,48) + Lo3p(i,58) + Lo3p(i,61) + Lo3p(i,69) +
3853	$ Lo3p(i,72) + Lo3p(i,76) + Lo3p(i,90)
3854
3855
3856	c CO
3857	c Prod: R1, R16, R38, R48, R50, R52, R55, R57, R58, R59, R73
3858	c Loss: R4, R70, R77, R78, R91
3859
3860	Pcotot(i) = Pco(i,1) + Pco(i,16) + Pco(i,38) + Pco(i,48) +
3861	$ Pco(i,50) + Pco(i,52) + Pco(i,55) + Pco(i,57) + Pco(i,58) +
3862	$ Pco(i,59) + Pco(i,73)
3863	Lcotot(i) = Lco(i,4) + Lco(i,70) + Lco(i,77) + Lco(i,78) +
3864	$ Lco(i,91)
3865
3866
3867	c H
3868	c Prod: R4, R8, R11, R14, R15, R27, R33, R76, R86, R87
3869	c Loss: R2, R13, R22, R43, R73, R74, R83
3870
3871	Phtot(i) = Ph(i,4) + Ph(i,8) + Ph(i,11) + Ph(i,14) + Ph(i,15) +
3872	$ Ph(i,27) + Ph(i,33) + Ph(i,76) + Ph(i,86) + Ph(i,87)
3873	Lhtot(i) = Lh(i,2) + Lh(i,13) + Lh(i,22) + Lh(i,43) + Lh(i,73) +
3874	$ Lh(i,74) + Lh(i,83)
3875
3876
3877	c OH
3878	c Prod: R3, R6, R8, R9, R14, R22, R24, R35, R39, R43
3879	c Loss: R4, R7, R11, R15, R18, R23, R33
3880
3881	Pohtot(i) = Poh(i,3) + Poh(i,6) + Poh(i,8) + Poh(i,9) +
3882	$ Poh(i,14) + Poh(i,22) + Poh(i,24) + Poh(i,35) + Poh(i,39) +
3883	$ Poh(i,43)
3884	Lohtot(i) = Loh(i,4) + Loh(i,7) + Loh(i,11) + Loh(i,15) +
3885	$ Loh(i,18) + Loh(i,23) + Loh(i,33)
3886
3887
3888	c HO2
3889	c Prod: R2, R18, R23
3890	c Loss: R3, R5, R7, R13, R24, R35, R39
3891
3892	Pho2tot(i) = Pho2(i,2) + Pho2(i,18) + Pho2(i,23)
3893	Lho2tot(i) = Lho2(i,3) + Lho2(i,5) + Lho2(i,7) + Lho2(i,13) +
3894	$ Lho2(i,24) + Lho2(i,35) + Lho2(i,39)
3895
3896
3897	c H2
3898	c Prod: R13
3899	c Loss: R14, R15, R27, R86
3900
3901	Ph2tot(i) = Ph2(i,13)
3902	Lh2tot(i) = Lh2(i,14) + Lh2(i,15) + Lh2(i,27) + Lh2(i,86)
3903
3904
3905	c H2O
3906	c Prod: R7, R15, R18
3907	c Loss: R8, R9
3908
3909	Ph2otot(i) = Ph2o(i,7) + Ph2o(i,15) + Ph2o(i,18)
3910	Lh2otot(i) = Lh2o(i,8) + Lh2o(i,9)
3911
3912
3913	c O(1d)
3914	c Prod: R16, R17, R26
3915	c Loss: R9, R14, R19, R20, R31
3916
3917	Po1dtot(i) = Po1d(i,16) + Po1d(i,17) + Po1d(i,26)
3918	Lo1dtot(i) = Lo1d(i,9) + Lo1d(i,14) + Lo1d(i,19) + Lo1d(i,20) +
3919	$ Lo1d(i,31)
3920
3921	c H2O2
3922	c Prod: R5
3923	c Loss: R6, R18
3924
3925	Ph2o2tot(i) = Ph2o2(i,5)
3926	Lh2o2tot(i) = Lh2o2(i,6) + Lh2o2(i,18)
3927
3928
3929	!Only if O3, N or ion chemistry requested
3930	if(chemthermod.ge.1) then
3931
3932	c O3
3933	c Prod: R21
3934	c Loss: R22, R23, R24, R25, R26, R34, R42
3935
3936	Po3tot(i) = Po3(i,21)
3937	Lo3tot(i) = Lo3(i,22) + Lo3(i,23) + Lo3(i,24) + Lo3(i,25) +
3938	$ Lo3(i,26) + Lo3(i,34) + Lo3(i,42)
3939
3940	endif
3941
3942
3943	!Only if N or ion chemistry requested
3944	if(chemthermod.ge.2) then
3945
3946	c N
3947	c Prod: R28, R29, R36, R37, R67, R69, R71, R75, R82, R85
3948	c Loss: R30, R32, R33, R34, R35, R45, R63, R66, R84, R88
3949
3950	Pntot(i) = Pn(i,28) + Pn(i,29) + Pn(i,36) + Pn(i,37) +Pn(i,67)+
3951	$ Pn(i,69) + Pn(i,71) + Pn(i,75) + Pn(i,82) + Pn(i,85)
3952	Lntot(i) = Ln(i,30) + Ln(i,32) + Ln(i,33) + Ln(i,34) +Ln(i,35)+
3953	$ Ln(i,45) + Ln(i,63) + Ln(i,66) + Ln(i,84) + Ln(i,88)
3954
3955
3956	c NO
3957	c Prod: R32, R33, R34, R35, R38, R41, R43, R44, R45, R63, R65
3958	c Loss: R29, R30, R39, R40, R42, R62, R64, R80
3959
3960	Pnotot(i) = Pno(i,32) + Pno(i,33) + Pno(i,34) + Pno(i,35) +
3961	$ Pno(i,38) + Pno(i,41) + Pno(i,43) + Pno(i,44) + Pno(i,45)+
3962	$ Pno(i,63) + Pno(i,65)
3963	Lnotot(i) = Lno(i,29) + Lno(i,30) + Lno(i,39) + Lno(i,40) +
3964	$ Lno(i,42) + Lno(i,62) + Lno(i,64) + Lno(i,80)
3965
3966
3967	c N2
3968	c Prod: R30, R68, R70, R72
3969	c Loss: R28, R65, R67, R81, R82, R89
3970
3971	Pn2tot(i) = Pn2(i,30) + Pn2(i,68) + Pn2(i,70) + Pn2(i,72)
3972	Ln2tot(i) = Ln2(i,28) + Ln2(i,65) + Ln2(i,67) + Ln2(i,81) +
3973	$ Ln2(i,82) + Ln2(i,89)
3974
3975
3976	c N(2d)
3977	c Prod: R28
3978	c Loss: R36, R37, R38
3979
3980	Pn2dtot(i) = Pn2d(i,28)
3981	Ln2dtot(i) = Ln2d(i,36) + Ln2d(i,37) + Ln2d(i,38)
3982
3983
3984	c NO2
3985	c Prod: R39, R40, R42
3986	c Loss: R41, R43, R44
3987
3988	Pno2tot(i) = Pno2(i,39) + Pno2(i,40) + Pno2(i,42)
3989	Lno2tot(i) = Lno2(i,41) + Lno2(i,43) + Lno2(i,44)
3990
3991	endif !Of chemthermod.ge.2
3992	!
3993
3994	!Only if ion chemistry requested
3995	if(chemthermod.eq.3) then
3996
3997	c CO2+
3998	c Prod: R51,R57, R68, R85, R92
3999	c Loss: R46, R47, R48, R55, R62, R63, R86
4000
4001	Pco2plustot(i) = Pco2plus(i,51) + Pco2plus(i,57) +
4002	$ Pco2plus(i,68) + Pco2plus(i,85) + Pco2plus(i,92)
4003	Lco2plustot(i) = Lco2plus(i,46) + Lco2plus(i,47) +
4004	$ Lco2plus(i,48) + Lco2plus(i,55) + Lco2plus(i,62) +
4005	$ Lco2plus(i,63) + Lco2plus(i,86)
4006
4007
4008	c O+
4009	c Prod: R47, R52, R58, R61, R72, R76, R90
4010	c Loss: 50,67,74
4011
4012	Poplustot(i) = Poplus(i,47) + Poplus(i,52) + Poplus(i,58) +
4013	$ Poplus(i,61) + Poplus(i,72) + Poplus(i,76) + Poplus(i,90)
4014	Loplustot(i) = Loplus(i,50) + Loplus(i,56) + Loplus(i,67) +
4015	$ Loplus(i,74)
4016
4017	c O2+
4018	c Prod: R46, R48, R50, R60, R93
4019	c Loss: R49, R64, R65, R66
4020
4021	Po2plustot(i) = Po2plus(i,46) + Po2plus(i,48) + Po2plus(i,50) +
4022	$ Po2plus(i,60) + Po2plus(i,93)
4023	Lo2plustot(i) = Lo2plus(i,49) + Lo2plus(i,64) + Lo2plus(i,65) +
4024	$ Lo2plus(i,66)
4025
4026
4027	c CO+
4028	c Prod: R53, R59, R63, R70, R77, R91
4029	c Loss: R57, R58, R73
4030
4031	Pcoplustot(i) = Pcoplus(i,53) + Pcoplus(i,56) + Pcoplus(i,59) +
4032	$ Pcoplus(i,63) + Pcoplus(i,70) + Pcoplus(i,77) +
4033	$ Pcoplus(i,91)
4034	Lcoplustot(i) = Lcoplus(i,57) + Lcoplus(i,58) + Lcoplus(i,73)
4035
4036
4037	c C+
4038	c Prod: R54, R78
4039	c Loss: R59
4040
4041	Pcplustot(i) = Pcplus(i,54) + Pcplus(i,78)
4042	Lcplustot(i) = Lcplus(i,59)
4043
4044
4045	c N+
4046	c Prod: R82, R84, R88
4047	c Loss: R85
4048
4049	Pnplustot(i) = Pnplus(i,82) + Pnplus(i,84) + Pnplus(i,88)
4050	Lnplustot(i) = Lnplus(i,85)
4051
4052
4053	c N2+
4054	c Prod: R81, R89
4055	c Loss: R68, R69, R70, R71, R72
4056
4057	Pn2plustot(i) = Pn2plus(i,81) + Pn2plus(i,89)
4058	Ln2plustot(i) = Ln2plus(i,68) + Ln2plus(i,69) + Ln2plus(i,70) +
4059	$ Ln2plus(i,71) + Ln2plus(i,72)
4060
4061
4062	c NO+
4063	c Prod: R62, R64, R65, R66, R67, R69, R80
4064	c Loss: R75
4065
4066	Pnoplustot(i) = Pnoplus(i,62) + Pnoplus(i,64) + Pnoplus(i,65) +
4067	$ Pnoplus(i,66) + Pnoplus(i,67) + Pnoplus(i,69) +
4068	$ Pnoplus(i,80)
4069	Lnoplustot(i) = Lnoplus(i,75)
4070
4071
4072	c H+
4073	c Prod: R73, R74, R83
4074	c Loss: R76
4075
4076	Phplustot(i) = Phplus(i,73) + Phplus(i,74) + Phplus(i,83)
4077	Lhplustot(i) = Lhplus(i,76)
4078
4079
4080	c HCO2+
4081	c Prod: R86
4082	c Loss: R87
4083
4084	Phco2plustot(i) = Phco2plus(i,86)
4085	Lhco2plustot(i) = Lhco2plus(i,87)
4086
4087	c e-
4088	c Prod: R51, R52, R53, R54, R60, R61, R77, R78, R80, R81, R82, R83, R84,
4089	c R88, R89, R90, R91, R92, R93
4090	c Loss: R49, R55, R71, R75, R87
4091
4092	Pelecttot(i) = Pelect(i,51) + Pelect(i,52) + Pelect(i,53) +
4093	$ Pelect(i,54) + Pelect(i,60) + Pelect(i,61) + Pelect(i,77)+
4094	$ Pelect(i,78) + Pelect(i,80) + Pelect(i,81) + Pelect(i,82)+
4095	$ Pelect(i,83) + Pelect(i,84) + Pelect(i,88) + Pelect(i,89)+
4096	$ Pelect(i,90) + Pelect(i,91) + Pelect(i,92) + Pelect(i,93)
4097	Lelecttot(i) = Lelect(i,49) + Lelect(i,55) + Lelect(i,71) +
4098	$ Lelect(i,75) + Lelect(i,87)
4099
4100	endif !Of chemthermod.eq.3
4101
4102	return
4103	c END
4104	end
4105
4106
4107
4108
4109	c**********************************************************************
4110	c**********************************************************************
4111
4112	subroutine EF_oscilacion
4113	& (ig,i,nlayer,paso,chemthermod,zenit, zx,
4114	& jdistot8, jdistot8_b,jion8,
4115	& tminaux,
4116	& co2xoutput, co2xinput,
4117	$ o2xoutput, o2xinput,
4118	$ o3pxoutput, o3pxinput,
4119	$ coxoutput, coxinput,
4120	$ h2xoutput, h2xinput,
4121	$ h2oxoutput, h2oxinput,
4122	$ h2o2xoutput, h2o2xinput,
4123	$ o3xoutput, o3xinput,
4124	$ nxoutput, nxinput,
4125	$ noxoutput, noxinput,
4126	$ n2xoutput, n2xinput,
4127	& o1dxoutput, o1dxinput,
4128	& ohxoutput, ohxinput,
4129	& ho2xoutput, ho2xinput,
4130	& hxoutput, hxinput,
4131	& n2dxoutput, n2dxinput,
4132	& no2xoutput, no2xinput,
4133	& co2plusxoutput, co2plusxinput,
4134	& o2plusxoutput, o2plusxinput,
4135	& coplusxoutput, coplusxinput,
4136	& oplusxoutput, oplusxinput,
4137	& cplusxoutput, cplusxinput,
4138	& noplusxoutput, noplusxinput,
4139	& n2plusxoutput, n2plusxinput,
4140	& hplusxoutput, hplusxinput,
4141	& nplusxoutput, nplusxinput,
4142	$ hco2plusxoutput,hco2plusxinput,
4143	& electxoutput, electxinput,
4144	& electxoutput_timemarching )
4145
4146
4147	c Calculates the concentrations of the fast species in PE. Includes a
4148	c procedure to avoid oscillations
4149
4150	c
4151
4152	c 2009 FGG Adaptation to GCM
4153	c jul 2008 MA Version en subrutina
4154	c nov 2007 MA Original. Evita oscilacion en EF
4155	c**********************************************************************
4156
4157	use iono_h
4158	use param_v4_h, only: nabs,
4159	. ch2, ch3, ch4, ch5, ch7,ch9,ch10,ch11,ch13,ch14,ch15,ch18,
4160	. ch19,ch20,ch21,ch22,ch23,ch24,ch30,ch31,ch32,ch33,ch34,
4161	. ch35,ch36,ch37,ch38,ch39,ch40,ch41,ch42,ch43,ch45,
4162	. ch46,ch47,ch48,ch49,ch50,ch55,ch56,ch57,ch58,ch59,ch62,
4163	. ch63,ch64,ch65,ch66,ch67,ch68,ch69,ch70,ch71,
4164	. ch72,ch73,ch74,ch75,ch76,ch85,ch86,ch87
4165
4166
4167	implicit none
4168
4169	c arguments
4170
4171	integer ig,nlayer
4172	integer i ! I. Layer
4173	integer paso ! I. paso temporal del timemarching, 1,numpasos
4174	integer chemthermod
4175	real*8 tminaux ! I.
4176	real zx(nlayer)
4177	real zenit
4178	real*8 jdistot8(nabs,nlayer) ! I.
4179	real*8 jdistot8_b(nabs,nlayer) ! I.
4180	real*8 jion8(nabs,nlayer,4)
4181
4182	real*8 co2xoutput,o2xoutput,o3pxoutput,coxoutput,h2xoutput
4183	real*8 h2o2xoutput,o3xoutput,h2oxoutput
4184	real*8 nxoutput,noxoutput,n2xoutput
4185	real*8 ho2xoutput, hxoutput, ohxoutput ! O.
4186	real*8 o1dxoutput, n2dxoutput, no2xoutput ! O.
4187
4188	real*8 co2xinput,o2xinput,o3pxinput,coxinput,h2xinput
4189	real*8 h2o2xinput,o3xinput,h2oxinput
4190	real*8 nxinput,noxinput,n2xinput
4191	real*8 ho2xinput, hxinput, ohxinput ! I.
4192	real*8 o1dxinput, n2dxinput, no2xinput ! I.
4193
4194	real*8 co2plusxoutput, coplusxoutput, oplusxoutput ! O.
4195	real*8 o2plusxoutput, cplusxoutput, noplusxoutput ! O.
4196	real*8 n2plusxoutput, hplusxoutput ! O.
4197	real*8 electxoutput, nplusxoutput, hco2plusxoutput ! O.
4198	real*8 co2plusxinput, coplusxinput, oplusxinput ! I.
4199	real*8 o2plusxinput, cplusxinput, noplusxinput ! I.
4200	real*8 n2plusxinput, hplusxinput ! I.
4201	real*8 electxinput, nplusxinput, hco2plusxinput ! I.
4202
4203	real*8 electxoutput_timemarching ! I.
4204
4205
4206
4207	c local variables
4208
4209	integer npasitos
4210	parameter (npasitos=20)
4211	real*8 electxoutput_neutr
4212
4213	real*8 ho2xoutput2, hxoutput2, ohxoutput2
4214	real*8 o1dxoutput2, n2dxoutput2, no2xoutput2
4215	real*8 co2plusxoutput2, coplusxoutput2, oplusxoutput2
4216	real*8 o2plusxoutput2, hplusxoutput2, nplusxoutput2
4217	real*8 cplusxoutput2, noplusxoutput2, n2plusxoutput2
4218	real*8 hco2plusxoutput2
4219	!real*8 electxoutput2
4220
4221	integer correc_oscil, ip
4222	real*8 avg_pares, avg_impar, dispersion
4223	real*8 dif_impar, dif_pares , dif_pares_impar
4224	real*8 ho2xpares(3), hxpares(3), ohxpares(3)
4225	real*8 o1dxpares(3), n2dxpares(3), no2xpares(3)
4226	real*8 co2plusxpares(3), coplusxpares(3), oplusxpares(3)
4227	real*8 o2plusxpares(3), hplusxpares(3), nplusxpares(3)
4228	real*8 cplusxpares(3), noplusxpares(3), n2plusxpares(3)
4229	real*8 hco2plusxpares(3)
4230	real*8 ho2ximpar(3), hximpar(3), ohximpar(3)
4231	real*8 o1dximpar(3), n2dximpar(3), no2ximpar(3)
4232	real*8 co2plusximpar(3), coplusximpar(3), oplusximpar(3)
4233	real*8 o2plusximpar(3), hplusximpar(3), nplusximpar(3)
4234	real*8 cplusximpar(3), noplusximpar(3), n2plusximpar(3)
4235	real*8 hco2plusximpar(3)
4236	real*8 auxp,auxl,ca, cb,cc, cd1,cd2,ce
4237	real*8 IonMostAbundant
4238
4239	real*8 cocimin
4240	parameter (cocimin=1.d-30)
4241
4242	integer cociopt
4243	parameter (cociopt=1)
4244
4245	c external functions
4246	c
4247	external cociente
4248	real*8 cociente
4249
4250	external ionsec_nplus
4251	real*8 ionsec_nplus
4252
4253	external ionsec_n2plus
4254	real*8 ionsec_n2plus
4255
4256	external ionsec_oplus
4257	real*8 ionsec_oplus
4258
4259	external ionsec_coplus
4260	real*8 ionsec_coplus
4261
4262	external ionsec_co2plus
4263	real*8 ionsec_co2plus
4264
4265	external ionsec_o2plus
4266	real*8 ionsec_o2plus
4267
4268	external avg
4269	real*8 avg
4270
4271	external dif
4272	real*8 dif
4273
4274	real*8 log1
4275	real*8 log2
4276	real*8 log3
4277
4278	ccccccccccccccc CODE STARTS
4279
4280	!!! Starts iteration to avoid oscilations
4281
4282	correc_oscil=0
4283
4284	do ip = 1, npasitos
4285
4286	if (ip.eq.1) then
4287	if (o1d_eq(i).eq.'Y') o1dxoutput = o1dxinput
4288	if (oh_eq(i).eq.'Y') ohxoutput = ohxinput
4289	if (ho2_eq(i).eq.'Y') ho2xoutput = ho2xinput
4290	if (h_eq(i).eq.'Y') hxoutput = hxinput
4291	!Only if N or ion chemistry requested
4292	if(chemthermod.ge.2) then
4293	if (n2d_eq(i).eq.'Y') n2dxoutput = n2dxinput
4294	if (no2_eq(i).eq.'Y') no2xoutput = no2xinput
4295	endif
4296	!
4297	!Only if ion chemistry requested
4298	if(chemthermod.ge.3) then
4299	if (n2plus_eq(i).eq.'Y') n2plusxoutput=n2plusxinput
4300	if (cplus_eq(i).eq.'Y') cplusxoutput=cplusxinput
4301	if (coplus_eq(i).eq.'Y') coplusxoutput=coplusxinput
4302	if (oplus_eq(i).eq.'Y') oplusxoutput=oplusxinput
4303	if (hplus_eq(i).eq.'Y') hplusxoutput=hplusxinput
4304	if (co2plus_eq(i).eq.'Y') co2plusxoutput=co2plusxinput
4305	if (noplus_eq(i).eq.'Y') noplusxoutput=noplusxinput
4306	if (o2plus_eq(i).eq.'Y') o2plusxoutput=o2plusxinput
4307	if (nplus_eq(i).eq.'Y') nplusxoutput=nplusxinput
4308	if (hco2plus_eq(i).eq.'Y') hco2plusxoutput=hco2plusxinput
4309	electxoutput = electxinput
4310	endif
4311	endif
4312
4313
4314	! 6 neutral , O1D, OH, HO2, H, N2D, NO2
4315
4316	!O1D
4317	if (o1d_eq(i).eq.'Y') then
4318	auxp = jdistot8_b(1,i) * co2xoutput
4319	& + jdistot8_b(2,i) * o2xoutput
4320	& + jdistot8_b(7,i) * o3xoutput
4321	auxl = ch9 * h2oxoutput
4322	& + ch14 * h2xoutput
4323	& + ch19 * co2xoutput
4324	& + ch20 * o2xoutput
4325	& + ch31 * n2xoutput
4326	o1dxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4327	end if
4328
4329	!OH
4330	if (oh_eq(i).eq.'Y') then
4331	auxp = ch3 * o3pxoutput * ho2xoutput
4332	& + 2.d0 * jdistot8(6,i) * h2o2xoutput
4333	& + jdistot8(4,i) * h2oxoutput
4334	& + 2.d0 * ch9 * o1dxoutput * h2oxoutput
4335	& + ch14 * o1dxoutput * h2xoutput
4336	& + ch22 * o3xoutput * hxoutput
4337	& + ch24 * o3xoutput * ho2xoutput
4338	& + ch35 * nxoutput * ho2xoutput
4339	& + ch39 * noxoutput * ho2xoutput
4340	& + ch43 * no2xoutput * hxoutput
4341	auxl = ch4 * coxoutput
4342	& + ch7 * ho2xoutput
4343	& + ch11 * o3pxoutput
4344	& + ch15 * h2xoutput
4345	& + ch18 * h2o2xoutput
4346	& + ch23 * o3xoutput
4347	& + ch33 * nxoutput
4348	ohxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4349
4350	end if
4351
4352
4353	!HO2
4354	if (ho2_eq(i).eq.'Y') then
4355	auxp = ch2 * hxoutput * o2xoutput * co2xoutput
4356	& + ch18 * ohxoutput * h2o2xoutput
4357	& + ch23 * ohxoutput * o3xoutput
4358	auxl = ch3 * o3pxoutput
4359	& + 2.d0 * ch5 * ho2xoutput
4360	& + ch7 * ohxoutput
4361	& + ch13 * hxoutput
4362	& + ch24 * o3xoutput
4363	& + ch35 * nxoutput
4364	& + ch39 * noxoutput
4365	ho2xoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4366	end if
4367
4368	!H
4369	if (h_eq(i).eq.'Y') then
4370	auxp = ch4 * coxoutput * ohxoutput
4371	& + jdistot8(4,i) * h2oxoutput
4372	& + ch11 * o3pxoutput * ohxoutput
4373	& + ch14 * o1dxoutput * h2xoutput
4374	& + ch15 * ohxoutput * h2xoutput
4375	& + 2.d0 * jdistot8(5,i) * h2xoutput
4376	& + ch33 * nxoutput * ohxoutput
4377	& + ch76 * hplusxoutput * o3pxoutput
4378	& + hxoutput * jion8(12,i,1)
4379	$ + ch86 * h2xoutput * co2plusxoutput
4380	$ + ch87 * hco2plusxoutput * electxoutput
4381	auxl = ch2 * o2xoutput * co2xoutput
4382	& + ch13 * ho2xoutput
4383	& + ch22 * o3xoutput
4384	& + ch43 * no2xoutput
4385	& + ch73 * coplusxoutput
4386	& + ch74 * oplusxoutput
4387	& + jion8(12,i,1)
4388	hxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4389	end if
4390
4391	!Only if N or ion chemistry requested
4392	if(chemthermod.ge.2) then
4393	!N2D
4394	if (n2d_eq(i).eq.'Y') then
4395	auxp = jdistot8(8,i) * n2xoutput
4396	auxl = ch36 * o3pxoutput
4397	& + ch37 * n2xoutput
4398	& + ch38 * co2xoutput
4399	n2dxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4400	end if
4401
4402	!NO2
4403	if (no2_eq(i).eq.'Y') then
4404	auxp = ch39 * noxoutput * ho2xoutput
4405	& + ch40 * o3pxoutput * noxoutput * co2xoutput
4406	& + ch42 * noxoutput * o3xoutput
4407	auxl = ch41 * o3pxoutput
4408	& + ch43 * hxoutput
4409	& + jdistot8(13,i)
4410	no2xoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4411	end if
4412
4413	endif !Of chemthermod.ge.2
4414
4415	! 9 ions
4416
4417	!Only if ion chemistry requested
4418	if(chemthermod.eq.3) then
4419	! N2+
4420	if (n2plus_eq(i).eq.'Y') then
4421	auxp = jion8(8,i,1)n2xoutput
4422	$ (1.d0+ionsec_n2plus(zenit,zx(i)))
4423	auxl = ch68*co2xoutput
4424	& + ch69*o3pxoutput
4425	& + ch70*coxoutput
4426	$ + ch71*electxoutput
4427	& + ch72*o3pxoutput
4428	n2plusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4429	endif
4430
4431	! C+
4432	if (cplus_eq(i).eq.'Y') then
4433	auxp = jion8(1,i,4)*co2xoutput
4434	& + jion8(11,i,2)*coxoutput
4435	auxl = ch59*co2xoutput
4436	cplusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4437	end if
4438
4439	! CO+
4440	if (coplus_eq(i).eq.'Y') then
4441	auxp = jion8(1,i,3)*co2xoutput
4442	$ + ch59cplusxoutput co2xoutput
4443	$ + ch63co2plusxoutputnxoutput
4444	$ + ch70n2plusxoutputcoxoutput
4445	$ + jion8(11,i,1)coxoutput
4446	$ (1.d0+ionsec_coplus(zenit,zx(i)))
4447	auxl = ch57*co2xoutput
4448	& + ch58*o3pxoutput
4449	& + ch73*hxoutput
4450	coplusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4451	end if
4452
4453	! CO2+
4454	if (co2plus_eq(i).eq.'Y') then
4455	auxp = jion8(1,i,1)co2xoutput
4456	$ (1.d0+ionsec_co2plus(zenit,zx(i)))
4457	& + ch57coplusxoutputco2xoutput
4458	$ + ch68n2plusxoutputco2xoutput
4459	& + ch85nplusxoutputco2xoutput
4460	auxl = ch46*o2xoutput
4461	& + ch47*o3pxoutput
4462	& + ch55*electxoutput
4463	$ + ch62*noxoutput
4464	& + ch63*nxoutput
4465	& + ch48*o3pxoutput
4466	$ + ch86*h2xoutput
4467	co2plusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4468	end if
4469
4470	!
4471	! [O+] [H+] are linked: 2 equations with 2 unknowns
4472	!
4473	! [H+] = (ca + cb [O+])/cd1
4474	! [O+] = (cc + ce [H+])/cd2
4475	!
4476	ca = ch73 * coplusxoutput * hxoutput
4477	& + jion8(12,i,1)*hxoutput
4478	cb = ch74 * hxoutput
4479	cd1 = ch76 * o3pxoutput
4480	cc = ch47co2plusxoutputo3pxoutput
4481	& + jion8(1,i,2)*co2xoutput
4482	$ + jion8(3,i,1)o3pxoutput
4483	$ (1.d0+ionsec_oplus(zenit,zx(i)))
4484	& + ch58coplusxoutputo3pxoutput
4485	$ + ch72n2plusxoutputo3pxoutput
4486	ce = ch76 * o3pxoutput
4487	cd2 = ch50co2xoutput + ch67n2xoutput + ch74*hxoutput
4488
4489	! O+
4490	if (oplus_eq(i).eq.'Y') then
4491	auxp = cccd1 + ceca
4492	auxl = cd1cd2 -cecb
4493	oplusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4494	end if
4495
4496	! H+
4497	if (hplus_eq(i).eq.'Y') then
4498	auxp = ca + cb * oplusxoutput
4499	auxl = cd1
4500	hplusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4501	endif
4502
4503	! O2+
4504	if (o2plus_eq(i).eq.'Y') then
4505	auxp = ch46co2plusxoutputo2xoutput
4506	$ + ch48co2plusxoutputo3pxoutput
4507	$ + ch50oplusxoutputco2xoutput
4508	$ + jion8(2,i,1)o2xoutput
4509	$ (1.d0+ionsec_o2plus(zenit,zx(i)))
4510	auxl = ch49electxoutput + ch64noxoutput
4511	$ + ch65n2xoutput + ch66nxoutput
4512	o2plusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4513	endif
4514
4515	! NO+
4516	if(noplus_eq(i).eq.'Y') then
4517	auxp = ch62coplusxoutputnoxoutput
4518	$ + ch64o2plusxoutputnoxoutput
4519	$ + ch65o2plusxoutputn2xoutput
4520	$ + ch66o2plusxoutputnxoutput
4521	$ + ch67oplusxoutputn2xoutput
4522	$ + ch69n2plusxoutputo3pxoutput
4523	$ + jion8(10,i,1)*noxoutput
4524	auxl = ch75*electxoutput
4525	noplusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4526	endif
4527
4528	! N+
4529	if(nplus_eq(i).eq.'Y') then
4530	auxp = jion8(8,i,2)*n2xoutput
4531	& + jion8(9,i,1)nxoutput
4532	$ (1.d0+ionsec_nplus(zenit,zx(i)))
4533	auxl = ch85*co2xoutput
4534	nplusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4535	endif
4536
4537	! HCO2+
4538	if(hco2plus_eq(i).eq.'Y') then
4539	auxp = ch86h2xoutputco2plusxoutput
4540	auxl = ch87*electxoutput
4541	hco2plusxoutput2 = cociente(auxp,auxl,cocimin,cociopt)
4542	endif
4543
4544	endif !Of chemthermod.eq.3
4545
4546	! Detection of oscilations and elimination
4547
4548	if (ip.eq.4 .or. ip.eq.14) then ! ***pares(1)
4549
4550	if (o1d_eq(i).eq.'Y') o1dxpares(1)=o1dxoutput2
4551	if (oh_eq(i).eq.'Y') ohxpares(1)=ohxoutput2
4552	if (ho2_eq(i).eq.'Y') ho2xpares(1)=ho2xoutput2
4553	if (h_eq(i).eq.'Y') hxpares(1)=hxoutput2
4554	!Only if N or ion chemistry requested
4555	if(chemthermod.ge.2) then
4556	if (n2d_eq(i).eq.'Y') n2dxpares(1)=n2dxoutput2
4557	if (no2_eq(i).eq.'Y') no2xpares(1)=no2xoutput2
4558	endif
4559	!
4560	!Only if ion chemistry requested
4561	if(chemthermod.eq.3) then
4562	if (n2plus_eq(i).eq.'Y') n2plusxpares(1)=n2plusxoutput2
4563	if (cplus_eq(i).eq.'Y') cplusxpares(1)=cplusxoutput2
4564	if (coplus_eq(i).eq.'Y') coplusxpares(1)=coplusxoutput2
4565	if (oplus_eq(i).eq.'Y') oplusxpares(1)=oplusxoutput2
4566	if (hplus_eq(i).eq.'Y') hplusxpares(1)=hplusxoutput2
4567	if (co2plus_eq(i).eq.'Y')co2plusxpares(1)=co2plusxoutput2
4568	if (noplus_eq(i).eq.'Y') noplusxpares(1)=noplusxoutput2
4569	if (o2plus_eq(i).eq.'Y') o2plusxpares(1)=o2plusxoutput2
4570	if (nplus_eq(i).eq.'Y') nplusxpares(1)=nplusxoutput2
4571	if (hco2plus_eq(i).eq.'Y')
4572	$ hco2plusxpares(1)=hco2plusxoutput2
4573	endif
4574
4575	elseif (ip.eq.6 .or. ip.eq.16) then ! ***pares(2)
4576
4577	if (o1d_eq(i).eq.'Y') o1dxpares(2)=o1dxoutput2
4578	if (oh_eq(i).eq.'Y') ohxpares(2)=ohxoutput2
4579	if (ho2_eq(i).eq.'Y') ho2xpares(2)=ho2xoutput2
4580	if (h_eq(i).eq.'Y') hxpares(2)=hxoutput2
4581	!Only if N or ion chemistry requested
4582	if(chemthermod.ge.2) then
4583	if (n2d_eq(i).eq.'Y') n2dxpares(2)=n2dxoutput2
4584	if (no2_eq(i).eq.'Y') no2xpares(2)=no2xoutput2
4585	endif
4586	!
4587	!Only if ion chemistry requested
4588	if(chemthermod.eq.3) then
4589	if (n2plus_eq(i).eq.'Y') n2plusxpares(2)=n2plusxoutput2
4590	if (cplus_eq(i).eq.'Y') cplusxpares(2)=cplusxoutput2
4591	if (coplus_eq(i).eq.'Y') coplusxpares(2)=coplusxoutput2
4592	if (oplus_eq(i).eq.'Y') oplusxpares(2)=oplusxoutput2
4593	if (hplus_eq(i).eq.'Y') hplusxpares(2)=hplusxoutput2
4594	if (co2plus_eq(i).eq.'Y')co2plusxpares(2)=co2plusxoutput2
4595	if (noplus_eq(i).eq.'Y') noplusxpares(2)=noplusxoutput2
4596	if (o2plus_eq(i).eq.'Y') o2plusxpares(2)=o2plusxoutput2
4597	if (nplus_eq(i).eq.'Y') nplusxpares(2)=nplusxoutput2
4598	if (hco2plus_eq(i).eq.'Y')
4599	$ hco2plusxpares(2)=hco2plusxoutput2
4600	endif
4601
4602	elseif (ip.eq.8 .or. ip.eq.18) then ! ***pares(3)
4603
4604	if (o1d_eq(i).eq.'Y') o1dxpares(3)=o1dxoutput2
4605	if (oh_eq(i).eq.'Y') ohxpares(3)=ohxoutput2
4606	if (ho2_eq(i).eq.'Y') ho2xpares(3)=ho2xoutput2
4607	if (h_eq(i).eq.'Y') hxpares(3)=hxoutput2
4608	!Only if N or ion chemistry requested
4609	if(chemthermod.ge.2) then
4610	if (n2d_eq(i).eq.'Y') n2dxpares(3)=n2dxoutput2
4611	if (no2_eq(i).eq.'Y') no2xpares(3)=no2xoutput2
4612	endif
4613	!
4614	!Only if ion chemistry requested
4615	if(chemthermod.eq.3) then
4616	if (n2plus_eq(i).eq.'Y') n2plusxpares(3)=n2plusxoutput2
4617	if (cplus_eq(i).eq.'Y') cplusxpares(3)=cplusxoutput2
4618	if (coplus_eq(i).eq.'Y') coplusxpares(3)=coplusxoutput2
4619	if (oplus_eq(i).eq.'Y') oplusxpares(3)=oplusxoutput2
4620	if (hplus_eq(i).eq.'Y') hplusxpares(3)=hplusxoutput2
4621	if (co2plus_eq(i).eq.'Y')co2plusxpares(3)=co2plusxoutput2
4622	if (noplus_eq(i).eq.'Y') noplusxpares(3)=noplusxoutput2
4623	if (o2plus_eq(i).eq.'Y') o2plusxpares(3)=o2plusxoutput2
4624	if (nplus_eq(i).eq.'Y') nplusxpares(3)=nplusxoutput2
4625	if (hco2plus_eq(i).eq.'Y')
4626	$ hco2plusxpares(3)=hco2plusxoutput2
4627	endif
4628
4629	elseif (ip.eq.5 .or. ip.eq.15) then ! ***impar(1)
4630
4631	if (o1d_eq(i).eq.'Y') o1dximpar(1)=o1dxoutput2
4632	if (oh_eq(i).eq.'Y') ohximpar(1)=ohxoutput2
4633	if (ho2_eq(i).eq.'Y') ho2ximpar(1)=ho2xoutput2
4634	if (h_eq(i).eq.'Y') hximpar(1)=hxoutput2
4635	!Only if N or ion chemistry requested
4636	if(chemthermod.ge.2) then
4637	if (n2d_eq(i).eq.'Y') n2dximpar(1)=n2dxoutput2
4638	if (no2_eq(i).eq.'Y') no2ximpar(1)=no2xoutput2
4639	endif
4640	!
4641	!Only if ion chemistry requested
4642	if(chemthermod.eq.3) then
4643	if (n2plus_eq(i).eq.'Y') n2plusximpar(1)=n2plusxoutput2
4644	if (cplus_eq(i).eq.'Y') cplusximpar(1)=cplusxoutput2
4645	if (coplus_eq(i).eq.'Y') coplusximpar(1)=coplusxoutput2
4646	if (oplus_eq(i).eq.'Y') oplusximpar(1)=oplusxoutput2
4647	if (hplus_eq(i).eq.'Y') hplusximpar(1)=hplusxoutput2
4648	if (co2plus_eq(i).eq.'Y')co2plusximpar(1)=co2plusxoutput2
4649	if (noplus_eq(i).eq.'Y') noplusximpar(1)=noplusxoutput2
4650	if (o2plus_eq(i).eq.'Y') o2plusximpar(1)=o2plusxoutput2
4651	if (nplus_eq(i).eq.'Y') nplusximpar(1)=nplusxoutput2
4652	if (hco2plus_eq(i).eq.'Y')
4653	$ hco2plusximpar(1)=hco2plusxoutput2
4654	endif
4655
4656	elseif (ip.eq.7 .or. ip.eq.17) then ! ***impar(2)
4657
4658	if (o1d_eq(i).eq.'Y') o1dximpar(2)=o1dxoutput2
4659	if (oh_eq(i).eq.'Y') ohximpar(2)=ohxoutput2
4660	if (ho2_eq(i).eq.'Y') ho2ximpar(2)=ho2xoutput2
4661	if (h_eq(i).eq.'Y') hximpar(2)=hxoutput2
4662	!Only if N or ion chemistry requested
4663	if(chemthermod.ge.2) then
4664	if (n2d_eq(i).eq.'Y') n2dximpar(2)=n2dxoutput2
4665	if (no2_eq(i).eq.'Y') no2ximpar(2)=no2xoutput2
4666	endif
4667	!
4668	!Only if ion chemistry requested
4669	if(chemthermod.eq.3) then
4670	if (n2plus_eq(i).eq.'Y') n2plusximpar(2)=n2plusxoutput2
4671	if (cplus_eq(i).eq.'Y') cplusximpar(2)=cplusxoutput2
4672	if (coplus_eq(i).eq.'Y') coplusximpar(2)=coplusxoutput2
4673	if (oplus_eq(i).eq.'Y') oplusximpar(2)=oplusxoutput2
4674	if (hplus_eq(i).eq.'Y') hplusximpar(2)=hplusxoutput2
4675	if (co2plus_eq(i).eq.'Y')co2plusximpar(2)=co2plusxoutput2
4676	if (noplus_eq(i).eq.'Y') noplusximpar(2)=noplusxoutput2
4677	if (o2plus_eq(i).eq.'Y') o2plusximpar(2)=o2plusxoutput2
4678	if (nplus_eq(i).eq.'Y') nplusximpar(2)=nplusxoutput2
4679	if (hco2plus_eq(i).eq.'Y')
4680	$ hco2plusximpar(2)=hco2plusxoutput2
4681	endif
4682
4683	elseif (ip.eq.9 .or. ip.eq.19) then ! ***impar(3)
4684
4685	if (o1d_eq(i).eq.'Y') o1dximpar(3)=o1dxoutput2
4686	if (oh_eq(i).eq.'Y') ohximpar(3)=ohxoutput2
4687	if (ho2_eq(i).eq.'Y') ho2ximpar(3)=ho2xoutput2
4688	if (h_eq(i).eq.'Y') hximpar(3)=hxoutput2
4689	!Only if N or ion chemistry requested
4690	if(chemthermod.ge.2) then
4691	if (n2d_eq(i).eq.'Y') n2dximpar(3)=n2dxoutput2
4692	if (no2_eq(i).eq.'Y') no2ximpar(3)=no2xoutput2
4693	endif
4694	!
4695	!Only if ion chemistry requested
4696	if(chemthermod.eq.3) then
4697	if (n2plus_eq(i).eq.'Y') n2plusximpar(3)=n2plusxoutput2
4698	if (cplus_eq(i).eq.'Y') cplusximpar(3)=cplusxoutput2
4699	if (coplus_eq(i).eq.'Y') coplusximpar(3)=coplusxoutput2
4700	if (oplus_eq(i).eq.'Y') oplusximpar(3)=oplusxoutput2
4701	if (hplus_eq(i).eq.'Y') hplusximpar(3)=hplusxoutput2
4702	if (co2plus_eq(i).eq.'Y')co2plusximpar(3)=co2plusxoutput2
4703	if (noplus_eq(i).eq.'Y') noplusximpar(3)=noplusxoutput2
4704	if (o2plus_eq(i).eq.'Y') o2plusximpar(3)=o2plusxoutput2
4705	if (nplus_eq(i).eq.'Y') nplusximpar(3)=nplusxoutput2
4706	if (hco2plus_eq(i).eq.'Y')
4707	$ hco2plusximpar(3)=hco2plusxoutput2
4708	endif
4709
4710
4711	if (o1d_eq(i).eq.'Y') then
4712	log1 = log10(o1dxpares(1))
4713	log2 = log10(o1dxpares(2))
4714	log3 = log10(o1dxpares(3))
4715	avg_pares = avg(log1,log2,log3)
4716	dif_pares = dif(log1,log2,log3,avg_pares)
4717	log1 = log10(o1dximpar(1))
4718	log2 = log10(o1dximpar(2))
4719	log3 = log10(o1dximpar(3))
4720	avg_impar = avg(log1,log2,log3)
4721	dif_impar = dif(log1,log2,log3,avg_impar)
4722	dispersion = dif_pares + dif_impar
4723	dif_pares_impar = abs(avg_pares-avg_impar)
4724	if (dif_pares_impar .gt. 5.d0*dispersion) then
4725	correc_oscil=correc_oscil+1
4726	o1dxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4727	endif
4728	endif
4729
4730	if (oh_eq(i).eq.'Y') then
4731	log1 = log10(ohxpares(1))
4732	log2 = log10(ohxpares(2))
4733	log3 = log10(ohxpares(3))
4734	avg_pares = avg(log1,log2,log3)
4735	dif_pares = dif(log1,log2,log3,avg_pares)
4736	log1 = log10(ohximpar(1))
4737	log2 = log10(ohximpar(2))
4738	log3 = log10(ohximpar(3))
4739	avg_impar = avg(log1,log2,log3)
4740	dif_impar = dif(log1,log2,log3,avg_impar)
4741	dispersion = dif_pares + dif_impar
4742	dif_pares_impar = abs(avg_pares-avg_impar)
4743	if (dif_pares_impar .gt. 5.*dispersion) then
4744	correc_oscil=correc_oscil+1
4745	ohxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4746	endif
4747
4748	endif
4749
4750	if (ho2_eq(i).eq.'Y') then
4751	log1 = log10(ho2xpares(1))
4752	log2 = log10(ho2xpares(2))
4753	log3 = log10(ho2xpares(3))
4754	avg_pares = avg(log1,log2,log3)
4755	dif_pares = dif(log1,log2,log3,avg_pares)
4756	log1 = log10(ho2ximpar(1))
4757	log2 = log10(ho2ximpar(2))
4758	log3 = log10(ho2ximpar(3))
4759	avg_impar = avg(log1,log2,log3)
4760	dif_impar = dif(log1,log2,log3,avg_impar)
4761	dispersion = dif_pares + dif_impar
4762	dif_pares_impar = abs(avg_pares-avg_impar)
4763	if (dif_pares_impar .gt. 5.*dispersion) then
4764	correc_oscil=correc_oscil+1
4765	ho2xoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4766	endif
4767	endif
4768
4769	if (h_eq(i).eq.'Y') then
4770	log1 = log10(hxpares(1))
4771	log2 = log10(hxpares(2))
4772	log3 = log10(hxpares(3))
4773	avg_pares = avg(log1,log2,log3)
4774	dif_pares = dif(log1,log2,log3,avg_pares)
4775	log1 = log10(hximpar(1))
4776	log2 = log10(hximpar(2))
4777	log3 = log10(hximpar(3))
4778	avg_impar = avg(log1,log2,log3)
4779	dif_impar = dif(log1,log2,log3,avg_impar)
4780	dispersion = dif_pares + dif_impar
4781	dif_pares_impar = abs(avg_pares-avg_impar)
4782	if (dif_pares_impar .gt. 5.*dispersion) then
4783	correc_oscil=correc_oscil+1
4784	hxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4785	endif
4786	endif
4787
4788	!Only if N or ion chemistry requested
4789	if(chemthermod.ge.2) then
4790	if (n2d_eq(i).eq.'Y') then
4791	log1 = log10(n2dxpares(1))
4792	log2 = log10(n2dxpares(2))
4793	log3 = log10(n2dxpares(3))
4794	avg_pares = avg(log1,log2,log3)
4795	dif_pares = dif(log1,log2,log3,avg_pares)
4796	log1 = log10(n2dximpar(1))
4797	log2 = log10(n2dximpar(2))
4798	log3 = log10(n2dximpar(3))
4799	avg_impar = avg(log1,log2,log3)
4800	dif_impar = dif(log1,log2,log3,avg_impar)
4801	dispersion = dif_pares + dif_impar
4802	dif_pares_impar = abs(avg_pares-avg_impar)
4803	if (dif_pares_impar .gt. 5.*dispersion) then
4804	correc_oscil=correc_oscil+1
4805	n2dxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4806	endif
4807	endif
4808
4809	if (no2_eq(i).eq.'Y') then
4810	log1 = log10(no2xpares(1))
4811	log2 = log10(no2xpares(2))
4812	log3 = log10(no2xpares(3))
4813	avg_pares = avg(log1,log2,log3)
4814	dif_pares = dif(log1,log2,log3,avg_pares)
4815	log1 = log10(no2ximpar(1))
4816	log2 = log10(no2ximpar(2))
4817	log3 = log10(no2ximpar(3))
4818	avg_impar = avg(log1,log2,log3)
4819	dif_impar = dif(log1,log2,log3,avg_impar)
4820	dispersion = dif_pares + dif_impar
4821	dif_pares_impar = abs(avg_pares-avg_impar)
4822	if (dif_pares_impar .gt. 5.*dispersion) then
4823	correc_oscil=correc_oscil+1
4824	no2xoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4825	endif
4826	endif
4827
4828	endif !Of chemthermod.ge.2
4829
4830	! IONS
4831	!Only if ion chemistry requested
4832	if(chemthermod.eq.3) then
4833	if (cplus_eq(i).eq.'Y') then
4834	log1 = log10(cplusxpares(1))
4835	log2 = log10(cplusxpares(2))
4836	log3 = log10(cplusxpares(3))
4837	avg_pares = avg(log1,log2,log3)
4838	dif_pares = dif(log1,log2,log3,avg_pares)
4839	log1 = log10(cplusximpar(1))
4840	log2 = log10(cplusximpar(2))
4841	log3 = log10(cplusximpar(3))
4842	avg_impar = avg(log1,log2,log3)
4843	dif_impar = dif(log1,log2,log3,avg_impar)
4844	dispersion = dif_pares + dif_impar
4845	dif_pares_impar = abs(avg_pares-avg_impar)
4846	if (dif_pares_impar .gt. 5.*dispersion) then
4847	correc_oscil=correc_oscil+1
4848	cplusxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4849	endif
4850	endif
4851
4852	if (coplus_eq(i).eq.'Y') then
4853	log1 = log10(coplusxpares(1))
4854	log2 = log10(coplusxpares(2))
4855	log3 = log10(coplusxpares(3))
4856	avg_pares = avg(log1,log2,log3)
4857	dif_pares = dif(log1,log2,log3,avg_pares)
4858	log1 = log10(coplusximpar(1))
4859	log2 = log10(coplusximpar(2))
4860	log3 = log10(coplusximpar(3))
4861	avg_impar = avg(log1,log2,log3)
4862	dif_impar = dif(log1,log2,log3,avg_impar)
4863	dispersion = dif_pares + dif_impar
4864	dif_pares_impar = abs(avg_pares-avg_impar)
4865	if (dif_pares_impar .gt. 5.*dispersion) then
4866	correc_oscil=correc_oscil+1
4867	coplusxoutput2=10.d0*(0.5(avg_pares+avg_impar))
4868	endif
4869	endif
4870
4871	if (oplus_eq(i).eq.'Y') then
4872	log1 = log10(oplusxpares(1))
4873	log2 = log10(oplusxpares(2))
4874	log3 = log10(oplusxpares(3))
4875	avg_pares = avg(log1,log2,log3)
4876	dif_pares = dif(log1,log2,log3,avg_pares)
4877	log1 = log10(oplusximpar(1))
4878	log2 = log10(oplusximpar(2))
4879	log3 = log10(oplusximpar(3))
4880	avg_impar = avg(log1,log2,log3)
4881	dif_impar = dif(log1,log2,log3,avg_impar)
4882	dispersion = dif_pares + dif_impar
4883	dif_pares_impar = abs(avg_pares-avg_impar)
4884	if (dif_pares_impar .gt. 5.*dispersion) then
4885	correc_oscil=correc_oscil+1
4886	oplusxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4887	endif
4888	endif
4889
4890	if (n2plus_eq(i).eq.'Y') then
4891	log1 = log10(n2plusxpares(1))
4892	log2 = log10(n2plusxpares(2))
4893	log3 = log10(n2plusxpares(3))
4894	avg_pares = avg(log1,log2,log3)
4895	dif_pares = dif(log1,log2,log3,avg_pares)
4896	log1 = log10(n2plusximpar(1))
4897	log2 = log10(n2plusximpar(2))
4898	log3 = log10(n2plusximpar(3))
4899	avg_impar = avg(log1,log2,log3)
4900	dif_impar = dif(log1,log2,log3,avg_impar)
4901	dispersion = dif_pares + dif_impar
4902	dif_pares_impar = abs(avg_pares-avg_impar)
4903	if (dif_pares_impar .gt. 5.*dispersion) then
4904	correc_oscil=correc_oscil+1
4905	n2plusxoutput2 =10.d0*(0.5(avg_pares+avg_impar))
4906	endif
4907	endif
4908
4909	if (hplus_eq(i).eq.'Y') then
4910	log1 = log10(hplusxpares(1))
4911	log2 = log10(hplusxpares(2))
4912	log3 = log10(hplusxpares(3))
4913	avg_pares = avg(log1,log2,log3)
4914	dif_pares = dif(log1,log2,log3,avg_pares)
4915	log1 = log10(hplusximpar(1))
4916	log2 = log10(hplusximpar(2))
4917	log3 = log10(hplusximpar(3))
4918	avg_impar = avg(log1,log2,log3)
4919	dif_impar = dif(log1,log2,log3,avg_impar)
4920	dispersion = dif_pares + dif_impar
4921	dif_pares_impar = abs(avg_pares-avg_impar)
4922	if (dif_pares_impar .gt. 5.*dispersion) then
4923	correc_oscil=correc_oscil+1
4924	hplusxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
4925	endif
4926	endif
4927
4928	if (co2plus_eq(i).eq.'Y') then
4929	log1 = log10(co2plusxpares(1))
4930	log2 = log10(co2plusxpares(2))
4931	log3 = log10(co2plusxpares(3))
4932	avg_pares = avg(log1,log2,log3)
4933	dif_pares = dif(log1,log2,log3,avg_pares)
4934	log1 = log10(co2plusximpar(1))
4935	log2 = log10(co2plusximpar(2))
4936	log3 = log10(co2plusximpar(3))
4937	avg_impar = avg(log1,log2,log3)
4938	dif_impar = dif(log1,log2,log3,avg_impar)
4939	dispersion = dif_pares + dif_impar
4940	dif_pares_impar = abs(avg_pares-avg_impar)
4941	if (dif_pares_impar .gt. 5.*dispersion) then
4942	correc_oscil=correc_oscil+1
4943	co2plusxoutput2=10.d0*(0.5(avg_pares+avg_impar))
4944	endif
4945	endif
4946
4947	if (o2plus_eq(i).eq.'Y') then
4948	log1 = log10(o2plusxpares(1))
4949	log2 = log10(o2plusxpares(2))
4950	log3 = log10(o2plusxpares(3))
4951	avg_pares = avg(log1,log2,log3)
4952	dif_pares = dif(log1,log2,log3,avg_pares)
4953	log1 = log10(o2plusximpar(1))
4954	log2 = log10(o2plusximpar(2))
4955	log3 = log10(o2plusximpar(3))
4956	avg_impar = avg(log1,log2,log3)
4957	dif_impar = dif(log1,log2,log3,avg_impar)
4958	dispersion = dif_pares + dif_impar
4959	dif_pares_impar = abs(avg_pares-avg_impar)
4960	if (dif_pares_impar .gt. 5.*dispersion) then
4961	correc_oscil=correc_oscil+1
4962	o2plusxoutput2 =10.d0*(0.5(avg_pares+avg_impar))
4963	endif
4964	endif
4965
4966	if (noplus_eq(i).eq.'Y') then
4967	log1 = log10(noplusxpares(1))
4968	log2 = log10(noplusxpares(2))
4969	log3 = log10(noplusxpares(3))
4970	avg_pares = avg(log1,log2,log3)
4971	dif_pares = dif(log1,log2,log3,avg_pares)
4972	log1 = log10(noplusximpar(1))
4973	log2 = log10(noplusximpar(2))
4974	log3 = log10(noplusximpar(3))
4975	avg_impar = avg(log1,log2,log3)
4976	dif_impar = dif(log1,log2,log3,avg_impar)
4977	dispersion = dif_pares + dif_impar
4978	dif_pares_impar = abs(avg_pares-avg_impar)
4979	if (dif_pares_impar .gt. 5.*dispersion) then
4980	correc_oscil=correc_oscil+1
4981	noplusxoutput2 =10.d0*(0.5(avg_pares+avg_impar))
4982	endif
4983	endif
4984
4985	if (nplus_eq(i).eq.'Y') then
4986	log1 = log10(nplusxpares(1))
4987	log2 = log10(nplusxpares(2))
4988	log3 = log10(nplusxpares(3))
4989	avg_pares = avg(log1,log2,log3)
4990	dif_pares = dif(log1,log2,log3,avg_pares)
4991	log1 = log10(nplusximpar(1))
4992	log2 = log10(nplusximpar(2))
4993	log3 = log10(nplusximpar(3))
4994	avg_impar = avg(log1,log2,log3)
4995	dif_impar = dif(log1,log2,log3,avg_impar)
4996	dispersion = dif_pares + dif_impar
4997	dif_pares_impar = abs(avg_pares-avg_impar)
4998	if (dif_pares_impar .gt. 5.*dispersion) then
4999	correc_oscil=correc_oscil+1
5000	nplusxoutput2 = 10.d0*(0.5(avg_pares+avg_impar))
5001	endif
5002	endif
5003
5004	if (hco2plus_eq(i).eq.'Y') then
5005	log1 = log10(hco2plusxpares(1))
5006	log2 = log10(hco2plusxpares(2))
5007	log3 = log10(hco2plusxpares(3))
5008	avg_pares = avg(log1,log2,log3)
5009	dif_pares = dif(log1,log2,log3,avg_pares)
5010	log1 = log10(hco2plusximpar(1))
5011	log2 = log10(hco2plusximpar(2))
5012	log3 = log10(hco2plusximpar(3))
5013	avg_impar = avg(log1,log2,log3)
5014	dif_impar = dif(log1,log2,log3,avg_impar)
5015	dispersion = dif_pares + dif_impar
5016	dif_pares_impar = abs(avg_pares-avg_impar)
5017	if (dif_pares_impar .gt. 5.*dispersion) then
5018	correc_oscil=correc_oscil+1
5019	hco2plusxoutput2=10.d0*(0.5(avg_pares+avg_impar))
5020	endif
5021	endif
5022
5023	endif !Of chemthermod.eq.3
5024
5025	endif
5026
5027
5028	! Preparation of next step
5029
5030	if (o1d_eq(i).eq.'Y') o1dxoutput = o1dxoutput2
5031	if (oh_eq(i).eq.'Y') ohxoutput = ohxoutput2
5032	if (ho2_eq(i).eq.'Y') ho2xoutput = ho2xoutput2
5033	if (h_eq(i).eq.'Y') hxoutput = hxoutput2
5034	!Only if N or ion chemistry requested
5035	if(chemthermod.ge.2) then
5036	if (n2d_eq(i).eq.'Y') n2dxoutput = n2dxoutput2
5037	if (no2_eq(i).eq.'Y') no2xoutput = no2xoutput2
5038	endif
5039	!
5040	!Only if ion chemistry requested
5041	if(chemthermod.eq.3) then
5042	if (n2plus_eq(i).eq.'Y') n2plusxoutput=n2plusxoutput2
5043	if (cplus_eq(i).eq.'Y') cplusxoutput=cplusxoutput2
5044	if (coplus_eq(i).eq.'Y') coplusxoutput=coplusxoutput2
5045	if (oplus_eq(i).eq.'Y') oplusxoutput=oplusxoutput2
5046	if (hplus_eq(i).eq.'Y') hplusxoutput=hplusxoutput2
5047	if (co2plus_eq(i).eq.'Y') co2plusxoutput=co2plusxoutput2
5048	if (noplus_eq(i).eq.'Y') noplusxoutput=noplusxoutput2
5049	if (o2plus_eq(i).eq.'Y') o2plusxoutput=o2plusxoutput2
5050	if (nplus_eq(i).eq.'Y') nplusxoutput=nplusxoutput2
5051	if (hco2plus_eq(i).eq.'Y') hco2plusxoutput=hco2plusxoutput2
5052	electxoutput = o2plusxoutput +
5053	@ co2plusxoutput +
5054	@ coplusxoutput +
5055	@ oplusxoutput +
5056	@ cplusxoutput +
5057	@ n2plusxoutput +
5058	@ nplusxoutput +
5059	@ noplusxoutput +
5060	@ hplusxoutput +
5061	$ hco2plusxoutput
5062
5063	electxoutput_neutr = electxoutput
5064
5065	IonMostAbundant = o2plusxoutput
5066	IonMostAbundant = max( co2plusxoutput, IonMostAbundant)
5067	IonMostAbundant = max( coplusxoutput, IonMostAbundant)
5068	IonMostAbundant = max( oplusxoutput, IonMostAbundant)
5069	IonMostAbundant = max( cplusxoutput, IonMostAbundant)
5070	IonMostAbundant = max( n2plusxoutput, IonMostAbundant)
5071	IonMostAbundant = max( noplusxoutput, IonMostAbundant)
5072	IonMostAbundant = max( nplusxoutput, IonMostAbundant)
5073	IonMostAbundant = max( hplusxoutput, IonMostAbundant)
5074	IonMostAbundant = max( hco2plusxoutput, IonMostAbundant)
5075	IonMostAbundant = IonMostAbundant / electxoutput
5076
5077	endif !Of chemthermod.eq.3
5078
5079
5080
5081	enddo
5082	!!! End of iteration
5083
5084
5085	return
5086	end
5087
5088	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5089	function avg(log1,log2,log3)
5090	implicit none
5091	real*8 avg
5092	real*8 log1,log2,log3
5093	avg = (log1+log2+log3)*0.333
5094	return
5095	end
5096	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5097	function dif(log1,log2,log3,avg)
5098	implicit none
5099	real*8 dif
5100	real*8 avg
5101	real*8 log1,log2,log3
5102	dif = (abs(log1-avg) +
5103	& abs(log2-avg) +
5104	& abs(log3-avg) ) * 0.333
5105	return
5106	end
5107	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5108
5109	c**********************************************************************
5110	c***********************************************************************
5111	function cociente (xnum, xdenom, minimo, option)
5112
5113	c Returns the ratio between XNUM and XDENOM avoiding null values
5114	c according to the action given by OPTION. Checks that the are not
5115	c negative values
5116
5117	c If XNUM = 0 -> cociente=0
5118	c If XDENOM < minimo, increases XDENOM in a factor=1e20 to avoid
5119	c overflow, and after the ratio XNUM/XDENOM the value is back to normal
5120	c if cociente < minimo -> cociente=minimo
5121	c If XDENOM = 0 then :
5122	c option = 0 .... warning message and stop
5123	c option = 1 .... " " and cociente=minimo
5124
5125	c MALV Jul-08 Original
5126	c***********************************************************************
5127
5128	implicit none
5129
5130	! Arguments
5131	real*8 cociente
5132	real*8 xnum
5133	real*8 xdenom
5134	real*8 minimo
5135	integer option
5136
5137	! Local variables
5138	real*8 factor
5139
5140	!!!!!!! Program starts
5141
5142	if (xnum .lt. 0.d0) then
5143	write (,) log10( xnum )
5144	STOP ' ERROR. Negative productions. XNUM=0.'
5145	elseif (xdenom .lt. 0.d0) then
5146	STOP ' ERROR. Negative losses. XDENOM=0.'
5147	endif
5148
5149	if (xnum .eq. 0.d0) then
5150	cociente = minimo
5151	return
5152	endif
5153
5154	if (xdenom .eq. 0.d0) then
5155	if (option .eq. 0) then
5156	STOP ' ERROR. xdenom=0. '
5157	elseif (option .eq. 1) then
5158	! write (,) 'WARNING !! xdenom=0 '
5159	! write (,) 'WARNING !! option=2 => cociente=minimo',
5160	! $ xdenom
5161	cociente = minimo
5162	return
5163	else
5164	STOP ' ERROR. option undefined in call to cociente'
5165	endif
5166	endif
5167
5168	if (xdenom .lt. minimo) then
5169	factor = xdenom * 1.d20
5170	cociente = xnum / factor * 1.d20
5171	else
5172	cociente = xnum / xdenom
5173	endif
5174
5175	if (cociente .lt. minimo) cociente = minimo
5176
5177	return
5178	c END
5179	end

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

Download in other formats:

Original Format