Context Navigation

paramfoto_compact.F @ 3567

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

Mars PCM:
Some tidying in aeronomars:

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

File size: 174.8 KB

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

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

Download in other formats:

Original Format