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

Last change on this file since 2589 was 2554, checked in by emillour, 4 years ago

Mars GCM:
Make NO and O2 Nightglow emissions computations more flexible (reaction index

no longer hard coded in photochemistry.F90).

FGG

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