Context Navigation

photochemistry.F90 @ 2160

Last change on this file since 2160 was 2158, checked in by emillour, 6 years ago

Mars GCM:

Updated chemical core to include ionospheric chemistry

FGG

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

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

Download in other formats:

Original Format