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

Last change on this file was 3726, checked in by emillour, 3 months ago

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