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

Last change on this file since 3162 was 3141, checked in by flefevre, 19 months ago

Correction of a duplicated ionic reaction (i026 and i030) + a bit of cosmetics
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: spare slot
1872
1873         i026(:) = 0.
1874
1875         nb_reaction_4 = nb_reaction_4 + 1
1876         v_4(:,nb_reaction_4) = i026(:)
1877
1878!---     i027+i028: HCO2+ + e -> H + O + CO
1879
1880!        UMIST
1881         !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
1882         !i028: 0.5 (HCO2+ + e-) -> O + CO
1883
1884         i027(:) = 8.1e-7*((300./t_elect(:))**0.64)
1885
1886         nb_reaction_4 = nb_reaction_4 + 1
1887         v_4(:,nb_reaction_4) = i027(:)
1888
1889         nb_reaction_4 = nb_reaction_4 + 1
1890         v_4(:,nb_reaction_4) = i027(:)
1891
1892!---     i029: HCO2+ + e -> OH + CO
1893
1894!        UMIST
1895
1896         i029(:) = 3.2e-7*((300./t_elect(:))**0.64)
1897
1898         nb_reaction_4 = nb_reaction_4 + 1
1899         v_4(:,nb_reaction_4) = i029(:)
1900
1901!---     i030: HCO2+ + e -> H + CO2
1902
1903!        UMIST
1904
1905         i030(:) = 6.0e-8*((300./t_elect(:))**0.64)
1906
1907         nb_reaction_4 = nb_reaction_4 + 1
1908         v_4(:,nb_reaction_4) = i030(:)
1909
1910!---     i031: HCO2+ + O -> HCO+ + O2
1911
1912!        UMIST
1913
1914         i031(:) = 1.e-9
1915
1916         nb_reaction_4 = nb_reaction_4 + 1
1917         v_4(:,nb_reaction_4) = i031(:)
1918
1919!---     i032: HCO2+ + CO -> HCO+ + CO2
1920
1921!        UMIST, from Prassad & Huntress 1980
1922
1923         i032(:) = 7.8e-10
1924
1925         nb_reaction_4 = nb_reaction_4 + 1
1926         v_4(:,nb_reaction_4) = i032(:)
1927
1928!---     i033: H+ + CO2 -> HCO+ + O
1929
1930!        UMIST, from Smith et al., Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
1931
1932         i033(:) = 3.5e-9
1933
1934         nb_reaction_4 = nb_reaction_4 + 1
1935         v_4(:,nb_reaction_4) = i033(:)
1936
1937!---     i034: CO2+ + H -> HCO+ + O
1938
1939!        Seen in Fox 2015, from Borodi et al., Int. J. Mass Spectrom. 280, 218-225, 2009
1940
1941         i034(:) = 4.5e-10
1942         nb_reaction_4 = nb_reaction_4 + 1
1943         v_4(:,nb_reaction_4) = i034(:)
1944
1945!---     i035: CO+ + H2 -> HCO+ + H
1946
1947         !UMIST, from Scott et al., J. Chem. Phys., 106, 3982-3987(1997)
1948
1949         i035(:) = 7.5e-10
1950
1951         nb_reaction_4 = nb_reaction_4 + 1
1952         v_4(:,nb_reaction_4) = i035(:)
1953
1954!---     i036: HCO+ + e- -> CO + H
1955
1956         !UMIST, from Mitchell, Phys. Rep., 186, 215 (1990)
1957
1958         i036(:) = 2.4e-7 *((300./t_elect(:))**0.69)
1959
1960         nb_reaction_4 = nb_reaction_4 + 1
1961         v_4(:,nb_reaction_4) = i036(:)
1962
1963!---     i037: CO2+ + H2O -> H2O+ + CO2
1964
1965         !UMIST, from Karpas, Z., Anicich, V.G., and Huntress, W.T., Chem. Phys. Lett., 59, 84 (1978)
1966
1967         i037(:) = 2.04e-9 *((300./t_elect(:))**0.5)
1968
1969         nb_reaction_4 = nb_reaction_4 + 1
1970         v_4(:,nb_reaction_4) = i037(:)
1971
1972!---     i038: CO+ + H2O -> H2O+ + CO
1973
1974         !UMIST, from Huntress, W.T., McEwan, M.J., Karpas, Z., and Anicich, V.G., Astrophys. J. Supp. Series, 44, 481 (1980)
1975
1976         i038(:) = 1.72e-9*((300./t_elect(:))**0.5)
1977
1978         nb_reaction_4 = nb_reaction_4 + 1
1979         v_4(:,nb_reaction_4) = i038(:)
1980
1981!---     i039: O+ + H2O -> H2O+ + O
1982
1983         !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)
1984
1985         i039(:) = 3.2e-9*((300./t_elect(:))**0.5)
1986
1987         nb_reaction_4 = nb_reaction_4 + 1
1988         v_4(:,nb_reaction_4) = i039(:)
1989
1990!---     i040: N2+ + H2O -> H2O+ + N2
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         i040(:) = 2.3e-9*((300./t_elect(:))**0.5)
1995
1996         nb_reaction_4 = nb_reaction_4 + 1
1997         v_4(:,nb_reaction_4) = i040(:)
1998
1999!---     i041: N+ + H2O -> H2O+ + N
2000
2001         !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)
2002
2003         i041(:) = 2.8e-9*((300./t_elect(:))**0.5)
2004
2005         nb_reaction_4 = nb_reaction_4 + 1
2006         v_4(:,nb_reaction_4) = i041(:)
2007
2008!---     i042: H+ + H2O -> H2O+ + H
2009
2010         !UMIST, from D. Smith, P. Spanel and C. A. Mayhew, Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
2011
2012         i042(:) = 6.9e-9*((300./t_elect(:))**0.5)
2013
2014         nb_reaction_4 = nb_reaction_4 + 1
2015         v_4(:,nb_reaction_4) = i042(:)
2016
2017!---     i043: H2O+ + O2 -> O2+ + H2O
2018
2019         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
2020
2021         i043(:) = 4.6e-10
2022
2023         nb_reaction_4 = nb_reaction_4 + 1
2024         v_4(:,nb_reaction_4) = i043(:)
2025
2026!---     i044: H2O+ + CO -> HCO+ + OH
2027
2028         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2029
2030         i044(:) = 5.0e-10
2031
2032         nb_reaction_4 = nb_reaction_4 + 1
2033         v_4(:,nb_reaction_4) = i044(:)
2034
2035!---     i045: H2O+ + O -> O2+ + H2
2036
2037         !UMIST, from Viggiano, A.A, Howarka, F., Albritton, D.L., Fehsenfeld, F.C., Adams, N.G., and Smith, D., Astrophys. J., 236, 492 (1980)
2038
2039         i045(:) = 4.0e-11
2040
2041         nb_reaction_4 = nb_reaction_4 + 1
2042         v_4(:,nb_reaction_4) = i045(:)
2043
2044!---     i046: H2O+ + NO -> NO+ + H2O
2045
2046         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
2047
2048         i046(:) = 2.7e-10
2049
2050         nb_reaction_4 = nb_reaction_4 + 1
2051         v_4(:,nb_reaction_4) = i046(:)
2052
2053!---     i047: H2O+ + e- -> H + H + O
2054
2055         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
2056         
2057         i047(:) = 3.05e-7*((300./t_elect(:))**0.5)
2058
2059         nb_reaction_4 = nb_reaction_4 + 1
2060         v_4(:,nb_reaction_4) = i047(:)
2061
2062!---     i048: H2O+ + e- -> H + OH
2063         
2064         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
2065
2066         i048(:) = 8.6e-8*((300./t_elect(:))**0.5)
2067
2068         nb_reaction_4 = nb_reaction_4 + 1
2069         v_4(:,nb_reaction_4) = i048(:)
2070
2071!---     i049: H2O+ + e- -> O + H2
2072
2073         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
2074
2075         i049(:) = 3.9e-8*((300./t_elect(:))**0.5)
2076
2077         nb_reaction_4 = nb_reaction_4 + 1
2078         v_4(:,nb_reaction_4) = i049(:)
2079
2080!---     i050: H2O+ + H2O -> H3O+ + OH
2081
2082         !UMIST, from Huntress, W.T. and Pinizzotto, R.F., J. Chem. Phys., 59, 4742 (1973)
2083
2084         i050(:) = 2.1e-9*((300./t_elect(:))**0.5)
2085
2086         nb_reaction_4 = nb_reaction_4 + 1
2087         v_4(:,nb_reaction_4) = i050(:)
2088
2089!---     i051: H2O+ + H2 -> H3O+ + H
2090
2091         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
2092
2093         i051(:) = 6.4e-10
2094
2095         nb_reaction_4 = nb_reaction_4 + 1
2096         v_4(:,nb_reaction_4) = i051(:)
2097
2098!---     i052: HCO+ + H2O -> H3O+ + CO
2099
2100         !UMIST, from Adams, N.G., Smith, D., and Grief, D., Int. J. Mass Spectrom. Ion Phys., 26, 405 (1978)
2101
2102         i052(:) = 2.5e-9*((300./t_elect(:))**0.5)
2103
2104         nb_reaction_4 = nb_reaction_4 + 1
2105         v_4(:,nb_reaction_4) = i052(:)
2106
2107!---     i053: H3O+ + e -> OH + H + H
2108
2109         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2110
2111         i053(:) = 3.05e-7*((300./t_elect(:))**0.5)
2112
2113         nb_reaction_4 = nb_reaction_4 + 1
2114         v_4(:,nb_reaction_4) = i053(:)
2115
2116!---     i054: H3O + + e -> H2O + H
2117
2118         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2119         
2120         i054(:) = 7.09e-8*((300./t_elect(:))**0.5)
2121
2122         nb_reaction_4 = nb_reaction_4 + 1
2123         v_4(:,nb_reaction_4) = i054(:)
2124
2125!---     i055: H3O+ + e -> OH + H2
2126
2127         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2128
2129         i055(:) = 5.37e-8*((300./t_elect(:))**0.5)
2130
2131         nb_reaction_4 = nb_reaction_4 + 1
2132         v_4(:,nb_reaction_4) = i055(:)
2133
2134!---     i056: H3O+ + e -> O + H2 + H
2135
2136         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
2137
2138         i056(:) = 5.6e-9*((300./t_elect(:))**0.5)
2139
2140         nb_reaction_4 = nb_reaction_4 + 1
2141         v_4(:,nb_reaction_4) = i056(:)
2142
2143         nb_reaction_4 = nb_reaction_4 + 1
2144         v_4(:,nb_reaction_4) = i056(:)
2145
2146!---     i057: O+ + H2 -> OH+ + H
2147
2148         !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)
2149
2150         i057(:) = 1.7e-9
2151
2152         nb_reaction_4 = nb_reaction_4 + 1
2153         v_4(:,nb_reaction_4) = i057(:)
2154
2155!---     i058: OH+ + O -> O2+ + H
2156
2157         !UMIST, from Prasad & Huntress, 1980, ApJS, 43, 1
2158
2159         i058(:) = 7.1e-10
2160
2161         nb_reaction_4 = nb_reaction_4 + 1
2162         v_4(:,nb_reaction_4) = i058(:)
2163
2164!---     i059: OH+ + CO2 -> HCO2+ + O
2165
2166         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2167
2168         i059(:) = 1.44e-9
2169
2170         nb_reaction_4 = nb_reaction_4 + 1
2171         v_4(:,nb_reaction_4) = i059(:)
2172
2173!---     i060: OH+ + CO -> HCO+ + O
2174
2175         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2176
2177         i060(:) = 1.05e-9
2178
2179         nb_reaction_4 = nb_reaction_4 + 1
2180         v_4(:,nb_reaction_4) = i060(:)
2181
2182!---     i061: OH+ + NO -> NO+ + OH (tasa de reacción UMIST 3.59e-10)
2183
2184         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2185
2186         i061(:) = 3.59e-10
2187
2188         nb_reaction_4 = nb_reaction_4 + 1
2189         v_4(:,nb_reaction_4) = i061(:)
2190
2191!---     i062: OH+ + H2 -> H2O+ + H (tasa de reacción UMIST 1.01e-9,
2192
2193         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2194
2195         i062(:) = 1.01e-9
2196
2197         nb_reaction_4 = nb_reaction_4 + 1
2198         v_4(:,nb_reaction_4) = i062(:)
2199
2200!---     i063: OH+ + O2 -> O2+ + OH (tasa de reacción UMIST 5.9e-10
2201
2202         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
2203
2204         i063(:) = 5.9e-10
2205
2206         nb_reaction_4 = nb_reaction_4 + 1
2207         v_4(:,nb_reaction_4) = i063(:)
2208
2209      end if   !ionchem
2210
2211!----------------------------------------------------------------------
2212!     heterogeneous chemistry
2213!----------------------------------------------------------------------
2214
2215      if (hetero_ice) then
2216
2217!        k = (surface*v*gamma)/4 (s-1)
2218!        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
2219 
2220!---     h001: ho2 + ice -> products
2221 
2222!        cooper and abbatt, 1996: gamma = 0.025
2223     
2224         gam = 0.025
2225         h001(:) = surfice1d(:)       &
2226                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
2227 
2228!        h002: oh + ice -> products
2229 
2230!        cooper and abbatt, 1996: gamma = 0.03
2231 
2232         gam = 0.03
2233         h002(:) = surfice1d(:)       &
2234                   *100.*sqrt(8.*8.31*t(:)/(17.e-3*pi))*gam/4.
2235
2236!---     h003: h2o2 + ice -> products
2237 
2238!        gamma = 0.    test value
2239 
2240         gam = 0.
2241         h003(:) = surfice1d(:)       &
2242                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
2243      else
2244         h001(:) = 0.
2245         h002(:) = 0.
2246         h003(:) = 0.
2247      end if
2248
2249      nb_phot = nb_phot + 1
2250      v_phot(:,nb_phot) = h001(:)
2251
2252      nb_phot = nb_phot + 1
2253      v_phot(:,nb_phot) = h002(:)
2254
2255      nb_phot = nb_phot + 1
2256      v_phot(:,nb_phot) = h003(:)
2257
2258      if (hetero_dust) then
2259 
2260!---     h004: ho2 + dust -> products
2261 
2262!        jacob, 2000: gamma = 0.2
2263!        see dereus et al., atm. chem. phys., 2005
2264 
2265         gam = 0.2
2266         h004(:) = surfdust1d(:)  &
2267                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
2268 
2269!---     h005: h2o2 + dust -> products
2270 
2271!        gamma = 5.e-4
2272!        see dereus et al., atm. chem. phys., 2005
2273 
2274         gam = 5.e-4
2275         h005(:) = surfdust1d(:)  &
2276                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
2277      else
2278         h004(:) = 0.
2279         h005(:) = 0.
2280      end if
2281 
2282      nb_phot = nb_phot + 1
2283      v_phot(:,nb_phot) = h004(:)
2284
2285      nb_phot = nb_phot + 1
2286      v_phot(:,nb_phot) = h005(:)
2287
2288end subroutine reactionrates
2289
2290!======================================================================
2291
2292 subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer,            &
2293                        nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
2294                        v_phot, v_3, v_4)
2295
2296!======================================================================
2297! filling of the jacobian matrix
2298!======================================================================
2299
2300use types_asis
2301
2302implicit none
2303
2304! input
2305
2306integer             :: ilev    ! level index
2307integer             :: nesp    ! number of species in the chemistry
2308integer, intent(in) :: nlayer  ! number of atmospheric layers
2309integer, intent(in) :: nb_reaction_3_max
2310                               ! number of quadratic reactions
2311integer, intent(in) :: nb_reaction_4_max
2312                               ! number of bimolecular reactions
2313integer, intent(in) :: nb_phot_max
2314                               ! number of processes treated numerically as photodissociations
2315
2316real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
2317real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
2318real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
2319real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
2320
2321! output
2322
2323real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
2324real (kind = 8), dimension(nesp), intent(out)      :: prod, loss
2325
2326! local
2327
2328integer :: iesp
2329integer :: ind_phot_2,ind_phot_4,ind_phot_6
2330integer :: ind_3_2,ind_3_4,ind_3_6
2331integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
2332integer :: iphot,i3,i4
2333
2334real(kind = 8) :: eps, eps_4  ! implicit/explicit coefficient
2335
2336! initialisations
2337
2338mat(:,:) = 0.
2339prod(:)  = 0.
2340loss(:)  = 0.
2341
2342! photodissociations
2343! or reactions a + c -> b + c
2344! or reactions a + ice -> b + c
2345
2346do iphot = 1,nb_phot_max
2347
2348  ind_phot_2 = indice_phot(iphot)%z2
2349  ind_phot_4 = indice_phot(iphot)%z4
2350  ind_phot_6 = indice_phot(iphot)%z6
2351
2352  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
2353  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
2354  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
2355
2356  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
2357  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
2358  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
2359
2360end do
2361
2362! reactions a + a -> b + c
2363
2364do i3 = 1,nb_reaction_3_max
2365
2366  ind_3_2 = indice_3(i3)%z2
2367  ind_3_4 = indice_3(i3)%z4
2368  ind_3_6 = indice_3(i3)%z6
2369
2370  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)
2371  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)
2372  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)
2373
2374  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
2375  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)
2376  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)
2377
2378end do
2379
2380! reactions a + b -> c + d
2381
2382eps = 1.d-10
2383
2384do i4 = 1,nb_reaction_4_max
2385
2386  ind_4_2 = indice_4(i4)%z2
2387  ind_4_4 = indice_4(i4)%z4
2388  ind_4_6 = indice_4(i4)%z6
2389  ind_4_8 = indice_4(i4)%z8
2390
2391  eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
2392  eps_4 = min(eps_4,1.0)
2393
2394  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)
2395  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)
2396  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)
2397  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)   
2398  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)
2399  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)
2400  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)
2401  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)
2402
2403
2404  loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
2405  loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
2406  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)
2407  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)
2408
2409end do
2410
2411end subroutine fill_matrix
2412
2413!================================================================
2414
2415 subroutine indice(nb_reaction_3_max, nb_reaction_4_max,                &
2416                   nb_phot_max, ionchem, deutchem, i_co2, i_co, i_o,    &
2417                   i_o1d, i_o2, i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2,    &
2418                   i_h2o, i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,     &
2419                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
2420                   i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
2421                   i_h2oplus, i_h3oplus, i_ohplus, i_elec,              &
2422                   i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2)
2423
2424!================================================================
2425! set the "indice" arrays used to fill the jacobian matrix      !
2426!----------------------------------------------------------------
2427! reaction                                   type               !
2428!----------------------------------------------------------------
2429! A + hv   --> B + C     photolysis          indice_phot        !
2430! A + B    --> C + D     bimolecular         indice_4           !
2431! A + A    --> B + C     quadratic           indice_3           !
2432! A + C    --> B + C     quenching           indice_phot        !
2433! A + ice  --> B + C     heterogeneous       indice_phot        !
2434!================================================================
2435
2436use types_asis
2437
2438implicit none
2439
2440! input
2441
2442integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
2443           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
2444           i_n, i_n2d, i_no, i_no2, i_n2,                   &
2445           i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
2446           i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, &
2447           i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec, &
2448           i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2
2449integer, intent(in) :: nb_reaction_3_max
2450                       ! number of quadratic reactions
2451integer, intent(in) :: nb_reaction_4_max
2452                       ! number of bimolecular reactions
2453integer, intent(in) :: nb_phot_max
2454                       ! number of processes treated numerically as photodissociations
2455logical, intent(in) :: ionchem
2456logical, intent(in) :: deutchem
2457
2458! local
2459
2460integer :: nb_phot, nb_reaction_3, nb_reaction_4
2461integer :: i_dummy
2462
2463allocate (indice_phot(nb_phot_max))
2464allocate (indice_3(nb_reaction_3_max))
2465allocate (indice_4(nb_reaction_4_max))
2466
2467i_dummy = 1
2468
2469nb_phot       = 0
2470nb_reaction_3 = 0
2471nb_reaction_4 = 0
2472
2473!===========================================================
2474!      O2 + hv -> O + O
2475!===========================================================
2476
2477nb_phot = nb_phot + 1
2478
2479indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy)
2480
2481!===========================================================
2482!      O2 + hv -> O + O(1D)
2483!===========================================================
2484
2485nb_phot = nb_phot + 1
2486
2487indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d)
2488
2489!===========================================================
2490!      CO2 + hv -> CO + O
2491!===========================================================
2492
2493nb_phot = nb_phot + 1
2494
2495indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o)
2496
2497!===========================================================
2498!      CO2 + hv -> CO + O(1D)
2499!===========================================================
2500
2501nb_phot = nb_phot + 1
2502
2503indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d)
2504
2505!===========================================================
2506!      O3 + hv -> O2 + O(1D)
2507!===========================================================
2508
2509nb_phot = nb_phot + 1
2510
2511indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d)
2512
2513!===========================================================
2514!      O3 + hv -> O2 + O
2515!===========================================================
2516
2517nb_phot = nb_phot + 1
2518
2519indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o)
2520
2521!===========================================================
2522!      H2O + hv -> H + OH
2523!===========================================================
2524
2525nb_phot = nb_phot + 1
2526
2527indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh)
2528
2529!===========================================================
2530!      H2O2 + hv -> OH + OH
2531!===========================================================
2532
2533nb_phot = nb_phot + 1
2534
2535indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy)
2536
2537!===========================================================
2538!      HO2 + hv -> OH + O
2539!===========================================================
2540
2541nb_phot = nb_phot + 1
2542
2543indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o)
2544
2545!===========================================================
2546!      H2 + hv -> H + H
2547!===========================================================
2548
2549nb_phot = nb_phot + 1
2550
2551indice_phot(nb_phot) = z3spec(1.0, i_h2, 1.0, i_h, 1.0, i_h)
2552
2553!===========================================================
2554!      NO + hv -> N + O
2555!===========================================================
2556
2557nb_phot = nb_phot + 1
2558
2559indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o)
2560
2561!===========================================================
2562!      NO2 + hv -> NO + O
2563!===========================================================
2564
2565nb_phot = nb_phot + 1
2566
2567indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o)
2568
2569!===========================================================
2570!      N2 + hv -> N + N
2571!===========================================================
2572
2573nb_phot = nb_phot + 1
2574
2575indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n)
2576
2577!Only if deuterium chemistry included
2578if (deutchem) then
2579!===========================================================
2580!      HDO + hv -> OD + H
2581!===========================================================
2582
2583   nb_phot = nb_phot + 1
2584
2585   indice_phot(nb_phot) = z3spec(1.0, i_hdo, 1.0, i_h, 1.0, i_od)
2586
2587!===========================================================
2588!      HDO + hv -> D + OH
2589!===========================================================
2590
2591   nb_phot = nb_phot + 1
2592
2593   indice_phot(nb_phot) = z3spec(1.0, i_hdo, 1.0, i_d, 1.0, i_oh)
2594   
2595endif !deutchem
2596
2597!Only if ion chemistry included
2598if (ionchem) then
2599
2600!===========================================================
2601!      CO2 + hv -> CO2+ + e-
2602!===========================================================
2603
2604   nb_phot = nb_phot + 1
2605
2606   indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
2607
2608!===========================================================
2609!      CO2 + hv -> O+ + CO + e-
2610!===========================================================
2611!We divide this reaction in two
2612
2613!0.5 CO2 + hv -> CO
2614   nb_phot = nb_phot + 1
2615
2616   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
2617
2618!0.5 CO2 + hv -> O+ + e-
2619   nb_phot = nb_phot + 1
2620
2621   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
2622
2623!===========================================================
2624!      CO2 + hv -> CO+ + O + e-
2625!===========================================================
2626!We divide this reaction in two
2627
2628!0.5 CO2 + hv -> O
2629   nb_phot = nb_phot + 1
2630
2631   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
2632
2633!0.5 CO2 + hv -> CO+ + e-
2634   nb_phot = nb_phot + 1
2635
2636   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
2637
2638!===========================================================
2639!      CO2 + hv -> C+ + O2 + e-
2640!===========================================================
2641!We divide this reaction in two
2642
2643!0.5 CO2 + hv -> O2
2644   nb_phot = nb_phot + 1
2645
2646   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
2647
2648!0.5 CO2 + hv -> C+ + e-
2649   nb_phot = nb_phot + 1
2650
2651   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
2652
2653!===========================================================
2654!      O2 + hv -> O2+ + e-
2655!===========================================================
2656
2657   nb_phot = nb_phot + 1
2658
2659   indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
2660
2661!===========================================================
2662!      O + hv -> O+ + e-
2663!===========================================================
2664
2665   nb_phot = nb_phot + 1
2666
2667   indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
2668
2669!===========================================================
2670!      NO + hv -> NO+ + e-
2671!===========================================================
2672
2673   nb_phot = nb_phot + 1
2674
2675   indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
2676
2677!===========================================================
2678!      CO + hv -> CO+ + e-
2679!===========================================================
2680
2681   nb_phot = nb_phot + 1
2682
2683   indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
2684
2685!===========================================================
2686!      CO + hv -> C+ + O + e-
2687!===========================================================
2688!We divide this reaction in two
2689
2690!0.5 CO + hv -> O
2691   nb_phot = nb_phot + 1
2692
2693   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
2694
2695!0.5 CO + hv -> C+ + e-
2696   nb_phot = nb_phot + 1
2697
2698   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
2699
2700!===========================================================
2701!      N2 + hv -> N2+ + e-
2702!===========================================================
2703
2704   nb_phot = nb_phot + 1
2705
2706   indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
2707
2708!===========================================================
2709!      N2 + hv -> N+ + N + e-
2710!===========================================================
2711!We divide this reaction in two
2712
2713!0.5 N2 + hv -> N
2714   nb_phot = nb_phot + 1
2715
2716   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
2717
2718!0.5 N2 + hv -> N+ + e-
2719   nb_phot = nb_phot + 1
2720
2721   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
2722
2723!===========================================================
2724!      N + hv -> N+ + e-
2725!===========================================================
2726
2727   nb_phot = nb_phot + 1
2728
2729   indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
2730
2731!===========================================================
2732!      H + hv -> H+ + e-
2733!===========================================================
2734
2735   nb_phot = nb_phot + 1
2736
2737   indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
2738
2739end if   !ionchem
2740
2741!===========================================================
2742!      a001 : O + O2 + CO2 -> O3 + CO2
2743!===========================================================
2744
2745nb_reaction_4 = nb_reaction_4 + 1
2746
2747indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
2748
2749!===========================================================
2750!      a002 : O + O + CO2 -> O2 + CO2
2751!===========================================================
2752
2753nb_reaction_3 = nb_reaction_3 + 1
2754
2755indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
2756
2757!===========================================================
2758!      a003 : O + O3 -> O2 + O2
2759!===========================================================
2760
2761nb_reaction_4 = nb_reaction_4 + 1
2762
2763indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
2764
2765!===========================================================
2766!      b001 : O(1D) + CO2 -> O + CO2
2767!===========================================================
2768
2769nb_phot = nb_phot + 1
2770
2771indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
2772
2773!===========================================================
2774!      b002 : O(1D) + H2O -> OH + OH
2775!===========================================================
2776
2777nb_reaction_4 = nb_reaction_4 + 1
2778
2779indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
2780
2781!===========================================================
2782!      b003 : O(1D) + H2 -> OH + H
2783!===========================================================
2784
2785nb_reaction_4 = nb_reaction_4 + 1
2786
2787indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
2788
2789!===========================================================
2790!      b004 : O(1D) + O2 -> O + O2
2791!===========================================================
2792
2793nb_phot = nb_phot + 1
2794
2795indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
2796
2797!===========================================================
2798!      b005 : O(1D) + O3 -> O2 + O2
2799!===========================================================
2800
2801nb_reaction_4 = nb_reaction_4 + 1
2802
2803indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
2804
2805!===========================================================
2806!      b006 : O(1D) + O3 -> O2 + O + O
2807!===========================================================
2808
2809nb_reaction_4 = nb_reaction_4 + 1
2810
2811indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
2812
2813
2814if(deutchem) then
2815!===========================================================
2816!      b010 : O(1D) + HDO -> OD + OH
2817!===========================================================
2818
2819   nb_reaction_4 = nb_reaction_4 + 1
2820
2821   indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hdo, 1.0, i_od, 1.0, i_oh)
2822
2823
2824!===========================================================
2825!      b011 : O(1D) + HD -> H + OD
2826!===========================================================
2827
2828   nb_reaction_4 = nb_reaction_4 + 1
2829
2830   indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hd, 1.0, i_h, 1.0, i_od)
2831
2832
2833!===========================================================
2834!      b012 : O(1D) + HD -> D + OH
2835!===========================================================
2836
2837   nb_reaction_4 = nb_reaction_4 + 1
2838
2839   indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hd, 1.0, i_d, 1.0, i_oh)
2840
2841endif !deutchem
2842
2843!===========================================================
2844!      c001 : O + HO2 -> OH + O2
2845!===========================================================
2846
2847nb_reaction_4 = nb_reaction_4 + 1
2848
2849indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
2850
2851!===========================================================
2852!      c002 : O + OH -> O2 + H
2853!===========================================================
2854
2855nb_reaction_4 = nb_reaction_4 + 1
2856
2857indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
2858
2859!===========================================================
2860!      c003 : H + O3 -> OH + O2
2861!===========================================================
2862
2863nb_reaction_4 = nb_reaction_4 + 1
2864
2865indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
2866
2867!===========================================================
2868!      c004 : H + HO2 -> OH + OH
2869!===========================================================
2870
2871nb_reaction_4 = nb_reaction_4 + 1
2872
2873indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
2874
2875!===========================================================
2876!      c005 : H + HO2 -> H2 + O2
2877!===========================================================
2878
2879nb_reaction_4 = nb_reaction_4 + 1
2880
2881indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
2882
2883!===========================================================
2884!      c006 : H + HO2 -> H2O + O
2885!===========================================================
2886
2887nb_reaction_4 = nb_reaction_4 + 1
2888
2889indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
2890
2891!===========================================================
2892!      c007 : OH + HO2 -> H2O + O2
2893!===========================================================
2894
2895nb_reaction_4 = nb_reaction_4 + 1
2896
2897indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
2898
2899!===========================================================
2900!      c008 : HO2 + HO2 -> H2O2 + O2
2901!===========================================================
2902
2903nb_reaction_3 = nb_reaction_3 + 1
2904
2905indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2906
2907!===========================================================
2908!      c009 : OH + H2O2 -> H2O + HO2
2909!===========================================================
2910
2911nb_reaction_4 = nb_reaction_4 + 1
2912
2913indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
2914
2915!===========================================================
2916!      c010 : OH + H2 -> H2O + H
2917!===========================================================
2918
2919nb_reaction_4 = nb_reaction_4 + 1
2920
2921indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
2922
2923!===========================================================
2924!      c011 : H + O2 + CO2 -> HO2 + CO2
2925!===========================================================
2926
2927nb_reaction_4 = nb_reaction_4 + 1
2928
2929indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
2930
2931!===========================================================
2932!      c012 : O + H2O2 -> OH + HO2
2933!===========================================================
2934
2935nb_reaction_4 = nb_reaction_4 + 1
2936
2937indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
2938
2939!===========================================================
2940!      c013 : OH + OH -> H2O + O
2941!===========================================================
2942
2943nb_reaction_3 = nb_reaction_3 + 1
2944
2945indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
2946
2947!===========================================================
2948!      c014 : OH + O3 -> HO2 + O2
2949!===========================================================
2950
2951nb_reaction_4 = nb_reaction_4 + 1
2952
2953indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
2954
2955!===========================================================
2956!      c015 : HO2 + O3 -> OH + O2 + O2
2957!===========================================================
2958
2959nb_reaction_4 = nb_reaction_4 + 1
2960
2961indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
2962
2963!===========================================================
2964!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
2965!===========================================================
2966
2967nb_reaction_3 = nb_reaction_3 + 1
2968
2969indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
2970
2971!===========================================================
2972!      c017 : OH + OH + CO2 -> H2O2 + CO2
2973!===========================================================
2974
2975nb_reaction_3 = nb_reaction_3 + 1
2976
2977indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
2978
2979!===========================================================
2980!      c018 : H + H + CO2 -> H2 + CO2
2981!===========================================================
2982
2983nb_reaction_3 = nb_reaction_3 + 1
2984
2985indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
2986
2987
2988!===========================================================
2989!      d001 : NO2 + O -> NO + O2
2990!===========================================================
2991
2992nb_reaction_4 = nb_reaction_4 + 1
2993
2994indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
2995
2996!===========================================================
2997!      d002 : NO + O3 -> NO2 + O2
2998!===========================================================
2999
3000nb_reaction_4 = nb_reaction_4 + 1
3001
3002indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
3003
3004!===========================================================
3005!      d003 : NO + HO2 -> NO2 + OH
3006!===========================================================
3007
3008nb_reaction_4 = nb_reaction_4 + 1
3009
3010indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
3011
3012!===========================================================
3013!      d004 : N + NO -> N2 + O
3014!===========================================================
3015
3016nb_reaction_4 = nb_reaction_4 + 1
3017
3018indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
3019
3020!===========================================================
3021!      d005 : N + O2 -> NO + O
3022!===========================================================
3023
3024nb_reaction_4 = nb_reaction_4 + 1
3025
3026indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
3027
3028!===========================================================
3029!      d006 : NO2 + H -> NO + OH
3030!===========================================================
3031
3032nb_reaction_4 = nb_reaction_4 + 1
3033
3034indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
3035
3036!===========================================================
3037!      d007 : N + O -> NO
3038!===========================================================
3039
3040nb_reaction_4 = nb_reaction_4 + 1
3041
3042indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
3043
3044!===========================================================
3045!      d008 : N + HO2 -> NO + OH
3046!===========================================================
3047
3048nb_reaction_4 = nb_reaction_4 + 1
3049
3050indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
3051
3052!===========================================================
3053!      d009 : N + OH -> NO + H
3054!===========================================================
3055
3056nb_reaction_4 = nb_reaction_4 + 1
3057
3058indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
3059
3060!===========================================================
3061!      d010 : N(2D) + O -> N + O
3062!===========================================================
3063
3064nb_phot = nb_phot + 1
3065
3066indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
3067
3068!===========================================================
3069!      d011 : N(2D) + N2 -> N + N2
3070!===========================================================
3071
3072nb_phot = nb_phot + 1
3073
3074indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
3075
3076!===========================================================
3077!      d012 : N(2D) + CO2 -> NO + CO
3078!===========================================================
3079
3080nb_reaction_4 = nb_reaction_4 + 1
3081
3082indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
3083
3084!===========================================================
3085!      e001 : CO + OH -> CO2 + H
3086!===========================================================
3087
3088nb_reaction_4 = nb_reaction_4 + 1
3089
3090indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
3091
3092!===========================================================
3093!      e002 : CO + O + M -> CO2 + M
3094!===========================================================
3095
3096nb_reaction_4 = nb_reaction_4 + 1
3097
3098indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
3099if(deutchem) then
3100!===========================================================
3101!      f001 : OD + OH -> HDO + O
3102!===========================================================
3103
3104   nb_reaction_4 = nb_reaction_4 + 1
3105
3106   indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_oh, 1.0, i_hdo, 1.0, i_o)
3107
3108!===========================================================
3109!      f002 : OD + H2 -> HDO + H
3110!===========================================================
3111
3112   nb_reaction_4 = nb_reaction_4 + 1
3113
3114   indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_h2, 1.0, i_hdo, 1.0, i_h)
3115
3116!===========================================================
3117!      f003 : OD + HO2 -> HDO + O2
3118!===========================================================
3119
3120   nb_reaction_4 = nb_reaction_4 + 1
3121
3122   indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_ho2, 1.0, i_hdo, 1.0, i_o2)
3123
3124!===========================================================
3125!      f004 : OD + H2O2 -> HDO + HO2
3126!===========================================================
3127   
3128   nb_reaction_4 = nb_reaction_4 + 1
3129
3130   indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_h2o2, 1.0, i_hdo, 1.0, i_ho2)
3131
3132!===========================================================
3133!      f005 : O + OD -> O2 + D
3134!===========================================================
3135
3136   nb_reaction_4 = nb_reaction_4 + 1
3137
3138   indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_od, 1.0, i_o2, 1.0, i_d)
3139
3140!===========================================================
3141!      f006 : OD + H2 -> H2O + D
3142!===========================================================
3143
3144   nb_reaction_4 = nb_reaction_4 + 1
3145
3146   indice_4(nb_reaction_4) = z4spec(1.0, i_h2, 1.0, i_od, 1.0, i_h2o, 1.0, i_d)
3147
3148!===========================================================
3149!      f007 : OD + H -> OH + D
3150!===========================================================
3151
3152   nb_reaction_4 = nb_reaction_4 + 1
3153
3154   indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_od, 1.0, i_oh, 1.0, i_d)
3155
3156!===========================================================
3157!      f008 : CO + OD -> CO2 + D
3158!===========================================================
3159
3160   nb_reaction_4 = nb_reaction_4 + 1
3161
3162   indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_od, 1.0, i_co2, 1.0, i_d)
3163
3164!===========================================================
3165!      f009 : O3 + D -> O2 + OD
3166!===========================================================
3167
3168   nb_reaction_4 = nb_reaction_4 + 1
3169
3170   indice_4(nb_reaction_4) = z4spec(1.0, i_o3, 1.0, i_d, 1.0, i_o2, 1.0, i_od)
3171
3172!===========================================================
3173!      f010 : HO2 + D -> OH + OD
3174!===========================================================
3175
3176   nb_reaction_4 = nb_reaction_4 + 1
3177
3178   indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_d, 1.0, i_oh, 1.0, i_od)
3179
3180!===========================================================
3181!      f011 : HO2 + D -> HDO + O
3182!===========================================================
3183
3184   nb_reaction_4 = nb_reaction_4 + 1
3185
3186   indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_d, 1.0, i_hdo, 1.0, i_o)
3187
3188!===========================================================
3189!      f012 : OH + D -> H + OD
3190!===========================================================
3191
3192   nb_reaction_4 = nb_reaction_4 + 1
3193
3194   indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_d, 1.0, i_h, 1.0, i_od)
3195
3196!===========================================================
3197!      f013 : H + D + CO2 -> HD + CO2
3198!===========================================================
3199
3200   nb_reaction_4 = nb_reaction_4 + 1
3201
3202   indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_d, 1.0, i_hd, 0.0, i_dummy)
3203
3204!===========================================================
3205!      f014 : D + HO2 -> HD + O2
3206!===========================================================
3207
3208   nb_reaction_4 = nb_reaction_4 + 1
3209
3210   indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_ho2, 1.0, i_hd, 1.0, i_o2)
3211
3212
3213!===========================================================
3214!      f015 : OH + HD -> HDO + H
3215!===========================================================
3216
3217   nb_reaction_4 = nb_reaction_4 + 1
3218
3219   indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hd, 1.0, i_hdo, 1.0, i_h)
3220
3221!===========================================================
3222!      f016 : OH + HD -> H2O + D
3223!===========================================================
3224
3225   nb_reaction_4 = nb_reaction_4 + 1
3226
3227   indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hd, 1.0, i_h2o, 1.0, i_d)
3228
3229!===========================================================
3230!      f017 : D + O2 + CO2 -> DO2 + CO2
3231!===========================================================
3232
3233   nb_reaction_4 = nb_reaction_4 + 1
3234
3235   indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_o2, 1.0, i_do2, 0.0, i_dummy)
3236
3237!===========================================================
3238!      f018 : OD + O3 -> DO2 + O2
3239!===========================================================
3240
3241   nb_reaction_4 = nb_reaction_4 + 1
3242
3243   indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_o3, 1.0, i_do2, 1.0, i_o2)
3244
3245!===========================================================
3246!      f019 : D + HO2 -> DO2 + H
3247!===========================================================
3248
3249   nb_reaction_4 = nb_reaction_4 + 1
3250
3251   indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_ho2, 1.0, i_do2, 1.0, i_h)
3252
3253!===========================================================
3254!      f020 : O + DO2 -> OD + O2
3255!===========================================================
3256
3257   nb_reaction_4 = nb_reaction_4 + 1
3258
3259   indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_do2, 1.0, i_od, 1.0, i_o2)
3260
3261!===========================================================
3262!      f021 : H + DO2 -> OH + OD
3263!===========================================================
3264   
3265   nb_reaction_4 = nb_reaction_4 + 1
3266
3267   indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_od, 1.0, i_oh)
3268
3269!===========================================================
3270!      f022 : H + DO2 -> HD + O2
3271!===========================================================
3272
3273   nb_reaction_4 = nb_reaction_4 + 1
3274
3275   indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_hd, 1.0, i_o2)
3276
3277!===========================================================
3278!      f023 : H + DO2 -> HDO + O
3279!===========================================================
3280
3281   nb_reaction_4 = nb_reaction_4 + 1
3282
3283   indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_hdo, 1.0, i_o)
3284
3285!===========================================================
3286!      f024 : H + DO2 -> HO2 + D
3287!===========================================================
3288
3289   nb_reaction_4 = nb_reaction_4 + 1
3290
3291   indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_ho2, 1.0, i_d)
3292
3293!===========================================================
3294!      f025 : OH + DO2 -> HDO + O2
3295!===========================================================
3296
3297   nb_reaction_4 = nb_reaction_4 + 1
3298
3299   indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_do2, 1.0, i_hdo, 1.0, i_o2)
3300
3301!===========================================================
3302!      f026 : DO2 + O3 -> OD + O2 + O2
3303!===========================================================
3304
3305   nb_reaction_4 = nb_reaction_4 + 1
3306
3307   indice_4(nb_reaction_4) = z4spec(1.0, i_do2, 1.0, i_o3, 1.0, i_od, 2.0, i_o2)
3308
3309!===========================================================
3310!      f027 : OD + OH + CO2 -> HDO2 + CO2
3311!===========================================================
3312
3313   nb_reaction_4 = nb_reaction_4 + 1
3314
3315   indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_oh, 1.0, i_hdo2, 0.0, i_dummy)
3316
3317!===========================================================
3318!      f028 : DO2 + HO2 -> HDO2 + O2
3319!===========================================================
3320
3321   nb_reaction_4 = nb_reaction_4 + 1
3322
3323   indice_4(nb_reaction_4) = z4spec(1.0, i_do2, 1.0, i_ho2, 1.0, i_hdo2, 1.0, i_o2)
3324
3325!===========================================================
3326!      f029 : O + HDO2 -> OD + HO2
3327!===========================================================
3328
3329   nb_reaction_4 = nb_reaction_4 + 1
3330
3331   indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_hdo2, 1.0, i_od, 1.0, i_ho2)
3332
3333!===========================================================
3334!      f030 : O + HDO2 -> OH + DO2
3335!===========================================================
3336
3337   nb_reaction_4 = nb_reaction_4 + 1
3338
3339   indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_hdo2, 1.0, i_oh, 1.0, i_do2)
3340
3341!===========================================================
3342!      f031 : OH + HDO2 -> HDO + HO2
3343!===========================================================
3344
3345   nb_reaction_4 = nb_reaction_4 + 1
3346
3347   indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hdo2, 1.0, i_hdo, 1.0, i_ho2)
3348
3349!===========================================================
3350!      f032 : OH + HDO2 -> H2O + DO2
3351!===========================================================
3352
3353   nb_reaction_4 = nb_reaction_4 + 1
3354
3355   indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hdo2, 1.0, i_h2o, 1.0, i_do2)
3356
3357endif !deutchem
3358
3359!Only if ion chemistry
3360if (ionchem) then
3361
3362!===========================================================
3363!      i001 : CO2+ + O2 -> O2+ + CO2
3364!===========================================================
3365
3366   nb_reaction_4 = nb_reaction_4 + 1
3367
3368   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
3369
3370!===========================================================
3371!      i002 : CO2+ + O -> O+ + CO2
3372!===========================================================
3373
3374   nb_reaction_4 = nb_reaction_4 + 1
3375
3376   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
3377
3378!===========================================================
3379!      i003 : CO2+ + O -> O2+ + CO
3380!===========================================================
3381
3382   nb_reaction_4 = nb_reaction_4 + 1
3383
3384   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
3385
3386!===========================================================
3387!      i004 : O2+ + e- -> O + O
3388!===========================================================
3389
3390   nb_reaction_4 = nb_reaction_4 + 1
3391
3392   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
3393
3394!===========================================================
3395!      i005 : O+ + CO2 -> O2+ + CO
3396!===========================================================
3397
3398   nb_reaction_4 = nb_reaction_4 + 1
3399
3400   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
3401
3402!===========================================================
3403!      i006 : CO2+ + e -> CO + O
3404!===========================================================
3405
3406   nb_reaction_4 = nb_reaction_4 + 1
3407
3408   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
3409
3410!===========================================================
3411!      i007 : CO2+ + NO -> NO+ + CO2
3412!===========================================================
3413
3414   nb_reaction_4 = nb_reaction_4 + 1
3415
3416   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
3417
3418!===========================================================
3419!      i008 : O2+ + NO -> NO+ + O2
3420!===========================================================
3421
3422   nb_reaction_4 = nb_reaction_4 + 1
3423
3424   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
3425
3426!===========================================================
3427!      i009 : O2+ + N2 -> NO+ + NO
3428!===========================================================
3429
3430   nb_reaction_4 = nb_reaction_4 + 1
3431
3432   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
3433
3434!===========================================================
3435!      i010 : O2+ + N -> NO+ + O
3436!===========================================================
3437
3438   nb_reaction_4 = nb_reaction_4 + 1
3439
3440   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
3441
3442!===========================================================
3443!      i011 : O+ + N2 -> NO+ + N
3444!===========================================================
3445
3446   nb_reaction_4 = nb_reaction_4 + 1
3447
3448   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
3449
3450!===========================================================
3451!      i012 : NO+ + e -> N + O
3452!===========================================================
3453
3454   nb_reaction_4 = nb_reaction_4 + 1
3455
3456   indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
3457
3458!===========================================================
3459!      i013 : CO+ + CO2 -> CO2+ + CO
3460!===========================================================
3461
3462   nb_reaction_4 = nb_reaction_4 + 1
3463
3464   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
3465
3466!===========================================================
3467!      i014 : CO+ + O -> O+ + CO
3468!===========================================================
3469
3470   nb_reaction_4 = nb_reaction_4 + 1
3471
3472   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
3473
3474!===========================================================
3475!      i015 : C+ + CO2 -> CO+ + CO
3476!===========================================================
3477
3478   nb_reaction_4 = nb_reaction_4 + 1
3479
3480   indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
3481
3482!===========================================================
3483!      i016 : N2+ + CO2 -> CO2+ + N2
3484!===========================================================
3485
3486   nb_reaction_4 = nb_reaction_4 + 1
3487
3488   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
3489
3490!===========================================================
3491!      i017 : N2+ + O -> NO+ + N
3492!===========================================================
3493
3494   nb_reaction_4 = nb_reaction_4 + 1
3495
3496   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
3497
3498!===========================================================
3499!      i018 : N2+ + CO -> CO+ + N2
3500!===========================================================
3501
3502   nb_reaction_4 = nb_reaction_4 + 1
3503
3504   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
3505
3506!===========================================================
3507!      i019 : N2+ + e -> N + N
3508!===========================================================
3509
3510   nb_reaction_4 = nb_reaction_4 + 1
3511
3512   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
3513
3514!===========================================================
3515!      i020 : N2+ + O -> O+ + N2
3516!===========================================================
3517
3518   nb_reaction_4 = nb_reaction_4 + 1
3519
3520   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
3521
3522!===========================================================
3523!      i021 : N+ + CO2 -> CO2+ + N
3524!===========================================================
3525
3526   nb_reaction_4 = nb_reaction_4 + 1
3527
3528   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
3529
3530!===========================================================
3531!      i022 : CO+ + H -> H+ + CO
3532!===========================================================
3533
3534   nb_reaction_4 = nb_reaction_4 + 1
3535
3536   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
3537
3538!===========================================================
3539!      i023 : O+ + H -> H+ + O
3540!===========================================================
3541
3542   nb_reaction_4 = nb_reaction_4 + 1
3543
3544   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
3545
3546!===========================================================
3547!      i024 : H+ + O -> O+ + H
3548!===========================================================
3549
3550   nb_reaction_4 = nb_reaction_4 + 1
3551
3552   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
3553
3554!===========================================================
3555!      i025 : CO2+ + H2 -> HCO2+ + H
3556!===========================================================
3557
3558   nb_reaction_4 = nb_reaction_4 + 1
3559
3560   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
3561
3562!===========================================================
3563!      i026 : spare slot (reaction rate set to zero)
3564!===========================================================
3565
3566   nb_reaction_4 = nb_reaction_4 + 1
3567
3568   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
3569
3570!===========================================================
3571!      i027 : HCO2+ + e -> H + O + CO
3572!===========================================================
3573!We divide this reaction in two
3574
3575!0.5HCO2+ + 0.5e -> H
3576
3577   nb_reaction_4 = nb_reaction_4 + 1
3578
3579   indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
3580
3581!0.5 HCO2+ + 0.5 e -> O + CO
3582
3583   nb_reaction_4 = nb_reaction_4 + 1
3584
3585   indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
3586
3587!===========================================================
3588!      i029 : HCO2+ + e -> OH + CO
3589!===========================================================
3590
3591   nb_reaction_4 = nb_reaction_4 + 1
3592
3593   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
3594
3595!===========================================================
3596!      i030 : HCO2+ + e -> H + CO2
3597!===========================================================
3598
3599   nb_reaction_4 = nb_reaction_4 + 1
3600
3601   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
3602
3603!===========================================================
3604!      i031 : HCO2+ + O -> HCO+ + O2
3605!===========================================================
3606
3607   nb_reaction_4 = nb_reaction_4 + 1
3608
3609   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_o, 1.0, i_hcoplus, 1.0, i_o2)
3610
3611!===========================================================
3612!      i032 : HCO2+ + CO -> HCO+ + CO2
3613!===========================================================
3614
3615   nb_reaction_4 = nb_reaction_4 + 1
3616   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_co2)
3617
3618!===========================================================
3619!      i033 : H+ + CO2 -> HCO+ + O
3620!===========================================================
3621
3622   nb_reaction_4 = nb_reaction_4 + 1
3623   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_co2, 1.0, i_hcoplus, 1.0, i_o)
3624
3625!===========================================================
3626!      i034 : CO2+ + H -> HCO+ + O
3627!===========================================================
3628
3629   nb_reaction_4 = nb_reaction_4 + 1
3630   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h, 1.0, i_hcoplus, 1.0, i_o)
3631
3632!===========================================================
3633!      i035 : CO+ + H2 -> HCO+ + H
3634!===========================================================
3635
3636   nb_reaction_4 = nb_reaction_4 + 1
3637   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2, 1.0, i_hcoplus, 1.0, i_h)
3638
3639
3640!===========================================================
3641!      i036 : HCO+ + e- -> CO + H
3642!===========================================================
3643
3644   nb_reaction_4 = nb_reaction_4 + 1
3645   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_elec, 1.0, i_co, 1.0, i_h)
3646
3647!===========================================================
3648!      i037 : CO2+ + H2O -> H2O+ + CO2
3649!===========================================================
3650
3651   nb_reaction_4 = nb_reaction_4 + 1
3652   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co2)
3653
3654!===========================================================
3655!      i038 : CO+ + H2O -> H2O+ + CO
3656!===========================================================
3657
3658   nb_reaction_4 = nb_reaction_4 + 1
3659   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co)
3660
3661!===========================================================
3662!      i039 : O+ + H2O -> H2O+ + O
3663!===========================================================
3664
3665   nb_reaction_4 = nb_reaction_4 + 1
3666   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_o)
3667
3668!===========================================================
3669!      i040 : N2+ + H2O -> H2O+ + N2
3670!===========================================================
3671
3672   nb_reaction_4 = nb_reaction_4 + 1
3673   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n2)
3674
3675!===========================================================
3676!      i041 : N+ + H2O -> H2O+ + N
3677!===========================================================
3678
3679   nb_reaction_4 = nb_reaction_4 + 1
3680   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n)
3681
3682!===========================================================
3683!      i042 : H+ + H2O -> H2O+ + H
3684!===========================================================
3685
3686   nb_reaction_4 = nb_reaction_4 + 1
3687   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_h)
3688
3689!===========================================================
3690!      i043 : H2O+ + O2 -> O2+ + H2O
3691!===========================================================
3692
3693   nb_reaction_4 = nb_reaction_4 + 1
3694   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_h2o)
3695
3696!===========================================================
3697!      i044 : H2O+ + CO -> HCO+ + OH
3698!===========================================================
3699
3700   nb_reaction_4 = nb_reaction_4 + 1
3701   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_oh)
3702
3703!===========================================================
3704!      i045 : H2O+ + O -> O2+ + H2
3705!===========================================================
3706
3707   nb_reaction_4 = nb_reaction_4 + 1
3708   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h2)
3709
3710!===========================================================
3711!      i046 : H2O+ + NO -> NO+ + H2O
3712!===========================================================
3713
3714   nb_reaction_4 = nb_reaction_4 + 1
3715   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_h2o)
3716
3717!===========================================================
3718!      i047 : H2O+ + e- -> H + H + O
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, 2.0, i_h, 1.0, i_o)
3723
3724!===========================================================
3725!      i048 : H2O+ + e- -> H + OH
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_h, 1.0, i_oh)
3730
3731!===========================================================
3732!      i049 : H2O+ + e- -> H2 + O
3733!===========================================================
3734
3735   nb_reaction_4 = nb_reaction_4 + 1
3736   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o)
3737
3738!===========================================================
3739!      i050 : H2O+ + H2O -> H3O+ + OH
3740!===========================================================
3741
3742   nb_reaction_4 = nb_reaction_4 + 1
3743   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_oh)
3744
3745!===========================================================
3746!      i051 : H2O+ + H2 -> H3O+ + H
3747!===========================================================
3748
3749   nb_reaction_4 = nb_reaction_4 + 1
3750   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2, 1.0, i_h3oplus, 1.0, i_h)
3751
3752!===========================================================
3753!      i052 : HCO+ + H2O -> H3O+ + CO
3754!===========================================================
3755
3756   nb_reaction_4 = nb_reaction_4 + 1
3757   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_co)
3758
3759!===========================================================
3760!      i053: H3O+ + e -> OH + H + 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_oh, 2.0, i_h)
3765
3766!===========================================================
3767!      i054: H3O+ + e -> H2O + H
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_h2o, 1.0, i_h)
3772
3773!===========================================================
3774!      i055: H3O+ + e -> HO + H2
3775!===========================================================
3776
3777   nb_reaction_4 = nb_reaction_4 + 1
3778   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 1.0, i_h2)
3779
3780!===========================================================
3781!      i056: H3O+ + e -> O + H2 + H
3782!===========================================================
3783!We divide this reaction in two
3784
3785!0.5H3O+ + 0.5e -> O
3786
3787   nb_reaction_4 = nb_reaction_4 + 1
3788   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_o, 0.0, i_dummy)
3789
3790!0.5H3O+ + 0.5e -> H2 + H
3791
3792   nb_reaction_4 = nb_reaction_4 + 1
3793   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_h2, 1.0, i_h)
3794
3795!===========================================================
3796!      i057: O+ + H2 -> OH+ + H
3797!===========================================================
3798
3799   nb_reaction_4 = nb_reaction_4 + 1
3800   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2, 1.0, i_ohplus, 1.0, i_h)
3801
3802!===========================================================
3803!      i058: OH+ + O -> O2+ + H
3804!===========================================================
3805
3806   nb_reaction_4 = nb_reaction_4 + 1
3807   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h)
3808
3809!===========================================================
3810!      i059: OH+ + CO2 -> HCO2+ + 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_co2, 1.0, i_hco2plus, 1.0, i_o)
3815
3816!===========================================================
3817!      i060: OH+ + CO -> HCO+ + O
3818!===========================================================
3819
3820   nb_reaction_4 = nb_reaction_4 + 1
3821   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_o)
3822
3823!===========================================================
3824!      i061: OH+ + NO -> NO+ + OH
3825!===========================================================
3826
3827   nb_reaction_4 = nb_reaction_4 + 1
3828   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_oh)
3829
3830!===========================================================
3831!      i062: OH+ + H2 -> H2O+ + H
3832!===========================================================
3833
3834   nb_reaction_4 = nb_reaction_4 + 1
3835   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_h2, 1.0, i_h2oplus, 1.0, i_h)
3836
3837!===========================================================
3838!      i063: OH+ + O2 -> O2+ + OH
3839!===========================================================
3840
3841   nb_reaction_4 = nb_reaction_4 + 1
3842   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_oh)
3843
3844end if    !ionchem
3845
3846!===========================================================
3847!      h001: HO2 + ice -> products
3848!            treated as
3849!            HO2 -> 0.5 H2O + 0.75 O2
3850!===========================================================
3851
3852nb_phot = nb_phot + 1
3853
3854indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
3855
3856!===========================================================
3857!      h002: OH + ice -> products
3858!            treated as
3859!            OH -> 0.5 H2O + 0.25 O2
3860!===========================================================
3861
3862nb_phot = nb_phot + 1
3863
3864indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2)
3865
3866!===========================================================
3867!      h003: H2O2 + ice -> products
3868!            treated as
3869!            H2O2 -> H2O + 0.5 O2
3870!===========================================================
3871
3872nb_phot = nb_phot + 1
3873
3874indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
3875
3876!===========================================================
3877!      h004: HO2 + dust -> products
3878!            treated as
3879!            HO2 -> 0.5 H2O + 0.75 O2
3880!===========================================================
3881
3882nb_phot = nb_phot + 1
3883
3884indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
3885
3886!===========================================================
3887!      h005: H2O2 + dust -> products
3888!            treated as
3889!            H2O2 -> H2O + 0.5 O2
3890!===========================================================
3891
3892nb_phot = nb_phot + 1
3893
3894indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
3895
3896!===========================================================
3897!  check dimensions
3898!===========================================================
3899
3900print*, 'nb_phot       = ', nb_phot
3901print*, 'nb_reaction_4 = ', nb_reaction_4
3902print*, 'nb_reaction_3 = ', nb_reaction_3
3903
3904if ((nb_phot /= nb_phot_max)             .or.  &
3905    (nb_reaction_3 /= nb_reaction_3_max) .or.  &
3906    (nb_reaction_4 /= nb_reaction_4_max)) then
3907   print*, 'wrong dimensions in indice'
3908   call abort_physic("indice","wrong array dimensions",1)
3909end if 
3910
3911end subroutine indice
3912
3913!*****************************************************************
3914
3915      subroutine gcmtochim(nlayer, ionchem, deutchem, nq, zycol,     &
3916                           lswitch, nesp,                            &
3917                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
3918                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
3919                           i_n, i_n2d, i_no, i_no2, i_n2,            &
3920                           i_co2plus, i_oplus, i_o2plus, i_noplus,   &
3921                           i_coplus, i_cplus, i_n2plus, i_nplus,     &
3922                           i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,&
3923                           i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, &
3924                           i_d, i_hd, i_do2, i_hdo2, dens, rm, c)
3925       
3926!*****************************************************************
3927
3928      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
3929     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,  &
3930     &                      igcm_ho2, igcm_h2o2, igcm_h2o_vap,           &
3931     &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2,&
3932     &                      igcm_co2plus, igcm_oplus, igcm_o2plus,       &
3933     &                      igcm_noplus, igcm_coplus, igcm_cplus,        &
3934     &                      igcm_n2plus, igcm_nplus, igcm_hplus,         &
3935     &                      igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,   &
3936     &                      igcm_h3oplus, igcm_ohplus, igcm_elec,        &
3937     &                      igcm_hdo_vap, igcm_od, igcm_d, igcm_hd,      &
3938     &                      igcm_do2, igcm_hdo2
3939
3940      implicit none
3941
3942      include "callkeys.h"
3943
3944!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3945!     input:
3946!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3947     
3948      integer, intent(in) :: nlayer ! number of atmospheric layers
3949      integer, intent(in) :: nq     ! number of tracers in the gcm
3950      logical, intent(in) :: ionchem
3951      logical, intent(in) :: deutchem
3952      integer :: nesp               ! number of species in the chemistry
3953      integer :: lswitch            ! interface level between chemistries
3954
3955      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,           &
3956                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                   &
3957                 i_n, i_n2d, i_no, i_no2, i_n2,                      &
3958                 i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus,   &
3959                 i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus,    &
3960                 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec,  &
3961                 i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2
3962
3963      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
3964      real :: dens(nlayer)          ! total number density (molecule.cm-3)
3965
3966!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3967!     output:
3968!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3969
3970      real, dimension(nlayer,nesp)            :: rm ! volume mixing ratios
3971      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3)
3972     
3973!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3974!     local:
3975!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3976
3977      integer      :: l, iesp
3978      logical,save :: firstcall = .true.
3979
3980!$OMP THREADPRIVATE(firstcall)
3981     
3982!     first call initializations
3983
3984      if (firstcall) then
3985
3986!       identify the indexes of the tracers we need
3987
3988         if (igcm_co2 == 0) then
3989            write(*,*) "gcmtochim: Error; no CO2 tracer !!!"
3990            call abort_physic("gcmtochim","missing co2 tracer",1)
3991         endif
3992         if (igcm_co == 0) then
3993            write(*,*) "gcmtochim: Error; no CO tracer !!!"
3994            call abort_physic("gcmtochim","missing co tracer",1)
3995         end if
3996         if (igcm_o == 0) then
3997            write(*,*) "gcmtochim: Error; no O tracer !!!"
3998            call abort_physic("gcmtochim","missing o tracer",1)
3999         end if
4000         if (igcm_o1d == 0) then
4001            write(*,*) "gcmtochim: Error; no O1D tracer !!!"
4002            call abort_physic("gcmtochim","missing o1d tracer",1)
4003         end if
4004         if (igcm_o2 == 0) then
4005            write(*,*) "gcmtochim: Error; no O2 tracer !!!"
4006            call abort_physic("gcmtochim","missing o2 tracer",1)
4007         end if
4008         if (igcm_o3 == 0) then
4009            write(*,*) "gcmtochim: Error; no O3 tracer !!!"
4010            call abort_physic("gcmtochim","missing o3 tracer",1)
4011         end if
4012         if (igcm_h == 0) then
4013            write(*,*) "gcmtochim: Error; no H tracer !!!"
4014            call abort_physic("gcmtochim","missing h tracer",1)
4015         end if
4016         if (igcm_h2 == 0) then
4017            write(*,*) "gcmtochim: Error; no H2 tracer !!!"
4018            call abort_physic("gcmtochim","missing h2 tracer",1)
4019         end if
4020         if (igcm_oh == 0) then
4021            write(*,*) "gcmtochim: Error; no OH tracer !!!"
4022            call abort_physic("gcmtochim","missing oh tracer",1)
4023         end if
4024         if (igcm_ho2 == 0) then
4025            write(*,*) "gcmtochim: Error; no HO2 tracer !!!"
4026            call abort_physic("gcmtochim","missing ho2 tracer",1)
4027         end if
4028         if (igcm_h2o2 == 0) then
4029            write(*,*) "gcmtochim: Error; no H2O2 tracer !!!"
4030            call abort_physic("gcmtochim","missing h2o2 tracer",1)
4031         end if
4032         if (igcm_n == 0) then
4033            write(*,*) "gcmtochim: Error; no N tracer !!!"
4034            call abort_physic("gcmtochim","missing n tracer",1)
4035         end if
4036         if (igcm_n2d == 0) then
4037            write(*,*) "gcmtochim: Error; no N2D tracer !!!"
4038            call abort_physic("gcmtochim","missing n2d tracer",1)
4039         end if
4040         if (igcm_no == 0) then
4041            write(*,*) "gcmtochim: Error; no NO tracer !!!"
4042            call abort_physic("gcmtochim","missing no tracer",1)
4043         end if
4044         if (igcm_no2 == 0) then
4045            write(*,*) "gcmtochim: Error; no NO2 tracer !!!"
4046            call abort_physic("gcmtochim","missing no2 tracer",1)
4047         end if
4048         if (igcm_n2 == 0) then
4049            write(*,*) "gcmtochim: Error; no N2 tracer !!!"
4050            call abort_physic("gcmtochim","missing n2 tracer",1)
4051         end if
4052         if (igcm_h2o_vap == 0) then
4053            write(*,*) "gcmtochim: Error; no water vapor tracer !!!"
4054            call abort_physic("gcmtochim","missing h2o_vap tracer",1)
4055         end if
4056         if(deutchem) then
4057            if (igcm_hdo_vap == 0) then
4058               write(*,*) "gcmtochim: Error; no HDO tracer !!!"
4059               call abort_physic("gcmtochim","missing hdo_vap tracer",1)
4060            end if
4061            if (igcm_od == 0) then
4062               write(*,*) "gcmtochim: Error; no OD tracer !!!"
4063               call abort_physic("gcmtochim","missing od tracer",1)
4064            end if
4065            if (igcm_d == 0) then
4066               write(*,*) "gcmtochim: Error; no D tracer !!!"
4067               call abort_physic("gcmtochim","missing d tracer",1)
4068            end if
4069            if (igcm_hd == 0) then
4070               write(*,*) "gcmtochim: Error; no HD tracer !!!"
4071               call abort_physic("gcmtochim","missing hd tracer",1)
4072            end if
4073            if (igcm_do2 == 0) then
4074               write(*,*) "gcmtochim: Error; no DO2 tracer !!!"
4075               call abort_physic("gcmtochim","missing do2 tracer",1)
4076            end if
4077            if (igcm_hdo2 == 0) then
4078               write(*,*) "gcmtochim: Error; no HDO2 tracer !!!"
4079               call abort_physic("gcmtochim","missing hdo2 tracer",1)
4080            end if
4081         endif !deutchem
4082         if (ionchem) then
4083            if (igcm_co2plus == 0) then
4084               write(*,*) "gcmtochim: Error; no CO2+ tracer !!!"
4085               call abort_physic("gcmtochim","missing co2plus tracer",1)
4086            end if
4087            if (igcm_oplus == 0) then
4088               write(*,*) "gcmtochim: Error; no O+ tracer !!!"
4089               call abort_physic("gcmtochim","missing oplus tracer",1)
4090            end if
4091            if (igcm_o2plus == 0) then
4092               write(*,*) "gcmtochim: Error; no O2+ tracer !!!"
4093               call abort_physic("gcmtochim","missing o2plus tracer",1)
4094            end if
4095            if (igcm_noplus == 0) then
4096               write(*,*) "gcmtochim: Error; no NO+ tracer !!!"
4097               call abort_physic("gcmtochim","missing noplus tracer",1)
4098            endif
4099            if (igcm_coplus == 0) then
4100               write(*,*) "gcmtochim: Error; no CO+ tracer !!!"
4101               call abort_physic("gcmtochim","missing coplus tracer",1)
4102            endif
4103            if (igcm_cplus == 0) then
4104               write(*,*) "gcmtochim: Error; no C+ tracer !!!"
4105               call abort_physic("gcmtochim","missing cplus tracer",1)
4106            endif
4107            if (igcm_n2plus == 0) then
4108               write(*,*) "gcmtochim: Error; no N2+ tracer !!!"
4109               call abort_physic("gcmtochim","missing n2plus tracer",1)
4110            endif
4111            if (igcm_nplus == 0) then
4112               write(*,*) "gcmtochim: Error; no N+ tracer !!!"
4113               call abort_physic("gcmtochim","missing nplus tracer",1)
4114            endif
4115            if (igcm_hplus == 0) then
4116               write(*,*) "gcmtochim: Error; no H+ tracer !!!"
4117               call abort_physic("gcmtochim","missing hplus tracer",1)
4118            endif
4119            if (igcm_hco2plus == 0) then
4120               write(*,*) "gcmtochim: Error; no HCO2+ tracer !!!"
4121               call abort_physic("gcmtochim","missing hco2plus tracer",1)
4122            endif
4123            if (igcm_hcoplus == 0) then
4124               write(*,*) "gcmtochim: Error; no HCO+ tracer !!!"
4125               call abort_physic("gcmtochim","missing hcoplus tracer",1)
4126            endif
4127            if (igcm_h2oplus == 0) then
4128               write(*,*) "gcmtochim: Error; no H2O+ tracer !!!"
4129               call abort_physic("gcmtochim","missing h2oplus tracer",1)
4130            endif
4131            if (igcm_h3oplus == 0) then
4132               write(*,*) "gcmtochim: Error; no H3O+ tracer !!!"
4133               call abort_physic("gcmtochim","missing h3oplus tracer",1)
4134            endif
4135            if (igcm_ohplus == 0) then
4136               write(*,*) "gcmtochim: Error; no OH+ tracer !!!"
4137               call abort_physic("gcmtochim","missing ohplus tracer",1)
4138            endif
4139            if (igcm_elec == 0) then
4140               write(*,*) "gcmtochim: Error; no e- tracer !!!"
4141               call abort_physic("gcmtochim","missing elec tracer",1)
4142            end if
4143         end if  ! ionchem
4144         firstcall = .false.
4145      end if ! of if (firstcall)
4146
4147!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4148!     initialise mixing ratios
4149!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4150
4151      do l = 1,nlayer
4152         rm(l,i_co2)  = zycol(l, igcm_co2)
4153         rm(l,i_co)   = zycol(l, igcm_co)
4154         rm(l,i_o)    = zycol(l, igcm_o)
4155         rm(l,i_o1d)  = zycol(l, igcm_o1d)
4156         rm(l,i_o2)   = zycol(l, igcm_o2)
4157         rm(l,i_o3)   = zycol(l, igcm_o3)
4158         rm(l,i_h)    = zycol(l, igcm_h)
4159         rm(l,i_h2)   = zycol(l, igcm_h2)
4160         rm(l,i_oh)   = zycol(l, igcm_oh)
4161         rm(l,i_ho2)  = zycol(l, igcm_ho2)
4162         rm(l,i_h2o2) = zycol(l, igcm_h2o2)
4163         rm(l,i_h2o)  = zycol(l, igcm_h2o_vap)
4164         rm(l,i_n)    = zycol(l, igcm_n)
4165         rm(l,i_n2d)  = zycol(l, igcm_n2d)
4166         rm(l,i_no)   = zycol(l, igcm_no)
4167         rm(l,i_no2)  = zycol(l, igcm_no2)
4168         rm(l,i_n2)   = zycol(l, igcm_n2)
4169      enddo
4170
4171      if (deutchem) then
4172         do l = 1,nlayer
4173            rm(l,i_hdo)  = zycol(l, igcm_hdo_vap)
4174            rm(l,i_od)   = zycol(l, igcm_od)
4175            rm(l,i_d)    = zycol(l, igcm_d)
4176            rm(l,i_hd)   = zycol(l, igcm_hd)
4177            rm(l,i_do2)  = zycol(l, igcm_do2)
4178            rm(l,i_hdo2)  = zycol(l, igcm_hdo2)
4179         end do
4180      endif
4181
4182      if (ionchem) then
4183         do l = 1,nlayer
4184            rm(l,i_co2plus)  = zycol(l, igcm_co2plus)
4185            rm(l,i_oplus)    = zycol(l, igcm_oplus)
4186            rm(l,i_o2plus)   = zycol(l, igcm_o2plus)
4187            rm(l,i_noplus)   = zycol(l, igcm_noplus)
4188            rm(l,i_coplus)   = zycol(l, igcm_coplus)
4189            rm(l,i_cplus)    = zycol(l, igcm_cplus)
4190            rm(l,i_n2plus)   = zycol(l, igcm_n2plus)
4191            rm(l,i_nplus)    = zycol(l, igcm_nplus)
4192            rm(l,i_hplus)    = zycol(l, igcm_hplus)
4193            rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
4194            rm(l,i_hcoplus)  = zycol(l, igcm_hcoplus)
4195            rm(l,i_h2oplus)  = zycol(l, igcm_h2oplus)
4196            rm(l,i_h3oplus)  = zycol(l, igcm_h3oplus)
4197            rm(l,i_ohplus)   = zycol(l, igcm_ohplus)
4198            rm(l,i_elec)     = zycol(l, igcm_elec)
4199         end do
4200      end if
4201
4202      where (rm(:,:) < 1.e-30)
4203         rm(:,:) = 0.
4204      end where
4205
4206!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4207!     initialise number densities
4208!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4209
4210      do iesp = 1,nesp
4211         do l = 1,nlayer
4212            c(l,iesp) = rm(l,iesp)*dens(l)
4213         end do
4214      end do
4215
4216      end subroutine gcmtochim
4217
4218!*****************************************************************
4219 
4220      subroutine chimtogcm(nlayer, ionchem, deutchem, nq, zycol,      &
4221                           lswitch, nesp,                             &
4222                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
4223                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
4224                           i_n, i_n2d, i_no, i_no2, i_n2,             &
4225                           i_co2plus, i_oplus, i_o2plus, i_noplus,    &
4226                           i_coplus, i_cplus, i_n2plus, i_nplus,      &
4227                           i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, &
4228                           i_h3oplus, i_ohplus, i_elec, i_hdo, i_od,  &
4229                           i_d, i_hd, i_do2, i_hdo2, dens, c)
4230 
4231!*****************************************************************
4232 
4233      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,          &
4234                            igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,   &
4235                            igcm_ho2, igcm_h2o2, igcm_h2o_vap,            &
4236                            igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2, &
4237                            igcm_co2plus, igcm_oplus, igcm_o2plus,        &
4238                            igcm_noplus, igcm_coplus, igcm_cplus,         &
4239                            igcm_n2plus, igcm_nplus, igcm_hplus,          &
4240                            igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
4241                            igcm_h3oplus, igcm_ohplus, igcm_elec,         &
4242                            igcm_hdo_vap, igcm_od, igcm_d, igcm_hd,       &
4243                            igcm_do2, igcm_hdo2
4244
4245      implicit none
4246
4247#include "callkeys.h"
4248
4249!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4250!     inputs:
4251!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4252 
4253      integer, intent(in) :: nlayer  ! number of atmospheric layers
4254      integer, intent(in) :: nq      ! number of tracers in the gcm
4255      logical, intent(in) :: ionchem
4256      logical, intent(in) :: deutchem
4257      integer :: nesp                ! number of species in the chemistry
4258      integer :: lswitch             ! interface level between chemistries
4259      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
4260                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
4261                 i_n, i_n2d, i_no, i_no2, i_n2,                  &
4262                 i_co2plus, i_oplus, i_o2plus, i_noplus,         &
4263                 i_coplus, i_cplus, i_n2plus, i_nplus,           &
4264                 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,      &
4265                 i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, i_d,  &
4266                 i_hd, i_do2, i_hdo2
4267
4268      real :: dens(nlayer)     ! total number density (molecule.cm-3)
4269      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3)
4270
4271!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4272!     output:
4273!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4274       
4275      real zycol(nlayer,nq)  ! volume mixing ratios in the gcm
4276
4277!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4278!     local:
4279!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4280
4281      integer l
4282     
4283!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4284!     save mixing ratios for the gcm
4285!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4286
4287      do l = 1,lswitch-1
4288         zycol(l, igcm_co2)     = c(l,i_co2)/dens(l)
4289         zycol(l, igcm_co)      = c(l,i_co)/dens(l)
4290         zycol(l, igcm_o)       = c(l,i_o)/dens(l)
4291         zycol(l, igcm_o1d)     = c(l,i_o1d)/dens(l)
4292         zycol(l, igcm_o2)      = c(l,i_o2)/dens(l)
4293         zycol(l, igcm_o3)      = c(l,i_o3)/dens(l)
4294         zycol(l, igcm_h)       = c(l,i_h)/dens(l) 
4295         zycol(l, igcm_h2)      = c(l,i_h2)/dens(l)
4296         zycol(l, igcm_oh)      = c(l,i_oh)/dens(l)
4297         zycol(l, igcm_ho2)     = c(l,i_ho2)/dens(l)
4298         zycol(l, igcm_h2o2)    = c(l,i_h2o2)/dens(l)
4299         zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
4300         zycol(l, igcm_n)       = c(l,i_n)/dens(l)
4301         zycol(l, igcm_n2d)     = c(l,i_n2d)/dens(l)
4302         zycol(l, igcm_no)      = c(l,i_no)/dens(l)
4303         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
4304         zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
4305      enddo
4306     
4307      if(deutchem) then
4308         do l=1,lswitch-1
4309            zycol(l, igcm_hdo_vap) = c(l,i_hdo)/dens(l)
4310            zycol(l, igcm_od)      = c(l,i_od)/dens(l)
4311            zycol(l, igcm_d)       = c(l,i_d)/dens(l)
4312            zycol(l, igcm_hd)      = c(l,i_hd)/dens(l)
4313            zycol(l, igcm_do2)     = c(l,i_do2)/dens(l)
4314            zycol(l, igcm_hdo2)    = c(l,i_hdo2)/dens(l)
4315         end do
4316      endif !deutchem
4317
4318      if (ionchem) then
4319         do l = 1,lswitch-1
4320            zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
4321            zycol(l, igcm_oplus)   = c(l,i_oplus)/dens(l)
4322            zycol(l, igcm_o2plus)  = c(l,i_o2plus)/dens(l)
4323            zycol(l, igcm_noplus)  = c(l,i_noplus)/dens(l)
4324            zycol(l, igcm_coplus)  = c(l,i_coplus)/dens(l)
4325            zycol(l, igcm_cplus)   = c(l,i_cplus)/dens(l)
4326            zycol(l, igcm_n2plus)  = c(l,i_n2plus)/dens(l)
4327            zycol(l, igcm_nplus)   = c(l,i_nplus)/dens(l)
4328            zycol(l, igcm_hplus)   = c(l,i_hplus)/dens(l)
4329            zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
4330            zycol(l, igcm_hcoplus) = c(l,i_hcoplus)/dens(l)
4331            zycol(l, igcm_h2oplus) = c(l,i_h2oplus)/dens(l)
4332            zycol(l, igcm_h3oplus) = c(l,i_h3oplus)/dens(l)
4333            zycol(l, igcm_ohplus)  = c(l,i_ohplus)/dens(l)
4334            zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
4335         end do
4336      end if
4337
4338      end subroutine chimtogcm
4339
4340end subroutine photochemistry
Note: See TracBrowser for help on using the repository browser.