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

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

Mars PCM:
More tidying in aeronomars:

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

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