Context Navigation

paramfoto_compact.F @ 3026

Last change on this file since 3026 was 2615, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars
File size: 175.4 KB

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

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

Download in other formats:

Original Format