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

Last change on this file since 2324 was 2321, checked in by emillour, 6 years ago

Mars GCM:

  • Finalized water ion chemisty (added H3O+ and OH+ ions). Added an example of relevent traceur.def file in deftank.
  • Added the computation of NO nightglow.

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