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

Last change on this file was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

  • remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
  • turn aeronomars routines to modules, for those which aren't in modules yet.

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