source: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90 @ 2505

Last change on this file since 2505 was 2490, checked in by flefevre, 4 years ago

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