Context Navigation

paramfoto_compact.F @ 1448

Last change on this file since 1448 was 1331, checked in by emillour, 10 years ago

Mars GCM:

Bug fix: Enforce recomputation of reaction rates at every call of routine prodsandlosses (in paramfoto_compact.F)

FGG

File size: 176.2 KB

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

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

Download in other formats:

Original Format