Context Navigation

paramfoto_compact.F @ 3567

Last change on this file since 3567 was 3464, checked in by emillour, 2 months ago

Mars PCM:
Some tidying in aeronomars:

make a jthermcalc_util.F to contain utility routines (used by jthermcal & jthermcalc_e107). Also add the possibility (turned off by default) in the interfast routine to do extra sanity checks.
turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107, paramphoto_compact and photochemistry into modules.

File size: 174.8 KB

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

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

Download in other formats:

Original Format