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

Last change on this file since 3026 was 2929, checked in by abierjon, 21 months ago

Mars GCM:
Replicate Venus commits r2927 and r2928 for Mars photochemistry :
Correction of the reaction N2+ + e- > N + N

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