Context Navigation

photochemistry.F90 @ 2236

Last change on this file since 2236 was 2170, checked in by flefevre, 5 years ago
la chimie des ions est désormais optionnelle (FGG) regroupement des options de la photochimie dans calchim_mod (FL) mise à jour cinétique CO + OH, inactive pour le moment (FL)
File size: 87.1 KB

Line
1	!****************************************************************
2	!
3	! Photochemical routine
4	!
5	! Author: Franck Lefevre
6	! ------
7	!
8	! Version: 27/04/2017
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, nb_reaction_3_max, &
16	nb_reaction_4_max, nb_phot_max, nphotion, &
17	jonline, ig, lswitch, zycol, sza, ptimestep, press, &
18	alt, temp, temp_elect, dens, zmmean, &
19	dist_sol, zday, &
20	surfdust1d, surfice1d, jo3, jh2o,tau, iter)
21
22	use param_v4_h, only: jion
23
24	implicit none
25
26	#include "callkeys.h"
27
28	!===================================================================
29	! inputs:
30	!===================================================================
31
32	integer, intent(in) :: nlayer ! number of atmospheric layers
33	integer, intent(in) :: nq ! number of tracers in traceur.def
34	integer, intent(in) :: nesp ! number of traceurs in chemistry
35	integer, intent(in) :: nb_reaction_3_max
36	! number of quadratic reactions
37	integer, intent(in) :: nb_reaction_4_max
38	! number of bimolecular reactions
39	integer, intent(in) :: nb_phot_max
40	! number of reactions treated numerically as photodissociations
41	integer, intent(in) :: nphotion
42	! number of photoionizations
43	logical, intent(in) :: ionchem! switch for ion chemistry
44	logical, intent(in) :: jonline! switch for on-line calculation of photolysis rates
45	integer :: ig ! grid point index
46
47	real :: sza ! solar zenith angle (deg)
48	real :: ptimestep ! physics timestep (s)
49	real :: press(nlayer) ! pressure (hpa)
50	real :: alt(nlayer) ! altitude (km)
51	real :: temp(nlayer) ! temperature (k)
52	real :: temp_elect(nlayer) ! electronic temperature (K)
53	real :: dens(nlayer) ! density (cm-3)
54	real :: zmmean(nlayer) ! mean molar mass (g/mole)
55	real :: dist_sol ! sun distance (au)
56	real :: zday ! date (time since Ls=0, in martian days)
57	real :: surfdust1d(nlayer) ! dust surface area (cm2/cm3)
58	real :: surfice1d(nlayer) ! ice surface area (cm2/cm3)
59	real :: tau ! dust optical depth
60
61	!===================================================================
62	! input/output:
63	!===================================================================
64
65	real :: zycol(nlayer,nq) ! chemical species volume mixing ratio
66
67	!===================================================================
68	! output:
69	!===================================================================
70
71	integer :: iter(nlayer) ! iteration counter
72	real :: jo3(nlayer) ! photodissociation rate o3 -> o1d
73	real :: jh2o(nlayer) ! photodissociation rate h2o -> h + oh
74
75	!===================================================================
76	! local:
77	!===================================================================
78
79	integer :: phychemrat ! (physical timestep)/(nominal chemical timestep)
80	integer :: j_o3_o1d, j_h2o, ilev, iesp
81	integer :: lswitch
82	logical, save :: firstcall = .true.
83	logical :: jionos ! switch for J parameterization
84
85	! tracer indexes in the chemistry:
86
87	integer,parameter :: i_co2 = 1
88	integer,parameter :: i_co = 2
89	integer,parameter :: i_o = 3
90	integer,parameter :: i_o1d = 4
91	integer,parameter :: i_o2 = 5
92	integer,parameter :: i_o3 = 6
93	integer,parameter :: i_h = 7
94	integer,parameter :: i_h2 = 8
95	integer,parameter :: i_oh = 9
96	integer,parameter :: i_ho2 = 10
97	integer,parameter :: i_h2o2 = 11
98	integer,parameter :: i_h2o = 12
99	integer,parameter :: i_n = 13
100	integer,parameter :: i_n2d = 14
101	integer,parameter :: i_no = 15
102	integer,parameter :: i_no2 = 16
103	integer,parameter :: i_n2 = 17
104	integer,parameter :: i_co2plus = 18
105	integer,parameter :: i_oplus = 19
106	integer,parameter :: i_o2plus = 20
107	integer,parameter :: i_noplus = 21
108	integer,parameter :: i_coplus = 22
109	integer,parameter :: i_cplus = 23
110	integer,parameter :: i_n2plus = 24
111	integer,parameter :: i_nplus = 25
112	integer,parameter :: i_hplus = 26
113	integer,parameter :: i_hco2plus= 27
114	integer,parameter :: i_elec = 28
115
116	integer :: ilay
117
118	real :: ctimestep ! standard timestep for the chemistry (s)
119	real :: dt_guess ! first-guess timestep (s)
120	real :: dt_corrected ! corrected timestep (s)
121	real :: time ! internal time (between 0 and ptimestep, in s)
122
123	real, dimension(nlayer,nesp) :: rm ! mixing ratios
124	real (kind = 8), dimension(nesp) :: cold ! number densities at previous timestep (molecule.cm-3)
125	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities at current timestep (molecule.cm-3)
126	real (kind = 8), dimension(nesp) :: cnew ! number densities at next timestep (molecule.cm-3)
127
128	! reaction rates
129
130	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
131	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
132	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
133	logical :: hetero_dust, hetero_ice
134
135	! matrix
136
137	real (kind = 8), dimension(nesp,nesp) :: mat, mat1
138	integer, dimension(nesp) :: indx
139	integer :: code
140
141	! production and loss terms (for first-guess solution only)
142
143	real (kind = 8), dimension(nesp) :: prod, loss
144
145	!===================================================================
146	! initialisation of the reaction indexes
147	!===================================================================
148
149	if (firstcall) then
150	print*,'photochemistry: initialize indexes'
151	call indice(nb_reaction_3_max,nb_reaction_4_max, &
152	nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
153	i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
154	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
155	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
156	i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
157	firstcall = .false.
158	end if
159
160	!===================================================================
161	! initialisation of mixing ratios and densities
162	!===================================================================
163
164	call gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp, &
165	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
166	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
167	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
168	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
169	i_n2plus, i_nplus, i_hplus, i_hco2plus, &
170	i_elec, dens, rm, c)
171
172	!===================================================================
173	! photolysis rates
174	!===================================================================
175
176	jionos = .true.
177
178	if (jonline) then
179	if (sza <= 113.) then ! day at 300 km
180	call photolysis_online(nlayer, nb_phot_max, alt, press, temp, zmmean, &
181	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
182	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
183	i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, &
184	tau, sza, dist_sol, v_phot)
185
186	if (jionos .and. ionchem) then
187	call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
188	do ilay=1,lswitch-1
189	call phdisrate(ig,nlayer,2,sza,ilay)
190	enddo
191	v_phot(:,14)=jion(1,:,1)
192	v_phot(:,15)=jion(1,:,2)
193	v_phot(:,16)=jion(1,:,2)
194	v_phot(:,17)=jion(1,:,3)
195	v_phot(:,18)=jion(1,:,3)
196	v_phot(:,19)=jion(1,:,4)
197	v_phot(:,20)=jion(1,:,4)
198	v_phot(:,21)=jion(2,:,1)
199	v_phot(:,22)=jion(3,:,1)
200	v_phot(:,23)=jion(10,:,1)
201	v_phot(:,24)=jion(11,:,1)
202	v_phot(:,25)=jion(11,:,2)
203	v_phot(:,26)=jion(11,:,2)
204	v_phot(:,27)=jion(8,:,1)
205	v_phot(:,28)=jion(8,:,2)
206	v_phot(:,29)=jion(8,:,2)
207	v_phot(:,30)=jion(9,:,1)
208	v_phot(:,31)=jion(12,:,1)
209	endif
210	! write(,)'photochemistry/205',c(:,i_co2),ig
211	! write(,)'photochemistry/206',v_phot(:,3),ig
212
213	else ! night
214	v_phot(:,:) = 0.
215	end if
216	!else if(jparam) then
217	! call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
218	! do ilay=1,lswitch-1
219	! call phdisrate(ig,nlayer,2,sza,ilay)
220	! enddo
221	! v_phot(:,1)=jdistot(2,:)
222	! v_phot(:,2)=jdistot_b(2,:)
223	! v_phot(:,3)=jdistot(1,:)
224	! v_phot(:,4)=jdistot_b(1,:)
225	! v_phot(:,5)=jdistot(7,:)
226	! v_phot(:,6)=jdistot_b(7,:)
227	! v_phot(:,7)=jdistot(4,:)
228	! v_phot(:,8)=jdistot(6,:)
229	! v_phot(:,10)=jdistot(5,:)
230	! v_phot(:,11)=jdistot(10,:)
231	! v_phot(:,12)=jdistot(13,:)
232	! v_phot(:,13)=jdistot(8,:)
233	! v_phot(:,14)=jion(1,:,1)
234	! v_phot(:,15)=jion(1,:,2)
235	! v_phot(:,16)=jion(1,:,2)
236	! v_phot(:,17)=jion(1,:,3)
237	! v_phot(:,18)=jion(1,:,3)
238	! v_phot(:,19)=jion(1,:,4)
239	! v_phot(:,20)=jion(1,:,4)
240	! v_phot(:,21)=jion(2,:,1)
241	! v_phot(:,22)=jion(3,:,1)
242	! v_phot(:,23)=jion(10,:,1)
243	! v_phot(:,24)=jion(11,:,1)
244	! v_phot(:,25)=jion(11,:,2)
245	! v_phot(:,26)=jion(11,:,2)
246	! v_phot(:,27)=jion(8,:,1)
247	! v_phot(:,28)=jion(8,:,2)
248	! v_phot(:,29)=jion(8,:,2)
249	! v_phot(:,30)=jion(9,:,1)
250	! v_phot(:,31)=jion(12,:,1)
251	else
252	tau = tau*7./press(1) ! dust in the lookup table is at 7 hpa
253	call photolysis(nlayer, nb_phot_max, lswitch, press, temp, sza, tau, &
254	zmmean, dist_sol,rm(:,i_co2), rm(:,i_o3), v_phot)
255	end if
256	! save o3 and h2o photolysis for output
257
258	j_o3_o1d = 5
259	jo3(:) = v_phot(:,j_o3_o1d)
260	j_h2o = 7
261	jh2o(:) = v_phot(:,j_h2o)
262
263	!===================================================================
264	! reaction rates
265	!===================================================================
266	! switches for heterogeneous chemistry
267	! hetero_ice : reactions on ice clouds
268	! hetero_dust : reactions on dust
269	!===================================================================
270
271	hetero_dust = .false.
272	hetero_ice = .true.
273
274	call reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
275	nb_phot_max, nphotion, lswitch, dens, c(:,i_co2), &
276	c(:,i_o2), c(:,i_o), c(:,i_n2), press, temp, &
277	temp_elect, hetero_dust, hetero_ice, &
278	surfdust1d, surfice1d, v_phot, v_3, v_4)
279
280	!===================================================================
281	! ctimestep : standard chemical timestep (s), defined as
282	! the fraction phychemrat of the physical timestep
283	!===================================================================
284
285	phychemrat = 1
286
287	ctimestep = ptimestep/real(phychemrat)
288
289	!print*, "ptimestep = ", ptimestep
290	!print*, "phychemrat = ", phychemrat
291	!print*, "ctimestep = ", ctimestep
292	!stop
293
294	!===================================================================
295	! loop over levels
296	!===================================================================
297
298	do ilev = 1,lswitch - 1
299
300	! initializations
301
302	time = 0.
303	iter(ilev) = 0
304	dt_guess = ctimestep
305	cold(:) = c(ilev,:)
306
307	! internal loop for the chemistry
308
309	do while (time < ptimestep)
310
311	iter(ilev) = iter(ilev) + 1 ! iteration counter
312
313	! first-guess: fill matrix
314
315	call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, &
316	nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
317	v_phot, v_3, v_4)
318
319	! adaptative evaluation of the sub time step
320
321	call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:), &
322	mat1, prod, loss, dens(ilev))
323
324	if (time + dt_corrected > ptimestep) then
325	dt_corrected = ptimestep - time
326	end if
327
328	! if (dt_corrected /= dt_guess) then ! the timestep has been modified
329
330	! form the matrix identity + mat*dt_corrected
331
332	mat(:,:) = mat1(:,:)*dt_corrected
333	do iesp = 1,nesp
334	mat(iesp,iesp) = 1. + mat(iesp,iesp)
335	end do
336
337	! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
338	cnew(:) = c(ilev,:)
339
340	#ifdef LAPACK
341	call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
342	#else
343	write(,) "photochemistry error, missing LAPACK routine dgesv"
344	stop
345	#endif
346
347	! end if
348
349	! eliminate small values
350
351	where (cnew(:)/dens(ilev) < 1.e-30)
352	cnew(:) = 0.
353	end where
354
355	! update concentrations
356
357	cold(:) = c(ilev,:)
358	c(ilev,:) = cnew(:)
359
360	! force charge neutrality (mod fgg, july 2019)
361
362	if (ionchem) then
363	if(c(ilev,i_elec).ne.c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+&
364	c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+&
365	c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)) then
366	c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
367	c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+ &
368	c(ilev,i_n2plus)+c(ilev,i_nplus)+c(ilev,i_hplus)+ &
369	c(ilev,i_hco2plus)
370	! write(,)'photochemistry/359'
371	! write(,)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
372	end if
373	end if
374	cnew(:) = 0.
375
376	! increment internal time
377
378	time = time + dt_corrected
379	dt_guess = dt_corrected ! first-guess timestep for next iteration
380
381	end do ! while (time < ptimestep)
382
383	end do ! ilev
384
385	!===================================================================
386	! save chemical species for the gcm
387	!===================================================================
388
389	call chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp, &
390	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
391	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
392	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
393	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
394	i_n2plus, i_nplus, i_hplus, i_hco2plus, &
395	i_elec, dens, c)
396	contains
397
398	!================================================================
399
400	subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
401	prod, loss, dens)
402
403	!================================================================
404	! iterative evaluation of the appropriate time step dtnew
405	! according to curvature criterion based on
406	! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
407	! with r = (tn - tn-1)/(tn+1 - tn)
408	!================================================================
409
410	implicit none
411
412	! input
413
414	integer :: nesp ! number of species in the chemistry
415
416	real :: dtold, ctimestep
417	real (kind = 8), dimension(nesp) :: cold, ccur
418	real (kind = 8), dimension(nesp,nesp) :: mat1
419	real (kind = 8), dimension(nesp) :: prod, loss
420	real :: dens
421
422	! output
423
424	real :: dtnew
425
426	! local
427
428	real (kind = 8), dimension(nesp) :: cnew
429	real (kind = 8), dimension(nesp,nesp) :: mat
430	real (kind = 8) :: atol, ratio, e, es, coef
431
432	integer :: code, iesp, iter
433	integer, dimension(nesp) :: indx
434
435	real :: dttest
436
437	! parameters
438
439	real (kind = 8), parameter :: dtmin = 10. ! minimum time step (s)
440	real (kind = 8), parameter :: vmrtol = 1.e-11 ! absolute tolerance on vmr
441	real (kind = 8), parameter :: rtol = 0.05 ! rtol recommended value : 0.1-0.02
442	integer, parameter :: niter = 3 ! number of iterations
443	real (kind = 8), parameter :: coefmax = 2.
444	real (kind = 8), parameter :: coefmin = 0.1
445	logical :: fast_guess = .true.
446
447	dttest = dtold ! dttest = dtold = dt_guess
448
449	atol = vmrtol*dens ! absolute tolerance in molecule.cm-3
450
451	do iter = 1,niter
452
453	if (fast_guess) then
454
455	! first guess : fast semi-implicit method
456
457	do iesp = 1, nesp
458	cnew(iesp) = (ccur(iesp) + prod(iesp)dttest)/(1. + loss(iesp)dttest)
459	end do
460
461	else
462
463	! first guess : form the matrix identity + mat*dt_guess
464
465	mat(:,:) = mat1(:,:)*dttest
466	do iesp = 1,nesp
467	mat(iesp,iesp) = 1. + mat(iesp,iesp)
468	end do
469
470	! form right-hand side (RHS) of the system
471
472	cnew(:) = ccur(:)
473
474	! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
475
476	#ifdef LAPACK
477	call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
478	#else
479	write(,) "photochemistry error, missing LAPACK routine dgesv"
480	stop
481	#endif
482
483	end if
484
485	! ratio old/new subtimestep
486
487	ratio = dtold/dttest
488
489	! e : local error indicatocitr
490
491	e = 0.
492
493	do iesp = 1,nesp
494	es = 2.abs((ratiocnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp)) &
495	/(1. + ratio)/max(ccur(iesp)*rtol,atol))
496
497	if (es > e) then
498	e = es
499	end if
500	end do
501
502	! timestep correction
503
504	coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
505
506	dttest = max(dtmin,dttest*coef)
507	dttest = min(ctimestep,dttest)
508
509	end do ! iter
510
511	! new timestep
512
513	dtnew = dttest
514
515	end subroutine define_dt
516
517	!======================================================================
518
519	subroutine reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
520	nb_phot_max, nphotion, lswitch, dens, co2, o2, o, &
521	n2, press, t, t_elect, hetero_dust, hetero_ice, &
522	surfdust1d, surfice1d, v_phot, v_3, v_4)
523
524	!================================================================
525	! compute reaction rates !
526	!----------------------------------------------------------------
527	! reaction type array !
528	!----------------------------------------------------------------
529	! A + B --> C + D bimolecular v_4 !
530	! A + A --> B + C quadratic v_3 !
531	! A + C --> B + C quenching v_phot !
532	! A + ice --> B + C heterogeneous v_phot !
533	!================================================================
534
535	use comcstfi_h
536	use photolysis_mod, only : nphot
537
538	implicit none
539
540	!----------------------------------------------------------------------
541	! input
542	!----------------------------------------------------------------------
543
544	integer, intent(in) :: nlayer ! number of atmospheric layers
545	integer, intent(in) :: nb_reaction_3_max ! number of quadratic reactions
546	integer, intent(in) :: nb_reaction_4_max ! number of bimolecular reactions
547	integer, intent(in) :: nb_phot_max ! number of reactions treated numerically as photodissociations
548	integer, intent(in) :: nphotion ! number of photoionizations
549	logical, intent(in) :: ionchem
550	integer :: lswitch ! interface level between lower
551	! atmosphere and thermosphere chemistries
552	real, dimension(nlayer) :: dens ! total number density (molecule.cm-3)
553	real, dimension(nlayer) :: press ! pressure (hPa)
554	real, dimension(nlayer) :: t ! temperature (K)
555	real, dimension(nlayer) :: t_elect ! electronic temperature (K)
556	real, dimension(nlayer) :: surfdust1d ! dust surface area (cm2.cm-3)
557	real, dimension(nlayer) :: surfice1d ! ice surface area (cm2.cm-3)
558	real (kind = 8), dimension(nlayer) :: co2 ! co2 number density (molecule.cm-3)
559	real (kind = 8), dimension(nlayer) :: o2 ! o2 number density (molecule.cm-3)
560	real (kind = 8), dimension(nlayer) :: o ! o number density (molecule.cm-3)
561	real (kind = 8), dimension(nlayer) :: n2 ! n2 number density (molecule.cm-3)
562	logical :: hetero_dust, hetero_ice ! switches for heterogeneous chemistry
563
564	!----------------------------------------------------------------------
565	! output
566	!----------------------------------------------------------------------
567
568	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
569	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
570	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
571
572	!----------------------------------------------------------------------
573	! local
574	!----------------------------------------------------------------------
575
576	integer :: ilev
577	integer :: nb_phot, nb_reaction_3, nb_reaction_4
578	real :: ak0, ak1, xpo, rate, rate1, rate2
579	real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
580	real, dimension(nlayer) :: deq
581	real, dimension(nlayer) :: a001, a002, a003, &
582	b001, b002, b003, b004, b005, b006, b007, &
583	b008, b009, &
584	c001, c002, c003, c004, c005, c006, c007, &
585	c008, c009, c010, c011, c012, c013, c014, &
586	c015, c016, c017, c018, &
587	d001, d002, d003, d004, d005, d006, d007, &
588	d008, d009, d010, d011, d012, &
589	e001, e002, &
590	i001, i002, i003, i004, i005, i006, &
591	i007, i008, i009, i010, i011, i012, &
592	i013, i014, i015, i016, i017, i018, i019, &
593	i020, i021, i022, i023, i024, i025, i026, &
594	i027, i028, i029, &
595	h001, h002, h003, h004, h005
596
597	!----------------------------------------------------------------------
598	! initialisation
599	!----------------------------------------------------------------------
600
601	nb_phot = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates
602	nb_reaction_3 = 0
603	nb_reaction_4 = 0
604
605	!----------------------------------------------------------------------
606	! oxygen reactions
607	!----------------------------------------------------------------------
608
609	!--- a001: o + o2 + co2 -> o3 + co2
610
611	! jpl 2003
612	!
613	! co2/n2 efficiency as a third body = 2.075
614	! from sehested et al., j. geophys. res., 100, 1995.
615
616	a001(:) = 2.0756.0e-34(t(:)/300.)*(-2.4)dens(:)
617
618	nb_reaction_4 = nb_reaction_4 + 1
619	v_4(:,nb_reaction_4) = a001(:)
620
621	!--- a002: o + o + co2 -> o2 + co2
622
623	! Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
624
625	! a002(:) = 2.55.2e-35exp(900./t(:))*dens(:)
626
627	! Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
628
629	! a002(:) = 1.2e-32(300./t(:))(2.0)dens(:) ! yung expression
630	a002(:) = 2.59.46e-34exp(485./t(:))*dens(:) ! nist expression
631
632	nb_reaction_3 = nb_reaction_3 + 1
633	v_3(:,nb_reaction_3) = a002(:)
634
635	!--- a003: o + o3 -> o2 + o2
636
637	! jpl 2003
638
639	a003(:) = 8.0e-12*exp(-2060./t(:))
640
641	nb_reaction_4 = nb_reaction_4 + 1
642	v_4(:,nb_reaction_4) = a003(:)
643
644	!----------------------------------------------------------------------
645	! o(1d) reactions
646	!----------------------------------------------------------------------
647
648	!--- b001: o(1d) + co2 -> o + co2
649
650	! jpl 2006
651
652	b001(:) = 7.5e-11*exp(115./t(:))
653
654	nb_phot = nb_phot + 1
655	v_phot(:,nb_phot) = b001(:)*co2(:)
656
657	!--- b002: o(1d) + h2o -> oh + oh
658
659	! jpl 2006
660
661	b002(:) = 1.63e-10*exp(60./t(:))
662
663	nb_reaction_4 = nb_reaction_4 + 1
664	v_4(:,nb_reaction_4) = b002(:)
665
666	!--- b003: o(1d) + h2 -> oh + h
667
668	! jpl 2011
669
670	b003(:) = 1.2e-10
671
672	nb_reaction_4 = nb_reaction_4 + 1
673	v_4(:,nb_reaction_4) = b003(:)
674
675	!--- b004: o(1d) + o2 -> o + o2
676
677	! jpl 2006
678
679	b004(:) = 3.3e-11*exp(55./t(:))
680
681	nb_phot = nb_phot + 1
682	v_phot(:,nb_phot) = b004(:)*o2(:)
683
684	!--- b005: o(1d) + o3 -> o2 + o2
685
686	! jpl 2003
687
688	b005(:) = 1.2e-10
689
690	nb_reaction_4 = nb_reaction_4 + 1
691	v_4(:,nb_reaction_4) = b005(:)
692
693	!--- b006: o(1d) + o3 -> o2 + o + o
694
695	! jpl 2003
696
697	b006(:) = 1.2e-10
698
699	nb_reaction_4 = nb_reaction_4 + 1
700	v_4(:,nb_reaction_4) = b006(:)
701
702	!--- b007: o(1d) + ch4 -> ch3 + oh
703
704	! jpl 2003
705
706	b007(:) = 1.5e-10*0.75
707
708	!--- b008: o(1d) + ch4 -> ch3o + h
709
710	! jpl 2003
711
712	b008(:) = 1.5e-10*0.20
713	!
714	!--- b009: o(1d) + ch4 -> ch2o + h2
715
716	! jpl 2003
717
718	b009(:) = 1.5e-10*0.05
719
720	!----------------------------------------------------------------------
721	! hydrogen reactions
722	!----------------------------------------------------------------------
723
724	!--- c001: o + ho2 -> oh + o2
725
726	! jpl 2003
727
728	c001(:) = 3.0e-11*exp(200./t(:))
729
730	nb_reaction_4 = nb_reaction_4 + 1
731	v_4(:,nb_reaction_4) = c001(:)
732
733	!--- c002: o + oh -> o2 + h
734
735	! jpl 2011
736
737	c002(:) = 1.8e-11*exp(180./t(:))
738
739	! robertson and smith, j. chem. phys. a 110, 6673, 2006
740
741	! c002(:) = 11.2e-11t(:)(-0.32)exp(177./t(:))
742
743	nb_reaction_4 = nb_reaction_4 + 1
744	v_4(:,nb_reaction_4) = c002(:)
745
746	!--- c003: h + o3 -> oh + o2
747
748	! jpl 2003
749
750	c003(:) = 1.4e-10*exp(-470./t(:))
751
752	nb_reaction_4 = nb_reaction_4 + 1
753	v_4(:,nb_reaction_4) = c003(:)
754
755	!--- c004: h + ho2 -> oh + oh
756
757	! jpl 2006
758
759	c004(:) = 7.2e-11
760
761	nb_reaction_4 = nb_reaction_4 + 1
762	v_4(:,nb_reaction_4) = c004(:)
763
764	!--- c005: h + ho2 -> h2 + o2
765
766	! jpl 2006
767
768	c005(:) = 6.9e-12
769
770	nb_reaction_4 = nb_reaction_4 + 1
771	v_4(:,nb_reaction_4) = c005(:)
772
773	!--- c006: h + ho2 -> h2o + o
774
775	! jpl 2006
776
777	c006(:) = 1.6e-12
778
779	nb_reaction_4 = nb_reaction_4 + 1
780	v_4(:,nb_reaction_4) = c006(:)
781
782	!--- c007: oh + ho2 -> h2o + o2
783
784	! jpl 2003
785
786	! canty et al., grl, 2006 suggest to increase this rate
787	! by 20%. not done here.
788
789	c007(:) = 4.8e-11*exp(250./t(:))
790
791	nb_reaction_4 = nb_reaction_4 + 1
792	v_4(:,nb_reaction_4) = c007(:)
793
794	!--- c008: ho2 + ho2 -> h2o2 + o2
795
796	! jpl 2015
797
798	! c008(:) = 3.0e-13*exp(460./t(:))
799
800	! christensen et al., grl, 13, 2002
801
802	c008(:) = 1.5e-12*exp(19./t(:))
803
804	nb_reaction_3 = nb_reaction_3 + 1
805	v_3(:,nb_reaction_3) = c008(:)
806
807	!--- c009: oh + h2o2 -> h2o + ho2
808
809	! jpl 2006
810
811	c009(:) = 1.8e-12
812
813	nb_reaction_4 = nb_reaction_4 + 1
814	v_4(:,nb_reaction_4) = c009(:)
815
816	!--- c010: oh + h2 -> h2o + h
817
818	! jpl 2006
819
820	c010(:) = 2.8e-12*exp(-1800./t(:))
821
822	nb_reaction_4 = nb_reaction_4 + 1
823	v_4(:,nb_reaction_4) = c010(:)
824
825	!--- c011: h + o2 + co2 -> ho2 + co2
826
827	! jpl 2011
828	! co2/n2 efficiency as a third body = 2.4
829	! from ashman and haynes, 27th symposium on combustion, 1998.
830
831	do ilev = 1,lswitch-1
832	ak0 = 2.44.4e-32(t(ilev)/300.)**(-1.3)
833	ak1 = 7.5e-11(t(ilev)/300.)*(0.2)
834
835	rate = (ak0dens(ilev))/(1. + ak0dens(ilev)/ak1)
836	xpo = 1./(1. + alog10((ak0dens(ilev))/ak1)*2)
837	c011(ilev) = rate0.6*xpo
838	end do
839
840	nb_reaction_4 = nb_reaction_4 + 1
841	v_4(:,nb_reaction_4) = c011(:)
842
843	!--- c012: o + h2o2 -> oh + ho2
844
845	! jpl 2003
846
847	c012(:) = 1.4e-12*exp(-2000./t(:))
848
849	nb_reaction_4 = nb_reaction_4 + 1
850	v_4(:,nb_reaction_4) = c012(:)
851
852	!--- c013: oh + oh -> h2o + o
853
854	! jpl 2006
855
856	c013(:) = 1.8e-12
857
858	nb_reaction_3 = nb_reaction_3 + 1
859	v_3(:,nb_reaction_3) = c013(:)
860
861	!--- c014: oh + o3 -> ho2 + o2
862
863	! jpl 2003
864
865	c014(:) = 1.7e-12*exp(-940./t(:))
866
867	nb_reaction_4 = nb_reaction_4 + 1
868	v_4(:,nb_reaction_4) = c014(:)
869
870	!--- c015: ho2 + o3 -> oh + o2 + o2
871
872	! jpl 2003
873
874	c015(:) = 1.0e-14*exp(-490./t(:))
875
876	nb_reaction_4 = nb_reaction_4 + 1
877	v_4(:,nb_reaction_4) = c015(:)
878
879	!--- c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
880
881	! jpl 2011
882
883	c016(:) = 2.52.1e-33exp(920./t(:))*dens(:)
884
885	nb_reaction_3 = nb_reaction_3 + 1
886	v_3(:,nb_reaction_3) = c016(:)
887
888	!--- c017: oh + oh + co2 -> h2o2 + co2
889
890	! jpl 2003
891
892	do ilev = 1,lswitch-1
893	ak0 = 2.56.9e-31(t(ilev)/300.)**(-1.0)
894	ak1 = 2.6e-11(t(ilev)/300.)*(0.0)
895
896	rate = (ak0dens(ilev))/(1. + ak0dens(ilev)/ak1)
897	xpo = 1./(1. + alog10((ak0dens(ilev))/ak1)*2)
898	c017(ilev) = rate0.6*xpo
899	end do
900
901	nb_reaction_3 = nb_reaction_3 + 1
902	v_3(:,nb_reaction_3) = c017(:)
903
904	!--- c018: h + h + co2 -> h2 + co2
905
906	! baulch et al., 2005
907
908	c018(:) = 2.51.8e-30(t(:)*(-1.0))dens(:)
909
910	nb_reaction_3 = nb_reaction_3 + 1
911	v_3(:,nb_reaction_3) = c018(:)
912
913	!----------------------------------------------------------------------
914	! nitrogen reactions
915	!----------------------------------------------------------------------
916
917	!--- d001: no2 + o -> no + o2
918
919	! jpl 2006
920
921	d001(:) = 5.1e-12*exp(210./t(:))
922
923	nb_reaction_4 = nb_reaction_4 + 1
924	v_4(:,nb_reaction_4) = d001(:)
925
926	!--- d002: no + o3 -> no2 + o2
927
928	! jpl 2006
929
930	d002(:) = 3.0e-12*exp(-1500./t(:))
931
932	nb_reaction_4 = nb_reaction_4 + 1
933	v_4(:,nb_reaction_4) = d002(:)
934
935	!--- d003: no + ho2 -> no2 + oh
936
937	! jpl 2011
938
939	d003(:) = 3.3e-12*exp(270./t(:))
940
941	nb_reaction_4 = nb_reaction_4 + 1
942	v_4(:,nb_reaction_4) = d003(:)
943
944	!--- d004: n + no -> n2 + o
945
946	! jpl 2011
947
948	d004(:) = 2.1e-11*exp(100./t(:))
949
950	nb_reaction_4 = nb_reaction_4 + 1
951	v_4(:,nb_reaction_4) = d004(:)
952
953	!--- d005: n + o2 -> no + o
954
955	! jpl 2011
956
957	d005(:) = 1.5e-11*exp(-3600./t(:))
958
959	nb_reaction_4 = nb_reaction_4 + 1
960	v_4(:,nb_reaction_4) = d005(:)
961
962	!--- d006: no2 + h -> no + oh
963
964	! jpl 2011
965
966	d006(:) = 4.0e-10*exp(-340./t(:))
967
968	nb_reaction_4 = nb_reaction_4 + 1
969	v_4(:,nb_reaction_4) = d006(:)
970
971	!--- d007: n + o -> no
972
973	d007(:) = 2.8e-17(300./t(:))*0.5
974
975	nb_reaction_4 = nb_reaction_4 + 1
976	v_4(:,nb_reaction_4) = d007(:)
977
978	!--- d008: n + ho2 -> no + oh
979
980	! brune et al., j. chem. phys., 87, 1983
981
982	d008(:) = 2.19e-11
983
984	nb_reaction_4 = nb_reaction_4 + 1
985	v_4(:,nb_reaction_4) = d008(:)
986
987	!--- d009: n + oh -> no + h
988
989	! atkinson et al., j. phys. chem. ref. data, 18, 881, 1989
990
991	d009(:) = 3.8e-11*exp(85./t(:))
992
993	nb_reaction_4 = nb_reaction_4 + 1
994	v_4(:,nb_reaction_4) = d009(:)
995
996	!--- d010: n(2d) + o -> n + o
997
998	! herron, j. phys. chem. ref. data, 1999
999
1000	d010(:) = 3.3e-12*exp(-260./t(:))
1001
1002	nb_phot = nb_phot + 1
1003	v_phot(:,nb_phot) = d010(:)*o(:)
1004
1005	!--- d011: n(2d) + n2 -> n + n2
1006
1007	! herron, j. phys. chem. ref. data, 1999
1008
1009	d011(:) = 1.7e-14
1010
1011	nb_phot = nb_phot + 1
1012	v_phot(:,nb_phot) = d011(:)*n2(:)
1013
1014	!--- d012: n(2d) + co2 -> no + co
1015
1016	! herron, j. phys. chem. ref. data, 1999
1017
1018	d012(:) = 3.6e-13
1019
1020	nb_reaction_4 = nb_reaction_4 + 1
1021	v_4(:,nb_reaction_4) = d012(:)
1022
1023	!----------------------------------------------------------------------
1024	! carbon reactions
1025	!----------------------------------------------------------------------
1026
1027	!--- e001: oh + co -> co2 + h
1028
1029	! jpl 2003
1030
1031	! e001(:) = 1.5e-13(1 + 0.6press(:)/1013.)
1032
1033	! mccabe et al., grl, 28, 3135, 2001
1034
1035	! e001(:) = 1.57e-13 + 3.54e-33*dens(:)
1036
1037	! jpl 2015
1038
1039	! do ilev = 1,lswitch-1
1040
1041	! branch 1 : oh + co -> h + co2
1042
1043	! rate1 = 1.5e-13(t(ilev)/300.)*(0.0)
1044
1045	! branch 2 : oh + co + m -> hoco + m
1046
1047	! ak0 = 5.9e-33(t(ilev)/300.)*(-1.0)
1048	! ak1 = 1.1e-12(t(ilev)/300.)*(1.3)
1049	! rate2 = (ak0dens(ilev))/(1. + ak0dens(ilev)/ak1)
1050	! xpo = 1./(1. + alog10((ak0dens(ilev))/ak1)*2)
1051
1052	! e001(ilev) = rate1 + rate20.6*xpo
1053	! end do
1054
1055	! joshi et al., 2006
1056
1057	do ilev = 1,lswitch-1
1058	k1a0 = 1.342.5dens(ilev) &
1059	1/(1/(3.62e-26t(ilev)*(-2.739)exp(-20./t(ilev))) &
1060	+ 1/(6.48e-33t(ilev)(0.14)exp(-57./t(ilev)))) ! typo in paper corrected
1061	k1b0 = 1.17e-19t(ilev)(2.053)exp(139./t(ilev)) &
1062	+ 9.56e-12t(ilev)(-0.664)exp(-167./t(ilev))
1063	k1ainf = 1.52e-17t(ilev)(1.858)exp(28.8/t(ilev)) &
1064	+ 4.78e-8t(ilev)(-1.851)exp(-318./t(ilev))
1065	x = k1a0/(k1ainf - k1b0)
1066	y = k1b0/(k1ainf - k1b0)
1067	fc = 0.628exp(-1223./t(ilev)) + (1. - 0.628)exp(-39./t(ilev)) &
1068	+ exp(-t(ilev)/255.)
1069	fx = fc(1./(1. + (alog(x))2)) ! typo in paper corrected
1070	k1a = k1a0((1. + y)/(1. + x))fx
1071	k1b = k1b0(1./(1.+x))fx
1072	e001(ilev) = k1a + k1b
1073	end do
1074
1075	nb_reaction_4 = nb_reaction_4 + 1
1076	v_4(:,nb_reaction_4) = e001(:)
1077
1078	!--- e002: o + co + m -> co2 + m
1079
1080	! tsang and hampson, 1986.
1081
1082	e002(:) = 2.56.5e-33exp(-2184./t(:))*dens(:)
1083
1084	nb_reaction_4 = nb_reaction_4 + 1
1085	v_4(:,nb_reaction_4) = e002(:)
1086
1087	!----------------------------------------------------------------------
1088	! ionospheric reactions
1089	! only if ionchem=true
1090	!----------------------------------------------------------------------
1091
1092	if (ionchem) then
1093
1094	!--- i001: co2+ + o2 -> o2+ + co2
1095
1096	! aninich, j. phys. chem. ref. data 1993
1097
1098	i001(:) = 5.5e-11(300./t_elect(:))*0.5
1099
1100	nb_reaction_4 = nb_reaction_4 + 1
1101	v_4(:,nb_reaction_4) = i001(:)
1102
1103	!--- i002: co2+ + o -> o+ + co2
1104
1105	! UMIST database
1106
1107	i002(:) = 9.6e-11
1108
1109	nb_reaction_4 = nb_reaction_4 + 1
1110	v_4(:,nb_reaction_4) = i002(:)
1111
1112	!--- i003: co2+ + o -> o2+ + co
1113
1114	! UMIST database
1115
1116	i003(:) = 1.64e-10
1117
1118	nb_reaction_4 = nb_reaction_4 + 1
1119	v_4(:,nb_reaction_4) = i003(:)
1120
1121	!--- i004: o2+ + e- -> o + o
1122
1123	! Alge et al., J. Phys. B, At. Mol. Phys. 1983
1124
1125	i004(:) = 2.0e-7(300./t_elect(:))*0.7
1126
1127	nb_reaction_4 = nb_reaction_4 + 1
1128	v_4(:,nb_reaction_4) = i004(:)
1129
1130	!--- i005: o+ + co2 -> o2+ + co
1131
1132	! UMIST database
1133
1134	i005(:) = 9.4e-10
1135
1136	nb_reaction_4 = nb_reaction_4 + 1
1137	v_4(:,nb_reaction_4) = i005(:)
1138
1139
1140	!--- i006: co2+ + e- -> co + o
1141
1142	! UMIST database
1143
1144	i006(:) = 3.8e-7(300./t_elect(:))*0.5
1145
1146	nb_reaction_4 = nb_reaction_4 + 1
1147	v_4(:,nb_reaction_4) = i006(:)
1148
1149
1150	!--- i007: co2+ + no -> no+ + co2
1151
1152	! UMIST database
1153
1154	i007(:) = 1.2e-10
1155
1156	nb_reaction_4 = nb_reaction_4 + 1
1157	v_4(:,nb_reaction_4) = i007(:)
1158
1159	!--- i008: o2+ + no -> no+ + o2
1160
1161	! UMIST database
1162
1163	i008(:) = 4.6e-10
1164
1165	nb_reaction_4 = nb_reaction_4 + 1
1166	v_4(:,nb_reaction_4) = i008(:)
1167
1168	!--- i009: o2+ + n2 -> no+ + no
1169
1170	! Fox & Sung 2001
1171
1172	i009(:) = 1.0e-15
1173
1174	nb_reaction_4 = nb_reaction_4 + 1
1175	v_4(:,nb_reaction_4) = i009(:)
1176
1177	!--- i010: o2+ + n -> no+ + o
1178
1179	! Fox & Sung 2001
1180
1181	i010(:) = 1.0e-10
1182
1183	nb_reaction_4 = nb_reaction_4 + 1
1184	v_4(:,nb_reaction_4) = i010(:)
1185
1186	!--- i011: o+ + n2 -> no+ + n
1187
1188	! Fox & Sung 2001
1189
1190	i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
1191
1192	nb_reaction_4 = nb_reaction_4 + 1
1193	v_4(:,nb_reaction_4) = i011(:)
1194
1195	!--- i012: no+ + e -> n + o
1196
1197	! UMIST database
1198
1199	i012(:) = 4.3e-7(300./t_elect(:))*0.37
1200
1201	nb_reaction_4 = nb_reaction_4 + 1
1202	v_4(:,nb_reaction_4) = i012(:)
1203
1204
1205	!--- i013: co+ + co2 -> co2+ + co
1206
1207	! UMIST database
1208
1209	i013(:) = 1.0e-9
1210
1211	nb_reaction_4 = nb_reaction_4 + 1
1212	v_4(:,nb_reaction_4) = i013(:)
1213
1214
1215	!--- i014: co+ + o -> o+ + co
1216
1217	! UMIST database
1218
1219	i014(:) = 1.4e-10
1220
1221	nb_reaction_4 = nb_reaction_4 + 1
1222	v_4(:,nb_reaction_4) = i014(:)
1223
1224	!--- i015: c+ + co2 -> co+ + co
1225
1226	! UMIST database
1227
1228	i015(:) = 1.1e-9
1229
1230	nb_reaction_4 = nb_reaction_4 + 1
1231	v_4(:,nb_reaction_4) = i015(:)
1232
1233
1234	!--- i016: N2+ + co2 -> co2+ + N2
1235
1236	! Fox & Song 2001
1237
1238	i016(:) = 9.0e-10(300./t_elect(:))*0.23
1239
1240	nb_reaction_4 = nb_reaction_4 + 1
1241	v_4(:,nb_reaction_4) = i016(:)
1242
1243
1244	!--- i017: N2+ + o -> no+ + N
1245
1246	! Fox & Song 2001
1247
1248	i017(:) = 1.33e-10(300./t_elect(:))*0.44
1249
1250	nb_reaction_4 = nb_reaction_4 + 1
1251	v_4(:,nb_reaction_4) = i017(:)
1252
1253	!--- i018: N2+ + co -> co+ + N2
1254
1255	! UMIST
1256
1257	i018(:) = 7.4e-11
1258
1259	nb_reaction_4 = nb_reaction_4 + 1
1260	v_4(:,nb_reaction_4) = i018(:)
1261
1262	!--- i019: N2+ + e -> N + N
1263
1264	! UMIST
1265
1266	i019(:) = 7.7e-7(300./t_elect(:))*0.3
1267
1268	nb_reaction_4 = nb_reaction_4 + 1
1269	v_4(:,nb_reaction_4) = i016(:)
1270
1271	!--- i020: N2+ + o -> o+ + N2
1272
1273	! Fox & Song 2001
1274
1275	i020(:) = 7.0e-12(300./t_elect(:))*0.23
1276
1277	nb_reaction_4 = nb_reaction_4 + 1
1278	v_4(:,nb_reaction_4) = i020(:)
1279
1280	!--- i021: N+ + co2 -> co2+ + N
1281
1282	! UMIST
1283
1284	i021(:) = 7.5e-10
1285
1286	nb_reaction_4 = nb_reaction_4 + 1
1287	v_4(:,nb_reaction_4) = i021(:)
1288
1289	!--- i022: CO+ + H -> H+ + CO
1290
1291	! Fox & Sung 2001
1292
1293	i022(:) = 4.0e-10
1294
1295	nb_reaction_4 = nb_reaction_4 + 1
1296	v_4(:,nb_reaction_4) = i022(:)
1297
1298	!--- i023: O+ + H -> H+ + O
1299
1300	! UMIST
1301
1302	i023(:) = 5.66e-10((t_elect(:)/300.)0.36)exp(8.6/t_elect(:))
1303
1304	nb_reaction_4 = nb_reaction_4 + 1
1305	v_4(:,nb_reaction_4) = i023(:)
1306
1307	!--- i024: H+ + O -> O+ + H
1308
1309	! UMIST
1310
1311	i024(:) = 6.86e-10((t_elect(:)/300.)0.26)exp(-224.3/t_elect(:))
1312
1313	nb_reaction_4 = nb_reaction_4 + 1
1314	v_4(:,nb_reaction_4) = i024(:)
1315
1316	!--- i025: CO+ + H2 -> HCO2+ + H
1317
1318	! UMIST
1319
1320	i025(:) = 9.5e-10
1321
1322	nb_reaction_4 = nb_reaction_4 + 1
1323	v_4(:,nb_reaction_4) = i025(:)
1324
1325	!--- i026: HCO2+ + e -> H + CO2
1326
1327	! UMIST
1328
1329	i026(:) = 1.75e-8((300./t_elect(:))*0.5)
1330
1331	nb_reaction_4 = nb_reaction_4 + 1
1332	v_4(:,nb_reaction_4) = i026(:)
1333
1334	!--- i027+i028: HCO2+ + e -> H + O + CO
1335
1336	! UMIST
1337	!Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
1338	!i028: 0.5 (HCO2+ + e-) -> O + CO
1339
1340	i027(:) = 2.38e-7((300./t_elect(:))*0.5)
1341
1342	nb_reaction_4 = nb_reaction_4 + 1
1343	v_4(:,nb_reaction_4) = i027(:)
1344
1345	nb_reaction_4 = nb_reaction_4 + 1
1346	v_4(:,nb_reaction_4) = i027(:)
1347
1348	!--- i029: HCO2+ + e -> H + CO2
1349
1350	! UMIST
1351
1352	i029(:) = 9.45e-8((300./t_elect(:))*0.5)
1353
1354	nb_reaction_4 = nb_reaction_4 + 1
1355	v_4(:,nb_reaction_4) = i029(:)
1356
1357	end if !ionchem
1358
1359	!----------------------------------------------------------------------
1360	! heterogeneous chemistry
1361	!----------------------------------------------------------------------
1362
1363	if (hetero_ice) then
1364
1365	! k = (surfacevgamma)/4 (s-1)
1366	! v = 100sqrt(8rt/(pim)) (cm s-1)
1367
1368	!--- h001: ho2 + ice -> products
1369
1370	! cooper and abbatt, 1996: gamma = 0.025
1371
1372	gam = 0.025
1373	h001(:) = surfice1d(:) &
1374	100.sqrt(8.8.31t(:)/(33.e-3pi))gam/4.
1375
1376	! h002: oh + ice -> products
1377
1378	! cooper and abbatt, 1996: gamma = 0.03
1379
1380	gam = 0.03
1381	h002(:) = surfice1d(:) &
1382	100.sqrt(8.8.31t(:)/(17.e-3pi))gam/4.
1383
1384	!--- h003: h2o2 + ice -> products
1385
1386	! gamma = 0. test value
1387
1388	gam = 0.
1389	h003(:) = surfice1d(:) &
1390	100.sqrt(8.8.31t(:)/(34.e-3pi))gam/4.
1391	else
1392	h001(:) = 0.
1393	h002(:) = 0.
1394	h003(:) = 0.
1395	end if
1396
1397	nb_phot = nb_phot + 1
1398	v_phot(:,nb_phot) = h001(:)
1399
1400	nb_phot = nb_phot + 1
1401	v_phot(:,nb_phot) = h002(:)
1402
1403	nb_phot = nb_phot + 1
1404	v_phot(:,nb_phot) = h003(:)
1405
1406	if (hetero_dust) then
1407
1408	!--- h004: ho2 + dust -> products
1409
1410	! jacob, 2000: gamma = 0.2
1411	! see dereus et al., atm. chem. phys., 2005
1412
1413	gam = 0.2
1414	h004(:) = surfdust1d(:) &
1415	100.sqrt(8.8.31t(:)/(33.e-3pi))gam/4.
1416
1417	!--- h005: h2o2 + dust -> products
1418
1419	! gamma = 5.e-4
1420	! see dereus et al., atm. chem. phys., 2005
1421
1422	gam = 5.e-4
1423	h005(:) = surfdust1d(:) &
1424	100.sqrt(8.8.31t(:)/(34.e-3pi))gam/4.
1425	else
1426	h004(:) = 0.
1427	h005(:) = 0.
1428	end if
1429
1430	nb_phot = nb_phot + 1
1431	v_phot(:,nb_phot) = h004(:)
1432
1433	nb_phot = nb_phot + 1
1434	v_phot(:,nb_phot) = h005(:)
1435
1436	end subroutine reactionrates
1437
1438	!======================================================================
1439
1440	subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, &
1441	nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
1442	v_phot, v_3, v_4)
1443
1444	!======================================================================
1445	! filling of the jacobian matrix
1446	!======================================================================
1447
1448	use types_asis
1449
1450	implicit none
1451
1452	! input
1453
1454	integer :: ilev ! level index
1455	integer :: nesp ! number of species in the chemistry
1456	integer, intent(in) :: nlayer ! number of atmospheric layers
1457	integer, intent(in) :: nb_reaction_3_max
1458	! number of quadratic reactions
1459	integer, intent(in) :: nb_reaction_4_max
1460	! number of bimolecular reactions
1461	integer, intent(in) :: nb_phot_max
1462	! number of processes treated numerically as photodissociations
1463
1464	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities
1465	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
1466	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
1467	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
1468
1469	! output
1470
1471	real (kind = 8), dimension(nesp,nesp), intent(out) :: mat ! matrix
1472	real (kind = 8), dimension(nesp), intent(out) :: prod, loss
1473
1474	! local
1475
1476	integer :: iesp
1477	integer :: ind_phot_2,ind_phot_4,ind_phot_6
1478	integer :: ind_3_2,ind_3_4,ind_3_6
1479	integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
1480	integer :: iphot,i3,i4
1481
1482	real(kind = 8) :: eps, eps_4 ! implicit/explicit coefficient
1483
1484	! initialisations
1485
1486	mat(:,:) = 0.
1487	prod(:) = 0.
1488	loss(:) = 0.
1489
1490	! photodissociations
1491	! or reactions a + c -> b + c
1492	! or reactions a + ice -> b + c
1493
1494	do iphot = 1,nb_phot_max
1495
1496	ind_phot_2 = indice_phot(iphot)%z2
1497	ind_phot_4 = indice_phot(iphot)%z4
1498	ind_phot_6 = indice_phot(iphot)%z6
1499
1500	mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
1501	mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
1502	mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
1503
1504	loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
1505	prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3v_phot(ilev,iphot)c(ilev,ind_phot_2)
1506	prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5v_phot(ilev,iphot)c(ilev,ind_phot_2)
1507
1508	end do
1509
1510	! reactions a + a -> b + c
1511
1512	do i3 = 1,nb_reaction_3_max
1513
1514	ind_3_2 = indice_3(i3)%z2
1515	ind_3_4 = indice_3(i3)%z4
1516	ind_3_6 = indice_3(i3)%z6
1517
1518	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)
1519	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)
1520	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)
1521
1522	loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1v_3(ilev,i3)c(ilev,ind_3_2)
1523	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)
1524	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)
1525
1526	end do
1527
1528	! reactions a + b -> c + d
1529
1530	eps = 1.d-10
1531
1532	do i4 = 1,nb_reaction_4_max
1533
1534	ind_4_2 = indice_4(i4)%z2
1535	ind_4_4 = indice_4(i4)%z4
1536	ind_4_6 = indice_4(i4)%z6
1537	ind_4_8 = indice_4(i4)%z8
1538
1539	eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
1540	eps_4 = min(eps_4,1.0)
1541
1542	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)
1543	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)
1544	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)
1545	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)
1546	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)
1547	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)
1548	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)
1549	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)
1550
1551
1552	loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1v_4(ilev,i4)c(ilev,ind_4_4)
1553	loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3v_4(ilev,i4)c(ilev,ind_4_2)
1554	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)
1555	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)
1556
1557	end do
1558
1559	end subroutine fill_matrix
1560
1561	!================================================================
1562
1563	subroutine indice(nb_reaction_3_max, nb_reaction_4_max, &
1564	nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
1565	i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
1566	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
1567	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
1568	i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
1569
1570	!================================================================
1571	! set the "indice" arrays used to fill the jacobian matrix !
1572	!----------------------------------------------------------------
1573	! reaction type !
1574	!----------------------------------------------------------------
1575	! A + hv --> B + C photolysis indice_phot !
1576	! A + B --> C + D bimolecular indice_4 !
1577	! A + A --> B + C quadratic indice_3 !
1578	! A + C --> B + C quenching indice_phot !
1579	! A + ice --> B + C heterogeneous indice_phot !
1580	!================================================================
1581
1582	use types_asis
1583
1584	implicit none
1585
1586	! input
1587
1588	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
1589	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
1590	i_n, i_n2d, i_no, i_no2, i_n2, &
1591	i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
1592	i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
1593	integer, intent(in) :: nb_reaction_3_max
1594	! number of quadratic reactions
1595	integer, intent(in) :: nb_reaction_4_max
1596	! number of bimolecular reactions
1597	integer, intent(in) :: nb_phot_max
1598	! number of processes treated numerically as photodissociations
1599	logical, intent(in) :: ionchem
1600
1601	! local
1602
1603	integer :: nb_phot, nb_reaction_3, nb_reaction_4
1604	integer :: i_dummy
1605
1606	allocate (indice_phot(nb_phot_max))
1607	allocate (indice_3(nb_reaction_3_max))
1608	allocate (indice_4(nb_reaction_4_max))
1609
1610	i_dummy = 1
1611
1612	nb_phot = 0
1613	nb_reaction_3 = 0
1614	nb_reaction_4 = 0
1615
1616	!===========================================================
1617	! O2 + hv -> O + O
1618	!===========================================================
1619
1620	nb_phot = nb_phot + 1
1621
1622	indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy)
1623
1624	!===========================================================
1625	! O2 + hv -> O + O(1D)
1626	!===========================================================
1627
1628	nb_phot = nb_phot + 1
1629
1630	indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d)
1631
1632	!===========================================================
1633	! CO2 + hv -> CO + O
1634	!===========================================================
1635
1636	nb_phot = nb_phot + 1
1637
1638	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o)
1639
1640	!===========================================================
1641	! CO2 + hv -> CO + O(1D)
1642	!===========================================================
1643
1644	nb_phot = nb_phot + 1
1645
1646	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d)
1647
1648	!===========================================================
1649	! O3 + hv -> O2 + O(1D)
1650	!===========================================================
1651
1652	nb_phot = nb_phot + 1
1653
1654	indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d)
1655
1656	!===========================================================
1657	! O3 + hv -> O2 + O
1658	!===========================================================
1659
1660	nb_phot = nb_phot + 1
1661
1662	indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o)
1663
1664	!===========================================================
1665	! H2O + hv -> H + OH
1666	!===========================================================
1667
1668	nb_phot = nb_phot + 1
1669
1670	indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh)
1671
1672	!===========================================================
1673	! H2O2 + hv -> OH + OH
1674	!===========================================================
1675
1676	nb_phot = nb_phot + 1
1677
1678	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy)
1679
1680	!===========================================================
1681	! HO2 + hv -> OH + O
1682	!===========================================================
1683
1684	nb_phot = nb_phot + 1
1685
1686	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o)
1687
1688	!===========================================================
1689	! H2 + hv -> H + H
1690	!===========================================================
1691
1692	nb_phot = nb_phot + 1
1693
1694	indice_phot(nb_phot) = z3spec(1.0, i_h2, 1.0, i_h, 1.0, i_h)
1695
1696	!===========================================================
1697	! NO + hv -> N + O
1698	!===========================================================
1699
1700	nb_phot = nb_phot + 1
1701
1702	indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o)
1703
1704	!===========================================================
1705	! NO2 + hv -> NO + O
1706	!===========================================================
1707
1708	nb_phot = nb_phot + 1
1709
1710	indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o)
1711
1712	!===========================================================
1713	! N2 + hv -> N + N
1714	!===========================================================
1715
1716	nb_phot = nb_phot + 1
1717
1718	indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n)
1719
1720	!Only if ion chemistry included
1721	if (ionchem) then
1722
1723	!===========================================================
1724	! CO2 + hv -> CO2+ + e-
1725	!===========================================================
1726
1727	nb_phot = nb_phot + 1
1728
1729	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
1730
1731	!===========================================================
1732	! CO2 + hv -> O+ + CO + e-
1733	!===========================================================
1734	!We divide this reaction in two
1735
1736	!0.5 CO2 + hv -> CO
1737	nb_phot = nb_phot + 1
1738
1739	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
1740
1741	!0.5 CO2 + hv -> O+ + e-
1742	nb_phot = nb_phot + 1
1743
1744	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
1745
1746	!===========================================================
1747	! CO2 + hv -> CO+ + O + e-
1748	!===========================================================
1749	!We divide this reaction in two
1750
1751	!0.5 CO2 + hv -> O
1752	nb_phot = nb_phot + 1
1753
1754	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
1755
1756	!0.5 CO2 + hv -> CO+ + e-
1757	nb_phot = nb_phot + 1
1758
1759	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
1760
1761	!===========================================================
1762	! CO2 + hv -> C+ + O2 + e-
1763	!===========================================================
1764	!We divide this reaction in two
1765
1766	!0.5 CO2 + hv -> O2
1767	nb_phot = nb_phot + 1
1768
1769	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
1770
1771	!0.5 CO2 + hv -> C+ + e-
1772	nb_phot = nb_phot + 1
1773
1774	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
1775
1776	!===========================================================
1777	! O2 + hv -> O2+ + e-
1778	!===========================================================
1779
1780	nb_phot = nb_phot + 1
1781
1782	indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
1783
1784	!===========================================================
1785	! O + hv -> O+ + e-
1786	!===========================================================
1787
1788	nb_phot = nb_phot + 1
1789
1790	indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
1791
1792	!===========================================================
1793	! NO + hv -> NO+ + e-
1794	!===========================================================
1795
1796	nb_phot = nb_phot + 1
1797
1798	indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
1799
1800	!===========================================================
1801	! CO + hv -> CO+ + e-
1802	!===========================================================
1803
1804	nb_phot = nb_phot + 1
1805
1806	indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
1807
1808	!===========================================================
1809	! CO + hv -> C+ + O + e-
1810	!===========================================================
1811	!We divide this reaction in two
1812
1813	!0.5 CO + hv -> O
1814	nb_phot = nb_phot + 1
1815
1816	indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
1817
1818	!0.5 CO + hv -> C+ + e-
1819	nb_phot = nb_phot + 1
1820
1821	indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
1822
1823	!===========================================================
1824	! N2 + hv -> N2+ + e-
1825	!===========================================================
1826
1827	nb_phot = nb_phot + 1
1828
1829	indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
1830
1831	!===========================================================
1832	! N2 + hv -> N+ + N + e-
1833	!===========================================================
1834	!We divide this reaction in two
1835
1836	!0.5 N2 + hv -> N
1837	nb_phot = nb_phot + 1
1838
1839	indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
1840
1841	!0.5 N2 + hv -> N+ + e-
1842	nb_phot = nb_phot + 1
1843
1844	indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
1845
1846	!===========================================================
1847	! N + hv -> N+ + e-
1848	!===========================================================
1849
1850	nb_phot = nb_phot + 1
1851
1852	indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
1853
1854	!===========================================================
1855	! H + hv -> H+ + e-
1856	!===========================================================
1857
1858	nb_phot = nb_phot + 1
1859
1860	indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
1861
1862	end if !ionchem
1863
1864	!===========================================================
1865	! a001 : O + O2 + CO2 -> O3 + CO2
1866	!===========================================================
1867
1868	nb_reaction_4 = nb_reaction_4 + 1
1869
1870	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
1871
1872	!===========================================================
1873	! a002 : O + O + CO2 -> O2 + CO2
1874	!===========================================================
1875
1876	nb_reaction_3 = nb_reaction_3 + 1
1877
1878	indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
1879
1880	!===========================================================
1881	! a003 : O + O3 -> O2 + O2
1882	!===========================================================
1883
1884	nb_reaction_4 = nb_reaction_4 + 1
1885
1886	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
1887
1888	!===========================================================
1889	! b001 : O(1D) + CO2 -> O + CO2
1890	!===========================================================
1891
1892	nb_phot = nb_phot + 1
1893
1894	indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
1895
1896	!===========================================================
1897	! b002 : O(1D) + H2O -> OH + OH
1898	!===========================================================
1899
1900	nb_reaction_4 = nb_reaction_4 + 1
1901
1902	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
1903
1904	!===========================================================
1905	! b003 : O(1D) + H2 -> OH + H
1906	!===========================================================
1907
1908	nb_reaction_4 = nb_reaction_4 + 1
1909
1910	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
1911
1912	!===========================================================
1913	! b004 : O(1D) + O2 -> O + O2
1914	!===========================================================
1915
1916	nb_phot = nb_phot + 1
1917
1918	indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
1919
1920	!===========================================================
1921	! b005 : O(1D) + O3 -> O2 + O2
1922	!===========================================================
1923
1924	nb_reaction_4 = nb_reaction_4 + 1
1925
1926	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
1927
1928	!===========================================================
1929	! b006 : O(1D) + O3 -> O2 + O + O
1930	!===========================================================
1931
1932	nb_reaction_4 = nb_reaction_4 + 1
1933
1934	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
1935
1936	!===========================================================
1937	! c001 : O + HO2 -> OH + O2
1938	!===========================================================
1939
1940	nb_reaction_4 = nb_reaction_4 + 1
1941
1942	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
1943
1944	!===========================================================
1945	! c002 : O + OH -> O2 + H
1946	!===========================================================
1947
1948	nb_reaction_4 = nb_reaction_4 + 1
1949
1950	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
1951
1952	!===========================================================
1953	! c003 : H + O3 -> OH + O2
1954	!===========================================================
1955
1956	nb_reaction_4 = nb_reaction_4 + 1
1957
1958	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
1959
1960	!===========================================================
1961	! c004 : H + HO2 -> OH + OH
1962	!===========================================================
1963
1964	nb_reaction_4 = nb_reaction_4 + 1
1965
1966	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
1967
1968	!===========================================================
1969	! c005 : H + HO2 -> H2 + O2
1970	!===========================================================
1971
1972	nb_reaction_4 = nb_reaction_4 + 1
1973
1974	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
1975
1976	!===========================================================
1977	! c006 : H + HO2 -> H2O + O
1978	!===========================================================
1979
1980	nb_reaction_4 = nb_reaction_4 + 1
1981
1982	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
1983
1984	!===========================================================
1985	! c007 : OH + HO2 -> H2O + O2
1986	!===========================================================
1987
1988	nb_reaction_4 = nb_reaction_4 + 1
1989
1990	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
1991
1992	!===========================================================
1993	! c008 : HO2 + HO2 -> H2O2 + O2
1994	!===========================================================
1995
1996	nb_reaction_3 = nb_reaction_3 + 1
1997
1998	indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
1999
2000	!===========================================================
2001	! c009 : OH + H2O2 -> H2O + HO2
2002	!===========================================================
2003
2004	nb_reaction_4 = nb_reaction_4 + 1
2005
2006	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
2007
2008	!===========================================================
2009	! c010 : OH + H2 -> H2O + H
2010	!===========================================================
2011
2012	nb_reaction_4 = nb_reaction_4 + 1
2013
2014	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
2015
2016	!===========================================================
2017	! c011 : H + O2 + CO2 -> HO2 + CO2
2018	!===========================================================
2019
2020	nb_reaction_4 = nb_reaction_4 + 1
2021
2022	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
2023
2024	!===========================================================
2025	! c012 : O + H2O2 -> OH + HO2
2026	!===========================================================
2027
2028	nb_reaction_4 = nb_reaction_4 + 1
2029
2030	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
2031
2032	!===========================================================
2033	! c013 : OH + OH -> H2O + O
2034	!===========================================================
2035
2036	nb_reaction_3 = nb_reaction_3 + 1
2037
2038	indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
2039
2040	!===========================================================
2041	! c014 : OH + O3 -> HO2 + O2
2042	!===========================================================
2043
2044	nb_reaction_4 = nb_reaction_4 + 1
2045
2046	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
2047
2048	!===========================================================
2049	! c015 : HO2 + O3 -> OH + O2 + O2
2050	!===========================================================
2051
2052	nb_reaction_4 = nb_reaction_4 + 1
2053
2054	indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
2055
2056	!===========================================================
2057	! c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
2058	!===========================================================
2059
2060	nb_reaction_3 = nb_reaction_3 + 1
2061
2062	indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2063
2064	!===========================================================
2065	! c017 : OH + OH + CO2 -> H2O2 + CO2
2066	!===========================================================
2067
2068	nb_reaction_3 = nb_reaction_3 + 1
2069
2070	indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
2071
2072	!===========================================================
2073	! c018 : H + H + CO2 -> H2 + CO2
2074	!===========================================================
2075
2076	nb_reaction_3 = nb_reaction_3 + 1
2077
2078	indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
2079
2080	!===========================================================
2081	! d001 : NO2 + O -> NO + O2
2082	!===========================================================
2083
2084	nb_reaction_4 = nb_reaction_4 + 1
2085
2086	indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
2087
2088	!===========================================================
2089	! d002 : NO + O3 -> NO2 + O2
2090	!===========================================================
2091
2092	nb_reaction_4 = nb_reaction_4 + 1
2093
2094	indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
2095
2096	!===========================================================
2097	! d003 : NO + HO2 -> NO2 + OH
2098	!===========================================================
2099
2100	nb_reaction_4 = nb_reaction_4 + 1
2101
2102	indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
2103
2104	!===========================================================
2105	! d004 : N + NO -> N2 + O
2106	!===========================================================
2107
2108	nb_reaction_4 = nb_reaction_4 + 1
2109
2110	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
2111
2112	!===========================================================
2113	! d005 : N + O2 -> NO + O
2114	!===========================================================
2115
2116	nb_reaction_4 = nb_reaction_4 + 1
2117
2118	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
2119
2120	!===========================================================
2121	! d006 : NO2 + H -> NO + OH
2122	!===========================================================
2123
2124	nb_reaction_4 = nb_reaction_4 + 1
2125
2126	indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
2127
2128	!===========================================================
2129	! d007 : N + O -> NO
2130	!===========================================================
2131
2132	nb_reaction_4 = nb_reaction_4 + 1
2133
2134	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
2135
2136	!===========================================================
2137	! d008 : N + HO2 -> NO + OH
2138	!===========================================================
2139
2140	nb_reaction_4 = nb_reaction_4 + 1
2141
2142	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
2143
2144	!===========================================================
2145	! d009 : N + OH -> NO + H
2146	!===========================================================
2147
2148	nb_reaction_4 = nb_reaction_4 + 1
2149
2150	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
2151
2152	!===========================================================
2153	! d010 : N(2D) + O -> N + O
2154	!===========================================================
2155
2156	nb_phot = nb_phot + 1
2157
2158	indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
2159
2160	!===========================================================
2161	! d011 : N(2D) + N2 -> N + N2
2162	!===========================================================
2163
2164	nb_phot = nb_phot + 1
2165
2166	indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
2167
2168	!===========================================================
2169	! d012 : N(2D) + CO2 -> NO + CO
2170	!===========================================================
2171
2172	nb_reaction_4 = nb_reaction_4 + 1
2173
2174	indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
2175
2176	!===========================================================
2177	! e001 : CO + OH -> CO2 + H
2178	!===========================================================
2179
2180	nb_reaction_4 = nb_reaction_4 + 1
2181
2182	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
2183
2184	!===========================================================
2185	! e002 : CO + O + M -> CO2 + M
2186	!===========================================================
2187
2188	nb_reaction_4 = nb_reaction_4 + 1
2189
2190	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
2191
2192	!Only if ion chemistry
2193	if (ionchem) then
2194
2195	!===========================================================
2196	! i001 : CO2+ + O2 -> O2+ + CO2
2197	!===========================================================
2198
2199	nb_reaction_4 = nb_reaction_4 + 1
2200
2201	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
2202
2203	!===========================================================
2204	! i002 : CO2+ + O -> O+ + CO2
2205	!===========================================================
2206
2207	nb_reaction_4 = nb_reaction_4 + 1
2208
2209	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
2210
2211	!===========================================================
2212	! i003 : CO2+ + O -> O2+ + CO
2213	!===========================================================
2214
2215	nb_reaction_4 = nb_reaction_4 + 1
2216
2217	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
2218
2219	!===========================================================
2220	! i004 : O2+ + e- -> O + O
2221	!===========================================================
2222
2223	nb_reaction_4 = nb_reaction_4 + 1
2224
2225	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
2226
2227	!===========================================================
2228	! i005 : O+ + CO2 -> O2+ + CO
2229	!===========================================================
2230
2231	nb_reaction_4 = nb_reaction_4 + 1
2232
2233	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
2234
2235	!===========================================================
2236	! i006 : CO2+ + e -> CO + O
2237	!===========================================================
2238
2239	nb_reaction_4 = nb_reaction_4 + 1
2240
2241	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
2242
2243	!===========================================================
2244	! i007 : CO2+ + NO -> NO+ + CO2
2245	!===========================================================
2246
2247	nb_reaction_4 = nb_reaction_4 + 1
2248
2249	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
2250
2251	!===========================================================
2252	! i008 : O2+ + NO -> NO+ + O2
2253	!===========================================================
2254
2255	nb_reaction_4 = nb_reaction_4 + 1
2256
2257	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
2258
2259	!===========================================================
2260	! i009 : O2+ + N2 -> NO+ + NO
2261	!===========================================================
2262
2263	nb_reaction_4 = nb_reaction_4 + 1
2264
2265	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
2266
2267	!===========================================================
2268	! i010 : O2+ + N -> NO+ + O
2269	!===========================================================
2270
2271	nb_reaction_4 = nb_reaction_4 + 1
2272
2273	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
2274
2275	!===========================================================
2276	! i011 : O+ + N2 -> NO+ + N
2277	!===========================================================
2278
2279	nb_reaction_4 = nb_reaction_4 + 1
2280
2281	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
2282
2283	!===========================================================
2284	! i012 : NO+ + e -> N + O
2285	!===========================================================
2286
2287	nb_reaction_4 = nb_reaction_4 + 1
2288
2289	indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
2290
2291	!===========================================================
2292	! i013 : CO+ + CO2 -> CO2+ + CO
2293	!===========================================================
2294
2295	nb_reaction_4 = nb_reaction_4 + 1
2296
2297	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
2298
2299	!===========================================================
2300	! i014 : CO+ + O -> O+ + CO
2301	!===========================================================
2302
2303	nb_reaction_4 = nb_reaction_4 + 1
2304
2305	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
2306
2307	!===========================================================
2308	! i015 : C+ + CO2 -> CO+ + CO
2309	!===========================================================
2310
2311	nb_reaction_4 = nb_reaction_4 + 1
2312
2313	indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
2314
2315	!===========================================================
2316	! i016 : N2+ + CO2 -> CO2+ + N2
2317	!===========================================================
2318
2319	nb_reaction_4 = nb_reaction_4 + 1
2320
2321	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
2322
2323	!===========================================================
2324	! i017 : N2+ + O -> NO+ + N
2325	!===========================================================
2326
2327	nb_reaction_4 = nb_reaction_4 + 1
2328
2329	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
2330
2331	!===========================================================
2332	! i018 : N2+ + CO -> CO+ + N2
2333	!===========================================================
2334
2335	nb_reaction_4 = nb_reaction_4 + 1
2336
2337	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
2338
2339	!===========================================================
2340	! i019 : N2+ + e -> N + N
2341	!===========================================================
2342
2343	nb_reaction_4 = nb_reaction_4 + 1
2344
2345	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
2346
2347	!===========================================================
2348	! i020 : N2+ + O -> O+ + N2
2349	!===========================================================
2350
2351	nb_reaction_4 = nb_reaction_4 + 1
2352
2353	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
2354
2355	!===========================================================
2356	! i021 : N+ + CO2 -> CO2+ + N
2357	!===========================================================
2358
2359	nb_reaction_4 = nb_reaction_4 + 1
2360
2361	indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
2362
2363	!===========================================================
2364	! i022 : CO+ + H -> H+ + CO
2365	!===========================================================
2366
2367	nb_reaction_4 = nb_reaction_4 + 1
2368
2369	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
2370
2371	!===========================================================
2372	! i023 : O+ + H -> H+ + O
2373	!===========================================================
2374
2375	nb_reaction_4 = nb_reaction_4 + 1
2376
2377	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
2378
2379	!===========================================================
2380	! i024 : H+ + O -> O+ + H
2381	!===========================================================
2382
2383	nb_reaction_4 = nb_reaction_4 + 1
2384
2385	indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
2386
2387	!===========================================================
2388	! i025 : CO2+ + H2 -> HCO2+ + H
2389	!===========================================================
2390
2391	nb_reaction_4 = nb_reaction_4 + 1
2392
2393	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
2394
2395	!===========================================================
2396	! i026 : HCO2+ + e -> H + CO2
2397	!===========================================================
2398
2399	nb_reaction_4 = nb_reaction_4 + 1
2400
2401	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
2402
2403	!===========================================================
2404	! i027 : HCO2+ + e -> H + O + CO
2405	!===========================================================
2406	!We divide this reaction in two
2407
2408	!0.5HCO2+ + 0.5e -> H
2409
2410	nb_reaction_4 = nb_reaction_4 + 1
2411
2412	indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
2413
2414	!0.5 HCO2+ + 0.5 e -> O + CO
2415
2416	nb_reaction_4 = nb_reaction_4 + 1
2417
2418	indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
2419
2420	!===========================================================
2421	! i029 : HCO2+ + e -> OH + CO
2422	!===========================================================
2423
2424	nb_reaction_4 = nb_reaction_4 + 1
2425
2426	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
2427
2428	end if !ionchem
2429
2430	!===========================================================
2431	! h001: HO2 + ice -> products
2432	! treated as
2433	! HO2 -> 0.5 H2O + 0.75 O2
2434	!===========================================================
2435
2436	nb_phot = nb_phot + 1
2437
2438	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
2439
2440	!===========================================================
2441	! h002: OH + ice -> products
2442	! treated as
2443	! OH -> 0.5 H2O + 0.25 O2
2444	!===========================================================
2445
2446	nb_phot = nb_phot + 1
2447
2448	indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2)
2449
2450	!===========================================================
2451	! h003: H2O2 + ice -> products
2452	! treated as
2453	! H2O2 -> H2O + 0.5 O2
2454	!===========================================================
2455
2456	nb_phot = nb_phot + 1
2457
2458	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
2459
2460	!===========================================================
2461	! h004: HO2 + dust -> products
2462	! treated as
2463	! HO2 -> 0.5 H2O + 0.75 O2
2464	!===========================================================
2465
2466	nb_phot = nb_phot + 1
2467
2468	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
2469
2470	!===========================================================
2471	! h005: H2O2 + dust -> products
2472	! treated as
2473	! H2O2 -> H2O + 0.5 O2
2474	!===========================================================
2475
2476	nb_phot = nb_phot + 1
2477
2478	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
2479
2480	!===========================================================
2481	! check dimensions
2482	!===========================================================
2483
2484	print*, 'nb_phot = ', nb_phot
2485	print*, 'nb_reaction_4 = ', nb_reaction_4
2486	print*, 'nb_reaction_3 = ', nb_reaction_3
2487
2488	if ((nb_phot /= nb_phot_max) .or. &
2489	(nb_reaction_3 /= nb_reaction_3_max) .or. &
2490	(nb_reaction_4 /= nb_reaction_4_max)) then
2491	print*, 'wrong dimensions in indice'
2492	stop
2493	end if
2494
2495	end subroutine indice
2496
2497	!*****************************************************************
2498
2499	subroutine gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,&
2500	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2501	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2502	i_n, i_n2d, i_no, i_no2, i_n2, &
2503	i_co2plus, i_oplus, i_o2plus, i_noplus, &
2504	i_coplus, i_cplus, i_n2plus, i_nplus, &
2505	i_hplus, i_hco2plus, i_elec, dens, rm, c)
2506
2507	!*****************************************************************
2508
2509	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, &
2510	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh, &
2511	& igcm_ho2, igcm_h2o2, igcm_h2o_vap, &
2512	& igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2,&
2513	& igcm_co2plus, igcm_oplus, igcm_o2plus, &
2514	& igcm_noplus, igcm_coplus, igcm_cplus, &
2515	& igcm_n2plus, igcm_nplus, igcm_hplus, &
2516	& igcm_hco2plus, igcm_elec
2517
2518	implicit none
2519
2520	#include "callkeys.h"
2521
2522	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2523	! input:
2524	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2525
2526	integer, intent(in) :: nlayer ! number of atmospheric layers
2527	integer, intent(in) :: nq ! number of tracers in the gcm
2528	logical, intent(in) :: ionchem
2529	integer :: nesp ! number of species in the chemistry
2530	integer :: lswitch ! interface level between chemistries
2531
2532	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2533	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2534	i_n, i_n2d, i_no, i_no2, i_n2, &
2535	i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
2536	i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
2537
2538	real :: zycol(nlayer,nq) ! volume mixing ratios in the gcm
2539	real :: dens(nlayer) ! total number density (molecule.cm-3)
2540
2541	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2542	! output:
2543	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2544
2545	real, dimension(nlayer,nesp) :: rm ! volume mixing ratios
2546	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
2547
2548	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2549	! local:
2550	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2551
2552	integer :: l, iesp
2553	logical,save :: firstcall = .true.
2554
2555	! first call initializations
2556
2557	if (firstcall) then
2558
2559	! identify the indexes of the tracers we need
2560
2561	if (igcm_co2 == 0) then
2562	write(,) "gcmtochim: Error; no CO2 tracer !!!"
2563	stop
2564	endif
2565	if (igcm_co == 0) then
2566	write(,) "gcmtochim: Error; no CO tracer !!!"
2567	stop
2568	end if
2569	if (igcm_o == 0) then
2570	write(,) "gcmtochim: Error; no O tracer !!!"
2571	stop
2572	end if
2573	if (igcm_o1d == 0) then
2574	write(,) "gcmtochim: Error; no O1D tracer !!!"
2575	stop
2576	end if
2577	if (igcm_o2 == 0) then
2578	write(,) "gcmtochim: Error; no O2 tracer !!!"
2579	stop
2580	end if
2581	if (igcm_o3 == 0) then
2582	write(,) "gcmtochim: Error; no O3 tracer !!!"
2583	stop
2584	end if
2585	if (igcm_h == 0) then
2586	write(,) "gcmtochim: Error; no H tracer !!!"
2587	stop
2588	end if
2589	if (igcm_h2 == 0) then
2590	write(,) "gcmtochim: Error; no H2 tracer !!!"
2591	stop
2592	end if
2593	if (igcm_oh == 0) then
2594	write(,) "gcmtochim: Error; no OH tracer !!!"
2595	stop
2596	end if
2597	if (igcm_ho2 == 0) then
2598	write(,) "gcmtochim: Error; no HO2 tracer !!!"
2599	stop
2600	end if
2601	if (igcm_h2o2 == 0) then
2602	write(,) "gcmtochim: Error; no H2O2 tracer !!!"
2603	stop
2604	end if
2605	if (igcm_n == 0) then
2606	write(,) "gcmtochim: Error; no N tracer !!!"
2607	stop
2608	end if
2609	if (igcm_n2d == 0) then
2610	write(,) "gcmtochim: Error; no N2D tracer !!!"
2611	stop
2612	end if
2613	if (igcm_no == 0) then
2614	write(,) "gcmtochim: Error; no NO tracer !!!"
2615	stop
2616	end if
2617	if (igcm_no2 == 0) then
2618	write(,) "gcmtochim: Error; no NO2 tracer !!!"
2619	stop
2620	end if
2621	if (igcm_n2 == 0) then
2622	write(,) "gcmtochim: Error; no N2 tracer !!!"
2623	stop
2624	end if
2625	if (igcm_h2o_vap == 0) then
2626	write(,) "gcmtochim: Error; no water vapor tracer !!!"
2627	stop
2628	end if
2629	if (ionchem) then
2630	if (igcm_co2plus == 0) then
2631	write(,) "gcmtochim: Error; no CO2+ tracer !!!"
2632	stop
2633	end if
2634	if (igcm_oplus == 0) then
2635	write(,) "gcmtochim: Error; no O+ tracer !!!"
2636	stop
2637	end if
2638	if (igcm_o2plus == 0) then
2639	write(,) "gcmtochim: Error; no O2+ tracer !!!"
2640	stop
2641	end if
2642	if (igcm_noplus == 0) then
2643	write(,) "gcmtochim: Error; no NO+ tracer !!!"
2644	stop
2645	endif
2646	if (igcm_coplus == 0) then
2647	write(,) "gcmtochim: Error; no CO+ tracer !!!"
2648	stop
2649	endif
2650	if (igcm_cplus == 0) then
2651	write(,) "gcmtochim: Error; no C+ tracer !!!"
2652	stop
2653	endif
2654	if (igcm_n2plus == 0) then
2655	write(,) "gcmtochim: Error; no N2+ tracer !!!"
2656	stop
2657	endif
2658	if (igcm_nplus == 0) then
2659	write(,) "gcmtochim: Error; no N+ tracer !!!"
2660	stop
2661	endif
2662	if (igcm_hplus == 0) then
2663	write(,) "gcmtochim: Error; no H+ tracer !!!"
2664	stop
2665	endif
2666	if (igcm_hco2plus == 0) then
2667	write(,) "gcmtochim: Error; no HCO2+ tracer !!!"
2668	stop
2669	endif
2670	if (igcm_elec == 0) then
2671	write(,) "gcmtochim: Error; no e- tracer !!!"
2672	stop
2673	end if
2674	end if ! ionchem
2675	firstcall = .false.
2676	end if ! of if (firstcall)
2677
2678	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2679	! initialise mixing ratios
2680	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2681
2682	do l = 1,nlayer
2683	rm(l,i_co2) = zycol(l, igcm_co2)
2684	rm(l,i_co) = zycol(l, igcm_co)
2685	rm(l,i_o) = zycol(l, igcm_o)
2686	rm(l,i_o1d) = zycol(l, igcm_o1d)
2687	rm(l,i_o2) = zycol(l, igcm_o2)
2688	rm(l,i_o3) = zycol(l, igcm_o3)
2689	rm(l,i_h) = zycol(l, igcm_h)
2690	rm(l,i_h2) = zycol(l, igcm_h2)
2691	rm(l,i_oh) = zycol(l, igcm_oh)
2692	rm(l,i_ho2) = zycol(l, igcm_ho2)
2693	rm(l,i_h2o2) = zycol(l, igcm_h2o2)
2694	rm(l,i_h2o) = zycol(l, igcm_h2o_vap)
2695	rm(l,i_n) = zycol(l, igcm_n)
2696	rm(l,i_n2d) = zycol(l, igcm_n2d)
2697	rm(l,i_no) = zycol(l, igcm_no)
2698	rm(l,i_no2) = zycol(l, igcm_no2)
2699	rm(l,i_n2) = zycol(l, igcm_n2)
2700	end do
2701
2702	if (ionchem) then
2703	do l = 1,nlayer
2704	rm(l,i_co2plus) = zycol(l, igcm_co2plus)
2705	rm(l,i_oplus) = zycol(l, igcm_oplus)
2706	rm(l,i_o2plus) = zycol(l, igcm_o2plus)
2707	rm(l,i_noplus) = zycol(l, igcm_noplus)
2708	rm(l,i_coplus) = zycol(l, igcm_coplus)
2709	rm(l,i_cplus) = zycol(l, igcm_cplus)
2710	rm(l,i_n2plus) = zycol(l, igcm_n2plus)
2711	rm(l,i_nplus) = zycol(l, igcm_nplus)
2712	rm(l,i_hplus) = zycol(l, igcm_hplus)
2713	rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
2714	rm(l,i_elec) = zycol(l, igcm_elec)
2715	end do
2716	end if
2717
2718	where (rm(:,:) < 1.e-30)
2719	rm(:,:) = 0.
2720	end where
2721
2722	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2723	! initialise number densities
2724	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2725
2726	do iesp = 1,nesp
2727	do l = 1,nlayer
2728	c(l,iesp) = rm(l,iesp)*dens(l)
2729	end do
2730	end do
2731
2732	end subroutine gcmtochim
2733
2734	!*****************************************************************
2735
2736	subroutine chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp, &
2737	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2738	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2739	i_n, i_n2d, i_no, i_no2, i_n2, &
2740	i_co2plus, i_oplus, i_o2plus, i_noplus, &
2741	i_coplus, i_cplus, i_n2plus, i_nplus, &
2742	i_hplus, i_hco2plus, i_elec, dens, c)
2743
2744	!*****************************************************************
2745
2746	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, &
2747	igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh, &
2748	igcm_ho2, igcm_h2o2, igcm_h2o_vap, &
2749	igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2, &
2750	igcm_co2plus, igcm_oplus, igcm_o2plus, &
2751	igcm_noplus, igcm_coplus, igcm_cplus, &
2752	igcm_n2plus, igcm_nplus, igcm_hplus, &
2753	igcm_hco2plus, igcm_elec
2754
2755	implicit none
2756
2757	#include "callkeys.h"
2758
2759	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2760	! inputs:
2761	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2762
2763	integer, intent(in) :: nlayer ! number of atmospheric layers
2764	integer, intent(in) :: nq ! number of tracers in the gcm
2765	logical, intent(in) :: ionchem
2766	integer :: nesp ! number of species in the chemistry
2767	integer :: lswitch ! interface level between chemistries
2768	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2769	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2770	i_n, i_n2d, i_no, i_no2, i_n2, &
2771	i_co2plus, i_oplus, i_o2plus, i_noplus, &
2772	i_coplus, i_cplus, i_n2plus, i_nplus, &
2773	i_hplus, i_hco2plus, i_elec
2774
2775	real :: dens(nlayer) ! total number density (molecule.cm-3)
2776	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
2777
2778	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2779	! output:
2780	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2781
2782	real zycol(nlayer,nq) ! volume mixing ratios in the gcm
2783
2784	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2785	! local:
2786	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2787
2788	integer l
2789
2790	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2791	! save mixing ratios for the gcm
2792	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2793
2794	do l = 1,lswitch-1
2795	zycol(l, igcm_co2) = c(l,i_co2)/dens(l)
2796	zycol(l, igcm_co) = c(l,i_co)/dens(l)
2797	zycol(l, igcm_o) = c(l,i_o)/dens(l)
2798	zycol(l, igcm_o1d) = c(l,i_o1d)/dens(l)
2799	zycol(l, igcm_o2) = c(l,i_o2)/dens(l)
2800	zycol(l, igcm_o3) = c(l,i_o3)/dens(l)
2801	zycol(l, igcm_h) = c(l,i_h)/dens(l)
2802	zycol(l, igcm_h2) = c(l,i_h2)/dens(l)
2803	zycol(l, igcm_oh) = c(l,i_oh)/dens(l)
2804	zycol(l, igcm_ho2) = c(l,i_ho2)/dens(l)
2805	zycol(l, igcm_h2o2) = c(l,i_h2o2)/dens(l)
2806	zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
2807	zycol(l, igcm_n) = c(l,i_n)/dens(l)
2808	zycol(l, igcm_n2d) = c(l,i_n2d)/dens(l)
2809	zycol(l, igcm_no) = c(l,i_no)/dens(l)
2810	zycol(l, igcm_no2) = c(l,i_no2)/dens(l)
2811	zycol(l, igcm_n2) = c(l,i_n2)/dens(l)
2812	end do
2813
2814	if (ionchem) then
2815	do l = 1,lswitch-1
2816	zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
2817	zycol(l, igcm_oplus) = c(l,i_oplus)/dens(l)
2818	zycol(l, igcm_o2plus) = c(l,i_o2plus)/dens(l)
2819	zycol(l, igcm_noplus) = c(l,i_noplus)/dens(l)
2820	zycol(l, igcm_coplus) = c(l,i_coplus)/dens(l)
2821	zycol(l, igcm_cplus) = c(l,i_cplus)/dens(l)
2822	zycol(l, igcm_n2plus) = c(l,i_n2plus)/dens(l)
2823	zycol(l, igcm_nplus) = c(l,i_nplus)/dens(l)
2824	zycol(l, igcm_hplus) = c(l,i_hplus)/dens(l)
2825	zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
2826	zycol(l, igcm_elec) = c(l,i_elec)/dens(l)
2827	end do
2828	end if
2829
2830	end subroutine chimtogcm
2831
2832	end subroutine photochemistry

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

Download in other formats:

Original Format