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

Last change on this file since 2814 was 2615, checked in by romain.vande, 3 years ago

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