Context Navigation

paramfoto_compact.F @ 757

Last change on this file since 757 was 690, checked in by acolaitis, 13 years ago
Code re-organization in diverse parts of the GCM code. These are NOT cosmetic changes, but are needed for compilation of the Mesoscale model in NESTED configuration
File size: 166.7 KB

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

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

Download in other formats:

Original Format