Context Navigation

paramfoto_compact.F @ 2030

Last change on this file since 2030 was 1888, checked in by emillour, 7 years ago

Mars GCM:
Updated the calculation of the dissociation and ionization
branching ratios, using the values in the Schunk and Nagy book.
New datafiles *branchingratio_schunkandnagy2000_param.dat must be loaded
for the "EUVDAT" subdirectory of the standard "datadir" directory.
Main changes are:

param_v4_h.F90 -> New definition for the O2 ionization branching ratio
param_read_e107.F -> Read the new files containing the S&N branching ratios
paramfoto_compact.F -> Mainly cleaning of the code and the comments.

Also correction of a bug affecting the calculation of CO losses

chemthermos.F90 -> Small modification to add the possibility of including

NO and O2 nightglow rates to the outputs

calchim.F90 and calchim_asis.F90 -> account for change in arguments in

calls to chemthermos
FGG

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

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

Download in other formats:

Original Format