Context Navigation

paramfoto_compact.F @ 2599

Last change on this file since 2599 was 2151, checked in by aslmd, 5 years ago
changed old functions dexp dlog in generic exp and log
File size: 175.3 KB

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

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

Download in other formats:

Original Format