Context Navigation

photochemistry.F90 @ 2275

Last change on this file since 2275 was 2273, checked in by emillour, 5 years ago
Mars GCM: A first step towards water ion chemistry. Add a new reaction and modify reaction rates for HCO2+ FGG
File size: 87.6 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, i030, &
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(:) = 8.1e-7((300./t_elect(:))*0.64)
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 -> OH + CO
1349
1350	! UMIST
1351
1352	i029(:) = 3.2e-7((300./t_elect(:))*0.64)
1353
1354	nb_reaction_4 = nb_reaction_4 + 1
1355	v_4(:,nb_reaction_4) = i029(:)
1356
1357	!--- i030: HCO2+ + e -> H + CO2
1358
1359	i030(:) = 6.0e-8((300./t_elect(:))*0.64)
1360	nb_reaction_4 = nb_reaction_4 + 1
1361	v_4(:,nb_reaction_4) = i030(:)
1362
1363	end if !ionchem
1364
1365	!----------------------------------------------------------------------
1366	! heterogeneous chemistry
1367	!----------------------------------------------------------------------
1368
1369	if (hetero_ice) then
1370
1371	! k = (surfacevgamma)/4 (s-1)
1372	! v = 100sqrt(8rt/(pim)) (cm s-1)
1373
1374	!--- h001: ho2 + ice -> products
1375
1376	! cooper and abbatt, 1996: gamma = 0.025
1377
1378	gam = 0.025
1379	h001(:) = surfice1d(:) &
1380	100.sqrt(8.8.31t(:)/(33.e-3pi))gam/4.
1381
1382	! h002: oh + ice -> products
1383
1384	! cooper and abbatt, 1996: gamma = 0.03
1385
1386	gam = 0.03
1387	h002(:) = surfice1d(:) &
1388	100.sqrt(8.8.31t(:)/(17.e-3pi))gam/4.
1389
1390	!--- h003: h2o2 + ice -> products
1391
1392	! gamma = 0. test value
1393
1394	gam = 0.
1395	h003(:) = surfice1d(:) &
1396	100.sqrt(8.8.31t(:)/(34.e-3pi))gam/4.
1397	else
1398	h001(:) = 0.
1399	h002(:) = 0.
1400	h003(:) = 0.
1401	end if
1402
1403	nb_phot = nb_phot + 1
1404	v_phot(:,nb_phot) = h001(:)
1405
1406	nb_phot = nb_phot + 1
1407	v_phot(:,nb_phot) = h002(:)
1408
1409	nb_phot = nb_phot + 1
1410	v_phot(:,nb_phot) = h003(:)
1411
1412	if (hetero_dust) then
1413
1414	!--- h004: ho2 + dust -> products
1415
1416	! jacob, 2000: gamma = 0.2
1417	! see dereus et al., atm. chem. phys., 2005
1418
1419	gam = 0.2
1420	h004(:) = surfdust1d(:) &
1421	100.sqrt(8.8.31t(:)/(33.e-3pi))gam/4.
1422
1423	!--- h005: h2o2 + dust -> products
1424
1425	! gamma = 5.e-4
1426	! see dereus et al., atm. chem. phys., 2005
1427
1428	gam = 5.e-4
1429	h005(:) = surfdust1d(:) &
1430	100.sqrt(8.8.31t(:)/(34.e-3pi))gam/4.
1431	else
1432	h004(:) = 0.
1433	h005(:) = 0.
1434	end if
1435
1436	nb_phot = nb_phot + 1
1437	v_phot(:,nb_phot) = h004(:)
1438
1439	nb_phot = nb_phot + 1
1440	v_phot(:,nb_phot) = h005(:)
1441
1442	end subroutine reactionrates
1443
1444	!======================================================================
1445
1446	subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, &
1447	nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
1448	v_phot, v_3, v_4)
1449
1450	!======================================================================
1451	! filling of the jacobian matrix
1452	!======================================================================
1453
1454	use types_asis
1455
1456	implicit none
1457
1458	! input
1459
1460	integer :: ilev ! level index
1461	integer :: nesp ! number of species in the chemistry
1462	integer, intent(in) :: nlayer ! number of atmospheric layers
1463	integer, intent(in) :: nb_reaction_3_max
1464	! number of quadratic reactions
1465	integer, intent(in) :: nb_reaction_4_max
1466	! number of bimolecular reactions
1467	integer, intent(in) :: nb_phot_max
1468	! number of processes treated numerically as photodissociations
1469
1470	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities
1471	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
1472	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
1473	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
1474
1475	! output
1476
1477	real (kind = 8), dimension(nesp,nesp), intent(out) :: mat ! matrix
1478	real (kind = 8), dimension(nesp), intent(out) :: prod, loss
1479
1480	! local
1481
1482	integer :: iesp
1483	integer :: ind_phot_2,ind_phot_4,ind_phot_6
1484	integer :: ind_3_2,ind_3_4,ind_3_6
1485	integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
1486	integer :: iphot,i3,i4
1487
1488	real(kind = 8) :: eps, eps_4 ! implicit/explicit coefficient
1489
1490	! initialisations
1491
1492	mat(:,:) = 0.
1493	prod(:) = 0.
1494	loss(:) = 0.
1495
1496	! photodissociations
1497	! or reactions a + c -> b + c
1498	! or reactions a + ice -> b + c
1499
1500	do iphot = 1,nb_phot_max
1501
1502	ind_phot_2 = indice_phot(iphot)%z2
1503	ind_phot_4 = indice_phot(iphot)%z4
1504	ind_phot_6 = indice_phot(iphot)%z6
1505
1506	mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
1507	mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
1508	mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
1509
1510	loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
1511	prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3v_phot(ilev,iphot)c(ilev,ind_phot_2)
1512	prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5v_phot(ilev,iphot)c(ilev,ind_phot_2)
1513
1514	end do
1515
1516	! reactions a + a -> b + c
1517
1518	do i3 = 1,nb_reaction_3_max
1519
1520	ind_3_2 = indice_3(i3)%z2
1521	ind_3_4 = indice_3(i3)%z4
1522	ind_3_6 = indice_3(i3)%z6
1523
1524	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)
1525	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)
1526	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)
1527
1528	loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1v_3(ilev,i3)c(ilev,ind_3_2)
1529	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)
1530	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)
1531
1532	end do
1533
1534	! reactions a + b -> c + d
1535
1536	eps = 1.d-10
1537
1538	do i4 = 1,nb_reaction_4_max
1539
1540	ind_4_2 = indice_4(i4)%z2
1541	ind_4_4 = indice_4(i4)%z4
1542	ind_4_6 = indice_4(i4)%z6
1543	ind_4_8 = indice_4(i4)%z8
1544
1545	eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
1546	eps_4 = min(eps_4,1.0)
1547
1548	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)
1549	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)
1550	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)
1551	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)
1552	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)
1553	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)
1554	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)
1555	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)
1556
1557
1558	loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1v_4(ilev,i4)c(ilev,ind_4_4)
1559	loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3v_4(ilev,i4)c(ilev,ind_4_2)
1560	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)
1561	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)
1562
1563	end do
1564
1565	end subroutine fill_matrix
1566
1567	!================================================================
1568
1569	subroutine indice(nb_reaction_3_max, nb_reaction_4_max, &
1570	nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
1571	i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
1572	i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, &
1573	i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
1574	i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
1575
1576	!================================================================
1577	! set the "indice" arrays used to fill the jacobian matrix !
1578	!----------------------------------------------------------------
1579	! reaction type !
1580	!----------------------------------------------------------------
1581	! A + hv --> B + C photolysis indice_phot !
1582	! A + B --> C + D bimolecular indice_4 !
1583	! A + A --> B + C quadratic indice_3 !
1584	! A + C --> B + C quenching indice_phot !
1585	! A + ice --> B + C heterogeneous indice_phot !
1586	!================================================================
1587
1588	use types_asis
1589
1590	implicit none
1591
1592	! input
1593
1594	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
1595	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
1596	i_n, i_n2d, i_no, i_no2, i_n2, &
1597	i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
1598	i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
1599	integer, intent(in) :: nb_reaction_3_max
1600	! number of quadratic reactions
1601	integer, intent(in) :: nb_reaction_4_max
1602	! number of bimolecular reactions
1603	integer, intent(in) :: nb_phot_max
1604	! number of processes treated numerically as photodissociations
1605	logical, intent(in) :: ionchem
1606
1607	! local
1608
1609	integer :: nb_phot, nb_reaction_3, nb_reaction_4
1610	integer :: i_dummy
1611
1612	allocate (indice_phot(nb_phot_max))
1613	allocate (indice_3(nb_reaction_3_max))
1614	allocate (indice_4(nb_reaction_4_max))
1615
1616	i_dummy = 1
1617
1618	nb_phot = 0
1619	nb_reaction_3 = 0
1620	nb_reaction_4 = 0
1621
1622	!===========================================================
1623	! O2 + hv -> O + O
1624	!===========================================================
1625
1626	nb_phot = nb_phot + 1
1627
1628	indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy)
1629
1630	!===========================================================
1631	! O2 + hv -> O + O(1D)
1632	!===========================================================
1633
1634	nb_phot = nb_phot + 1
1635
1636	indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d)
1637
1638	!===========================================================
1639	! CO2 + hv -> CO + O
1640	!===========================================================
1641
1642	nb_phot = nb_phot + 1
1643
1644	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o)
1645
1646	!===========================================================
1647	! CO2 + hv -> CO + O(1D)
1648	!===========================================================
1649
1650	nb_phot = nb_phot + 1
1651
1652	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d)
1653
1654	!===========================================================
1655	! O3 + hv -> O2 + O(1D)
1656	!===========================================================
1657
1658	nb_phot = nb_phot + 1
1659
1660	indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d)
1661
1662	!===========================================================
1663	! O3 + hv -> O2 + O
1664	!===========================================================
1665
1666	nb_phot = nb_phot + 1
1667
1668	indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o)
1669
1670	!===========================================================
1671	! H2O + hv -> H + OH
1672	!===========================================================
1673
1674	nb_phot = nb_phot + 1
1675
1676	indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh)
1677
1678	!===========================================================
1679	! H2O2 + hv -> OH + OH
1680	!===========================================================
1681
1682	nb_phot = nb_phot + 1
1683
1684	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy)
1685
1686	!===========================================================
1687	! HO2 + hv -> OH + O
1688	!===========================================================
1689
1690	nb_phot = nb_phot + 1
1691
1692	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o)
1693
1694	!===========================================================
1695	! H2 + hv -> H + H
1696	!===========================================================
1697
1698	nb_phot = nb_phot + 1
1699
1700	indice_phot(nb_phot) = z3spec(1.0, i_h2, 1.0, i_h, 1.0, i_h)
1701
1702	!===========================================================
1703	! NO + hv -> N + O
1704	!===========================================================
1705
1706	nb_phot = nb_phot + 1
1707
1708	indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o)
1709
1710	!===========================================================
1711	! NO2 + hv -> NO + O
1712	!===========================================================
1713
1714	nb_phot = nb_phot + 1
1715
1716	indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o)
1717
1718	!===========================================================
1719	! N2 + hv -> N + N
1720	!===========================================================
1721
1722	nb_phot = nb_phot + 1
1723
1724	indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n)
1725
1726	!Only if ion chemistry included
1727	if (ionchem) then
1728
1729	!===========================================================
1730	! CO2 + hv -> CO2+ + e-
1731	!===========================================================
1732
1733	nb_phot = nb_phot + 1
1734
1735	indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
1736
1737	!===========================================================
1738	! CO2 + hv -> O+ + CO + e-
1739	!===========================================================
1740	!We divide this reaction in two
1741
1742	!0.5 CO2 + hv -> CO
1743	nb_phot = nb_phot + 1
1744
1745	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
1746
1747	!0.5 CO2 + hv -> O+ + e-
1748	nb_phot = nb_phot + 1
1749
1750	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
1751
1752	!===========================================================
1753	! CO2 + hv -> CO+ + O + e-
1754	!===========================================================
1755	!We divide this reaction in two
1756
1757	!0.5 CO2 + hv -> O
1758	nb_phot = nb_phot + 1
1759
1760	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
1761
1762	!0.5 CO2 + hv -> CO+ + e-
1763	nb_phot = nb_phot + 1
1764
1765	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
1766
1767	!===========================================================
1768	! CO2 + hv -> C+ + O2 + e-
1769	!===========================================================
1770	!We divide this reaction in two
1771
1772	!0.5 CO2 + hv -> O2
1773	nb_phot = nb_phot + 1
1774
1775	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
1776
1777	!0.5 CO2 + hv -> C+ + e-
1778	nb_phot = nb_phot + 1
1779
1780	indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
1781
1782	!===========================================================
1783	! O2 + hv -> O2+ + e-
1784	!===========================================================
1785
1786	nb_phot = nb_phot + 1
1787
1788	indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
1789
1790	!===========================================================
1791	! O + hv -> O+ + e-
1792	!===========================================================
1793
1794	nb_phot = nb_phot + 1
1795
1796	indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
1797
1798	!===========================================================
1799	! NO + hv -> NO+ + e-
1800	!===========================================================
1801
1802	nb_phot = nb_phot + 1
1803
1804	indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
1805
1806	!===========================================================
1807	! CO + hv -> CO+ + e-
1808	!===========================================================
1809
1810	nb_phot = nb_phot + 1
1811
1812	indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
1813
1814	!===========================================================
1815	! CO + hv -> C+ + O + e-
1816	!===========================================================
1817	!We divide this reaction in two
1818
1819	!0.5 CO + hv -> O
1820	nb_phot = nb_phot + 1
1821
1822	indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
1823
1824	!0.5 CO + hv -> C+ + e-
1825	nb_phot = nb_phot + 1
1826
1827	indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
1828
1829	!===========================================================
1830	! N2 + hv -> N2+ + e-
1831	!===========================================================
1832
1833	nb_phot = nb_phot + 1
1834
1835	indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
1836
1837	!===========================================================
1838	! N2 + hv -> N+ + N + e-
1839	!===========================================================
1840	!We divide this reaction in two
1841
1842	!0.5 N2 + hv -> N
1843	nb_phot = nb_phot + 1
1844
1845	indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
1846
1847	!0.5 N2 + hv -> N+ + e-
1848	nb_phot = nb_phot + 1
1849
1850	indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
1851
1852	!===========================================================
1853	! N + hv -> N+ + e-
1854	!===========================================================
1855
1856	nb_phot = nb_phot + 1
1857
1858	indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
1859
1860	!===========================================================
1861	! H + hv -> H+ + e-
1862	!===========================================================
1863
1864	nb_phot = nb_phot + 1
1865
1866	indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
1867
1868	end if !ionchem
1869
1870	!===========================================================
1871	! a001 : O + O2 + CO2 -> O3 + CO2
1872	!===========================================================
1873
1874	nb_reaction_4 = nb_reaction_4 + 1
1875
1876	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
1877
1878	!===========================================================
1879	! a002 : O + O + CO2 -> O2 + CO2
1880	!===========================================================
1881
1882	nb_reaction_3 = nb_reaction_3 + 1
1883
1884	indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
1885
1886	!===========================================================
1887	! a003 : O + O3 -> O2 + O2
1888	!===========================================================
1889
1890	nb_reaction_4 = nb_reaction_4 + 1
1891
1892	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
1893
1894	!===========================================================
1895	! b001 : O(1D) + CO2 -> O + CO2
1896	!===========================================================
1897
1898	nb_phot = nb_phot + 1
1899
1900	indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
1901
1902	!===========================================================
1903	! b002 : O(1D) + H2O -> OH + OH
1904	!===========================================================
1905
1906	nb_reaction_4 = nb_reaction_4 + 1
1907
1908	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
1909
1910	!===========================================================
1911	! b003 : O(1D) + H2 -> OH + H
1912	!===========================================================
1913
1914	nb_reaction_4 = nb_reaction_4 + 1
1915
1916	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
1917
1918	!===========================================================
1919	! b004 : O(1D) + O2 -> O + O2
1920	!===========================================================
1921
1922	nb_phot = nb_phot + 1
1923
1924	indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
1925
1926	!===========================================================
1927	! b005 : O(1D) + O3 -> O2 + O2
1928	!===========================================================
1929
1930	nb_reaction_4 = nb_reaction_4 + 1
1931
1932	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
1933
1934	!===========================================================
1935	! b006 : O(1D) + O3 -> O2 + O + O
1936	!===========================================================
1937
1938	nb_reaction_4 = nb_reaction_4 + 1
1939
1940	indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
1941
1942	!===========================================================
1943	! c001 : O + HO2 -> OH + O2
1944	!===========================================================
1945
1946	nb_reaction_4 = nb_reaction_4 + 1
1947
1948	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
1949
1950	!===========================================================
1951	! c002 : O + OH -> O2 + H
1952	!===========================================================
1953
1954	nb_reaction_4 = nb_reaction_4 + 1
1955
1956	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
1957
1958	!===========================================================
1959	! c003 : H + O3 -> OH + O2
1960	!===========================================================
1961
1962	nb_reaction_4 = nb_reaction_4 + 1
1963
1964	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
1965
1966	!===========================================================
1967	! c004 : H + HO2 -> OH + OH
1968	!===========================================================
1969
1970	nb_reaction_4 = nb_reaction_4 + 1
1971
1972	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
1973
1974	!===========================================================
1975	! c005 : H + HO2 -> H2 + O2
1976	!===========================================================
1977
1978	nb_reaction_4 = nb_reaction_4 + 1
1979
1980	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
1981
1982	!===========================================================
1983	! c006 : H + HO2 -> H2O + O
1984	!===========================================================
1985
1986	nb_reaction_4 = nb_reaction_4 + 1
1987
1988	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
1989
1990	!===========================================================
1991	! c007 : OH + HO2 -> H2O + O2
1992	!===========================================================
1993
1994	nb_reaction_4 = nb_reaction_4 + 1
1995
1996	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
1997
1998	!===========================================================
1999	! c008 : HO2 + HO2 -> H2O2 + O2
2000	!===========================================================
2001
2002	nb_reaction_3 = nb_reaction_3 + 1
2003
2004	indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2005
2006	!===========================================================
2007	! c009 : OH + H2O2 -> H2O + HO2
2008	!===========================================================
2009
2010	nb_reaction_4 = nb_reaction_4 + 1
2011
2012	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
2013
2014	!===========================================================
2015	! c010 : OH + H2 -> H2O + H
2016	!===========================================================
2017
2018	nb_reaction_4 = nb_reaction_4 + 1
2019
2020	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
2021
2022	!===========================================================
2023	! c011 : H + O2 + CO2 -> HO2 + CO2
2024	!===========================================================
2025
2026	nb_reaction_4 = nb_reaction_4 + 1
2027
2028	indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
2029
2030	!===========================================================
2031	! c012 : O + H2O2 -> OH + HO2
2032	!===========================================================
2033
2034	nb_reaction_4 = nb_reaction_4 + 1
2035
2036	indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
2037
2038	!===========================================================
2039	! c013 : OH + OH -> H2O + O
2040	!===========================================================
2041
2042	nb_reaction_3 = nb_reaction_3 + 1
2043
2044	indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
2045
2046	!===========================================================
2047	! c014 : OH + O3 -> HO2 + O2
2048	!===========================================================
2049
2050	nb_reaction_4 = nb_reaction_4 + 1
2051
2052	indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
2053
2054	!===========================================================
2055	! c015 : HO2 + O3 -> OH + O2 + O2
2056	!===========================================================
2057
2058	nb_reaction_4 = nb_reaction_4 + 1
2059
2060	indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
2061
2062	!===========================================================
2063	! c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
2064	!===========================================================
2065
2066	nb_reaction_3 = nb_reaction_3 + 1
2067
2068	indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2069
2070	!===========================================================
2071	! c017 : OH + OH + CO2 -> H2O2 + CO2
2072	!===========================================================
2073
2074	nb_reaction_3 = nb_reaction_3 + 1
2075
2076	indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
2077
2078	!===========================================================
2079	! c018 : H + H + CO2 -> H2 + CO2
2080	!===========================================================
2081
2082	nb_reaction_3 = nb_reaction_3 + 1
2083
2084	indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
2085
2086	!===========================================================
2087	! d001 : NO2 + O -> NO + O2
2088	!===========================================================
2089
2090	nb_reaction_4 = nb_reaction_4 + 1
2091
2092	indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
2093
2094	!===========================================================
2095	! d002 : NO + O3 -> NO2 + O2
2096	!===========================================================
2097
2098	nb_reaction_4 = nb_reaction_4 + 1
2099
2100	indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
2101
2102	!===========================================================
2103	! d003 : NO + HO2 -> NO2 + OH
2104	!===========================================================
2105
2106	nb_reaction_4 = nb_reaction_4 + 1
2107
2108	indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
2109
2110	!===========================================================
2111	! d004 : N + NO -> N2 + O
2112	!===========================================================
2113
2114	nb_reaction_4 = nb_reaction_4 + 1
2115
2116	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
2117
2118	!===========================================================
2119	! d005 : N + O2 -> NO + O
2120	!===========================================================
2121
2122	nb_reaction_4 = nb_reaction_4 + 1
2123
2124	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
2125
2126	!===========================================================
2127	! d006 : NO2 + H -> NO + OH
2128	!===========================================================
2129
2130	nb_reaction_4 = nb_reaction_4 + 1
2131
2132	indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
2133
2134	!===========================================================
2135	! d007 : N + O -> NO
2136	!===========================================================
2137
2138	nb_reaction_4 = nb_reaction_4 + 1
2139
2140	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
2141
2142	!===========================================================
2143	! d008 : N + HO2 -> NO + OH
2144	!===========================================================
2145
2146	nb_reaction_4 = nb_reaction_4 + 1
2147
2148	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
2149
2150	!===========================================================
2151	! d009 : N + OH -> NO + H
2152	!===========================================================
2153
2154	nb_reaction_4 = nb_reaction_4 + 1
2155
2156	indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
2157
2158	!===========================================================
2159	! d010 : N(2D) + O -> N + O
2160	!===========================================================
2161
2162	nb_phot = nb_phot + 1
2163
2164	indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
2165
2166	!===========================================================
2167	! d011 : N(2D) + N2 -> N + N2
2168	!===========================================================
2169
2170	nb_phot = nb_phot + 1
2171
2172	indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
2173
2174	!===========================================================
2175	! d012 : N(2D) + CO2 -> NO + CO
2176	!===========================================================
2177
2178	nb_reaction_4 = nb_reaction_4 + 1
2179
2180	indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
2181
2182	!===========================================================
2183	! e001 : CO + OH -> CO2 + H
2184	!===========================================================
2185
2186	nb_reaction_4 = nb_reaction_4 + 1
2187
2188	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
2189
2190	!===========================================================
2191	! e002 : CO + O + M -> CO2 + M
2192	!===========================================================
2193
2194	nb_reaction_4 = nb_reaction_4 + 1
2195
2196	indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
2197
2198	!Only if ion chemistry
2199	if (ionchem) then
2200
2201	!===========================================================
2202	! i001 : CO2+ + O2 -> O2+ + CO2
2203	!===========================================================
2204
2205	nb_reaction_4 = nb_reaction_4 + 1
2206
2207	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
2208
2209	!===========================================================
2210	! i002 : CO2+ + O -> O+ + CO2
2211	!===========================================================
2212
2213	nb_reaction_4 = nb_reaction_4 + 1
2214
2215	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
2216
2217	!===========================================================
2218	! i003 : CO2+ + O -> O2+ + CO
2219	!===========================================================
2220
2221	nb_reaction_4 = nb_reaction_4 + 1
2222
2223	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
2224
2225	!===========================================================
2226	! i004 : O2+ + e- -> O + O
2227	!===========================================================
2228
2229	nb_reaction_4 = nb_reaction_4 + 1
2230
2231	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
2232
2233	!===========================================================
2234	! i005 : O+ + CO2 -> O2+ + CO
2235	!===========================================================
2236
2237	nb_reaction_4 = nb_reaction_4 + 1
2238
2239	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
2240
2241	!===========================================================
2242	! i006 : CO2+ + e -> CO + O
2243	!===========================================================
2244
2245	nb_reaction_4 = nb_reaction_4 + 1
2246
2247	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
2248
2249	!===========================================================
2250	! i007 : CO2+ + NO -> NO+ + CO2
2251	!===========================================================
2252
2253	nb_reaction_4 = nb_reaction_4 + 1
2254
2255	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
2256
2257	!===========================================================
2258	! i008 : O2+ + NO -> NO+ + O2
2259	!===========================================================
2260
2261	nb_reaction_4 = nb_reaction_4 + 1
2262
2263	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
2264
2265	!===========================================================
2266	! i009 : O2+ + N2 -> NO+ + NO
2267	!===========================================================
2268
2269	nb_reaction_4 = nb_reaction_4 + 1
2270
2271	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
2272
2273	!===========================================================
2274	! i010 : O2+ + N -> NO+ + O
2275	!===========================================================
2276
2277	nb_reaction_4 = nb_reaction_4 + 1
2278
2279	indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
2280
2281	!===========================================================
2282	! i011 : O+ + N2 -> NO+ + N
2283	!===========================================================
2284
2285	nb_reaction_4 = nb_reaction_4 + 1
2286
2287	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
2288
2289	!===========================================================
2290	! i012 : NO+ + e -> N + O
2291	!===========================================================
2292
2293	nb_reaction_4 = nb_reaction_4 + 1
2294
2295	indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
2296
2297	!===========================================================
2298	! i013 : CO+ + CO2 -> CO2+ + CO
2299	!===========================================================
2300
2301	nb_reaction_4 = nb_reaction_4 + 1
2302
2303	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
2304
2305	!===========================================================
2306	! i014 : CO+ + O -> O+ + CO
2307	!===========================================================
2308
2309	nb_reaction_4 = nb_reaction_4 + 1
2310
2311	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
2312
2313	!===========================================================
2314	! i015 : C+ + CO2 -> CO+ + CO
2315	!===========================================================
2316
2317	nb_reaction_4 = nb_reaction_4 + 1
2318
2319	indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
2320
2321	!===========================================================
2322	! i016 : N2+ + CO2 -> CO2+ + N2
2323	!===========================================================
2324
2325	nb_reaction_4 = nb_reaction_4 + 1
2326
2327	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
2328
2329	!===========================================================
2330	! i017 : N2+ + O -> NO+ + N
2331	!===========================================================
2332
2333	nb_reaction_4 = nb_reaction_4 + 1
2334
2335	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
2336
2337	!===========================================================
2338	! i018 : N2+ + CO -> CO+ + N2
2339	!===========================================================
2340
2341	nb_reaction_4 = nb_reaction_4 + 1
2342
2343	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
2344
2345	!===========================================================
2346	! i019 : N2+ + e -> N + N
2347	!===========================================================
2348
2349	nb_reaction_4 = nb_reaction_4 + 1
2350
2351	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
2352
2353	!===========================================================
2354	! i020 : N2+ + O -> O+ + N2
2355	!===========================================================
2356
2357	nb_reaction_4 = nb_reaction_4 + 1
2358
2359	indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
2360
2361	!===========================================================
2362	! i021 : N+ + CO2 -> CO2+ + N
2363	!===========================================================
2364
2365	nb_reaction_4 = nb_reaction_4 + 1
2366
2367	indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
2368
2369	!===========================================================
2370	! i022 : CO+ + H -> H+ + CO
2371	!===========================================================
2372
2373	nb_reaction_4 = nb_reaction_4 + 1
2374
2375	indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
2376
2377	!===========================================================
2378	! i023 : O+ + H -> H+ + O
2379	!===========================================================
2380
2381	nb_reaction_4 = nb_reaction_4 + 1
2382
2383	indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
2384
2385	!===========================================================
2386	! i024 : H+ + O -> O+ + H
2387	!===========================================================
2388
2389	nb_reaction_4 = nb_reaction_4 + 1
2390
2391	indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
2392
2393	!===========================================================
2394	! i025 : CO2+ + H2 -> HCO2+ + H
2395	!===========================================================
2396
2397	nb_reaction_4 = nb_reaction_4 + 1
2398
2399	indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
2400
2401	!===========================================================
2402	! i026 : HCO2+ + e -> H + CO2
2403	!===========================================================
2404
2405	nb_reaction_4 = nb_reaction_4 + 1
2406
2407	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
2408
2409	!===========================================================
2410	! i027 : HCO2+ + e -> H + O + CO
2411	!===========================================================
2412	!We divide this reaction in two
2413
2414	!0.5HCO2+ + 0.5e -> H
2415
2416	nb_reaction_4 = nb_reaction_4 + 1
2417
2418	indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
2419
2420	!0.5 HCO2+ + 0.5 e -> O + CO
2421
2422	nb_reaction_4 = nb_reaction_4 + 1
2423
2424	indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
2425
2426	!===========================================================
2427	! i029 : HCO2+ + e -> OH + CO
2428	!===========================================================
2429
2430	nb_reaction_4 = nb_reaction_4 + 1
2431
2432	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
2433
2434
2435	!===========================================================
2436	! i030 : HCO2+ + e -> H + CO2
2437	!===========================================================
2438
2439	nb_reaction_4 = nb_reaction_4 + 1
2440
2441	indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
2442
2443	end if !ionchem
2444
2445	!===========================================================
2446	! h001: HO2 + ice -> products
2447	! treated as
2448	! HO2 -> 0.5 H2O + 0.75 O2
2449	!===========================================================
2450
2451	nb_phot = nb_phot + 1
2452
2453	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
2454
2455	!===========================================================
2456	! h002: OH + ice -> products
2457	! treated as
2458	! OH -> 0.5 H2O + 0.25 O2
2459	!===========================================================
2460
2461	nb_phot = nb_phot + 1
2462
2463	indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2)
2464
2465	!===========================================================
2466	! h003: H2O2 + ice -> products
2467	! treated as
2468	! H2O2 -> H2O + 0.5 O2
2469	!===========================================================
2470
2471	nb_phot = nb_phot + 1
2472
2473	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
2474
2475	!===========================================================
2476	! h004: HO2 + dust -> products
2477	! treated as
2478	! HO2 -> 0.5 H2O + 0.75 O2
2479	!===========================================================
2480
2481	nb_phot = nb_phot + 1
2482
2483	indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
2484
2485	!===========================================================
2486	! h005: H2O2 + dust -> products
2487	! treated as
2488	! H2O2 -> H2O + 0.5 O2
2489	!===========================================================
2490
2491	nb_phot = nb_phot + 1
2492
2493	indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
2494
2495	!===========================================================
2496	! check dimensions
2497	!===========================================================
2498
2499	print*, 'nb_phot = ', nb_phot
2500	print*, 'nb_reaction_4 = ', nb_reaction_4
2501	print*, 'nb_reaction_3 = ', nb_reaction_3
2502
2503	if ((nb_phot /= nb_phot_max) .or. &
2504	(nb_reaction_3 /= nb_reaction_3_max) .or. &
2505	(nb_reaction_4 /= nb_reaction_4_max)) then
2506	print*, 'wrong dimensions in indice'
2507	stop
2508	end if
2509
2510	end subroutine indice
2511
2512	!*****************************************************************
2513
2514	subroutine gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,&
2515	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2516	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2517	i_n, i_n2d, i_no, i_no2, i_n2, &
2518	i_co2plus, i_oplus, i_o2plus, i_noplus, &
2519	i_coplus, i_cplus, i_n2plus, i_nplus, &
2520	i_hplus, i_hco2plus, i_elec, dens, rm, c)
2521
2522	!*****************************************************************
2523
2524	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, &
2525	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh, &
2526	& igcm_ho2, igcm_h2o2, igcm_h2o_vap, &
2527	& igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2,&
2528	& igcm_co2plus, igcm_oplus, igcm_o2plus, &
2529	& igcm_noplus, igcm_coplus, igcm_cplus, &
2530	& igcm_n2plus, igcm_nplus, igcm_hplus, &
2531	& igcm_hco2plus, igcm_elec
2532
2533	implicit none
2534
2535	#include "callkeys.h"
2536
2537	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2538	! input:
2539	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2540
2541	integer, intent(in) :: nlayer ! number of atmospheric layers
2542	integer, intent(in) :: nq ! number of tracers in the gcm
2543	logical, intent(in) :: ionchem
2544	integer :: nesp ! number of species in the chemistry
2545	integer :: lswitch ! interface level between chemistries
2546
2547	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2548	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2549	i_n, i_n2d, i_no, i_no2, i_n2, &
2550	i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
2551	i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
2552
2553	real :: zycol(nlayer,nq) ! volume mixing ratios in the gcm
2554	real :: dens(nlayer) ! total number density (molecule.cm-3)
2555
2556	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2557	! output:
2558	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2559
2560	real, dimension(nlayer,nesp) :: rm ! volume mixing ratios
2561	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
2562
2563	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2564	! local:
2565	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2566
2567	integer :: l, iesp
2568	logical,save :: firstcall = .true.
2569
2570	! first call initializations
2571
2572	if (firstcall) then
2573
2574	! identify the indexes of the tracers we need
2575
2576	if (igcm_co2 == 0) then
2577	write(,) "gcmtochim: Error; no CO2 tracer !!!"
2578	stop
2579	endif
2580	if (igcm_co == 0) then
2581	write(,) "gcmtochim: Error; no CO tracer !!!"
2582	stop
2583	end if
2584	if (igcm_o == 0) then
2585	write(,) "gcmtochim: Error; no O tracer !!!"
2586	stop
2587	end if
2588	if (igcm_o1d == 0) then
2589	write(,) "gcmtochim: Error; no O1D tracer !!!"
2590	stop
2591	end if
2592	if (igcm_o2 == 0) then
2593	write(,) "gcmtochim: Error; no O2 tracer !!!"
2594	stop
2595	end if
2596	if (igcm_o3 == 0) then
2597	write(,) "gcmtochim: Error; no O3 tracer !!!"
2598	stop
2599	end if
2600	if (igcm_h == 0) then
2601	write(,) "gcmtochim: Error; no H tracer !!!"
2602	stop
2603	end if
2604	if (igcm_h2 == 0) then
2605	write(,) "gcmtochim: Error; no H2 tracer !!!"
2606	stop
2607	end if
2608	if (igcm_oh == 0) then
2609	write(,) "gcmtochim: Error; no OH tracer !!!"
2610	stop
2611	end if
2612	if (igcm_ho2 == 0) then
2613	write(,) "gcmtochim: Error; no HO2 tracer !!!"
2614	stop
2615	end if
2616	if (igcm_h2o2 == 0) then
2617	write(,) "gcmtochim: Error; no H2O2 tracer !!!"
2618	stop
2619	end if
2620	if (igcm_n == 0) then
2621	write(,) "gcmtochim: Error; no N tracer !!!"
2622	stop
2623	end if
2624	if (igcm_n2d == 0) then
2625	write(,) "gcmtochim: Error; no N2D tracer !!!"
2626	stop
2627	end if
2628	if (igcm_no == 0) then
2629	write(,) "gcmtochim: Error; no NO tracer !!!"
2630	stop
2631	end if
2632	if (igcm_no2 == 0) then
2633	write(,) "gcmtochim: Error; no NO2 tracer !!!"
2634	stop
2635	end if
2636	if (igcm_n2 == 0) then
2637	write(,) "gcmtochim: Error; no N2 tracer !!!"
2638	stop
2639	end if
2640	if (igcm_h2o_vap == 0) then
2641	write(,) "gcmtochim: Error; no water vapor tracer !!!"
2642	stop
2643	end if
2644	if (ionchem) then
2645	if (igcm_co2plus == 0) then
2646	write(,) "gcmtochim: Error; no CO2+ tracer !!!"
2647	stop
2648	end if
2649	if (igcm_oplus == 0) then
2650	write(,) "gcmtochim: Error; no O+ tracer !!!"
2651	stop
2652	end if
2653	if (igcm_o2plus == 0) then
2654	write(,) "gcmtochim: Error; no O2+ tracer !!!"
2655	stop
2656	end if
2657	if (igcm_noplus == 0) then
2658	write(,) "gcmtochim: Error; no NO+ tracer !!!"
2659	stop
2660	endif
2661	if (igcm_coplus == 0) then
2662	write(,) "gcmtochim: Error; no CO+ tracer !!!"
2663	stop
2664	endif
2665	if (igcm_cplus == 0) then
2666	write(,) "gcmtochim: Error; no C+ tracer !!!"
2667	stop
2668	endif
2669	if (igcm_n2plus == 0) then
2670	write(,) "gcmtochim: Error; no N2+ tracer !!!"
2671	stop
2672	endif
2673	if (igcm_nplus == 0) then
2674	write(,) "gcmtochim: Error; no N+ tracer !!!"
2675	stop
2676	endif
2677	if (igcm_hplus == 0) then
2678	write(,) "gcmtochim: Error; no H+ tracer !!!"
2679	stop
2680	endif
2681	if (igcm_hco2plus == 0) then
2682	write(,) "gcmtochim: Error; no HCO2+ tracer !!!"
2683	stop
2684	endif
2685	if (igcm_elec == 0) then
2686	write(,) "gcmtochim: Error; no e- tracer !!!"
2687	stop
2688	end if
2689	end if ! ionchem
2690	firstcall = .false.
2691	end if ! of if (firstcall)
2692
2693	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2694	! initialise mixing ratios
2695	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2696
2697	do l = 1,nlayer
2698	rm(l,i_co2) = zycol(l, igcm_co2)
2699	rm(l,i_co) = zycol(l, igcm_co)
2700	rm(l,i_o) = zycol(l, igcm_o)
2701	rm(l,i_o1d) = zycol(l, igcm_o1d)
2702	rm(l,i_o2) = zycol(l, igcm_o2)
2703	rm(l,i_o3) = zycol(l, igcm_o3)
2704	rm(l,i_h) = zycol(l, igcm_h)
2705	rm(l,i_h2) = zycol(l, igcm_h2)
2706	rm(l,i_oh) = zycol(l, igcm_oh)
2707	rm(l,i_ho2) = zycol(l, igcm_ho2)
2708	rm(l,i_h2o2) = zycol(l, igcm_h2o2)
2709	rm(l,i_h2o) = zycol(l, igcm_h2o_vap)
2710	rm(l,i_n) = zycol(l, igcm_n)
2711	rm(l,i_n2d) = zycol(l, igcm_n2d)
2712	rm(l,i_no) = zycol(l, igcm_no)
2713	rm(l,i_no2) = zycol(l, igcm_no2)
2714	rm(l,i_n2) = zycol(l, igcm_n2)
2715	end do
2716
2717	if (ionchem) then
2718	do l = 1,nlayer
2719	rm(l,i_co2plus) = zycol(l, igcm_co2plus)
2720	rm(l,i_oplus) = zycol(l, igcm_oplus)
2721	rm(l,i_o2plus) = zycol(l, igcm_o2plus)
2722	rm(l,i_noplus) = zycol(l, igcm_noplus)
2723	rm(l,i_coplus) = zycol(l, igcm_coplus)
2724	rm(l,i_cplus) = zycol(l, igcm_cplus)
2725	rm(l,i_n2plus) = zycol(l, igcm_n2plus)
2726	rm(l,i_nplus) = zycol(l, igcm_nplus)
2727	rm(l,i_hplus) = zycol(l, igcm_hplus)
2728	rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
2729	rm(l,i_elec) = zycol(l, igcm_elec)
2730	end do
2731	end if
2732
2733	where (rm(:,:) < 1.e-30)
2734	rm(:,:) = 0.
2735	end where
2736
2737	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2738	! initialise number densities
2739	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2740
2741	do iesp = 1,nesp
2742	do l = 1,nlayer
2743	c(l,iesp) = rm(l,iesp)*dens(l)
2744	end do
2745	end do
2746
2747	end subroutine gcmtochim
2748
2749	!*****************************************************************
2750
2751	subroutine chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp, &
2752	i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2753	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2754	i_n, i_n2d, i_no, i_no2, i_n2, &
2755	i_co2plus, i_oplus, i_o2plus, i_noplus, &
2756	i_coplus, i_cplus, i_n2plus, i_nplus, &
2757	i_hplus, i_hco2plus, i_elec, dens, c)
2758
2759	!*****************************************************************
2760
2761	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, &
2762	igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh, &
2763	igcm_ho2, igcm_h2o2, igcm_h2o_vap, &
2764	igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2, &
2765	igcm_co2plus, igcm_oplus, igcm_o2plus, &
2766	igcm_noplus, igcm_coplus, igcm_cplus, &
2767	igcm_n2plus, igcm_nplus, igcm_hplus, &
2768	igcm_hco2plus, igcm_elec
2769
2770	implicit none
2771
2772	#include "callkeys.h"
2773
2774	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2775	! inputs:
2776	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2777
2778	integer, intent(in) :: nlayer ! number of atmospheric layers
2779	integer, intent(in) :: nq ! number of tracers in the gcm
2780	logical, intent(in) :: ionchem
2781	integer :: nesp ! number of species in the chemistry
2782	integer :: lswitch ! interface level between chemistries
2783	integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
2784	i_h2, i_oh, i_ho2, i_h2o2, i_h2o, &
2785	i_n, i_n2d, i_no, i_no2, i_n2, &
2786	i_co2plus, i_oplus, i_o2plus, i_noplus, &
2787	i_coplus, i_cplus, i_n2plus, i_nplus, &
2788	i_hplus, i_hco2plus, i_elec
2789
2790	real :: dens(nlayer) ! total number density (molecule.cm-3)
2791	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
2792
2793	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2794	! output:
2795	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2796
2797	real zycol(nlayer,nq) ! volume mixing ratios in the gcm
2798
2799	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2800	! local:
2801	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2802
2803	integer l
2804
2805	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2806	! save mixing ratios for the gcm
2807	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2808
2809	do l = 1,lswitch-1
2810	zycol(l, igcm_co2) = c(l,i_co2)/dens(l)
2811	zycol(l, igcm_co) = c(l,i_co)/dens(l)
2812	zycol(l, igcm_o) = c(l,i_o)/dens(l)
2813	zycol(l, igcm_o1d) = c(l,i_o1d)/dens(l)
2814	zycol(l, igcm_o2) = c(l,i_o2)/dens(l)
2815	zycol(l, igcm_o3) = c(l,i_o3)/dens(l)
2816	zycol(l, igcm_h) = c(l,i_h)/dens(l)
2817	zycol(l, igcm_h2) = c(l,i_h2)/dens(l)
2818	zycol(l, igcm_oh) = c(l,i_oh)/dens(l)
2819	zycol(l, igcm_ho2) = c(l,i_ho2)/dens(l)
2820	zycol(l, igcm_h2o2) = c(l,i_h2o2)/dens(l)
2821	zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
2822	zycol(l, igcm_n) = c(l,i_n)/dens(l)
2823	zycol(l, igcm_n2d) = c(l,i_n2d)/dens(l)
2824	zycol(l, igcm_no) = c(l,i_no)/dens(l)
2825	zycol(l, igcm_no2) = c(l,i_no2)/dens(l)
2826	zycol(l, igcm_n2) = c(l,i_n2)/dens(l)
2827	end do
2828
2829	if (ionchem) then
2830	do l = 1,lswitch-1
2831	zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
2832	zycol(l, igcm_oplus) = c(l,i_oplus)/dens(l)
2833	zycol(l, igcm_o2plus) = c(l,i_o2plus)/dens(l)
2834	zycol(l, igcm_noplus) = c(l,i_noplus)/dens(l)
2835	zycol(l, igcm_coplus) = c(l,i_coplus)/dens(l)
2836	zycol(l, igcm_cplus) = c(l,i_cplus)/dens(l)
2837	zycol(l, igcm_n2plus) = c(l,i_n2plus)/dens(l)
2838	zycol(l, igcm_nplus) = c(l,i_nplus)/dens(l)
2839	zycol(l, igcm_hplus) = c(l,i_hplus)/dens(l)
2840	zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
2841	zycol(l, igcm_elec) = c(l,i_elec)/dens(l)
2842	end do
2843	end if
2844
2845	end subroutine chimtogcm
2846
2847	end subroutine photochemistry

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

Download in other formats:

Original Format