Context Navigation

paramfoto_compact.F @ 1198

Last change on this file since 1198 was 1013, checked in by emillour, 11 years ago

Mars GCM:

Upgrade on the thermospheric photochemical reaction rates. These are now read in from a file "chemthermos_reactionrates.def".

FGG + JYC

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

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

Download in other formats:

Original Format