Context Navigation

photochemistry.F90 @ 2883

Last change on this file since 2883 was 2833, checked in by flefevre, 2 years ago
Kinetics data: JPL 2019 update
File size: 137.5 KB

Line
1	!****************************************************************
2	!
3	! Photochemical routine
4	!
5	! Authors: Franck Lefevre, Francisco Gonzalez-Galindo
6	! -------
7	!
8	! Version: 14/11/2020
9	!
10	! ASIS scheme : for details on the method see
11	! Cariolle et al., Geosci. Model Dev., 10, 1467-1485, 2017.
12	!
13	!*****************************************************************
14
15	subroutine photochemistry(nlayer, nq, nesp, ionchem, deutchem, &
16	nb_reaction_3_max, nb_reaction_4_max, nphot, &
17	nb_phot_max, nphotion, &
18	jonline, ig, lswitch, zycol, sza, ptimestep, press, &
19	alt, temp, temp_elect, dens, zmmean, &
20	dist_sol, zday, &
21	surfdust1d, surfice1d, jo3, jh2o,em_no,em_o2, &
22	tau, iter)
23
24	use param_v4_h, only: jion
25
26	implicit none
27
28	include "callkeys.h"
29
30	!===================================================================
31	! inputs:
32	!===================================================================
33
34	integer, intent(in) :: nlayer ! number of atmospheric layers
35	integer, intent(in) :: nq ! number of tracers in traceur.def
36	integer, intent(in) :: nesp ! number of traceurs in chemistry
37	integer, intent(in) :: nb_reaction_3_max
38	! number of quadratic reactions
39	integer, intent(in) :: nb_reaction_4_max
40	! number of bimolecular reactions
41	integer, intent(in) :: nphot
42	! number of photodissociations
43	integer, intent(in) :: nb_phot_max
44	! number of reactions treated numerically as photodissociations
45	integer, intent(in) :: nphotion
46	! number of photoionizations
47	logical, intent(in) :: ionchem! switch for ion chemistry
48	logical, intent(in) :: deutchem
49	! switch for deuterium chemistry
50	logical, intent(in) :: jonline! switch for on-line calculation of photolysis rates
51	integer :: ig ! grid point index
52
53	real :: sza ! solar zenith angle (deg)
54	real :: ptimestep ! physics timestep (s)
55	real :: press(nlayer) ! pressure (hpa)
56	real :: alt(nlayer) ! altitude (km)
57	real :: temp(nlayer) ! temperature (k)
58	real :: temp_elect(nlayer) ! electronic temperature (K)
59	real :: dens(nlayer) ! density (cm-3)
60	real :: zmmean(nlayer) ! mean molar mass (g/mole)
61	real :: dist_sol ! sun distance (au)
62	real :: zday ! date (time since Ls=0, in martian days)
63	real :: surfdust1d(nlayer) ! dust surface area (cm2/cm3)
64	real :: surfice1d(nlayer) ! ice surface area (cm2/cm3)
65	real :: tau ! dust optical depth
66
67	!===================================================================
68	! input/output:
69	!===================================================================
70
71	real :: zycol(nlayer,nq) ! chemical species volume mixing ratio
72
73	!===================================================================
74	! output:
75	!===================================================================
76
77	integer :: iter(nlayer) ! iteration counter
78	real :: jo3(nlayer) ! photodissociation rate o3 -> o1d
79	real :: jh2o(nlayer) ! photodissociation rate h2o -> h + oh
80	real :: em_no(nlayer) ! NO nightglow volume emission rate
81	real :: em_o2(nlayer) ! O2 nightglow volume emission rate
82
83	!===================================================================
84	! local:
85	!===================================================================
86
87	integer :: phychemrat ! (physical timestep)/(nominal chemical timestep)
88	integer :: j_o3_o1d, j_h2o, ilev, iesp
89	integer :: lswitch
90	logical, save :: firstcall = .true.
91
92	!$OMP THREADPRIVATE(firstcall)
93
94	logical :: jionos ! switch for J parameterization
95
96	! tracer indexes in the chemistry:
97	! Species always considered in the chemistry
98	integer,parameter :: i_co2 = 1
99	integer,parameter :: i_co = 2
100	integer,parameter :: i_o = 3
101	integer,parameter :: i_o1d = 4
102	integer,parameter :: i_o2 = 5
103	integer,parameter :: i_o3 = 6
104	integer,parameter :: i_h = 7
105	integer,parameter :: i_h2 = 8
106	integer,parameter :: i_oh = 9
107	integer,parameter :: i_ho2 = 10
108	integer,parameter :: i_h2o2 = 11
109	integer,parameter :: i_h2o = 12
110	integer,parameter :: i_n = 13
111	integer,parameter :: i_n2d = 14
112	integer,parameter :: i_no = 15
113	integer,parameter :: i_no2 = 16
114	integer,parameter :: i_n2 = 17
115	!Species that may or not be included
116	!Their indices will be defined, if species are to be considered, in firstcall
117	integer,save :: i_co2plus = 0
118	integer,save :: i_oplus = 0
119	integer,save :: i_o2plus = 0
120	integer,save :: i_noplus = 0
121	integer,save :: i_coplus = 0
122	integer,save :: i_cplus = 0
123	integer,save :: i_n2plus = 0
124	integer,save :: i_nplus = 0
125	integer,save :: i_hplus = 0
126	integer,save :: i_hco2plus = 0
127	integer,save :: i_hcoplus = 0
128	integer,save :: i_h2oplus = 0
129	integer,save :: i_h3oplus = 0
130	integer,save :: i_ohplus = 0
131	integer,save :: i_elec = 0
132	integer,save :: i_hdo = 0
133	integer,save :: i_od = 0
134	integer,save :: i_d = 0
135	integer,save :: i_hd = 0
136	integer,save :: i_do2 = 0
137	integer,save :: i_hdo2 = 0
138
139	!$OMP THREADPRIVATE(i_co2plus,i_oplus,i_o2plus,i_noplus,i_coplus,i_cplus,i_n2plus,i_nplus,i_hplus ,i_hco2plus,i_hcoplus,i_h2oplus,i_h3oplus,i_ohplus,i_elec,i_hdo,i_od,i_d,i_hd,i_do2,i_hdo2)
140
141	!integer,parameter :: i_co2plus = 18
142	!integer,parameter :: i_oplus = 19
143	!integer,parameter :: i_o2plus = 20
144	!integer,parameter :: i_noplus = 21
145	!integer,parameter :: i_coplus = 22
146	!integer,parameter :: i_cplus = 23
147	!integer,parameter :: i_n2plus = 24
148	!integer,parameter :: i_nplus = 25
149	!integer,parameter :: i_hplus = 26
150	!integer,parameter :: i_hco2plus= 27
151	!integer,parameter :: i_hcoplus = 28
152	!integer,parameter :: i_h2oplus = 29
153	!integer,parameter :: i_h3oplus = 30
154	!integer,parameter :: i_ohplus = 31
155	!integer,parameter :: i_elec = 32
156	!integer,parameter :: i_hdo = 33
157	!integer,parameter :: i_od = 34
158	!integer,parameter :: i_d = 35
159	!integer,parameter :: i_hd = 36
160	!integer,parameter :: i_do2 = 37
161	!integer,parameter :: i_hdo2 = 38
162
163	integer :: i_last
164
165	!Tracer indexes for photionization coeffs
166
167	integer,parameter :: induv_co2 = 1
168	integer,parameter :: induv_o2 = 2
169	integer,parameter :: induv_o = 3
170	integer,parameter :: induv_h2o = 4
171	integer,parameter :: induv_h2 = 5
172	integer,parameter :: induv_h2o2= 6
173	integer,parameter :: induv_o3 = 7
174	integer,parameter :: induv_n2 = 8
175	integer,parameter :: induv_n = 9
176	integer,parameter :: induv_no = 10
177	integer,parameter :: induv_co = 11
178	integer,parameter :: induv_h = 12
179	integer,parameter :: induv_no2 = 13
180
181	integer :: ilay
182	integer :: ind_norec
183	integer :: ind_orec
184
185	real :: ctimestep ! standard timestep for the chemistry (s)
186	real :: dt_guess ! first-guess timestep (s)
187	real :: dt_corrected ! corrected timestep (s)
188	real :: time ! internal time (between 0 and ptimestep, in s)
189
190	real, dimension(nlayer,nesp) :: rm ! mixing ratios
191	real (kind = 8), dimension(nesp) :: cold ! number densities at previous timestep (molecule.cm-3)
192	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities at current timestep (molecule.cm-3)
193	real (kind = 8), dimension(nesp) :: cnew ! number densities at next timestep (molecule.cm-3)
194
195	! reaction rates
196
197	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
198	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
199	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
200	logical :: hetero_dust, hetero_ice
201
202	! matrix
203
204	real (kind = 8), dimension(nesp,nesp) :: mat, mat1
205	integer, dimension(nesp) :: indx
206	integer :: code
207
208	! production and loss terms (for first-guess solution only)
209
210	real (kind = 8), dimension(nesp) :: prod, loss
211
212	!===================================================================
213	! initialisation of the reaction indexes
214	!===================================================================
215
216	if (firstcall) then
217	!Indices for species that may or not be considered
218	i_last = i_n2
219	if(ionchem) then
220	i_co2plus = i_last + 1
221	i_oplus = i_last + 2
222	i_o2plus = i_last + 3
223	i_noplus = i_last + 4
224	i_coplus = i_last + 5
225	i_cplus = i_last + 6
226	i_n2plus = i_last + 7
227	i_nplus = i_last + 8
228	i_hplus = i_last + 9
229	i_hco2plus = i_last + 10
230	i_hcoplus = i_last + 11
231	i_h2oplus = i_last + 12
232	i_h3oplus = i_last + 13
233	i_ohplus = i_last + 14
234	i_elec = i_last + 15
235	!Update last used index
236	i_last = i_elec
237	endif
238	if(deutchem) then
239	i_hdo = i_last + 1
240	i_od = i_last + 2
241	i_d = i_last + 3
242	i_hd = i_last + 4
243	i_do2 = i_last + 5
244	i_hdo2 = i_last + 6
245	i_last = i_hdo2
246	endif
247	print*,'photochemistry: check number of indices',i_last,nesp
248	print*,'photochemistry: initialize indexes'
249	call indice(nb_reaction_3_max,nb_reaction_4_max, &
250	nb_phot_max, ionchem, deutchem, i_co2, i_co, i_o, &
251	i_o1d, i_o2, i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, &
252	i_h2o, i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
253	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
254	i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus, &
255	i_h2oplus, i_h3oplus, i_ohplus, i_elec, &
256	i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2)
257	firstcall = .false.
258	end if
259
260	!===================================================================
261	! initialisation of mixing ratios and densities
262	!===================================================================
263
264	call gcmtochim(nlayer, ionchem, deutchem, nq, zycol, lswitch, &
265	nesp, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
266	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
267	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
268	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
269	i_n2plus, i_nplus, i_hplus, i_hco2plus, &
270	i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, &
271	i_elec, i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2, &
272	dens, rm, c)
273
274	!===================================================================
275	! photolysis rates
276	!===================================================================
277
278	jionos = .true.
279
280	if (jonline) then
281	if (sza <= 113.) then ! day at 300 km
282	call photolysis_online(nlayer, deutchem, nb_phot_max, alt, &
283	press, temp, zmmean, &
284	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
285	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
286	i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, &
287	tau, sza, dist_sol, v_phot)
288
289	!Calculation of photoionization rates, if needed
290	if (jionos .and. ionchem) then
291	call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
292	do ilay=1,lswitch-1
293	call phdisrate(ig,nlayer,2,sza,ilay)
294	end do
295	!CO2 photoionization
296	v_phot(:,nphot+ 1) = jion(induv_co2,:,1)
297	v_phot(:,nphot+ 2) = jion(induv_co2,:,2)
298	v_phot(:,nphot+ 3) = jion(induv_co2,:,2)
299	v_phot(:,nphot+ 4) = jion(induv_co2,:,3)
300	v_phot(:,nphot+ 5) = jion(induv_co2,:,3)
301	v_phot(:,nphot+ 6) = jion(induv_co2,:,4)
302	v_phot(:,nphot+ 7) = jion(induv_co2,:,4)
303	!O2 photoionization
304	v_phot(:,nphot+ 8) = jion(induv_o2,:,1)
305	!O photoionization
306	v_phot(:,nphot+ 9) = jion(induv_o,:,1)
307	!NO photoionization
308	v_phot(:,nphot+10) = jion(induv_no,:,1)
309	!CO photoionization
310	v_phot(:,nphot+11) = jion(induv_co,:,1)
311	v_phot(:,nphot+12) = jion(induv_co,:,2)
312	v_phot(:,nphot+13) = jion(induv_co,:,2)
313	!N2 photoionization
314	v_phot(:,nphot+14) = jion(induv_n2,:,1)
315	v_phot(:,nphot+15) = jion(induv_n2,:,2)
316	v_phot(:,nphot+16) = jion(induv_n2,:,2)
317	!N photoionization
318	v_phot(:,nphot+17) = jion(induv_n,:,1)
319	!H photoionization
320	v_phot(:,nphot+18) = jion(induv_h,:,1)
321	end if
322
323	else ! night
324	v_phot(:,:) = 0.
325	end if
326	!else if(jparam) then
327	! call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
328	! do ilay=1,lswitch-1
329	! call phdisrate(ig,nlayer,2,sza,ilay)
330	! enddo
331	! v_phot(:,1)=jdistot(2,:)
332	! v_phot(:,2)=jdistot_b(2,:)
333	! v_phot(:,3)=jdistot(1,:)
334	! v_phot(:,4)=jdistot_b(1,:)
335	! v_phot(:,5)=jdistot(7,:)
336	! v_phot(:,6)=jdistot_b(7,:)
337	! v_phot(:,7)=jdistot(4,:)
338	! v_phot(:,8)=jdistot(6,:)
339	! v_phot(:,10)=jdistot(5,:)
340	! v_phot(:,11)=jdistot(10,:)
341	! v_phot(:,12)=jdistot(13,:)
342	! v_phot(:,13)=jdistot(8,:)
343	! v_phot(:,14)=jion(1,:,1)
344	! v_phot(:,15)=jion(1,:,2)
345	! v_phot(:,16)=jion(1,:,2)
346	! v_phot(:,17)=jion(1,:,3)
347	! v_phot(:,18)=jion(1,:,3)
348	! v_phot(:,19)=jion(1,:,4)
349	! v_phot(:,20)=jion(1,:,4)
350	! v_phot(:,21)=jion(2,:,1)
351	! v_phot(:,22)=jion(3,:,1)
352	! v_phot(:,23)=jion(10,:,1)
353	! v_phot(:,24)=jion(11,:,1)
354	! v_phot(:,25)=jion(11,:,2)
355	! v_phot(:,26)=jion(11,:,2)
356	! v_phot(:,27)=jion(8,:,1)
357	! v_phot(:,28)=jion(8,:,2)
358	! v_phot(:,29)=jion(8,:,2)
359	! v_phot(:,30)=jion(9,:,1)
360	! v_phot(:,31)=jion(12,:,1)
361	else
362	tau = tau*7./press(1) ! dust in the lookup table is at 7 hpa
363	call photolysis(nlayer, nb_phot_max, lswitch, press, temp, sza, tau, &
364	zmmean, dist_sol,rm(:,i_co2), rm(:,i_o3), v_phot)
365	end if
366	! save o3 and h2o photolysis for output
367
368	j_o3_o1d = 5
369	jo3(:) = v_phot(:,j_o3_o1d)
370	j_h2o = 7
371	jh2o(:) = v_phot(:,j_h2o)
372
373	!===================================================================
374	! reaction rates
375	!===================================================================
376	! switches for heterogeneous chemistry
377	! hetero_ice : reactions on ice clouds
378	! hetero_dust : reactions on dust
379	!===================================================================
380
381	hetero_dust = .false.
382	hetero_ice = .true.
383
384	call reactionrates(nlayer, ionchem, deutchem, &
385	nb_reaction_3_max, nb_reaction_4_max, &
386	nb_phot_max, nphotion, lswitch, dens, c(:,i_co2), &
387	c(:,i_o2), c(:,i_o), c(:,i_n2), press, temp, &
388	temp_elect, hetero_dust, hetero_ice, &
389	surfdust1d, surfice1d, v_phot, v_3, v_4, ind_norec, &
390	ind_orec)
391
392	!===================================================================
393	! ctimestep : standard chemical timestep (s), defined as
394	! the fraction phychemrat of the physical timestep
395	!===================================================================
396
397	phychemrat = 1
398
399	ctimestep = ptimestep/real(phychemrat)
400
401	!print*, "ptimestep = ", ptimestep
402	!print*, "phychemrat = ", phychemrat
403	!print*, "ctimestep = ", ctimestep
404	!stop
405
406	!===================================================================
407	! loop over levels
408	!===================================================================
409
410	do ilev = 1,lswitch - 1
411
412	! initializations
413
414	time = 0.
415	iter(ilev) = 0
416	dt_guess = ctimestep
417	cold(:) = c(ilev,:)
418
419	! internal loop for the chemistry
420
421	do while (time < ptimestep)
422
423	iter(ilev) = iter(ilev) + 1 ! iteration counter
424
425	! first-guess: fill matrix
426
427	call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, &
428	nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
429	v_phot, v_3, v_4)
430
431	! adaptative evaluation of the sub time step
432
433	call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:), &
434	mat1, prod, loss, dens(ilev))
435
436	if (time + dt_corrected > ptimestep) then
437	dt_corrected = ptimestep - time
438	end if
439
440	! if (dt_corrected /= dt_guess) then ! the timestep has been modified
441
442	! form the matrix identity + mat*dt_corrected
443
444	mat(:,:) = mat1(:,:)*dt_corrected
445	do iesp = 1,nesp
446	mat(iesp,iesp) = 1. + mat(iesp,iesp)
447	end do
448
449	! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
450	cnew(:) = c(ilev,:)
451
452	#ifdef LAPACK
453	call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
454	#else
455	write(,) "photochemistry error, missing LAPACK routine dgesv"
456	call abort_physic("photochemistry","missing LAPACK routine",1)
457	#endif
458
459	! end if
460
461	! eliminate small values
462
463	where (cnew(:)/dens(ilev) < 1.e-30)
464	cnew(:) = 0.
465	end where
466
467	! update concentrations
468
469	cold(:) = c(ilev,:)
470	c(ilev,:) = cnew(:)
471
472	! force charge neutrality (mod fgg, july 2019)
473
474	if (ionchem) then
475	if(c(ilev,i_elec).ne.c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+&
476	c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+&
477	c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)+ &
478	c(ilev,i_hcoplus)+c(ilev,i_h2oplus)+c(ilev,i_h3oplus)+ &
479	c(ilev,i_ohplus)) then
480	c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
481	c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+ &
482	c(ilev,i_n2plus)+c(ilev,i_nplus)+c(ilev,i_hplus)+ &
483	c(ilev,i_hco2plus)+c(ilev,i_hcoplus)+c(ilev,i_h2oplus)+ &
484	c(ilev,i_h3oplus)+c(ilev,i_ohplus)
485	! write(,)'photochemistry/359'
486	! write(,)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
487	end if
488	end if
489
490	cnew(:) = 0.
491
492	! increment internal time
493
494	time = time + dt_corrected
495	dt_guess = dt_corrected ! first-guess timestep for next iteration
496
497	end do ! while (time < ptimestep)
498
499	end do ! ilev
500
501	!add calculation of NO and O2 nightglows
502
503	em_no(:)=c(:,i_o)c(:,i_n)v_4(:,ind_norec) !2.8e-17(300./temp(:)))*0.5
504	em_o2(:)=0.75c(:,i_o)c(:,i_o)c(:,i_co2)v_3(:,ind_orec) !2.59.46e-34exp(485./temp(:))*dens(:)
505
506	!===================================================================
507	! save chemical species for the gcm
508	!===================================================================
509
510	call chimtogcm(nlayer, ionchem, deutchem, nq, zycol, lswitch, &
511	nesp, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
512	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
513	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
514	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
515	i_n2plus, i_nplus, i_hplus, i_hco2plus, &
516	i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, &
517	i_elec, i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2, &
518	dens, c)
519	contains
520
521	!================================================================
522
523	subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
524	prod, loss, dens)
525
526	!================================================================
527	! iterative evaluation of the appropriate time step dtnew
528	! according to curvature criterion based on
529	! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
530	! with r = (tn - tn-1)/(tn+1 - tn)
531	!================================================================
532
533	implicit none
534
535	! input
536
537	integer :: nesp ! number of species in the chemistry
538
539	real :: dtold, ctimestep
540	real (kind = 8), dimension(nesp) :: cold, ccur
541	real (kind = 8), dimension(nesp,nesp) :: mat1
542	real (kind = 8), dimension(nesp) :: prod, loss
543	real :: dens
544
545	! output
546
547	real :: dtnew
548
549	! local
550
551	real (kind = 8), dimension(nesp) :: cnew
552	real (kind = 8), dimension(nesp,nesp) :: mat
553	real (kind = 8) :: atol, ratio, e, es, coef
554
555	integer :: code, iesp, iter
556	integer, dimension(nesp) :: indx
557
558	real :: dttest
559
560	! parameters
561
562	real (kind = 8), parameter :: dtmin = 10. ! minimum time step (s)
563	real (kind = 8), parameter :: vmrtol = 1.e-11 ! absolute tolerance on vmr
564	real (kind = 8), parameter :: rtol = 0.05 ! rtol recommended value : 0.1-0.02
565	integer, parameter :: niter = 3 ! number of iterations
566	real (kind = 8), parameter :: coefmax = 2.
567	real (kind = 8), parameter :: coefmin = 0.1
568	logical :: fast_guess = .true.
569
570	dttest = dtold ! dttest = dtold = dt_guess
571
572	atol = vmrtol*dens ! absolute tolerance in molecule.cm-3
573
574	do iter = 1,niter
575
576	if (fast_guess) then
577
578	! first guess : fast semi-implicit method
579
580	do iesp = 1, nesp
581	cnew(iesp) = (ccur(iesp) + prod(iesp)dttest)/(1. + loss(iesp)dttest)
582	end do
583
584	else
585
586	! first guess : form the matrix identity + mat*dt_guess
587
588	mat(:,:) = mat1(:,:)*dttest
589	do iesp = 1,nesp
590	mat(iesp,iesp) = 1. + mat(iesp,iesp)
591	end do
592
593	! form right-hand side (RHS) of the system
594
595	cnew(:) = ccur(:)
596
597	! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
598
599	#ifdef LAPACK
600	call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
601	#else
602	write(,) "photochemistry error, missing LAPACK routine dgesv"
603	call abort_physic("photochemistry","missing LAPACK routine",1)
604	#endif
605
606	end if
607
608	! ratio old/new subtimestep
609
610	ratio = dtold/dttest
611
612	! e : local error indicatocitr
613
614	e = 0.
615
616	do iesp = 1,nesp
617	es = 2.abs((ratiocnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp)) &
618	/(1. + ratio)/max(ccur(iesp)*rtol,atol))
619
620	if (es > e) then
621	e = es
622	end if
623	end do
624
625	! timestep correction
626
627	coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
628
629	dttest = max(dtmin,dttest*coef)
630	dttest = min(ctimestep,dttest)
631
632	end do ! iter
633
634	! new timestep
635
636	dtnew = dttest
637
638	end subroutine define_dt
639
640	!======================================================================
641
642	subroutine reactionrates(nlayer, ionchem, deutchem , &
643	nb_reaction_3_max, nb_reaction_4_max, &
644	nb_phot_max, nphotion, lswitch, dens, co2, o2, o, &
645	n2, press, t, t_elect, hetero_dust, hetero_ice, &
646	surfdust1d, surfice1d, v_phot, v_3, v_4,ind_norec, &
647	ind_orec)
648
649	!================================================================
650	! compute reaction rates !
651	!----------------------------------------------------------------
652	! reaction type array !
653	!----------------------------------------------------------------
654	! A + B --> C + D bimolecular v_4 !
655	! A + A --> B + C quadratic v_3 !
656	! A + C --> B + C quenching v_phot !
657	! A + ice --> B + C heterogeneous v_phot !
658	!================================================================
659
660	use comcstfi_h
661	use photolysis_mod, only : nphot
662
663	implicit none
664
665	!----------------------------------------------------------------------
666	! input
667	!----------------------------------------------------------------------
668
669	integer, intent(in) :: nlayer ! number of atmospheric layers
670	integer, intent(in) :: nb_reaction_3_max ! number of quadratic reactions
671	integer, intent(in) :: nb_reaction_4_max ! number of bimolecular reactions
672	integer, intent(in) :: nb_phot_max ! number of reactions treated numerically as photodissociations
673	integer, intent(in) :: nphotion ! number of photoionizations
674	logical, intent(in) :: ionchem
675	logical, intent(in) :: deutchem
676	integer :: lswitch ! interface level between lower
677	! atmosphere and thermosphere chemistries
678	real, dimension(nlayer) :: dens ! total number density (molecule.cm-3)
679	real, dimension(nlayer) :: press ! pressure (hPa)
680	real, dimension(nlayer) :: t ! temperature (K)
681	real, dimension(nlayer) :: t_elect ! electronic temperature (K)
682	real, dimension(nlayer) :: surfdust1d ! dust surface area (cm2.cm-3)
683	real, dimension(nlayer) :: surfice1d ! ice surface area (cm2.cm-3)
684	real (kind = 8), dimension(nlayer) :: co2 ! co2 number density (molecule.cm-3)
685	real (kind = 8), dimension(nlayer) :: o2 ! o2 number density (molecule.cm-3)
686	real (kind = 8), dimension(nlayer) :: o ! o number density (molecule.cm-3)
687	real (kind = 8), dimension(nlayer) :: n2 ! n2 number density (molecule.cm-3)
688	logical :: hetero_dust, hetero_ice ! switches for heterogeneous chemistry
689
690	!----------------------------------------------------------------------
691	! output
692	!----------------------------------------------------------------------
693
694	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
695	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
696	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
697	integer :: ind_norec
698	integer :: ind_orec
699
700	!----------------------------------------------------------------------
701	! local
702	!----------------------------------------------------------------------
703
704	integer :: ilev
705	integer :: nb_phot, nb_reaction_3, nb_reaction_4
706	real :: ak0, ak1, xpo, rate, rate1, rate2
707	real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
708	real :: k0, kinf, kf, kint, kca
709	real, dimension(nlayer) :: deq
710	real, dimension(nlayer) :: a001, a002, a003, &
711	b001, b002, b003, b004, b005, b006, b007, &
712	b008, b009, b010, b011, b012, &
713	c001, c002, c003, c004, c005, c006, c007, &
714	c008, c009, c010, c011, c012, c013, c014, &
715	c015, c016, c017, c018, &
716	d001, d002, d003, d004, d005, d006, d007, &
717	d008, d009, d010, d011, d012, &
718	e001, e002, &
719	f001, f002, f003, f004, f005, f006, f007, &
720	f008, f009, f010, f011, f012, f013, f014, &
721	f015, f016, f017, f018, f019, f020, f021, &
722	f022, f023, f024, f025, f026, f027, f028, &
723	f029, f030, f031, f032, &
724	i001, i002, i003, i004, i005, i006, &
725	i007, i008, i009, i010, i011, i012, &
726	i013, i014, i015, i016, i017, i018, i019, &
727	i020, i021, i022, i023, i024, i025, i026, &
728	i027, i028, i029, i030, i031, i032, i033, &
729	i034, i035, i036, i037, i038, i039, i040, &
730	i041, i042, i043, i044, i045, i046, i047, &
731	i048, i049, i050, i051, i052, i053, i054, &
732	i055, i056, i057, i058, i059, i060, i061, &
733	i062, i063, &
734	h001, h002, h003, h004, h005
735
736	!----------------------------------------------------------------------
737	! initialisation
738	!----------------------------------------------------------------------
739
740	nb_phot = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates
741	nb_reaction_3 = 0
742	nb_reaction_4 = 0
743
744	!----------------------------------------------------------------------
745	! oxygen reactions
746	!----------------------------------------------------------------------
747
748	!--- a001: o + o2 + co2 -> o3 + co2
749
750	! jpl 2003
751	!
752	! co2/n2 efficiency as a third body = 2.075
753	! from sehested et al., j. geophys. res., 100, 1995.
754
755	! a001(:) = 2.0756.0e-34(t(:)/300.)*(-2.4)dens(:)
756
757	! jpl 2019
758
759	a001(:) = 2.0756.1e-34(t(:)/298.)*(-2.4)dens(:)
760
761	nb_reaction_4 = nb_reaction_4 + 1
762	v_4(:,nb_reaction_4) = a001(:)
763
764	!--- a002: o + o + co2 -> o2 + co2
765
766	! Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
767
768	! a002(:) = 2.55.2e-35exp(900./t(:))*dens(:)
769
770	! Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
771
772	! a002(:) = 1.2e-32(300./t(:))(2.0)dens(:) ! yung expression
773	a002(:) = 2.59.46e-34exp(485./t(:))*dens(:) ! nist expression
774
775	nb_reaction_3 = nb_reaction_3 + 1
776	v_3(:,nb_reaction_3) = a002(:)
777	ind_orec=nb_reaction_3
778
779	!--- a003: o + o3 -> o2 + o2
780
781	! jpl 2003
782
783	a003(:) = 8.0e-12*exp(-2060./t(:))
784
785	nb_reaction_4 = nb_reaction_4 + 1
786	v_4(:,nb_reaction_4) = a003(:)
787
788	!----------------------------------------------------------------------
789	! o(1d) reactions
790	!----------------------------------------------------------------------
791
792	!--- b001: o(1d) + co2 -> o + co2
793
794	! jpl 2006
795
796	b001(:) = 7.5e-11*exp(115./t(:))
797
798	nb_phot = nb_phot + 1
799	v_phot(:,nb_phot) = b001(:)*co2(:)
800
801	!--- b002: o(1d) + h2o -> oh + oh
802
803	! jpl 2006
804
805	b002(:) = 1.63e-10*exp(60./t(:))
806
807	nb_reaction_4 = nb_reaction_4 + 1
808	v_4(:,nb_reaction_4) = b002(:)
809
810	!--- b003: o(1d) + h2 -> oh + h
811
812	! jpl 2011
813
814	b003(:) = 1.2e-10
815
816	nb_reaction_4 = nb_reaction_4 + 1
817	v_4(:,nb_reaction_4) = b003(:)
818
819	!--- b004: o(1d) + o2 -> o + o2
820
821	! jpl 2006
822
823	b004(:) = 3.3e-11*exp(55./t(:))
824
825	nb_phot = nb_phot + 1
826	v_phot(:,nb_phot) = b004(:)*o2(:)
827
828	!--- b005: o(1d) + o3 -> o2 + o2
829
830	! jpl 2003
831
832	b005(:) = 1.2e-10
833
834	nb_reaction_4 = nb_reaction_4 + 1
835	v_4(:,nb_reaction_4) = b005(:)
836
837	!--- b006: o(1d) + o3 -> o2 + o + o
838
839	! jpl 2003
840
841	b006(:) = 1.2e-10
842
843	nb_reaction_4 = nb_reaction_4 + 1
844	v_4(:,nb_reaction_4) = b006(:)
845
846	!--- b007: o(1d) + ch4 -> ch3 + oh
847
848	! jpl 2003
849
850	b007(:) = 1.5e-10*0.75
851
852	!--- b008: o(1d) + ch4 -> ch3o + h
853
854	! jpl 2003
855
856	b008(:) = 1.5e-10*0.20
857	!
858	!--- b009: o(1d) + ch4 -> ch2o + h2
859
860	! jpl 2003
861
862	b009(:) = 1.5e-10*0.05
863
864	if(deutchem) then
865	!
866	!--- b010: o(1d) + hdo -> od + oh
867	b010(:) = b002(:)
868
869	nb_reaction_4 = nb_reaction_4 + 1
870	v_4(:,nb_reaction_4) = b010(:)
871
872	!
873	!--- b011: o(1d) + hd -> h + od
874
875	!Laurent et al., 1995
876
877	b011(:) = 1.3e-10
878
879	nb_reaction_4 = nb_reaction_4 + 1
880	v_4(:,nb_reaction_4) = b011(:)
881
882	!
883	!--- b012: o(1d) + hd -> d + oh
884
885	!Laurent et al., 1995
886
887	b012(:) = 1.0e-10
888
889	nb_reaction_4 = nb_reaction_4 + 1
890	v_4(:,nb_reaction_4) = b012(:)
891
892	endif !deutchem
893	!----------------------------------------------------------------------
894	! hydrogen reactions
895	!----------------------------------------------------------------------
896
897	!--- c001: o + ho2 -> oh + o2
898
899	! jpl 2003
900
901	c001(:) = 3.0e-11*exp(200./t(:))
902
903	nb_reaction_4 = nb_reaction_4 + 1
904	v_4(:,nb_reaction_4) = c001(:)
905
906	!--- c002: o + oh -> o2 + h
907
908	! jpl 2011
909
910	c002(:) = 1.8e-11*exp(180./t(:))
911
912	! c002(:) = c002(:)*1.12 ! FL li et al. 2007
913
914	! robertson and smith, j. chem. phys. a 110, 6673, 2006
915
916	! c002(:) = 11.2e-11t(:)(-0.32)exp(177./t(:))
917
918	nb_reaction_4 = nb_reaction_4 + 1
919	v_4(:,nb_reaction_4) = c002(:)
920
921	!--- c003: h + o3 -> oh + o2
922
923	! jpl 2003
924
925	c003(:) = 1.4e-10*exp(-470./t(:))
926
927	nb_reaction_4 = nb_reaction_4 + 1
928	v_4(:,nb_reaction_4) = c003(:)
929
930	!--- c004: h + ho2 -> oh + oh
931
932	! jpl 2006
933
934	c004(:) = 7.2e-11
935
936	nb_reaction_4 = nb_reaction_4 + 1
937	v_4(:,nb_reaction_4) = c004(:)
938
939	!--- c005: h + ho2 -> h2 + o2
940
941	! jpl 2006
942
943	c005(:) = 6.9e-12
944
945	nb_reaction_4 = nb_reaction_4 + 1
946	v_4(:,nb_reaction_4) = c005(:)
947
948	!--- c006: h + ho2 -> h2o + o
949
950	! jpl 2006
951
952	c006(:) = 1.6e-12
953
954	nb_reaction_4 = nb_reaction_4 + 1
955	v_4(:,nb_reaction_4) = c006(:)
956
957	!--- c007: oh + ho2 -> h2o + o2
958
959	! jpl 2003
960
961	! canty et al., grl, 2006 suggest to increase this rate
962	! by 20%. not done here.
963
964	c007(:) = 4.8e-11*exp(250./t(:))
965
966	! c007(:) = c007(:)*0.9 ! FL li et al. 2007
967
968	nb_reaction_4 = nb_reaction_4 + 1
969	v_4(:,nb_reaction_4) = c007(:)
970
971	!--- c008: ho2 + ho2 -> h2o2 + o2
972
973	! jpl 2015
974
975	c008(:) = 3.0e-13*exp(460./t(:))
976
977	! christensen et al., grl, 13, 2002
978
979	! c008(:) = 1.5e-12*exp(19./t(:))
980
981	nb_reaction_3 = nb_reaction_3 + 1
982	v_3(:,nb_reaction_3) = c008(:)
983
984	!--- c009: oh + h2o2 -> h2o + ho2
985
986	! jpl 2006
987
988	c009(:) = 1.8e-12
989
990	nb_reaction_4 = nb_reaction_4 + 1
991	v_4(:,nb_reaction_4) = c009(:)
992
993	!--- c010: oh + h2 -> h2o + h
994
995	! jpl 2006
996
997	c010(:) = 2.8e-12*exp(-1800./t(:))
998
999	nb_reaction_4 = nb_reaction_4 + 1
1000	v_4(:,nb_reaction_4) = c010(:)
1001
1002	!--- c011: h + o2 + co2 -> ho2 + co2
1003
1004	! jpl 2011
1005	! co2/n2 efficiency as a third body = 2.4
1006	! from ashman and haynes, 27th symposium on combustion, 1998.
1007
1008	do ilev = 1,lswitch-1
1009	! ak0 = 3.12.44.4e-32(t(ilev)/300.)*(-1.3) ! FL li et al 2017
1010	! ak0 = 2.44.4e-32(t(ilev)/300.)**(-1.3)
1011	! ak1 = 7.5e-11(t(ilev)/300.)*(0.2)
1012
1013	! jpl 2019
1014
1015	ak0 = 2.45.3e-32(t(ilev)/298.)**(-1.8)
1016	ak1 = 9.5e-11(t(ilev)/298.)*(0.4)
1017
1018	rate = (ak0dens(ilev))/(1. + ak0dens(ilev)/ak1)
1019	xpo = 1./(1. + alog10((ak0dens(ilev))/ak1)*2)
1020	c011(ilev) = rate0.6*xpo
1021	end do
1022
1023	nb_reaction_4 = nb_reaction_4 + 1
1024	v_4(:,nb_reaction_4) = c011(:)
1025
1026	!--- c012: o + h2o2 -> oh + ho2
1027
1028	! jpl 2003
1029
1030	c012(:) = 1.4e-12*exp(-2000./t(:))
1031
1032	nb_reaction_4 = nb_reaction_4 + 1
1033	v_4(:,nb_reaction_4) = c012(:)
1034
1035	!--- c013: oh + oh -> h2o + o
1036
1037	! jpl 2006
1038
1039	c013(:) = 1.8e-12
1040
1041	nb_reaction_3 = nb_reaction_3 + 1
1042	v_3(:,nb_reaction_3) = c013(:)
1043
1044	!--- c014: oh + o3 -> ho2 + o2
1045
1046	! jpl 2003
1047
1048	c014(:) = 1.7e-12*exp(-940./t(:))
1049
1050	nb_reaction_4 = nb_reaction_4 + 1
1051	v_4(:,nb_reaction_4) = c014(:)
1052
1053	!--- c015: ho2 + o3 -> oh + o2 + o2
1054
1055	! jpl 2003
1056
1057	c015(:) = 1.0e-14*exp(-490./t(:))
1058
1059	nb_reaction_4 = nb_reaction_4 + 1
1060	v_4(:,nb_reaction_4) = c015(:)
1061
1062	!--- c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
1063
1064	! jpl 2011
1065
1066	c016(:) = 2.52.1e-33exp(920./t(:))*dens(:)
1067
1068	nb_reaction_3 = nb_reaction_3 + 1
1069	v_3(:,nb_reaction_3) = c016(:)
1070
1071	!--- c017: oh + oh + co2 -> h2o2 + co2
1072
1073	! jpl 2003
1074
1075	do ilev = 1,lswitch-1
1076	ak0 = 2.56.9e-31(t(ilev)/300.)**(-1.0)
1077	ak1 = 2.6e-11(t(ilev)/300.)*(0.0)
1078
1079	rate = (ak0dens(ilev))/(1. + ak0dens(ilev)/ak1)
1080	xpo = 1./(1. + alog10((ak0dens(ilev))/ak1)*2)
1081	c017(ilev) = rate0.6*xpo
1082	end do
1083
1084	nb_reaction_3 = nb_reaction_3 + 1
1085	v_3(:,nb_reaction_3) = c017(:)
1086
1087	!--- c018: h + h + co2 -> h2 + co2
1088
1089	! baulch et al., 2005
1090
1091	c018(:) = 2.51.8e-30(t(:)*(-1.0))dens(:)
1092
1093	nb_reaction_3 = nb_reaction_3 + 1
1094	v_3(:,nb_reaction_3) = c018(:)
1095
1096
1097	!----------------------------------------------------------------------
1098	! nitrogen reactions
1099	!----------------------------------------------------------------------
1100
1101	!--- d001: no2 + o -> no + o2
1102
1103	! jpl 2006
1104
1105	! d001(:) = 5.1e-12*exp(210./t(:))
1106
1107	! jpl 2019
1108
1109	! For the sake of simplicity, it is assumed that the association yield (kf)
1110	! gives the same product as the chemical activation yield (kca).
1111	! Thus the only products are no + o2. There is no production of no3.
1112
1113	do ilev = 1,lswitch-1
1114
1115	! association
1116
1117	k0 = 2.53.4e-31(298./t(ilev))**(1.6)
1118	kinf = 2.3e-11(298./t(ilev))*(0.2)
1119
1120	kf = (kinfk0dens(ilev)/(kinf + k0*dens(ilev))) &
1121	0.6(1. + (log10(k0dens(ilev)/kinf))2.)(-1.0)
1122
1123	! chemical activation
1124
1125	kint = 5.3e-12*exp(200./t(ilev))
1126
1127	kca = kint*(1. - kf/kinf)
1128
1129	! total : association + chemical activation
1130
1131	d001(ilev) = kf + kca
1132
1133	end do
1134
1135	nb_reaction_4 = nb_reaction_4 + 1
1136	v_4(:,nb_reaction_4) = d001(:)
1137
1138	!--- d002: no + o3 -> no2 + o2
1139
1140	! jpl 2006
1141
1142	d002(:) = 3.0e-12*exp(-1500./t(:))
1143
1144	nb_reaction_4 = nb_reaction_4 + 1
1145	v_4(:,nb_reaction_4) = d002(:)
1146
1147	!--- d003: no + ho2 -> no2 + oh
1148
1149	! jpl 2011
1150
1151	! d003(:) = 3.3e-12*exp(270./t(:))
1152
1153	! jpl 2019
1154
1155	d003(:) = 3.44e-12*exp(260./t(:))
1156
1157	nb_reaction_4 = nb_reaction_4 + 1
1158	v_4(:,nb_reaction_4) = d003(:)
1159
1160	!--- d004: n + no -> n2 + o
1161
1162	! jpl 2011
1163
1164	d004(:) = 2.1e-11*exp(100./t(:))
1165
1166	nb_reaction_4 = nb_reaction_4 + 1
1167	v_4(:,nb_reaction_4) = d004(:)
1168
1169	!--- d005: n + o2 -> no + o
1170
1171	! jpl 2011
1172
1173	! d005(:) = 1.5e-11*exp(-3600./t(:))
1174
1175	! jpl 2019
1176
1177	d005(:) = 3.3e-12*exp(-3150./t(:))
1178
1179	nb_reaction_4 = nb_reaction_4 + 1
1180	v_4(:,nb_reaction_4) = d005(:)
1181
1182	!--- d006: no2 + h -> no + oh
1183
1184	! jpl 2011
1185
1186	! d006(:) = 4.0e-10*exp(-340./t(:))
1187
1188	! jpl 2019
1189
1190	d006(:) = 1.35e-10
1191
1192	nb_reaction_4 = nb_reaction_4 + 1
1193	v_4(:,nb_reaction_4) = d006(:)
1194
1195	!--- d007: n + o -> no
1196
1197	d007(:) = 2.8e-17(300./t(:))*0.5
1198
1199	nb_reaction_4 = nb_reaction_4 + 1
1200	v_4(:,nb_reaction_4) = d007(:)
1201	ind_norec = nb_reaction_4
1202
1203	!--- d008: n + ho2 -> no + oh
1204
1205	! brune et al., j. chem. phys., 87, 1983
1206
1207	d008(:) = 2.19e-11
1208
1209	nb_reaction_4 = nb_reaction_4 + 1
1210	v_4(:,nb_reaction_4) = d008(:)
1211
1212	!--- d009: n + oh -> no + h
1213
1214	! atkinson et al., j. phys. chem. ref. data, 18, 881, 1989
1215
1216	d009(:) = 3.8e-11*exp(85./t(:))
1217
1218	nb_reaction_4 = nb_reaction_4 + 1
1219	v_4(:,nb_reaction_4) = d009(:)
1220
1221	!--- d010: n(2d) + o -> n + o
1222
1223	! herron, j. phys. chem. ref. data, 1999
1224
1225	d010(:) = 3.3e-12*exp(-260./t(:))
1226
1227	nb_phot = nb_phot + 1
1228	v_phot(:,nb_phot) = d010(:)*o(:)
1229
1230	!--- d011: n(2d) + n2 -> n + n2
1231
1232	! herron, j. phys. chem. ref. data, 1999
1233
1234	d011(:) = 1.7e-14
1235
1236	nb_phot = nb_phot + 1
1237	v_phot(:,nb_phot) = d011(:)*n2(:)
1238
1239	!--- d012: n(2d) + co2 -> no + co
1240
1241	! herron, j. phys. chem. ref. data, 1999
1242
1243	d012(:) = 3.6e-13
1244
1245	nb_reaction_4 = nb_reaction_4 + 1
1246	v_4(:,nb_reaction_4) = d012(:)
1247
1248	!----------------------------------------------------------------------
1249	! carbon reactions
1250	!----------------------------------------------------------------------
1251
1252	!--- e001: oh + co -> co2 + h
1253
1254	! jpl 2015
1255	!
1256	! do ilev = 1,lswitch-1
1257	!
1258	! branch 1 : oh + co -> h + co2
1259	!
1260	! rate1 = 1.5e-13(t(ilev)/300.)*(0.0)
1261	!
1262	! branch 2 : oh + co + m -> hoco + m
1263	!
1264	! ak0 = 5.9e-33(t(ilev)/300.)*(-1.0)
1265	! ak1 = 1.1e-12(t(ilev)/300.)*(1.3)
1266	! rate2 = (ak0dens(ilev))/(1. + ak0dens(ilev)/ak1)
1267	! xpo = 1./(1. + alog10((ak0dens(ilev))/ak1)*2)
1268	!
1269	! e001(ilev) = rate1 + rate20.6*xpo
1270	! end do
1271
1272	! jpl 2019
1273
1274	! For the sake of simplicity, it is assumed that the association yield (kf)
1275	! gives the same product as the chemical activation yield (kca).
1276	! Thus the only products are h + co2. There is no production of hoco.
1277
1278	do ilev = 1,lswitch-1
1279
1280	! association
1281
1282	k0 = 2.56.9e-33(298./t(ilev))**(2.1)
1283	kinf = 1.1e-12(298./t(ilev))*(-1.3)
1284
1285	kf = (kinfk0dens(ilev)/(kinf + k0*dens(ilev))) &
1286	0.6(1. + (log10(k0dens(ilev)/kinf))2.)(-1.0)
1287
1288	! chemical activation
1289
1290	kint = 1.85e-13*exp(-65./t(ilev))
1291
1292	kca = kint*(1. - kf/kinf)
1293
1294	! total : association + chemical activation
1295
1296	e001(ilev) = kf + kca
1297
1298	end do
1299
1300	nb_reaction_4 = nb_reaction_4 + 1
1301	v_4(:,nb_reaction_4) = e001(:)
1302
1303	!--- e002: o + co + m -> co2 + m
1304
1305	! tsang and hampson, 1986.
1306
1307	e002(:) = 2.56.5e-33exp(-2184./t(:))*dens(:)
1308
1309	nb_reaction_4 = nb_reaction_4 + 1
1310	v_4(:,nb_reaction_4) = e002(:)
1311
1312
1313	!----------------------------------------------------------------------
1314	! deuterium reactions
1315	! only if deutchem = true
1316	!----------------------------------------------------------------------
1317
1318	if(deutchem) then
1319
1320	!--- f001: od + oh -> hdo + o
1321
1322	! Bedjanian, Y.; Le Bras, G.; and Poulet, G.,
1323	!Kinetic Study of OH + OH and OD + OD Reactions,
1324	!J. Phys. Chem. A, 103, pp. 7017 - 7025, 1999
1325
1326	f001(:) = 1.5e-12
1327
1328	nb_reaction_4 = nb_reaction_4 + 1
1329	v_4(:,nb_reaction_4) = f001(:)
1330
1331	!--- f002: od + h2 -> HDO + H
1332
1333	!Talukdar, R.K. et al.,
1334	!Kinetics of hydroxyl radical reactions with isotopically labeled hydrogen,
1335	!J. Phys. Chem., 100, pp. 3037 - 3043, 1996
1336
1337	f002(:) = 7.41e-15
1338
1339	nb_reaction_4 = nb_reaction_4 + 1
1340	v_4(:,nb_reaction_4) = f002(:)
1341
1342	!--- f003: od + ho2 -> hdo + o2
1343
1344	f003(:) = c007(:)
1345
1346	nb_reaction_4 = nb_reaction_4 + 1
1347	v_4(:,nb_reaction_4) = f003(:)
1348
1349	!--- f004: od + h2o2 -> hdo + ho2
1350
1351	!Vaghjiani, G.L. & Ravishankara, A.R.,
1352	!Reactions of OH and OD with H2O2 and D2O2,
1353	!J. Phys. Chem., 93, pp. 7833 - 7837, 1989;
1354
1355	f004(:) = 1.79e-12
1356
1357	nb_reaction_4 = nb_reaction_4 + 1
1358	v_4(:,nb_reaction_4) = f004(:)
1359
1360	!--- f005: o + od -> o2 + d
1361
1362	!Following Yung+1998, rate equal to that of O + OH -> O2 + H (c002)
1363
1364	f005(:) = c002(:)
1365
1366	nb_reaction_4 = nb_reaction_4 + 1
1367	v_4(:,nb_reaction_4) = f005(:)
1368
1369	!--- f006: od + h2 -> h2o + d
1370
1371	!Following Yung+1998, rate equal to that of OH + H2 -> H2O + H (c010)
1372
1373	f006(:) = c010(:)
1374
1375	nb_reaction_4 = nb_reaction_4 + 1
1376	v_4(:,nb_reaction_4) = f006(:)
1377
1378	!--- f007: od + h -> oh + d
1379
1380	!Rate following Yung+1988 and the rate of the inverse reaction (f012) from Atahan et al. J. Chem. Phys. 2005
1381
1382	f007(:) = 1.61e-10((298./t(:))0.32)exp(-16./t(:))
1383
1384	nb_reaction_4 = nb_reaction_4 + 1
1385	v_4(:,nb_reaction_4) = f007(:)
1386
1387	!--- f008: co + od -> co2 + d
1388
1389	!Following Yung+1988 rate equal to that of reaction CO + OH -> CO2 + H
1390
1391	f008(:) = e001(:)
1392
1393	nb_reaction_4 = nb_reaction_4 + 1
1394	v_4(:,nb_reaction_4) = f008(:)
1395
1396	!--- f009: o3 + d -> o2 + od
1397
1398	!Rate from NIST, Yu & Varandas, J. Chem. Soc. Faraday Trans., 93, 2651-2656, 1997
1399
1400	f009(:) = 7.41e-11*exp(-379./t(:))
1401
1402	nb_reaction_4 = nb_reaction_4 + 1
1403	v_4(:,nb_reaction_4) = f009(:)
1404
1405	!--- f010: HO2 + D -> OH + OD
1406
1407	!Following Yung+1988, rate equal to 0.71 times the rate for reaction HO2 + H -> OH + OH (c004)
1408
1409	f010(:) = 0.71*c004(:)
1410
1411	nb_reaction_4 = nb_reaction_4 + 1
1412	v_4(:,nb_reaction_4) = f010(:)
1413
1414	!--- f011: HO2 + D -> HDO + O
1415
1416	!Following Yung+1988, rate equal to 0.71 times the rate for reaction HO2 + H -> H2O + O (c006)
1417
1418	f011(:) = 0.71*c006(:)
1419
1420	nb_reaction_4 = nb_reaction_4 + 1
1421	v_4(:,nb_reaction_4) = f011(:)
1422
1423	!--- f012: OH + D -> OD + H
1424
1425	!Rate from NIST, Atahan et al. J. Chem. Phys. 2005
1426
1427	f012(:) = 1.16e-10((298./t(:))0.32)exp(-16./t(:))
1428
1429	nb_reaction_4 = nb_reaction_4 + 1
1430	v_4(:,nb_reaction_4) = f012(:)
1431
1432	!--- f013: h + d + co2 -> hd + co2
1433
1434	!According to Yung et al. 1988, rate equal to that of H + H + CO2
1435	!(reaction c018). Source: baulch et al., 2005
1436
1437	f013(:) = c018(:)
1438
1439	nb_reaction_4 = nb_reaction_4 + 1
1440	v_4(:,nb_reaction_4) = f013(:)
1441
1442	!--- f014: D + HO2 -> HD + O2
1443
1444	!According to Yung et al., rate equal to 0.71 times the rate of
1445	! H + HO2 -> H2 + O2 (reaction c005, source JPL 2019)
1446
1447	f014(:) = 0.71*c005(:)
1448
1449	nb_reaction_4 = nb_reaction_4 + 1
1450	v_4(:,nb_reaction_4) = f014(:)
1451
1452	!--- f015: OH + HD -> HDO + H
1453
1454	!Talukdar et al., Kinetics of hydroxyl radical reactions with
1455	!isotopically labeled hydrogen, J. Phys. Chem. 100, 3037-3043, 1996
1456
1457	f015(:) = 8.5e-13*exp(-2130./t(:))
1458
1459	nb_reaction_4 = nb_reaction_4 + 1
1460	v_4(:,nb_reaction_4) = f015(:)
1461
1462	!--- f016: OH + HD -> H2O + D
1463
1464	!Talukdar et al., 1996
1465
1466	f016(:) = 4.15e-12*exp(-2130./t(:))
1467
1468	nb_reaction_4 = nb_reaction_4 + 1
1469	v_4(:,nb_reaction_4) = f016(:)
1470
1471	!--- f017: D + O2 + CO2 -> DO2 + CO2
1472
1473	!Breshears et al., Room-temperature rate constants for the reaction
1474	!D + O2 + M -> DO2 + M, M=Ar, D2, CO2 and F2. J. Chem. Soc. Faraday
1475	!Trans., 87, 2337-2355 (1991)
1476
1477	f017(:) = 1.59e-31
1478
1479	nb_reaction_4 = nb_reaction_4 + 1
1480	v_4(:,nb_reaction_4) = f017(:)
1481
1482	!--- f018: OD + O3 -> DO2 + O2
1483
1484	!According to Yung+1988, rate equal to that of reaccion
1485	!OH + O3 -> HO2 + O2, (reaction c014)
1486
1487	f018(:) = c014(:)
1488
1489	nb_reaction_4 = nb_reaction_4 + 1
1490	v_4(:,nb_reaction_4) = f018(:)
1491
1492	!--- f019: D + HO2 -> DO2 + H
1493
1494	!Yung et al., 1988
1495
1496	f019(:) = 1.0e-10
1497
1498	nb_reaction_4 = nb_reaction_4 + 1
1499	v_4(:,nb_reaction_4) = f019(:)
1500
1501	!--- f020: O + DO2 -> OD + O2
1502
1503	!According to Yung+1988, rate equal to that of O + HO2 -> OH + O2
1504	! -> reaction c001
1505
1506	f020(:) = c001(:)
1507
1508	nb_reaction_4 = nb_reaction_4 + 1
1509	v_4(:,nb_reaction_4) = f020(:)
1510
1511	!--- f021: H + DO2 -> OH + OD
1512
1513	!According to Yung+1988, rate equal to that of H + HO2 -> OH + OH
1514	! -> reaction c004
1515
1516	f021(:) = c004(:)
1517
1518	nb_reaction_4 = nb_reaction_4 + 1
1519	v_4(:,nb_reaction_4) = f021(:)
1520
1521	!--- f022: H + DO2 -> HD + O2
1522
1523	!According to Yung+1988, rate equal to that of H + HO2 -> H2 + O2
1524	! -> reaction c005
1525
1526	f022(:) = c005(:)
1527
1528	nb_reaction_4 = nb_reaction_4 + 1
1529	v_4(:,nb_reaction_4) = f022(:)
1530
1531	!--- f023: H + DO2 -> HDO + O
1532
1533	!According to Yung+1988, rate equal to that of H + HO2 -> H2O + O
1534	! -> reaction c006
1535
1536	f023(:) = c006(:)
1537
1538	nb_reaction_4 = nb_reaction_4 + 1
1539	v_4(:,nb_reaction_4) = f023(:)
1540
1541	!--- f024: H + DO2 -> HO2 + D
1542
1543	!Yung+1988
1544
1545	f024(:) = 1.85e-10*exp(-890./t(:))
1546
1547	nb_reaction_4 = nb_reaction_4 + 1
1548	v_4(:,nb_reaction_4) = f024(:)
1549
1550	!--- f025: OH + DO2 -> HDO + O2
1551
1552	!According to Yung+1988, rate equal to that of OH + HO2 -> H2O + O2
1553	! -> reaction c007
1554
1555	f025(:) = c007(:)
1556
1557	nb_reaction_4 = nb_reaction_4 + 1
1558	v_4(:,nb_reaction_4) = f025(:)
1559
1560	!--- f026: DO2 + O3 -> OD + O2 + O2
1561
1562	!According to Yung+1988, rate equal to that of the reaction
1563	! HO2 + O3 -> OH + O2 + O2 -> reaction c015
1564
1565	f026(:) = c015(:)
1566
1567	nb_reaction_4 = nb_reaction_4 + 1
1568	v_4(:,nb_reaction_4) = f026(:)
1569
1570	!--- f027: OD + OH + CO2 -> HDO2 + CO2
1571
1572	!According to Yung+1988, rate equal to that of the reaction
1573	! OH + OH + CO2 -> H2O2 + CO2 (reaction c017)
1574
1575	f027(:) = c017(:)
1576
1577	nb_reaction_4 = nb_reaction_4 + 1
1578	v_4(:,nb_reaction_4) = f027(:)
1579
1580	!--- f028: DO2 + HO2 -> HDO2 + O2
1581
1582	!According to Yung+1988, rate equal to that of HO2 + HO2 -> H2O2 + O2
1583	! (reaction c008)
1584
1585	f028(:) = c008(:)
1586
1587	nb_reaction_4 = nb_reaction_4 + 1
1588	v_4(:,nb_reaction_4) = f028(:)
1589
1590	!--- f029: O + HDO2 -> OD + HO2
1591
1592	!According to Yung+1988, rate half that of O + H2O2 -> OH + HO2
1593	! (reaction c012)
1594
1595	f029(:) = 0.5*c012(:)
1596
1597	nb_reaction_4 = nb_reaction_4 + 1
1598	v_4(:,nb_reaction_4) = f029(:)
1599
1600	!--- f030: O + HDO2 -> OH + DO2
1601
1602	!According to Yung+1988, rate half that of O + H2O2 -> OH + HO2
1603	! (reaction c012)
1604
1605	f030(:) = 0.5*c012(:)
1606
1607	nb_reaction_4 = nb_reaction_4 + 1
1608	v_4(:,nb_reaction_4) = f030(:)
1609
1610	!--- f031: OH + HDO2 -> HDO + HO2
1611
1612	!According to Yung+1988, rate half that of OH + H2O2 -> H2O + HO2
1613	! (reaction c009)
1614
1615	f031(:) = 0.5*c009(:)
1616
1617	nb_reaction_4 = nb_reaction_4 + 1
1618	v_4(:,nb_reaction_4) = f031(:)
1619
1620
1621	!--- f032: OH + HDO2 -> H2O + DO2
1622
1623	!According to Yung+1988, rate half that of OH + H2O2 -> H2O + HO2
1624	! (reaction c009)
1625
1626	f032(:) = 0.5*c009(:)
1627
1628	nb_reaction_4 = nb_reaction_4 + 1
1629	v_4(:,nb_reaction_4) = f032(:)
1630
1631	endif !deutchem
1632
1633	!----------------------------------------------------------------------
1634	! ionospheric reactions
1635	! only if ionchem=true
1636	!----------------------------------------------------------------------
1637
1638	if (ionchem) then
1639
1640	!--- i001: co2+ + o2 -> o2+ + co2
1641
1642	! aninich, j. phys. chem. ref. data 1993
1643
1644	i001(:) = 5.5e-11(300./t_elect(:))*0.5
1645
1646	nb_reaction_4 = nb_reaction_4 + 1
1647	v_4(:,nb_reaction_4) = i001(:)
1648
1649	!--- i002: co2+ + o -> o+ + co2
1650
1651	! UMIST database
1652
1653	i002(:) = 9.6e-11
1654
1655	nb_reaction_4 = nb_reaction_4 + 1
1656	v_4(:,nb_reaction_4) = i002(:)
1657
1658	!--- i003: co2+ + o -> o2+ + co
1659
1660	! UMIST database
1661
1662	i003(:) = 1.64e-10
1663
1664	nb_reaction_4 = nb_reaction_4 + 1
1665	v_4(:,nb_reaction_4) = i003(:)
1666
1667	!--- i004: o2+ + e- -> o + o
1668
1669	! Alge et al., J. Phys. B, At. Mol. Phys. 1983
1670
1671	i004(:) = 2.0e-7(300./t_elect(:))*0.7
1672
1673	nb_reaction_4 = nb_reaction_4 + 1
1674	v_4(:,nb_reaction_4) = i004(:)
1675
1676	!--- i005: o+ + co2 -> o2+ + co
1677
1678	! UMIST database
1679
1680	i005(:) = 9.4e-10
1681
1682	nb_reaction_4 = nb_reaction_4 + 1
1683	v_4(:,nb_reaction_4) = i005(:)
1684
1685
1686	!--- i006: co2+ + e- -> co + o
1687
1688	! UMIST database
1689
1690	i006(:) = 3.8e-7(300./t_elect(:))*0.5
1691
1692	nb_reaction_4 = nb_reaction_4 + 1
1693	v_4(:,nb_reaction_4) = i006(:)
1694
1695
1696	!--- i007: co2+ + no -> no+ + co2
1697
1698	! UMIST database
1699
1700	i007(:) = 1.2e-10
1701
1702	nb_reaction_4 = nb_reaction_4 + 1
1703	v_4(:,nb_reaction_4) = i007(:)
1704
1705	!--- i008: o2+ + no -> no+ + o2
1706
1707	! UMIST database
1708
1709	i008(:) = 4.6e-10
1710
1711	nb_reaction_4 = nb_reaction_4 + 1
1712	v_4(:,nb_reaction_4) = i008(:)
1713
1714	!--- i009: o2+ + n2 -> no+ + no
1715
1716	! Fox & Sung 2001
1717
1718	i009(:) = 1.0e-15
1719
1720	nb_reaction_4 = nb_reaction_4 + 1
1721	v_4(:,nb_reaction_4) = i009(:)
1722
1723	!--- i010: o2+ + n -> no+ + o
1724
1725	! Fox & Sung 2001
1726
1727	i010(:) = 1.0e-10
1728
1729	nb_reaction_4 = nb_reaction_4 + 1
1730	v_4(:,nb_reaction_4) = i010(:)
1731
1732	!--- i011: o+ + n2 -> no+ + n
1733
1734	! Fox & Sung 2001
1735
1736	i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
1737
1738	nb_reaction_4 = nb_reaction_4 + 1
1739	v_4(:,nb_reaction_4) = i011(:)
1740
1741	!--- i012: no+ + e -> n + o
1742
1743	! UMIST database
1744
1745	i012(:) = 4.3e-7(300./t_elect(:))*0.37
1746
1747	nb_reaction_4 = nb_reaction_4 + 1
1748	v_4(:,nb_reaction_4) = i012(:)
1749
1750
1751	!--- i013: co+ + co2 -> co2+ + co
1752
1753	! UMIST database
1754
1755	i013(:) = 1.0e-9
1756
1757	nb_reaction_4 = nb_reaction_4 + 1
1758	v_4(:,nb_reaction_4) = i013(:)
1759
1760
1761	!--- i014: co+ + o -> o+ + co
1762
1763	! UMIST database
1764
1765	i014(:) = 1.4e-10
1766
1767	nb_reaction_4 = nb_reaction_4 + 1
1768	v_4(:,nb_reaction_4) = i014(:)
1769
1770	!--- i015: c+ + co2 -> co+ + co
1771
1772	! UMIST database
1773
1774	i015(:) = 1.1e-9
1775
1776	nb_reaction_4 = nb_reaction_4 + 1
1777	v_4(:,nb_reaction_4) = i015(:)
1778
1779
1780	!--- i016: N2+ + co2 -> co2+ + N2
1781
1782	! Fox & Song 2001
1783
1784	i016(:) = 9.0e-10(300./t_elect(:))*0.23
1785
1786	nb_reaction_4 = nb_reaction_4 + 1
1787	v_4(:,nb_reaction_4) = i016(:)
1788
1789
1790	!--- i017: N2+ + o -> no+ + N
1791
1792	! Fox & Song 2001
1793
1794	i017(:) = 1.33e-10(300./t_elect(:))*0.44
1795
1796	nb_reaction_4 = nb_reaction_4 + 1
1797	v_4(:,nb_reaction_4) = i017(:)
1798
1799	!--- i018: N2+ + co -> co+ + N2
1800
1801	! UMIST
1802
1803	i018(:) = 7.4e-11
1804
1805	nb_reaction_4 = nb_reaction_4 + 1
1806	v_4(:,nb_reaction_4) = i018(:)
1807
1808	!--- i019: N2+ + e -> N + N
1809
1810	! UMIST
1811
1812	i019(:) = 7.7e-7(300./t_elect(:))*0.3
1813
1814	nb_reaction_4 = nb_reaction_4 + 1
1815	v_4(:,nb_reaction_4) = i016(:)
1816
1817	!--- i020: N2+ + o -> o+ + N2
1818
1819	! Fox & Song 2001
1820
1821	i020(:) = 7.0e-12(300./t_elect(:))*0.23
1822
1823	nb_reaction_4 = nb_reaction_4 + 1
1824	v_4(:,nb_reaction_4) = i020(:)
1825
1826	!--- i021: N+ + co2 -> co2+ + N
1827
1828	! UMIST
1829
1830	i021(:) = 7.5e-10
1831
1832	nb_reaction_4 = nb_reaction_4 + 1
1833	v_4(:,nb_reaction_4) = i021(:)
1834
1835	!--- i022: CO+ + H -> H+ + CO
1836
1837	! Fox & Sung 2001
1838
1839	i022(:) = 4.0e-10
1840
1841	nb_reaction_4 = nb_reaction_4 + 1
1842	v_4(:,nb_reaction_4) = i022(:)
1843
1844	!--- i023: O+ + H -> H+ + O
1845
1846	! UMIST
1847
1848	i023(:) = 5.66e-10((t_elect(:)/300.)0.36)exp(8.6/t_elect(:))
1849
1850	nb_reaction_4 = nb_reaction_4 + 1
1851	v_4(:,nb_reaction_4) = i023(:)
1852
1853	!--- i024: H+ + O -> O+ + H
1854
1855	! UMIST
1856
1857	i024(:) = 6.86e-10((t_elect(:)/300.)0.26)exp(-224.3/t_elect(:))
1858
1859	nb_reaction_4 = nb_reaction_4 + 1
1860	v_4(:,nb_reaction_4) = i024(:)
1861
1862	!--- i025: CO+ + H2 -> HCO2+ + H
1863
1864	! UMIST
1865
1866	i025(:) = 9.5e-10
1867
1868	nb_reaction_4 = nb_reaction_4 + 1
1869	v_4(:,nb_reaction_4) = i025(:)
1870
1871	!--- i026: HCO2+ + e -> H + CO2
1872
1873	! UMIST
1874
1875	i026(:) = 1.75e-8((300./t_elect(:))*0.5)
1876
1877	nb_reaction_4 = nb_reaction_4 + 1
1878	v_4(:,nb_reaction_4) = i026(:)
1879
1880	!--- i027+i028: HCO2+ + e -> H + O + CO
1881
1882	! UMIST
1883	!Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
1884	!i028: 0.5 (HCO2+ + e-) -> O + CO
1885
1886	i027(:) = 8.1e-7((300./t_elect(:))*0.64)
1887
1888	nb_reaction_4 = nb_reaction_4 + 1
1889	v_4(:,nb_reaction_4) = i027(:)
1890
1891	nb_reaction_4 = nb_reaction_4 + 1
1892	v_4(:,nb_reaction_4) = i027(:)
1893
1894	!--- i029: HCO2+ + e -> OH + CO
1895
1896	! UMIST
1897
1898	i029(:) = 3.2e-7((300./t_elect(:))*0.64)
1899
1900	nb_reaction_4 = nb_reaction_4 + 1
1901	v_4(:,nb_reaction_4) = i029(:)
1902
1903	!--- i030: HCO2+ + e -> H + CO2
1904
1905	i030(:) = 6.0e-8((300./t_elect(:))*0.64)
1906	nb_reaction_4 = nb_reaction_4 + 1
1907	v_4(:,nb_reaction_4) = i030(:)
1908
1909	!--- i031: HCO2+ + O -> HCO+ + O2
1910
1911	! UMIST
1912
1913	i031(:) = 1.e-9
1914	nb_reaction_4 = nb_reaction_4 + 1
1915	v_4(:,nb_reaction_4) = i031(:)
1916
1917	!--- i032: HCO2+ + CO -> HCO+ + CO2
1918
1919	! UMIST, from Prassad & Huntress 1980
1920
1921	i032(:) = 7.8e-10
1922	nb_reaction_4 = nb_reaction_4 + 1
1923	v_4(:,nb_reaction_4) = i032(:)
1924
1925	!--- i033: H+ + CO2 -> HCO+ + O
1926
1927	! UMIST, from Smith et al., Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
1928
1929	i033(:) = 3.5e-9
1930	nb_reaction_4 = nb_reaction_4 + 1
1931	v_4(:,nb_reaction_4) = i033(:)
1932
1933
1934	!--- i034: CO2+ + H -> HCO+ + O
1935
1936	! Seen in Fox 2015, from Borodi et al., Int. J. Mass Spectrom. 280, 218-225, 2009
1937
1938	i034(:) = 4.5e-10
1939	nb_reaction_4 = nb_reaction_4 + 1
1940	v_4(:,nb_reaction_4) = i034(:)
1941
1942	!--- i035: CO+ + H2 -> HCO+ + H
1943
1944	!UMIST, from Scott et al., J. Chem. Phys., 106, 3982-3987(1997)
1945
1946	i035(:) = 7.5e-10
1947	nb_reaction_4 = nb_reaction_4 + 1
1948	v_4(:,nb_reaction_4) = i035(:)
1949
1950	!--- i036: HCO+ + e- -> CO + H
1951
1952	!UMIST, from Mitchell, Phys. Rep., 186, 215 (1990)
1953
1954	i036(:) = 2.4e-7 ((300./t_elect(:))*0.69)
1955	nb_reaction_4 = nb_reaction_4 + 1
1956	v_4(:,nb_reaction_4) = i036(:)
1957
1958	!--- i037: CO2+ + H2O -> H2O+ + CO2
1959
1960	!UMIST, from Karpas, Z., Anicich, V.G., and Huntress, W.T., Chem. Phys. Lett., 59, 84 (1978)
1961
1962	i037(:) = 2.04e-9 ((300./t_elect(:))*0.5)
1963	nb_reaction_4 = nb_reaction_4 + 1
1964	v_4(:,nb_reaction_4) = i037(:)
1965
1966	!--- i038: CO+ + H2O -> H2O+ + CO
1967
1968	!UMIST, from Huntress, W.T., McEwan, M.J., Karpas, Z., and Anicich, V.G., Astrophys. J. Supp. Series, 44, 481 (1980)
1969
1970	i038(:) = 1.72e-9((300./t_elect(:))*0.5)
1971	nb_reaction_4 = nb_reaction_4 + 1
1972	v_4(:,nb_reaction_4) = i038(:)
1973
1974	!--- i039: O+ + H2O -> H2O+ + O
1975
1976	!UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
1977
1978	i039(:) = 3.2e-9((300./t_elect(:))*0.5)
1979	nb_reaction_4 = nb_reaction_4 + 1
1980	v_4(:,nb_reaction_4) = i039(:)
1981
1982	!--- i040: N2+ + H2O -> H2O+ + N2
1983
1984	!UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
1985
1986	i040(:) = 2.3e-9((300./t_elect(:))*0.5)
1987	nb_reaction_4 = nb_reaction_4 + 1
1988	v_4(:,nb_reaction_4) = i040(:)
1989
1990	!--- i041: N+ + H2O -> H2O+ + N
1991
1992	!UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
1993
1994	i041(:) = 2.8e-9((300./t_elect(:))*0.5)
1995	nb_reaction_4 = nb_reaction_4 + 1
1996	v_4(:,nb_reaction_4) = i041(:)
1997
1998
1999	!--- i042: H+ + H2O -> H2O+ + H
2000
2001	!UMIST, from D. Smith, P. Spanel and C. A. Mayhew, Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
2002
2003	i042(:) = 6.9e-9((300./t_elect(:))*0.5)
2004	nb_reaction_4 = nb_reaction_4 + 1
2005	v_4(:,nb_reaction_4) = i042(:)
2006
2007	!--- i043: H2O+ + O2 -> O2+ + H2O
2008
2009	!UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
2010
2011	i043(:) = 4.6e-10
2012	nb_reaction_4 = nb_reaction_4 + 1
2013	v_4(:,nb_reaction_4) = i043(:)
2014
2015	!--- i044: H2O+ + CO -> HCO+ + OH
2016
2017	!UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2018
2019	i044(:) = 5.0e-10
2020	nb_reaction_4 = nb_reaction_4 + 1
2021	v_4(:,nb_reaction_4) = i044(:)
2022
2023	!--- i045: H2O+ + O -> O2+ + H2
2024
2025	!UMIST, from Viggiano, A.A, Howarka, F., Albritton, D.L., Fehsenfeld, F.C., Adams, N.G., and Smith, D., Astrophys. J., 236, 492 (1980)
2026
2027	i045(:) = 4.0e-11
2028	nb_reaction_4 = nb_reaction_4 + 1
2029	v_4(:,nb_reaction_4) = i045(:)
2030
2031	!--- i046: H2O+ + NO -> NO+ + H2O
2032
2033	!UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
2034
2035	i046(:) = 2.7e-10
2036	nb_reaction_4 = nb_reaction_4 + 1
2037	v_4(:,nb_reaction_4) = i046(:)
2038
2039	!--- i047: H2O+ + e- -> H + H + O
2040
2041	!UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
2042
2043	i047(:) = 3.05e-7((300./t_elect(:))*0.5)
2044	nb_reaction_4 = nb_reaction_4 + 1
2045	v_4(:,nb_reaction_4) = i047(:)
2046
2047	!--- i048: H2O+ + e- -> H + OH
2048
2049	!UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
2050
2051	i048(:) = 8.6e-8((300./t_elect(:))*0.5)
2052	nb_reaction_4 = nb_reaction_4 + 1
2053	v_4(:,nb_reaction_4) = i048(:)
2054
2055	!--- i049: H2O+ + e- -> O + H2
2056
2057	!UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
2058
2059	i049(:) = 3.9e-8((300./t_elect(:))*0.5)
2060	nb_reaction_4 = nb_reaction_4 + 1
2061	v_4(:,nb_reaction_4) = i049(:)
2062
2063	!--- i050: H2O+ + H2O -> H3O+ + OH
2064
2065	!UMIST, from Huntress, W.T. and Pinizzotto, R.F., J. Chem. Phys., 59, 4742 (1973)
2066
2067	i050(:) = 2.1e-9((300./t_elect(:))*0.5)
2068	nb_reaction_4 = nb_reaction_4 + 1
2069	v_4(:,nb_reaction_4) = i050(:)
2070
2071
2072	!--- i051: H2O+ + H2 -> H3O+ + H
2073
2074	!UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
2075
2076	i051(:) = 6.4e-10
2077	nb_reaction_4 = nb_reaction_4 + 1
2078	v_4(:,nb_reaction_4) = i051(:)
2079
2080	!--- i052: HCO+ + H2O -> H3O+ + CO
2081
2082	!UMIST, from Adams, N.G., Smith, D., and Grief, D., Int. J. Mass Spectrom. Ion Phys., 26, 405 (1978)
2083
2084	i052(:) = 2.5e-9((300./t_elect(:))*0.5)
2085	nb_reaction_4 = nb_reaction_4 + 1
2086	v_4(:,nb_reaction_4) = i052(:)
2087
2088	!--- i053: H3O+ + e -> OH + H + H
2089
2090	!UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2091
2092	i053(:) = 3.05e-7((300./t_elect(:))*0.5)
2093	nb_reaction_4 = nb_reaction_4 + 1
2094	v_4(:,nb_reaction_4) = i053(:)
2095
2096	!--- i054: H3O + + e -> H2O + H
2097
2098	!UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2099
2100	i054(:) = 7.09e-8((300./t_elect(:))*0.5)
2101	nb_reaction_4 = nb_reaction_4 + 1
2102	v_4(:,nb_reaction_4) = i054(:)
2103
2104	!--- i055: H3O+ + e -> OH + H2
2105
2106	!UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2107
2108	i055(:) = 5.37e-8((300./t_elect(:))*0.5)
2109	nb_reaction_4 = nb_reaction_4 + 1
2110	v_4(:,nb_reaction_4) = i055(:)
2111
2112	!--- i056: H3O+ + e -> O + H2 + H
2113
2114	!UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2115
2116	i056(:) = 5.6e-9((300./t_elect(:))*0.5)
2117	nb_reaction_4 = nb_reaction_4 + 1
2118	v_4(:,nb_reaction_4) = i056(:)
2119
2120	nb_reaction_4 = nb_reaction_4 + 1
2121	v_4(:,nb_reaction_4) = i056(:)
2122
2123	!--- i057: O+ + H2 -> OH+ + H
2124
2125	!UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
2126
2127	i057(:) = 1.7e-9
2128	nb_reaction_4 = nb_reaction_4 + 1
2129	v_4(:,nb_reaction_4) = i057(:)
2130
2131	!--- i058: OH+ + O -> O2+ + H
2132
2133	!UMIST, from Prasad & Huntress, 1980, ApJS, 43, 1
2134
2135	i058(:) = 7.1e-10
2136	nb_reaction_4 = nb_reaction_4 + 1
2137	v_4(:,nb_reaction_4) = i058(:)
2138
2139	!--- i059: OH+ + CO2 -> HCO2+ + O
2140
2141	!UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2142
2143	i059(:) = 1.44e-9
2144	nb_reaction_4 = nb_reaction_4 + 1
2145	v_4(:,nb_reaction_4) = i059(:)
2146
2147	!--- i060: OH+ + CO -> HCO+ + O
2148
2149	!UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2150
2151	i060(:) = 1.05e-9
2152	nb_reaction_4 = nb_reaction_4 + 1
2153	v_4(:,nb_reaction_4) = i060(:)
2154
2155	!--- i061: OH+ + NO -> NO+ + OH (tasa de reacciÃ³n UMIST 3.59e-10)
2156
2157	!UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2158
2159	i061(:) = 3.59e-10
2160	nb_reaction_4 = nb_reaction_4 + 1
2161	v_4(:,nb_reaction_4) = i061(:)
2162
2163	!--- i062: OH+ + H2 -> H2O+ + H (tasa de reacciÃ³n UMIST 1.01e-9,
2164
2165	!UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2166
2167	i062(:) = 1.01e-9
2168	nb_reaction_4 = nb_reaction_4 + 1
2169	v_4(:,nb_reaction_4) = i062(:)
2170
2171	!--- i063: OH+ + O2 -> O2+ + OH (tasa de reacciÃ³n UMIST 5.9e-10
2172
2173	!UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2174
2175	i063(:) = 5.9e-10
2176	nb_reaction_4 = nb_reaction_4 + 1
2177	v_4(:,nb_reaction_4) = i063(:)
2178
2179	end if !ionchem
2180
2181	!----------------------------------------------------------------------
2182	! heterogeneous chemistry
2183	!----------------------------------------------------------------------
2184
2185	if (hetero_ice) then
2186
2187	! k = (surfacevgamma)/4 (s-1)
2188	! v = 100sqrt(8rt/(pim)) (cm s-1)
2189
2190	!--- h001: ho2 + ice -> products
2191
2192	! cooper and abbatt, 1996: gamma = 0.025
2193
2194	gam = 0.025
2195	h001(:) = surfice1d(:) &
2196	100.sqrt(8.8.31t(:)/(33.e-3pi))gam/4.
2197
2198	! h002: oh + ice -> products
2199
2200	! cooper and abbatt, 1996: gamma = 0.03
2201
2202	gam = 0.03
2203	h002(:) = surfice1d(:) &
2204	100.sqrt(8.8.31t(:)/(17.e-3pi))gam/4.
2205
2206	!--- h003: h2o2 + ice -> products
2207
2208	! gamma = 0. test value
2209
2210	gam = 0.
2211	h003(:) = surfice1d(:) &
2212	100.sqrt(8.8.31t(:)/(34.e-3pi))gam/4.
2213	else
2214	h001(:) = 0.
2215	h002(:) = 0.
2216	h003(:) = 0.
2217	end if
2218
2219	nb_phot = nb_phot + 1
2220	v_phot(:,nb_phot) = h001(:)
2221
2222	nb_phot = nb_phot + 1
2223	v_phot(:,nb_phot) = h002(:)
2224
2225	nb_phot = nb_phot + 1
2226	v_phot(:,nb_phot) = h003(:)
2227
2228	if (hetero_dust) then
2229
2230	!--- h004: ho2 + dust -> products
2231
2232	! jacob, 2000: gamma = 0.2
2233	! see dereus et al., atm. chem. phys., 2005
2234
2235	gam = 0.2
2236	h004(:) = surfdust1d(:) &
2237	100.sqrt(8.8.31t(:)/(33.e-3pi))gam/4.
2238
2239	!--- h005: h2o2 + dust -> products
2240
2241	! gamma = 5.e-4
2242	! see dereus et al., atm. chem. phys., 2005
2243
2244	gam = 5.e-4
2245	h005(:) = surfdust1d(:) &
2246	100.sqrt(8.8.31t(:)/(34.e-3pi))gam/4.
2247	else
2248	h004(:) = 0.
2249	h005(:) = 0.
2250	end if
2251
2252	nb_phot = nb_phot + 1
2253	v_phot(:,nb_phot) = h004(:)
2254
2255	nb_phot = nb_phot + 1
2256	v_phot(:,nb_phot) = h005(:)
2257
2258	end subroutine reactionrates
2259
2260	!======================================================================
2261
2262	subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, &
2263	nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
2264	v_phot, v_3, v_4)
2265
2266	!======================================================================
2267	! filling of the jacobian matrix
2268	!======================================================================
2269
2270	use types_asis
2271
2272	implicit none
2273
2274	! input
2275
2276	integer :: ilev ! level index
2277	integer :: nesp ! number of species in the chemistry
2278	integer, intent(in) :: nlayer ! number of atmospheric layers
2279	integer, intent(in) :: nb_reaction_3_max
2280	! number of quadratic reactions
2281	integer, intent(in) :: nb_reaction_4_max
2282	! number of bimolecular reactions
2283	integer, intent(in) :: nb_phot_max
2284	! number of processes treated numerically as photodissociations
2285
2286	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities
2287	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
2288	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
2289	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
2290
2291	! output
2292
2293	real (kind = 8), dimension(nesp,nesp), intent(out) :: mat ! matrix
2294	real (kind = 8), dimension(nesp), intent(out) :: prod, loss
2295
2296	! local
2297
2298	integer :: iesp
2299	integer :: ind_phot_2,ind_phot_4,ind_phot_6
2300	integer :: ind_3_2,ind_3_4,ind_3_6
2301	integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
2302	integer :: iphot,i3,i4
2303
2304	real(kind = 8) :: eps, eps_4 ! implicit/explicit coefficient
2305
2306	! initialisations
2307
2308	mat(:,:) = 0.
2309	prod(:) = 0.
2310	loss(:) = 0.
2311
2312	! photodissociations
2313	! or reactions a + c -> b + c
2314	! or reactions a + ice -> b + c
2315
2316	do iphot = 1,nb_phot_max
2317
2318	ind_phot_2 = indice_phot(iphot)%z2
2319	ind_phot_4 = indice_phot(iphot)%z4
2320	ind_phot_6 = indice_phot(iphot)%z6
2321
2322	mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
2323	mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
2324	mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
2325
2326	loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
2327	prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3v_phot(ilev,iphot)c(ilev,ind_phot_2)
2328	prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5v_phot(ilev,iphot)c(ilev,ind_phot_2)
2329
2330	end do
2331
2332	! reactions a + a -> b + c
2333
2334	do i3 = 1,nb_reaction_3_max
2335
2336	ind_3_2 = indice_3(i3)%z2
2337	ind_3_4 = indice_3(i3)%z4
2338	ind_3_6 = indice_3(i3)%z6
2339
2340	mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1v_3(ilev,i3)c(ilev,ind_3_2)
2341	mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3v_3(ilev,i3)c(ilev,ind_3_2)
2342	mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5v_3(ilev,i3)c(ilev,ind_3_2)
2343
2344	loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1v_3(ilev,i3)c(ilev,ind_3_2)
2345	prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3v_3(ilev,i3)c(ilev,ind_3_2)*c(ilev,ind_3_2)
2346	prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5v_3(ilev,i3)c(ilev,ind_3_2)*c(ilev,ind_3_2)
2347
2348	end do
2349
2350	! reactions a + b -> c + d
2351
2352	eps = 1.d-10
2353
2354	do i4 = 1,nb_reaction_4_max
2355
2356	ind_4_2 = indice_4(i4)%z2
2357	ind_4_4 = indice_4(i4)%z4
2358	ind_4_6 = indice_4(i4)%z6
2359	ind_4_8 = indice_4(i4)%z8
2360
2361	eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
2362	eps_4 = min(eps_4,1.0)
2363
2364	mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
2365	mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
2366	mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
2367	mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
2368	mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
2369	mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
2370	mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
2371	mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
2372
2373
2374	loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1v_4(ilev,i4)c(ilev,ind_4_4)
2375	loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3v_4(ilev,i4)c(ilev,ind_4_2)
2376	prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5v_4(ilev,i4)c(ilev,ind_4_2)*c(ilev,ind_4_4)
2377	prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7v_4(ilev,i4)c(ilev,ind_4_2)*c(ilev,ind_4_4)
2378
2379	end do
2380
2381	end subroutine fill_matrix
2382
2383	!================================================================
2384
2385	subroutine indice(nb_reaction_3_max, nb_reaction_4_max, &
2386	nb_phot_max, ionchem, deutchem, i_co2, i_co, i_o, &
2387	i_o1d, i_o2, i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, &
2388	i_h2o, i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
2389	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
2390	i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus, &
2391	i_h2oplus, i_h3oplus, i_ohplus, i_elec, &
2392	i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2)
2393
2394	!================================================================
2395	! set the "indice" arrays used to fill the jacobian matrix !
2396	!----------------------------------------------------------------
2397	! reaction type !
2398	!----------------------------------------------------------------
2399	! A + hv --> B + C photolysis indice_phot !
2400	! A + B --> C + D bimolecular indice_4 !
2401	! A + A --> B + C quadratic indice_3 !
2402	! A + C --> B + C quenching indice_phot !
2403	! A + ice --> B + C heterogeneous indice_phot !
2404	!================================================================
2405
2406	use types_asis
2407
2408	implicit none
2409
2410	! input
2411
2412	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2413	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2414	i_n, i_n2d, i_no, i_no2, i_n2, &
2415	i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
2416	i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, &
2417	i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec, &
2418	i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2
2419	integer, intent(in) :: nb_reaction_3_max
2420	! number of quadratic reactions
2421	integer, intent(in) :: nb_reaction_4_max
2422	! number of bimolecular reactions
2423	integer, intent(in) :: nb_phot_max
2424	! number of processes treated numerically as photodissociations
2425	logical, intent(in) :: ionchem
2426	logical, intent(in) :: deutchem
2427
2428	! local
2429
2430	integer :: nb_phot, nb_reaction_3, nb_reaction_4
2431	integer :: i_dummy
2432
2433	allocate (indice_phot(nb_phot_max))
2434	allocate (indice_3(nb_reaction_3_max))
2435	allocate (indice_4(nb_reaction_4_max))
2436
2437	i_dummy = 1
2438
2439	nb_phot = 0
2440	nb_reaction_3 = 0
2441	nb_reaction_4 = 0
2442
2443	!===========================================================
2444	! O2 + hv -> O + O
2445	!===========================================================
2446
2447	nb_phot = nb_phot + 1
2448
2449	indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy)
2450
2451	!===========================================================
2452	! O2 + hv -> O + O(1D)
2453	!===========================================================
2454
2455	nb_phot = nb_phot + 1
2456
2457	indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d)
2458
2459	!===========================================================
2460	! CO2 + hv -> CO + O
2461	!===========================================================
2462
2463	nb_phot = nb_phot + 1
2464
2465	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o)
2466
2467	!===========================================================
2468	! CO2 + hv -> CO + O(1D)
2469	!===========================================================
2470
2471	nb_phot = nb_phot + 1
2472
2473	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d)
2474
2475	!===========================================================
2476	! O3 + hv -> O2 + O(1D)
2477	!===========================================================
2478
2479	nb_phot = nb_phot + 1
2480
2481	indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d)
2482
2483	!===========================================================
2484	! O3 + hv -> O2 + O
2485	!===========================================================
2486
2487	nb_phot = nb_phot + 1
2488
2489	indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o)
2490
2491	!===========================================================
2492	! H2O + hv -> H + OH
2493	!===========================================================
2494
2495	nb_phot = nb_phot + 1
2496
2497	indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh)
2498
2499	!===========================================================
2500	! H2O2 + hv -> OH + OH
2501	!===========================================================
2502
2503	nb_phot = nb_phot + 1
2504
2505	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy)
2506
2507	!===========================================================
2508	! HO2 + hv -> OH + O
2509	!===========================================================
2510
2511	nb_phot = nb_phot + 1
2512
2513	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o)
2514
2515	!===========================================================
2516	! H2 + hv -> H + H
2517	!===========================================================
2518
2519	nb_phot = nb_phot + 1
2520
2521	indice_phot(nb_phot) = z3spec(1.0, i_h2, 1.0, i_h, 1.0, i_h)
2522
2523	!===========================================================
2524	! NO + hv -> N + O
2525	!===========================================================
2526
2527	nb_phot = nb_phot + 1
2528
2529	indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o)
2530
2531	!===========================================================
2532	! NO2 + hv -> NO + O
2533	!===========================================================
2534
2535	nb_phot = nb_phot + 1
2536
2537	indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o)
2538
2539	!===========================================================
2540	! N2 + hv -> N + N
2541	!===========================================================
2542
2543	nb_phot = nb_phot + 1
2544
2545	indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n)
2546
2547	!Only if deuterium chemistry included
2548	if (deutchem) then
2549	!===========================================================
2550	! HDO + hv -> OD + H
2551	!===========================================================
2552
2553	nb_phot = nb_phot + 1
2554
2555	indice_phot(nb_phot) = z3spec(1.0, i_hdo, 1.0, i_h, 1.0, i_od)
2556
2557	!===========================================================
2558	! HDO + hv -> D + OH
2559	!===========================================================
2560
2561	nb_phot = nb_phot + 1
2562
2563	indice_phot(nb_phot) = z3spec(1.0, i_hdo, 1.0, i_d, 1.0, i_oh)
2564
2565	endif !deutchem
2566
2567	!Only if ion chemistry included
2568	if (ionchem) then
2569
2570	!===========================================================
2571	! CO2 + hv -> CO2+ + e-
2572	!===========================================================
2573
2574	nb_phot = nb_phot + 1
2575
2576	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
2577
2578	!===========================================================
2579	! CO2 + hv -> O+ + CO + e-
2580	!===========================================================
2581	!We divide this reaction in two
2582
2583	!0.5 CO2 + hv -> CO
2584	nb_phot = nb_phot + 1
2585
2586	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
2587
2588	!0.5 CO2 + hv -> O+ + e-
2589	nb_phot = nb_phot + 1
2590
2591	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
2592
2593	!===========================================================
2594	! CO2 + hv -> CO+ + O + e-
2595	!===========================================================
2596	!We divide this reaction in two
2597
2598	!0.5 CO2 + hv -> O
2599	nb_phot = nb_phot + 1
2600
2601	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
2602
2603	!0.5 CO2 + hv -> CO+ + e-
2604	nb_phot = nb_phot + 1
2605
2606	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
2607
2608	!===========================================================
2609	! CO2 + hv -> C+ + O2 + e-
2610	!===========================================================
2611	!We divide this reaction in two
2612
2613	!0.5 CO2 + hv -> O2
2614	nb_phot = nb_phot + 1
2615
2616	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
2617
2618	!0.5 CO2 + hv -> C+ + e-
2619	nb_phot = nb_phot + 1
2620
2621	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
2622
2623	!===========================================================
2624	! O2 + hv -> O2+ + e-
2625	!===========================================================
2626
2627	nb_phot = nb_phot + 1
2628
2629	indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
2630
2631	!===========================================================
2632	! O + hv -> O+ + e-
2633	!===========================================================
2634
2635	nb_phot = nb_phot + 1
2636
2637	indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
2638
2639	!===========================================================
2640	! NO + hv -> NO+ + e-
2641	!===========================================================
2642
2643	nb_phot = nb_phot + 1
2644
2645	indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
2646
2647	!===========================================================
2648	! CO + hv -> CO+ + e-
2649	!===========================================================
2650
2651	nb_phot = nb_phot + 1
2652
2653	indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
2654
2655	!===========================================================
2656	! CO + hv -> C+ + O + e-
2657	!===========================================================
2658	!We divide this reaction in two
2659
2660	!0.5 CO + hv -> O
2661	nb_phot = nb_phot + 1
2662
2663	indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
2664
2665	!0.5 CO + hv -> C+ + e-
2666	nb_phot = nb_phot + 1
2667
2668	indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
2669
2670	!===========================================================
2671	! N2 + hv -> N2+ + e-
2672	!===========================================================
2673
2674	nb_phot = nb_phot + 1
2675
2676	indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
2677
2678	!===========================================================
2679	! N2 + hv -> N+ + N + e-
2680	!===========================================================
2681	!We divide this reaction in two
2682
2683	!0.5 N2 + hv -> N
2684	nb_phot = nb_phot + 1
2685
2686	indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
2687
2688	!0.5 N2 + hv -> N+ + e-
2689	nb_phot = nb_phot + 1
2690
2691	indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
2692
2693	!===========================================================
2694	! N + hv -> N+ + e-
2695	!===========================================================
2696
2697	nb_phot = nb_phot + 1
2698
2699	indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
2700
2701	!===========================================================
2702	! H + hv -> H+ + e-
2703	!===========================================================
2704
2705	nb_phot = nb_phot + 1
2706
2707	indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
2708
2709	end if !ionchem
2710
2711	!===========================================================
2712	! a001 : O + O2 + CO2 -> O3 + CO2
2713	!===========================================================
2714
2715	nb_reaction_4 = nb_reaction_4 + 1
2716
2717	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
2718
2719	!===========================================================
2720	! a002 : O + O + CO2 -> O2 + CO2
2721	!===========================================================
2722
2723	nb_reaction_3 = nb_reaction_3 + 1
2724
2725	indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
2726
2727	!===========================================================
2728	! a003 : O + O3 -> O2 + O2
2729	!===========================================================
2730
2731	nb_reaction_4 = nb_reaction_4 + 1
2732
2733	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
2734
2735	!===========================================================
2736	! b001 : O(1D) + CO2 -> O + CO2
2737	!===========================================================
2738
2739	nb_phot = nb_phot + 1
2740
2741	indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
2742
2743	!===========================================================
2744	! b002 : O(1D) + H2O -> OH + OH
2745	!===========================================================
2746
2747	nb_reaction_4 = nb_reaction_4 + 1
2748
2749	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
2750
2751	!===========================================================
2752	! b003 : O(1D) + H2 -> OH + H
2753	!===========================================================
2754
2755	nb_reaction_4 = nb_reaction_4 + 1
2756
2757	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
2758
2759	!===========================================================
2760	! b004 : O(1D) + O2 -> O + O2
2761	!===========================================================
2762
2763	nb_phot = nb_phot + 1
2764
2765	indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
2766
2767	!===========================================================
2768	! b005 : O(1D) + O3 -> O2 + O2
2769	!===========================================================
2770
2771	nb_reaction_4 = nb_reaction_4 + 1
2772
2773	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
2774
2775	!===========================================================
2776	! b006 : O(1D) + O3 -> O2 + O + O
2777	!===========================================================
2778
2779	nb_reaction_4 = nb_reaction_4 + 1
2780
2781	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
2782
2783
2784	if(deutchem) then
2785	!===========================================================
2786	! b010 : O(1D) + HDO -> OD + OH
2787	!===========================================================
2788
2789	nb_reaction_4 = nb_reaction_4 + 1
2790
2791	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hdo, 1.0, i_od, 1.0, i_oh)
2792
2793
2794	!===========================================================
2795	! b011 : O(1D) + HD -> H + OD
2796	!===========================================================
2797
2798	nb_reaction_4 = nb_reaction_4 + 1
2799
2800	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hd, 1.0, i_h, 1.0, i_od)
2801
2802
2803	!===========================================================
2804	! b012 : O(1D) + HD -> D + OH
2805	!===========================================================
2806
2807	nb_reaction_4 = nb_reaction_4 + 1
2808
2809	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hd, 1.0, i_d, 1.0, i_oh)
2810
2811	endif !deutchem
2812
2813	!===========================================================
2814	! c001 : O + HO2 -> OH + O2
2815	!===========================================================
2816
2817	nb_reaction_4 = nb_reaction_4 + 1
2818
2819	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
2820
2821	!===========================================================
2822	! c002 : O + OH -> O2 + H
2823	!===========================================================
2824
2825	nb_reaction_4 = nb_reaction_4 + 1
2826
2827	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
2828
2829	!===========================================================
2830	! c003 : H + O3 -> OH + O2
2831	!===========================================================
2832
2833	nb_reaction_4 = nb_reaction_4 + 1
2834
2835	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
2836
2837	!===========================================================
2838	! c004 : H + HO2 -> OH + OH
2839	!===========================================================
2840
2841	nb_reaction_4 = nb_reaction_4 + 1
2842
2843	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
2844
2845	!===========================================================
2846	! c005 : H + HO2 -> H2 + O2
2847	!===========================================================
2848
2849	nb_reaction_4 = nb_reaction_4 + 1
2850
2851	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
2852
2853	!===========================================================
2854	! c006 : H + HO2 -> H2O + O
2855	!===========================================================
2856
2857	nb_reaction_4 = nb_reaction_4 + 1
2858
2859	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
2860
2861	!===========================================================
2862	! c007 : OH + HO2 -> H2O + O2
2863	!===========================================================
2864
2865	nb_reaction_4 = nb_reaction_4 + 1
2866
2867	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
2868
2869	!===========================================================
2870	! c008 : HO2 + HO2 -> H2O2 + O2
2871	!===========================================================
2872
2873	nb_reaction_3 = nb_reaction_3 + 1
2874
2875	indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2876
2877	!===========================================================
2878	! c009 : OH + H2O2 -> H2O + HO2
2879	!===========================================================
2880
2881	nb_reaction_4 = nb_reaction_4 + 1
2882
2883	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
2884
2885	!===========================================================
2886	! c010 : OH + H2 -> H2O + H
2887	!===========================================================
2888
2889	nb_reaction_4 = nb_reaction_4 + 1
2890
2891	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
2892
2893	!===========================================================
2894	! c011 : H + O2 + CO2 -> HO2 + CO2
2895	!===========================================================
2896
2897	nb_reaction_4 = nb_reaction_4 + 1
2898
2899	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
2900
2901	!===========================================================
2902	! c012 : O + H2O2 -> OH + HO2
2903	!===========================================================
2904
2905	nb_reaction_4 = nb_reaction_4 + 1
2906
2907	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
2908
2909	!===========================================================
2910	! c013 : OH + OH -> H2O + O
2911	!===========================================================
2912
2913	nb_reaction_3 = nb_reaction_3 + 1
2914
2915	indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
2916
2917	!===========================================================
2918	! c014 : OH + O3 -> HO2 + O2
2919	!===========================================================
2920
2921	nb_reaction_4 = nb_reaction_4 + 1
2922
2923	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
2924
2925	!===========================================================
2926	! c015 : HO2 + O3 -> OH + O2 + O2
2927	!===========================================================
2928
2929	nb_reaction_4 = nb_reaction_4 + 1
2930
2931	indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
2932
2933	!===========================================================
2934	! c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
2935	!===========================================================
2936
2937	nb_reaction_3 = nb_reaction_3 + 1
2938
2939	indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2940
2941	!===========================================================
2942	! c017 : OH + OH + CO2 -> H2O2 + CO2
2943	!===========================================================
2944
2945	nb_reaction_3 = nb_reaction_3 + 1
2946
2947	indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
2948
2949	!===========================================================
2950	! c018 : H + H + CO2 -> H2 + CO2
2951	!===========================================================
2952
2953	nb_reaction_3 = nb_reaction_3 + 1
2954
2955	indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
2956
2957
2958	!===========================================================
2959	! d001 : NO2 + O -> NO + O2
2960	!===========================================================
2961
2962	nb_reaction_4 = nb_reaction_4 + 1
2963
2964	indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
2965
2966	!===========================================================
2967	! d002 : NO + O3 -> NO2 + O2
2968	!===========================================================
2969
2970	nb_reaction_4 = nb_reaction_4 + 1
2971
2972	indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
2973
2974	!===========================================================
2975	! d003 : NO + HO2 -> NO2 + OH
2976	!===========================================================
2977
2978	nb_reaction_4 = nb_reaction_4 + 1
2979
2980	indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
2981
2982	!===========================================================
2983	! d004 : N + NO -> N2 + O
2984	!===========================================================
2985
2986	nb_reaction_4 = nb_reaction_4 + 1
2987
2988	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
2989
2990	!===========================================================
2991	! d005 : N + O2 -> NO + O
2992	!===========================================================
2993
2994	nb_reaction_4 = nb_reaction_4 + 1
2995
2996	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
2997
2998	!===========================================================
2999	! d006 : NO2 + H -> NO + OH
3000	!===========================================================
3001
3002	nb_reaction_4 = nb_reaction_4 + 1
3003
3004	indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
3005
3006	!===========================================================
3007	! d007 : N + O -> NO
3008	!===========================================================
3009
3010	nb_reaction_4 = nb_reaction_4 + 1
3011
3012	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
3013
3014	!===========================================================
3015	! d008 : N + HO2 -> NO + OH
3016	!===========================================================
3017
3018	nb_reaction_4 = nb_reaction_4 + 1
3019
3020	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
3021
3022	!===========================================================
3023	! d009 : N + OH -> NO + H
3024	!===========================================================
3025
3026	nb_reaction_4 = nb_reaction_4 + 1
3027
3028	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
3029
3030	!===========================================================
3031	! d010 : N(2D) + O -> N + O
3032	!===========================================================
3033
3034	nb_phot = nb_phot + 1
3035
3036	indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
3037
3038	!===========================================================
3039	! d011 : N(2D) + N2 -> N + N2
3040	!===========================================================
3041
3042	nb_phot = nb_phot + 1
3043
3044	indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
3045
3046	!===========================================================
3047	! d012 : N(2D) + CO2 -> NO + CO
3048	!===========================================================
3049
3050	nb_reaction_4 = nb_reaction_4 + 1
3051
3052	indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
3053
3054	!===========================================================
3055	! e001 : CO + OH -> CO2 + H
3056	!===========================================================
3057
3058	nb_reaction_4 = nb_reaction_4 + 1
3059
3060	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
3061
3062	!===========================================================
3063	! e002 : CO + O + M -> CO2 + M
3064	!===========================================================
3065
3066	nb_reaction_4 = nb_reaction_4 + 1
3067
3068	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
3069	if(deutchem) then
3070	!===========================================================
3071	! f001 : OD + OH -> HDO + O
3072	!===========================================================
3073
3074	nb_reaction_4 = nb_reaction_4 + 1
3075
3076	indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_oh, 1.0, i_hdo, 1.0, i_o)
3077
3078	!===========================================================
3079	! f002 : OD + H2 -> HDO + H
3080	!===========================================================
3081
3082	nb_reaction_4 = nb_reaction_4 + 1
3083
3084	indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_h2, 1.0, i_hdo, 1.0, i_h)
3085
3086	!===========================================================
3087	! f003 : OD + HO2 -> HDO + O2
3088	!===========================================================
3089
3090	nb_reaction_4 = nb_reaction_4 + 1
3091
3092	indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_ho2, 1.0, i_hdo, 1.0, i_o2)
3093
3094	!===========================================================
3095	! f004 : OD + H2O2 -> HDO + HO2
3096	!===========================================================
3097
3098	nb_reaction_4 = nb_reaction_4 + 1
3099
3100	indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_h2o2, 1.0, i_hdo, 1.0, i_ho2)
3101
3102	!===========================================================
3103	! f005 : O + OD -> O2 + D
3104	!===========================================================
3105
3106	nb_reaction_4 = nb_reaction_4 + 1
3107
3108	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_od, 1.0, i_o2, 1.0, i_d)
3109
3110	!===========================================================
3111	! f006 : OD + H2 -> H2O + D
3112	!===========================================================
3113
3114	nb_reaction_4 = nb_reaction_4 + 1
3115
3116	indice_4(nb_reaction_4) = z4spec(1.0, i_h2, 1.0, i_od, 1.0, i_h2o, 1.0, i_d)
3117
3118	!===========================================================
3119	! f007 : OD + H -> OH + D
3120	!===========================================================
3121
3122	nb_reaction_4 = nb_reaction_4 + 1
3123
3124	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_od, 1.0, i_oh, 1.0, i_d)
3125
3126	!===========================================================
3127	! f008 : CO + OD -> CO2 + D
3128	!===========================================================
3129
3130	nb_reaction_4 = nb_reaction_4 + 1
3131
3132	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_od, 1.0, i_co2, 1.0, i_d)
3133
3134	!===========================================================
3135	! f009 : O3 + D -> O2 + OD
3136	!===========================================================
3137
3138	nb_reaction_4 = nb_reaction_4 + 1
3139
3140	indice_4(nb_reaction_4) = z4spec(1.0, i_o3, 1.0, i_d, 1.0, i_o2, 1.0, i_od)
3141
3142	!===========================================================
3143	! f010 : HO2 + D -> OH + OD
3144	!===========================================================
3145
3146	nb_reaction_4 = nb_reaction_4 + 1
3147
3148	indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_d, 1.0, i_oh, 1.0, i_od)
3149
3150	!===========================================================
3151	! f011 : HO2 + D -> HDO + O
3152	!===========================================================
3153
3154	nb_reaction_4 = nb_reaction_4 + 1
3155
3156	indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_d, 1.0, i_hdo, 1.0, i_o)
3157
3158	!===========================================================
3159	! f012 : OH + D -> H + OD
3160	!===========================================================
3161
3162	nb_reaction_4 = nb_reaction_4 + 1
3163
3164	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_d, 1.0, i_h, 1.0, i_od)
3165
3166	!===========================================================
3167	! f013 : H + D + CO2 -> HD + CO2
3168	!===========================================================
3169
3170	nb_reaction_4 = nb_reaction_4 + 1
3171
3172	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_d, 1.0, i_hd, 0.0, i_dummy)
3173
3174	!===========================================================
3175	! f014 : D + HO2 -> HD + O2
3176	!===========================================================
3177
3178	nb_reaction_4 = nb_reaction_4 + 1
3179
3180	indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_ho2, 1.0, i_hd, 1.0, i_o2)
3181
3182
3183	!===========================================================
3184	! f015 : OH + HD -> HDO + H
3185	!===========================================================
3186
3187	nb_reaction_4 = nb_reaction_4 + 1
3188
3189	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hd, 1.0, i_hdo, 1.0, i_h)
3190
3191
3192	!===========================================================
3193	! f016 : OH + HD -> H2O + D
3194	!===========================================================
3195
3196	nb_reaction_4 = nb_reaction_4 + 1
3197
3198	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hd, 1.0, i_h2o, 1.0, i_d)
3199
3200
3201	!===========================================================
3202	! f017 : D + O2 + CO2 -> DO2 + CO2
3203	!===========================================================
3204
3205	nb_reaction_4 = nb_reaction_4 + 1
3206
3207	indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_o2, 1.0, i_do2, 0.0, i_dummy)
3208
3209
3210	!===========================================================
3211	! f018 : OD + O3 -> DO2 + O2
3212	!===========================================================
3213
3214	nb_reaction_4 = nb_reaction_4 + 1
3215
3216	indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_o3, 1.0, i_do2, 1.0, i_o2)
3217
3218
3219	!===========================================================
3220	! f019 : D + HO2 -> DO2 + H
3221	!===========================================================
3222
3223	nb_reaction_4 = nb_reaction_4 + 1
3224
3225	indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_ho2, 1.0, i_do2, 1.0, i_h)
3226
3227
3228	!===========================================================
3229	! f020 : O + DO2 -> OD + O2
3230	!===========================================================
3231
3232	nb_reaction_4 = nb_reaction_4 + 1
3233
3234	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_do2, 1.0, i_od, 1.0, i_o2)
3235
3236
3237	!===========================================================
3238	! f021 : H + DO2 -> OH + OD
3239	!===========================================================
3240
3241	nb_reaction_4 = nb_reaction_4 + 1
3242
3243	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_od, 1.0, i_oh)
3244
3245
3246	!===========================================================
3247	! f022 : H + DO2 -> HD + O2
3248	!===========================================================
3249
3250	nb_reaction_4 = nb_reaction_4 + 1
3251
3252	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_hd, 1.0, i_o2)
3253
3254
3255	!===========================================================
3256	! f023 : H + DO2 -> HDO + O
3257	!===========================================================
3258
3259	nb_reaction_4 = nb_reaction_4 + 1
3260
3261	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_hdo, 1.0, i_o)
3262
3263
3264	!===========================================================
3265	! f024 : H + DO2 -> HO2 + D
3266	!===========================================================
3267
3268	nb_reaction_4 = nb_reaction_4 + 1
3269
3270	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_ho2, 1.0, i_d)
3271
3272
3273	!===========================================================
3274	! f025 : OH + DO2 -> HDO + O2
3275	!===========================================================
3276
3277	nb_reaction_4 = nb_reaction_4 + 1
3278
3279	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_do2, 1.0, i_hdo, 1.0, i_o2)
3280
3281
3282	!===========================================================
3283	! f026 : DO2 + O3 -> OD + O2 + O2
3284	!===========================================================
3285
3286	nb_reaction_4 = nb_reaction_4 + 1
3287
3288	indice_4(nb_reaction_4) = z4spec(1.0, i_do2, 1.0, i_o3, 1.0, i_od, 2.0, i_o2)
3289
3290
3291	!===========================================================
3292	! f027 : OD + OH + CO2 -> HDO2 + CO2
3293	!===========================================================
3294
3295	nb_reaction_4 = nb_reaction_4 + 1
3296
3297	indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_oh, 1.0, i_hdo2, 0.0, i_dummy)
3298
3299
3300	!===========================================================
3301	! f028 : DO2 + HO2 -> HDO2 + O2
3302	!===========================================================
3303
3304	nb_reaction_4 = nb_reaction_4 + 1
3305
3306	indice_4(nb_reaction_4) = z4spec(1.0, i_do2, 1.0, i_ho2, 1.0, i_hdo2, 1.0, i_o2)
3307
3308	!===========================================================
3309	! f029 : O + HDO2 -> OD + HO2
3310	!===========================================================
3311
3312	nb_reaction_4 = nb_reaction_4 + 1
3313
3314	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_hdo2, 1.0, i_od, 1.0, i_ho2)
3315
3316
3317	!===========================================================
3318	! f030 : O + HDO2 -> OH + DO2
3319	!===========================================================
3320
3321	nb_reaction_4 = nb_reaction_4 + 1
3322
3323	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_hdo2, 1.0, i_oh, 1.0, i_do2)
3324
3325
3326	!===========================================================
3327	! f031 : OH + HDO2 -> HDO + HO2
3328	!===========================================================
3329
3330	nb_reaction_4 = nb_reaction_4 + 1
3331
3332	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hdo2, 1.0, i_hdo, 1.0, i_ho2)
3333
3334
3335	!===========================================================
3336	! f032 : OH + HDO2 -> H2O + DO2
3337	!===========================================================
3338
3339	nb_reaction_4 = nb_reaction_4 + 1
3340
3341	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hdo2, 1.0, i_h2o, 1.0, i_do2)
3342
3343
3344	endif !deutchem
3345
3346	!Only if ion chemistry
3347	if (ionchem) then
3348
3349	!===========================================================
3350	! i001 : CO2+ + O2 -> O2+ + CO2
3351	!===========================================================
3352
3353	nb_reaction_4 = nb_reaction_4 + 1
3354
3355	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
3356
3357	!===========================================================
3358	! i002 : CO2+ + O -> O+ + CO2
3359	!===========================================================
3360
3361	nb_reaction_4 = nb_reaction_4 + 1
3362
3363	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
3364
3365	!===========================================================
3366	! i003 : CO2+ + O -> O2+ + CO
3367	!===========================================================
3368
3369	nb_reaction_4 = nb_reaction_4 + 1
3370
3371	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
3372
3373	!===========================================================
3374	! i004 : O2+ + e- -> O + O
3375	!===========================================================
3376
3377	nb_reaction_4 = nb_reaction_4 + 1
3378
3379	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
3380
3381	!===========================================================
3382	! i005 : O+ + CO2 -> O2+ + CO
3383	!===========================================================
3384
3385	nb_reaction_4 = nb_reaction_4 + 1
3386
3387	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
3388
3389	!===========================================================
3390	! i006 : CO2+ + e -> CO + O
3391	!===========================================================
3392
3393	nb_reaction_4 = nb_reaction_4 + 1
3394
3395	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
3396
3397	!===========================================================
3398	! i007 : CO2+ + NO -> NO+ + CO2
3399	!===========================================================
3400
3401	nb_reaction_4 = nb_reaction_4 + 1
3402
3403	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
3404
3405	!===========================================================
3406	! i008 : O2+ + NO -> NO+ + O2
3407	!===========================================================
3408
3409	nb_reaction_4 = nb_reaction_4 + 1
3410
3411	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
3412
3413	!===========================================================
3414	! i009 : O2+ + N2 -> NO+ + NO
3415	!===========================================================
3416
3417	nb_reaction_4 = nb_reaction_4 + 1
3418
3419	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
3420
3421	!===========================================================
3422	! i010 : O2+ + N -> NO+ + O
3423	!===========================================================
3424
3425	nb_reaction_4 = nb_reaction_4 + 1
3426
3427	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
3428
3429	!===========================================================
3430	! i011 : O+ + N2 -> NO+ + N
3431	!===========================================================
3432
3433	nb_reaction_4 = nb_reaction_4 + 1
3434
3435	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
3436
3437	!===========================================================
3438	! i012 : NO+ + e -> N + O
3439	!===========================================================
3440
3441	nb_reaction_4 = nb_reaction_4 + 1
3442
3443	indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
3444
3445	!===========================================================
3446	! i013 : CO+ + CO2 -> CO2+ + CO
3447	!===========================================================
3448
3449	nb_reaction_4 = nb_reaction_4 + 1
3450
3451	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
3452
3453	!===========================================================
3454	! i014 : CO+ + O -> O+ + CO
3455	!===========================================================
3456
3457	nb_reaction_4 = nb_reaction_4 + 1
3458
3459	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
3460
3461	!===========================================================
3462	! i015 : C+ + CO2 -> CO+ + CO
3463	!===========================================================
3464
3465	nb_reaction_4 = nb_reaction_4 + 1
3466
3467	indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
3468
3469	!===========================================================
3470	! i016 : N2+ + CO2 -> CO2+ + N2
3471	!===========================================================
3472
3473	nb_reaction_4 = nb_reaction_4 + 1
3474
3475	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
3476
3477	!===========================================================
3478	! i017 : N2+ + O -> NO+ + N
3479	!===========================================================
3480
3481	nb_reaction_4 = nb_reaction_4 + 1
3482
3483	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
3484
3485	!===========================================================
3486	! i018 : N2+ + CO -> CO+ + N2
3487	!===========================================================
3488
3489	nb_reaction_4 = nb_reaction_4 + 1
3490
3491	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
3492
3493	!===========================================================
3494	! i019 : N2+ + e -> N + N
3495	!===========================================================
3496
3497	nb_reaction_4 = nb_reaction_4 + 1
3498
3499	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
3500
3501	!===========================================================
3502	! i020 : N2+ + O -> O+ + N2
3503	!===========================================================
3504
3505	nb_reaction_4 = nb_reaction_4 + 1
3506
3507	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
3508
3509	!===========================================================
3510	! i021 : N+ + CO2 -> CO2+ + N
3511	!===========================================================
3512
3513	nb_reaction_4 = nb_reaction_4 + 1
3514
3515	indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
3516
3517	!===========================================================
3518	! i022 : CO+ + H -> H+ + CO
3519	!===========================================================
3520
3521	nb_reaction_4 = nb_reaction_4 + 1
3522
3523	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
3524
3525	!===========================================================
3526	! i023 : O+ + H -> H+ + O
3527	!===========================================================
3528
3529	nb_reaction_4 = nb_reaction_4 + 1
3530
3531	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
3532
3533	!===========================================================
3534	! i024 : H+ + O -> O+ + H
3535	!===========================================================
3536
3537	nb_reaction_4 = nb_reaction_4 + 1
3538
3539	indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
3540
3541	!===========================================================
3542	! i025 : CO2+ + H2 -> HCO2+ + H
3543	!===========================================================
3544
3545	nb_reaction_4 = nb_reaction_4 + 1
3546
3547	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
3548
3549	!===========================================================
3550	! i026 : HCO2+ + e -> H + CO2
3551	!===========================================================
3552
3553	nb_reaction_4 = nb_reaction_4 + 1
3554
3555	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
3556
3557	!===========================================================
3558	! i027 : HCO2+ + e -> H + O + CO
3559	!===========================================================
3560	!We divide this reaction in two
3561
3562	!0.5HCO2+ + 0.5e -> H
3563
3564	nb_reaction_4 = nb_reaction_4 + 1
3565
3566	indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
3567
3568	!0.5 HCO2+ + 0.5 e -> O + CO
3569
3570	nb_reaction_4 = nb_reaction_4 + 1
3571
3572	indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
3573
3574	!===========================================================
3575	! i029 : HCO2+ + e -> OH + CO
3576	!===========================================================
3577
3578	nb_reaction_4 = nb_reaction_4 + 1
3579
3580	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
3581
3582
3583	!===========================================================
3584	! i030 : HCO2+ + e -> H + CO2
3585	!===========================================================
3586
3587	nb_reaction_4 = nb_reaction_4 + 1
3588
3589	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
3590
3591
3592	!===========================================================
3593	! i031 : HCO2+ + O -> HCO+ + O2
3594	!===========================================================
3595
3596	nb_reaction_4 = nb_reaction_4 + 1
3597
3598	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_o, 1.0, i_hcoplus, 1.0, i_o2)
3599
3600
3601	!===========================================================
3602	! i032 : HCO2+ + CO -> HCO+ + CO2
3603	!===========================================================
3604
3605	nb_reaction_4 = nb_reaction_4 + 1
3606	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_co2)
3607
3608
3609	!===========================================================
3610	! i033 : H+ + CO2 -> HCO+ + O
3611	!===========================================================
3612
3613	nb_reaction_4 = nb_reaction_4 + 1
3614	indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_co2, 1.0, i_hcoplus, 1.0, i_o)
3615
3616
3617	!===========================================================
3618	! i034 : CO2+ + H -> HCO+ + O
3619	!===========================================================
3620
3621	nb_reaction_4 = nb_reaction_4 + 1
3622	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h, 1.0, i_hcoplus, 1.0, i_o)
3623
3624
3625	!===========================================================
3626	! i035 : CO+ + H2 -> HCO+ + H
3627	!===========================================================
3628
3629	nb_reaction_4 = nb_reaction_4 + 1
3630	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2, 1.0, i_hcoplus, 1.0, i_h)
3631
3632
3633	!===========================================================
3634	! i036 : HCO+ + e- -> CO + H
3635	!===========================================================
3636
3637	nb_reaction_4 = nb_reaction_4 + 1
3638	indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_elec, 1.0, i_co, 1.0, i_h)
3639
3640	!===========================================================
3641	! i037 : CO2+ + H2O -> H2O+ + CO2
3642	!===========================================================
3643
3644	nb_reaction_4 = nb_reaction_4 + 1
3645	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co2)
3646
3647	!===========================================================
3648	! i038 : CO+ + H2O -> H2O+ + CO
3649	!===========================================================
3650
3651	nb_reaction_4 = nb_reaction_4 + 1
3652	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co)
3653
3654	!===========================================================
3655	! i039 : O+ + H2O -> H2O+ + O
3656	!===========================================================
3657
3658	nb_reaction_4 = nb_reaction_4 + 1
3659	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_o)
3660
3661	!===========================================================
3662	! i040 : N2+ + H2O -> H2O+ + N2
3663	!===========================================================
3664
3665	nb_reaction_4 = nb_reaction_4 + 1
3666	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n2)
3667
3668	!===========================================================
3669	! i041 : N+ + H2O -> H2O+ + N
3670	!===========================================================
3671
3672	nb_reaction_4 = nb_reaction_4 + 1
3673	indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n)
3674
3675	!===========================================================
3676	! i042 : H+ + H2O -> H2O+ + H
3677	!===========================================================
3678
3679	nb_reaction_4 = nb_reaction_4 + 1
3680	indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_h)
3681
3682	!===========================================================
3683	! i043 : H2O+ + O2 -> O2+ + H2O
3684	!===========================================================
3685
3686	nb_reaction_4 = nb_reaction_4 + 1
3687	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_h2o)
3688
3689	!===========================================================
3690	! i044 : H2O+ + CO -> HCO+ + OH
3691	!===========================================================
3692
3693	nb_reaction_4 = nb_reaction_4 + 1
3694	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_oh)
3695
3696	!===========================================================
3697	! i045 : H2O+ + O -> O2+ + H2
3698	!===========================================================
3699
3700	nb_reaction_4 = nb_reaction_4 + 1
3701	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h2)
3702
3703	!===========================================================
3704	! i046 : H2O+ + NO -> NO+ + H2O
3705	!===========================================================
3706
3707	nb_reaction_4 = nb_reaction_4 + 1
3708	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_h2o)
3709
3710	!===========================================================
3711	! i047 : H2O+ + e- -> H + H + O
3712	!===========================================================
3713
3714	nb_reaction_4 = nb_reaction_4 + 1
3715	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 2.0, i_h, 1.0, i_o)
3716
3717	!===========================================================
3718	! i048 : H2O+ + e- -> H + OH
3719	!===========================================================
3720
3721	nb_reaction_4 = nb_reaction_4 + 1
3722	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h, 1.0, i_oh)
3723
3724	!===========================================================
3725	! i049 : H2O+ + e- -> H2 + O
3726	!===========================================================
3727
3728	nb_reaction_4 = nb_reaction_4 + 1
3729	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o)
3730
3731	!===========================================================
3732	! i050 : H2O+ + H2O -> H3O+ + OH
3733	!===========================================================
3734
3735	nb_reaction_4 = nb_reaction_4 + 1
3736	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_oh)
3737
3738	!===========================================================
3739	! i051 : H2O+ + H2 -> H3O+ + H
3740	!===========================================================
3741
3742	nb_reaction_4 = nb_reaction_4 + 1
3743	indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2, 1.0, i_h3oplus, 1.0, i_h)
3744
3745	!===========================================================
3746	! i052 : HCO+ + H2O -> H3O+ + CO
3747	!===========================================================
3748
3749	nb_reaction_4 = nb_reaction_4 + 1
3750	indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_co)
3751
3752	!===========================================================
3753	! i053: H3O+ + e -> OH + H + H
3754	!===========================================================
3755
3756	nb_reaction_4 = nb_reaction_4 + 1
3757	indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 2.0, i_h)
3758
3759	!===========================================================
3760	! i054: H3O+ + e -> H2O + H
3761	!===========================================================
3762
3763	nb_reaction_4 = nb_reaction_4 + 1
3764	indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_h2o, 1.0, i_h)
3765
3766	!===========================================================
3767	! i055: H3O+ + e -> HO + H2
3768	!===========================================================
3769
3770	nb_reaction_4 = nb_reaction_4 + 1
3771	indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 1.0, i_h2)
3772
3773	!===========================================================
3774	! i056: H3O+ + e -> O + H2 + H
3775	!===========================================================
3776	!We divide this reaction in two
3777
3778	!0.5H3O+ + 0.5e -> O
3779
3780	nb_reaction_4 = nb_reaction_4 + 1
3781	indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_o, 0.0, i_dummy)
3782
3783	!0.5H3O+ + 0.5e -> H2 + H
3784
3785	nb_reaction_4 = nb_reaction_4 + 1
3786	indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_h2, 1.0, i_h)
3787
3788	!===========================================================
3789	! i057: O+ + H2 -> OH+ + H
3790	!===========================================================
3791
3792	nb_reaction_4 = nb_reaction_4 + 1
3793	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2, 1.0, i_ohplus, 1.0, i_h)
3794
3795	!===========================================================
3796	! i058: OH+ + O -> O2+ + H
3797	!===========================================================
3798
3799	nb_reaction_4 = nb_reaction_4 + 1
3800	indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h)
3801
3802	!===========================================================
3803	! i059: OH+ + CO2 -> HCO2+ + O
3804	!===========================================================
3805
3806	nb_reaction_4 = nb_reaction_4 + 1
3807	indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co2, 1.0, i_hco2plus, 1.0, i_o)
3808
3809	!===========================================================
3810	! i060: OH+ + CO -> HCO+ + O
3811	!===========================================================
3812
3813	nb_reaction_4 = nb_reaction_4 + 1
3814	indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_o)
3815
3816	!===========================================================
3817	! i061: OH+ + NO -> NO+ + OH
3818	!===========================================================
3819
3820	nb_reaction_4 = nb_reaction_4 + 1
3821	indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_oh)
3822
3823	!===========================================================
3824	! i062: OH+ + H2 -> H2O+ + H
3825	!===========================================================
3826
3827	nb_reaction_4 = nb_reaction_4 + 1
3828	indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_h2, 1.0, i_h2oplus, 1.0, i_h)
3829
3830	!===========================================================
3831	! i063: OH+ + O2 -> O2+ + OH
3832	!===========================================================
3833
3834	nb_reaction_4 = nb_reaction_4 + 1
3835	indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_oh)
3836
3837	end if !ionchem
3838
3839	!===========================================================
3840	! h001: HO2 + ice -> products
3841	! treated as
3842	! HO2 -> 0.5 H2O + 0.75 O2
3843	!===========================================================
3844
3845	nb_phot = nb_phot + 1
3846
3847	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
3848
3849	!===========================================================
3850	! h002: OH + ice -> products
3851	! treated as
3852	! OH -> 0.5 H2O + 0.25 O2
3853	!===========================================================
3854
3855	nb_phot = nb_phot + 1
3856
3857	indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2)
3858
3859	!===========================================================
3860	! h003: H2O2 + ice -> products
3861	! treated as
3862	! H2O2 -> H2O + 0.5 O2
3863	!===========================================================
3864
3865	nb_phot = nb_phot + 1
3866
3867	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
3868
3869	!===========================================================
3870	! h004: HO2 + dust -> products
3871	! treated as
3872	! HO2 -> 0.5 H2O + 0.75 O2
3873	!===========================================================
3874
3875	nb_phot = nb_phot + 1
3876
3877	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
3878
3879	!===========================================================
3880	! h005: H2O2 + dust -> products
3881	! treated as
3882	! H2O2 -> H2O + 0.5 O2
3883	!===========================================================
3884
3885	nb_phot = nb_phot + 1
3886
3887	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
3888
3889	!===========================================================
3890	! check dimensions
3891	!===========================================================
3892
3893	print*, 'nb_phot = ', nb_phot
3894	print*, 'nb_reaction_4 = ', nb_reaction_4
3895	print*, 'nb_reaction_3 = ', nb_reaction_3
3896
3897	if ((nb_phot /= nb_phot_max) .or. &
3898	(nb_reaction_3 /= nb_reaction_3_max) .or. &
3899	(nb_reaction_4 /= nb_reaction_4_max)) then
3900	print*, 'wrong dimensions in indice'
3901	call abort_physic("indice","wrong array dimensions",1)
3902	end if
3903
3904	end subroutine indice
3905
3906	!*****************************************************************
3907
3908	subroutine gcmtochim(nlayer, ionchem, deutchem, nq, zycol, &
3909	lswitch, nesp, &
3910	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
3911	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
3912	i_n, i_n2d, i_no, i_no2, i_n2, &
3913	i_co2plus, i_oplus, i_o2plus, i_noplus, &
3914	i_coplus, i_cplus, i_n2plus, i_nplus, &
3915	i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,&
3916	i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, &
3917	i_d, i_hd, i_do2, i_hdo2, dens, rm, c)
3918
3919	!*****************************************************************
3920
3921	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, &
3922	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh, &
3923	& igcm_ho2, igcm_h2o2, igcm_h2o_vap, &
3924	& igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2,&
3925	& igcm_co2plus, igcm_oplus, igcm_o2plus, &
3926	& igcm_noplus, igcm_coplus, igcm_cplus, &
3927	& igcm_n2plus, igcm_nplus, igcm_hplus, &
3928	& igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, &
3929	& igcm_h3oplus, igcm_ohplus, igcm_elec, &
3930	& igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, &
3931	& igcm_do2, igcm_hdo2
3932
3933	implicit none
3934
3935	include "callkeys.h"
3936
3937	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3938	! input:
3939	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3940
3941	integer, intent(in) :: nlayer ! number of atmospheric layers
3942	integer, intent(in) :: nq ! number of tracers in the gcm
3943	logical, intent(in) :: ionchem
3944	logical, intent(in) :: deutchem
3945	integer :: nesp ! number of species in the chemistry
3946	integer :: lswitch ! interface level between chemistries
3947
3948	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
3949	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
3950	i_n, i_n2d, i_no, i_no2, i_n2, &
3951	i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
3952	i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, &
3953	i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec, &
3954	i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2
3955
3956	real :: zycol(nlayer,nq) ! volume mixing ratios in the gcm
3957	real :: dens(nlayer) ! total number density (molecule.cm-3)
3958
3959	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3960	! output:
3961	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3962
3963	real, dimension(nlayer,nesp) :: rm ! volume mixing ratios
3964	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
3965
3966	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3967	! local:
3968	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3969
3970	integer :: l, iesp
3971	logical,save :: firstcall = .true.
3972
3973	!$OMP THREADPRIVATE(firstcall)
3974
3975	! first call initializations
3976
3977	if (firstcall) then
3978
3979	! identify the indexes of the tracers we need
3980
3981	if (igcm_co2 == 0) then
3982	write(,) "gcmtochim: Error; no CO2 tracer !!!"
3983	call abort_physic("gcmtochim","missing co2 tracer",1)
3984	endif
3985	if (igcm_co == 0) then
3986	write(,) "gcmtochim: Error; no CO tracer !!!"
3987	call abort_physic("gcmtochim","missing co tracer",1)
3988	end if
3989	if (igcm_o == 0) then
3990	write(,) "gcmtochim: Error; no O tracer !!!"
3991	call abort_physic("gcmtochim","missing o tracer",1)
3992	end if
3993	if (igcm_o1d == 0) then
3994	write(,) "gcmtochim: Error; no O1D tracer !!!"
3995	call abort_physic("gcmtochim","missing o1d tracer",1)
3996	end if
3997	if (igcm_o2 == 0) then
3998	write(,) "gcmtochim: Error; no O2 tracer !!!"
3999	call abort_physic("gcmtochim","missing o2 tracer",1)
4000	end if
4001	if (igcm_o3 == 0) then
4002	write(,) "gcmtochim: Error; no O3 tracer !!!"
4003	call abort_physic("gcmtochim","missing o3 tracer",1)
4004	end if
4005	if (igcm_h == 0) then
4006	write(,) "gcmtochim: Error; no H tracer !!!"
4007	call abort_physic("gcmtochim","missing h tracer",1)
4008	end if
4009	if (igcm_h2 == 0) then
4010	write(,) "gcmtochim: Error; no H2 tracer !!!"
4011	call abort_physic("gcmtochim","missing h2 tracer",1)
4012	end if
4013	if (igcm_oh == 0) then
4014	write(,) "gcmtochim: Error; no OH tracer !!!"
4015	call abort_physic("gcmtochim","missing oh tracer",1)
4016	end if
4017	if (igcm_ho2 == 0) then
4018	write(,) "gcmtochim: Error; no HO2 tracer !!!"
4019	call abort_physic("gcmtochim","missing ho2 tracer",1)
4020	end if
4021	if (igcm_h2o2 == 0) then
4022	write(,) "gcmtochim: Error; no H2O2 tracer !!!"
4023	call abort_physic("gcmtochim","missing h2o2 tracer",1)
4024	end if
4025	if (igcm_n == 0) then
4026	write(,) "gcmtochim: Error; no N tracer !!!"
4027	call abort_physic("gcmtochim","missing n tracer",1)
4028	end if
4029	if (igcm_n2d == 0) then
4030	write(,) "gcmtochim: Error; no N2D tracer !!!"
4031	call abort_physic("gcmtochim","missing n2d tracer",1)
4032	end if
4033	if (igcm_no == 0) then
4034	write(,) "gcmtochim: Error; no NO tracer !!!"
4035	call abort_physic("gcmtochim","missing no tracer",1)
4036	end if
4037	if (igcm_no2 == 0) then
4038	write(,) "gcmtochim: Error; no NO2 tracer !!!"
4039	call abort_physic("gcmtochim","missing no2 tracer",1)
4040	end if
4041	if (igcm_n2 == 0) then
4042	write(,) "gcmtochim: Error; no N2 tracer !!!"
4043	call abort_physic("gcmtochim","missing n2 tracer",1)
4044	end if
4045	if (igcm_h2o_vap == 0) then
4046	write(,) "gcmtochim: Error; no water vapor tracer !!!"
4047	call abort_physic("gcmtochim","missing h2o_vap tracer",1)
4048	end if
4049	if(deutchem) then
4050	if (igcm_hdo_vap == 0) then
4051	write(,) "gcmtochim: Error; no HDO tracer !!!"
4052	call abort_physic("gcmtochim","missing hdo_vap tracer",1)
4053	end if
4054	if (igcm_od == 0) then
4055	write(,) "gcmtochim: Error; no OD tracer !!!"
4056	call abort_physic("gcmtochim","missing od tracer",1)
4057	end if
4058	if (igcm_d == 0) then
4059	write(,) "gcmtochim: Error; no D tracer !!!"
4060	call abort_physic("gcmtochim","missing d tracer",1)
4061	end if
4062	if (igcm_hd == 0) then
4063	write(,) "gcmtochim: Error; no HD tracer !!!"
4064	call abort_physic("gcmtochim","missing hd tracer",1)
4065	end if
4066	if (igcm_do2 == 0) then
4067	write(,) "gcmtochim: Error; no DO2 tracer !!!"
4068	call abort_physic("gcmtochim","missing do2 tracer",1)
4069	end if
4070	if (igcm_hdo2 == 0) then
4071	write(,) "gcmtochim: Error; no HDO2 tracer !!!"
4072	call abort_physic("gcmtochim","missing hdo2 tracer",1)
4073	end if
4074	endif !deutchem
4075	if (ionchem) then
4076	if (igcm_co2plus == 0) then
4077	write(,) "gcmtochim: Error; no CO2+ tracer !!!"
4078	call abort_physic("gcmtochim","missing co2plus tracer",1)
4079	end if
4080	if (igcm_oplus == 0) then
4081	write(,) "gcmtochim: Error; no O+ tracer !!!"
4082	call abort_physic("gcmtochim","missing oplus tracer",1)
4083	end if
4084	if (igcm_o2plus == 0) then
4085	write(,) "gcmtochim: Error; no O2+ tracer !!!"
4086	call abort_physic("gcmtochim","missing o2plus tracer",1)
4087	end if
4088	if (igcm_noplus == 0) then
4089	write(,) "gcmtochim: Error; no NO+ tracer !!!"
4090	call abort_physic("gcmtochim","missing noplus tracer",1)
4091	endif
4092	if (igcm_coplus == 0) then
4093	write(,) "gcmtochim: Error; no CO+ tracer !!!"
4094	call abort_physic("gcmtochim","missing coplus tracer",1)
4095	endif
4096	if (igcm_cplus == 0) then
4097	write(,) "gcmtochim: Error; no C+ tracer !!!"
4098	call abort_physic("gcmtochim","missing cplus tracer",1)
4099	endif
4100	if (igcm_n2plus == 0) then
4101	write(,) "gcmtochim: Error; no N2+ tracer !!!"
4102	call abort_physic("gcmtochim","missing n2plus tracer",1)
4103	endif
4104	if (igcm_nplus == 0) then
4105	write(,) "gcmtochim: Error; no N+ tracer !!!"
4106	call abort_physic("gcmtochim","missing nplus tracer",1)
4107	endif
4108	if (igcm_hplus == 0) then
4109	write(,) "gcmtochim: Error; no H+ tracer !!!"
4110	call abort_physic("gcmtochim","missing hplus tracer",1)
4111	endif
4112	if (igcm_hco2plus == 0) then
4113	write(,) "gcmtochim: Error; no HCO2+ tracer !!!"
4114	call abort_physic("gcmtochim","missing hco2plus tracer",1)
4115	endif
4116	if (igcm_hcoplus == 0) then
4117	write(,) "gcmtochim: Error; no HCO+ tracer !!!"
4118	call abort_physic("gcmtochim","missing hcoplus tracer",1)
4119	endif
4120	if (igcm_h2oplus == 0) then
4121	write(,) "gcmtochim: Error; no H2O+ tracer !!!"
4122	call abort_physic("gcmtochim","missing h2oplus tracer",1)
4123	endif
4124	if (igcm_h3oplus == 0) then
4125	write(,) "gcmtochim: Error; no H3O+ tracer !!!"
4126	call abort_physic("gcmtochim","missing h3oplus tracer",1)
4127	endif
4128	if (igcm_ohplus == 0) then
4129	write(,) "gcmtochim: Error; no OH+ tracer !!!"
4130	call abort_physic("gcmtochim","missing ohplus tracer",1)
4131	endif
4132	if (igcm_elec == 0) then
4133	write(,) "gcmtochim: Error; no e- tracer !!!"
4134	call abort_physic("gcmtochim","missing elec tracer",1)
4135	end if
4136	end if ! ionchem
4137	firstcall = .false.
4138	end if ! of if (firstcall)
4139
4140	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4141	! initialise mixing ratios
4142	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4143
4144	do l = 1,nlayer
4145	rm(l,i_co2) = zycol(l, igcm_co2)
4146	rm(l,i_co) = zycol(l, igcm_co)
4147	rm(l,i_o) = zycol(l, igcm_o)
4148	rm(l,i_o1d) = zycol(l, igcm_o1d)
4149	rm(l,i_o2) = zycol(l, igcm_o2)
4150	rm(l,i_o3) = zycol(l, igcm_o3)
4151	rm(l,i_h) = zycol(l, igcm_h)
4152	rm(l,i_h2) = zycol(l, igcm_h2)
4153	rm(l,i_oh) = zycol(l, igcm_oh)
4154	rm(l,i_ho2) = zycol(l, igcm_ho2)
4155	rm(l,i_h2o2) = zycol(l, igcm_h2o2)
4156	rm(l,i_h2o) = zycol(l, igcm_h2o_vap)
4157	rm(l,i_n) = zycol(l, igcm_n)
4158	rm(l,i_n2d) = zycol(l, igcm_n2d)
4159	rm(l,i_no) = zycol(l, igcm_no)
4160	rm(l,i_no2) = zycol(l, igcm_no2)
4161	rm(l,i_n2) = zycol(l, igcm_n2)
4162	enddo
4163
4164	if (deutchem) then
4165	do l = 1,nlayer
4166	rm(l,i_hdo) = zycol(l, igcm_hdo_vap)
4167	rm(l,i_od) = zycol(l, igcm_od)
4168	rm(l,i_d) = zycol(l, igcm_d)
4169	rm(l,i_hd) = zycol(l, igcm_hd)
4170	rm(l,i_do2) = zycol(l, igcm_do2)
4171	rm(l,i_hdo2) = zycol(l, igcm_hdo2)
4172	end do
4173	endif
4174
4175	if (ionchem) then
4176	do l = 1,nlayer
4177	rm(l,i_co2plus) = zycol(l, igcm_co2plus)
4178	rm(l,i_oplus) = zycol(l, igcm_oplus)
4179	rm(l,i_o2plus) = zycol(l, igcm_o2plus)
4180	rm(l,i_noplus) = zycol(l, igcm_noplus)
4181	rm(l,i_coplus) = zycol(l, igcm_coplus)
4182	rm(l,i_cplus) = zycol(l, igcm_cplus)
4183	rm(l,i_n2plus) = zycol(l, igcm_n2plus)
4184	rm(l,i_nplus) = zycol(l, igcm_nplus)
4185	rm(l,i_hplus) = zycol(l, igcm_hplus)
4186	rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
4187	rm(l,i_hcoplus) = zycol(l, igcm_hcoplus)
4188	rm(l,i_h2oplus) = zycol(l, igcm_h2oplus)
4189	rm(l,i_h3oplus) = zycol(l, igcm_h3oplus)
4190	rm(l,i_ohplus) = zycol(l, igcm_ohplus)
4191	rm(l,i_elec) = zycol(l, igcm_elec)
4192	end do
4193	end if
4194
4195	where (rm(:,:) < 1.e-30)
4196	rm(:,:) = 0.
4197	end where
4198
4199	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4200	! initialise number densities
4201	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4202
4203	do iesp = 1,nesp
4204	do l = 1,nlayer
4205	c(l,iesp) = rm(l,iesp)*dens(l)
4206	end do
4207	end do
4208
4209	end subroutine gcmtochim
4210
4211	!*****************************************************************
4212
4213	subroutine chimtogcm(nlayer, ionchem, deutchem, nq, zycol, &
4214	lswitch, nesp, &
4215	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
4216	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
4217	i_n, i_n2d, i_no, i_no2, i_n2, &
4218	i_co2plus, i_oplus, i_o2plus, i_noplus, &
4219	i_coplus, i_cplus, i_n2plus, i_nplus, &
4220	i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, &
4221	i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, &
4222	i_d, i_hd, i_do2, i_hdo2, dens, c)
4223
4224	!*****************************************************************
4225
4226	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, &
4227	igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh, &
4228	igcm_ho2, igcm_h2o2, igcm_h2o_vap, &
4229	igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2, &
4230	igcm_co2plus, igcm_oplus, igcm_o2plus, &
4231	igcm_noplus, igcm_coplus, igcm_cplus, &
4232	igcm_n2plus, igcm_nplus, igcm_hplus, &
4233	igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, &
4234	igcm_h3oplus, igcm_ohplus, igcm_elec, &
4235	igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, &
4236	igcm_do2, igcm_hdo2
4237
4238	implicit none
4239
4240	#include "callkeys.h"
4241
4242	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4243	! inputs:
4244	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4245
4246	integer, intent(in) :: nlayer ! number of atmospheric layers
4247	integer, intent(in) :: nq ! number of tracers in the gcm
4248	logical, intent(in) :: ionchem
4249	logical, intent(in) :: deutchem
4250	integer :: nesp ! number of species in the chemistry
4251	integer :: lswitch ! interface level between chemistries
4252	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
4253	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
4254	i_n, i_n2d, i_no, i_no2, i_n2, &
4255	i_co2plus, i_oplus, i_o2plus, i_noplus, &
4256	i_coplus, i_cplus, i_n2plus, i_nplus, &
4257	i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, &
4258	i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, i_d, &
4259	i_hd, i_do2, i_hdo2
4260
4261	real :: dens(nlayer) ! total number density (molecule.cm-3)
4262	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
4263
4264	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4265	! output:
4266	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4267
4268	real zycol(nlayer,nq) ! volume mixing ratios in the gcm
4269
4270	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4271	! local:
4272	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4273
4274	integer l
4275
4276	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4277	! save mixing ratios for the gcm
4278	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4279
4280	do l = 1,lswitch-1
4281	zycol(l, igcm_co2) = c(l,i_co2)/dens(l)
4282	zycol(l, igcm_co) = c(l,i_co)/dens(l)
4283	zycol(l, igcm_o) = c(l,i_o)/dens(l)
4284	zycol(l, igcm_o1d) = c(l,i_o1d)/dens(l)
4285	zycol(l, igcm_o2) = c(l,i_o2)/dens(l)
4286	zycol(l, igcm_o3) = c(l,i_o3)/dens(l)
4287	zycol(l, igcm_h) = c(l,i_h)/dens(l)
4288	zycol(l, igcm_h2) = c(l,i_h2)/dens(l)
4289	zycol(l, igcm_oh) = c(l,i_oh)/dens(l)
4290	zycol(l, igcm_ho2) = c(l,i_ho2)/dens(l)
4291	zycol(l, igcm_h2o2) = c(l,i_h2o2)/dens(l)
4292	zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
4293	zycol(l, igcm_n) = c(l,i_n)/dens(l)
4294	zycol(l, igcm_n2d) = c(l,i_n2d)/dens(l)
4295	zycol(l, igcm_no) = c(l,i_no)/dens(l)
4296	zycol(l, igcm_no2) = c(l,i_no2)/dens(l)
4297	zycol(l, igcm_n2) = c(l,i_n2)/dens(l)
4298	enddo
4299
4300	if(deutchem) then
4301	do l=1,lswitch-1
4302	zycol(l, igcm_hdo_vap) = c(l,i_hdo)/dens(l)
4303	zycol(l, igcm_od) = c(l,i_od)/dens(l)
4304	zycol(l, igcm_d) = c(l,i_d)/dens(l)
4305	zycol(l, igcm_hd) = c(l,i_hd)/dens(l)
4306	zycol(l, igcm_do2) = c(l,i_do2)/dens(l)
4307	zycol(l, igcm_hdo2) = c(l,i_hdo2)/dens(l)
4308	end do
4309	endif !deutchem
4310
4311	if (ionchem) then
4312	do l = 1,lswitch-1
4313	zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
4314	zycol(l, igcm_oplus) = c(l,i_oplus)/dens(l)
4315	zycol(l, igcm_o2plus) = c(l,i_o2plus)/dens(l)
4316	zycol(l, igcm_noplus) = c(l,i_noplus)/dens(l)
4317	zycol(l, igcm_coplus) = c(l,i_coplus)/dens(l)
4318	zycol(l, igcm_cplus) = c(l,i_cplus)/dens(l)
4319	zycol(l, igcm_n2plus) = c(l,i_n2plus)/dens(l)
4320	zycol(l, igcm_nplus) = c(l,i_nplus)/dens(l)
4321	zycol(l, igcm_hplus) = c(l,i_hplus)/dens(l)
4322	zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
4323	zycol(l, igcm_hcoplus) = c(l,i_hcoplus)/dens(l)
4324	zycol(l, igcm_h2oplus) = c(l,i_h2oplus)/dens(l)
4325	zycol(l, igcm_h3oplus) = c(l,i_h3oplus)/dens(l)
4326	zycol(l, igcm_ohplus) = c(l,i_ohplus)/dens(l)
4327	zycol(l, igcm_elec) = c(l,i_elec)/dens(l)
4328	end do
4329	end if
4330
4331	end subroutine chimtogcm
4332
4333	end subroutine photochemistry

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

Download in other formats:

Original Format