source: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F @ 1747

Last change on this file since 1747 was 1660, checked in by emillour, 8 years ago

Mars GCM:

  • Added possibility to run with an Helium "he" tracer (to be initialized at constant value of 3.6e-7 kg/kg_air, i.e. the 4ppm of Krasnopolsky 1996 EUVE satellite, using newstart).
  • corrected case for CH4 in aeronomars/photochemistry.F (was using undefined indexes when there is no CH4).
  • updated aki/cpi coefficients for Argon used to compute mean atmospheric Cp and thermal conductivity in aeronomars/concentrations.F

JYC+EM
File size: 50.2 KB
Line 
1c*****************************************************************
2c
3c     Photochemical routine
4c
5c     Author: Franck Lefevre
6c     ------
7c
8c     Version: 11/12/2013
9c
10c*****************************************************************
11
12      subroutine photochemistry(nlayer, nq,
13     $                  lswitch, zycol, sza, ptimestep,
14     $                  press,temp, dens, dist_sol, surfdust1d,
15     $                  surfice1d, jo3, tau)
16
17      implicit none
18
19#include "chimiedata.h"
20#include "callkeys.h"
21
22cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
23c     inputs:
24cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
25
26      integer, intent(in) :: nlayer ! number of atmospheric layers
27      integer, intent(in) :: nq     ! number of tracers
28     
29      real :: sza                   ! solar zenith angle (deg)
30      real :: ptimestep             ! physics timestep (s)
31      real :: press(nlayer)         ! pressure (hpa)
32      real :: temp(nlayer)          ! temperature (k)
33      real :: dens(nlayer)          ! density (cm-3)
34      real :: dist_sol              ! sun distance (au)
35      real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
36      real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
37      real :: tau                   ! optical depth at 7 hpa
38
39cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
40c     input/output:
41cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
42     
43      real :: zycol(nlayer,nq)      ! chemical species volume mixing ratio
44
45cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
46c     output:
47cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
48     
49      real :: jo3(nlayer)  ! photodissociation rate o3 -> o1d
50
51cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
52c     local:
53cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
54
55      integer, parameter :: nesp = 16 ! number of species in the chemistry
56
57      integer :: phychemrat       ! ratio physics timestep/chemistry timestep
58      integer :: istep
59      integer :: i_co2,i_o3,j_o3_o1d,lev
60      integer :: lswitch
61
62      real :: stimestep           ! standard timestep for the chemistry (s)
63      real :: ctimestep           ! real timestep for the chemistry (s)
64      real :: rm(nlayer,nesp)     ! species volume mixing ratio
65      real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
66      real :: rmco2(nlayer)       ! co2 mixing ratio
67      real :: rmo3(nlayer)        ! ozone mixing ratio
68
69c     reaction rates
70
71      real :: a001(nlayer), a002(nlayer), a003(nlayer)
72      real :: b001(nlayer), b002(nlayer), b003(nlayer),
73     $        b004(nlayer), b005(nlayer), b006(nlayer),
74     $        b007(nlayer), b008(nlayer), b009(nlayer)
75      real :: c001(nlayer), c002(nlayer), c003(nlayer),
76     $        c004(nlayer), c005(nlayer), c006(nlayer),
77     $        c007(nlayer), c008(nlayer), c009(nlayer),
78     $        c010(nlayer), c011(nlayer), c012(nlayer),
79     $        c013(nlayer), c014(nlayer), c015(nlayer),
80     $        c016(nlayer), c017(nlayer), c018(nlayer)
81      real :: d001(nlayer), d002(nlayer), d003(nlayer)
82      real :: e001(nlayer), e002(nlayer), e003(nlayer)
83      real :: h001(nlayer), h002(nlayer), h003(nlayer),
84     $        h004(nlayer), h005(nlayer)
85      real :: t001(nlayer), t002(nlayer), t003(nlayer)
86
87cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
88c     stimestep  : standard timestep for the chemistry (s)       c
89c     ctimestep  : real timestep for the chemistry (s)           c
90c     phychemrat : ratio physical/chemical timestep              c
91cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
92
93      stimestep = 600. ! standard value : 10 mn
94
95      phychemrat = nint(ptimestep/stimestep)
96
97      ctimestep = ptimestep/real(phychemrat)
98
99cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
100c     initialisation of chemical species and families            c
101cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
102
103      call gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
104
105cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
106c     compute reaction rates                                     c
107cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
108
109      call chemrates(nlayer,
110     $               lswitch, dens, press, temp,
111     $               surfdust1d, surfice1d,
112     $               a001, a002, a003,
113     $               b001, b002, b003, b004, b005, b006,
114     $               b007, b008, b009,
115     $               c001, c002, c003, c004, c005, c006,
116     $               c007, c008, c009, c010, c011, c012,
117     $               c013, c014, c015, c016, c017, c018,
118     $               d001, d002, d003,
119     $               e001, e002, e003,
120     $               h001, h002, h003, h004, h005,
121     $               t001, t002, t003, tau)
122
123cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
124c     interpolation of photolysis rates in the lookup table      c
125cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
126
127      i_co2 = 1
128      i_o3  = 6
129
130      do lev = 1,lswitch-1
131         rmco2(lev) = rm(lev,i_co2)
132         rmo3(lev)  = rm(lev, i_o3)
133      end do
134
135      call photolysis(nlayer, lswitch, press, temp, sza, tau, dist_sol,
136     $                rmco2, rmo3, j)
137
138      j_o3_o1d = 5
139
140      do lev = 1,lswitch-1
141         jo3(lev) = j(lev,j_o3_o1d)
142      end do
143
144cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
145c     loop over time within the physical timestep                c
146cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
147
148      do istep = 1,phychemrat
149         call chimie(nlayer, lswitch, nesp, rm, j, dens, ctimestep,
150     $               press, temp, sza, dist_sol,
151     $               a001, a002, a003,
152     $               b001, b002, b003, b004, b005, b006,
153     $               b007, b008, b009,
154     $               c001, c002, c003, c004, c005, c006,
155     $               c007, c008, c009, c010, c011, c012,
156     $               c013, c014, c015, c016, c017, c018,
157     $               d001, d002, d003,
158     $               e001, e002, e003,
159     $               h001, h002, h003, h004, h005,
160     $               t001, t002, t003)
161      end do
162
163cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
164c     save chemical species for the gcm                          c
165cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
166
167      call chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
168
169      contains
170
171
172c*****************************************************************
173
174      subroutine chimie(nlayer,
175     $                  lswitch, nesp, rm, j, dens, dt,
176     $                  press, t, sza, dist_sol,
177     $                  a001, a002, a003,
178     $                  b001, b002, b003, b004, b005, b006,
179     $                  b007, b008, b009,
180     $                  c001, c002, c003, c004, c005, c006,
181     $                  c007, c008, c009, c010, c011, c012,
182     $                  c013, c014, c015, c016, c017, c018,
183     $                  d001, d002, d003,
184     $                  e001, e002, e003,
185     $                  h001, h002, h003, h004, h005,
186     $                  t001, t002, t003)
187
188c*****************************************************************
189
190      use tracer_mod, only: igcm_ch4
191      implicit none
192
193#include "chimiedata.h"
194#include "callkeys.h"
195
196cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
197c     inputs:
198cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
199
200      integer, intent(in) :: nlayer ! number of atmospheric layers
201      integer :: lswitch            ! interface level between chemistries
202      integer :: nesp               ! number of species
203
204      real :: dens(nlayer)          ! density (cm-3)
205      real :: dt                    ! chemistry timestep (s)
206      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1)
207      real :: press(nlayer)         ! pressure (hpa)
208      real :: t(nlayer)             ! temperature (k)
209      real :: sza                   ! solar zenith angle (deg)
210      real :: dist_sol              ! sun distance (au)
211
212c     reaction rates
213
214      real :: a001(nlayer), a002(nlayer), a003(nlayer)
215      real :: b001(nlayer), b002(nlayer), b003(nlayer),
216     $        b004(nlayer), b005(nlayer), b006(nlayer),
217     $        b007(nlayer), b008(nlayer), b009(nlayer)
218      real :: c001(nlayer), c002(nlayer), c003(nlayer),
219     $        c004(nlayer), c005(nlayer), c006(nlayer),
220     $        c007(nlayer), c008(nlayer), c009(nlayer),
221     $        c010(nlayer), c011(nlayer), c012(nlayer),
222     $        c013(nlayer), c014(nlayer), c015(nlayer),
223     $        c016(nlayer), c017(nlayer), c018(nlayer)
224      real :: d001(nlayer), d002(nlayer), d003(nlayer)
225      real :: e001(nlayer), e002(nlayer), e003(nlayer)
226      real :: h001(nlayer), h002(nlayer), h003(nlayer),
227     $        h004(nlayer), h005(nlayer)
228      real :: t001(nlayer), t002(nlayer), t003(nlayer)
229
230cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
231c     input/output:
232cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
233
234      real :: rm(nlayer,nesp)  ! volume mixing ratios
235
236cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
237c     local:
238cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
239
240      real, parameter :: hetero_ice  = 1. ! switch for het. chem. on ice clouds
241      real, parameter :: hetero_dust = 0. ! switch for het. chem. on dust
242
243      integer :: iesp, iter, l
244      integer :: i_co2, i_co, i_o2, i_h2, i_h2o, i_h2o2, i_hox, i_ox,
245     $           i_o1d, i_o, i_o3, i_h, i_oh, i_ho2, i_ch4, i_n2
246      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
247     $           j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
248     $           j_ho2, j_no, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
249     $           j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
250     $           j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
251     $           j_hocl, j_clo, j_so2, j_so, j_h2s, j_so3,
252     $           j_hno3, j_hno4
253
254      integer, parameter :: niter = 5        ! iterations in the chemical scheme
255
256      real :: cc0(nlayer,nesp)        ! initial number density (cm-3)
257      real :: cc(nlayer,nesp)         ! final number density (cm-3)
258      real :: nox(nlayer)             ! nox number density (cm-3)
259      real :: no(nlayer)              ! no  number density (cm-3)
260      real :: no2(nlayer)             ! no2 number density (cm-3)
261      real :: production(nlayer,nesp) ! production rate
262      real :: loss(nlayer,nesp)       ! loss rate
263
264      real :: ro_o3, rh_ho2, roh_ho2, rno2_no
265
266cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
267c     tracer numbering in the chemistry
268cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
269
270      i_co2    =  1
271      i_co     =  2
272      i_o      =  3
273      i_o1d    =  4
274      i_o2     =  5
275      i_o3     =  6
276      i_h      =  7
277      i_h2     =  8
278      i_oh     =  9
279      i_ho2    = 10
280      i_h2o2   = 11
281      i_ch4    = 12
282      i_h2o    = 13
283      i_n2     = 14
284      i_hox    = 15
285      i_ox     = 16
286
287cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
288c     numbering of photolysis rates
289cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
290
291      j_o2_o         =  1      ! o2 + hv     -> o + o
292      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
293      j_co2_o        =  3      ! co2 + hv    -> co + o
294      j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
295      j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
296      j_o3_o         =  6      ! o3 + hv     -> o2 + o
297      j_h2o          =  7      ! h2o + hv    -> h + oh
298      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
299      j_ho2          =  9      ! ho2 + hv    -> oh + o
300      j_no           =  10     ! no + hv     -> n + o
301      j_no2          =  11     ! no2 + hv    -> no + o
302      j_hno3         =  12     ! hno3 + hv   -> oh + no2
303      j_hno4         =  13     ! hno4 + hv   -> ho2 + no2
304
305cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
306c     volume mixing ratio -> density conversion
307cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
308
309      do iesp = 1,nesp
310         do l = 1,lswitch-1
311            cc0(l,iesp) = rm(l,iesp)*dens(l)
312            cc(l,iesp)  = cc0(l,iesp)
313         end do
314      end do
315
316cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
317c     co2 and nox number densities (cm-3)   
318c     
319c     nox mixing ratio: 6.e-10 (yung and demore, 1999)
320cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
321
322      do l = 1,lswitch-1
323         nox(l) = 6.e-10*dens(l)
324      end do
325
326cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
327c     loop over iterations in the chemical scheme
328cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
329
330      do iter = 1,niter
331
332cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
333c        nox species
334cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
335c        no2/no
336cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
337
338      do l = 1,lswitch-1
339
340         rno2_no = (d002(l)*cc(l,i_o3) + d003(l)*cc(l,i_ho2))
341     $            /(j(l,j_no2) +
342     $              d001(l)*max(cc(l,i_o),1.e-30*dens(l)))
343
344cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
345c        no
346cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
347
348         no(l) = nox(l)/(1. + rno2_no)
349
350cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
351c        no2
352cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
353
354         no2(l) = no(l)*rno2_no
355
356      end do
357
358cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
359c        hox species
360cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
361c        photochemical equilibrium for oh and ho2
362cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
363c        h/ho2
364cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
365
366      do l = 1,lswitch-1
367
368         rh_ho2 = (c001(l)*cc(l,i_o)
369     $           + c004(l)*cc(l,i_h)
370     $           + c005(l)*cc(l,i_h)
371     $           + c006(l)*cc(l,i_h)
372     $           + c007(l)*cc(l,i_oh)
373     $           + 2.*c008(l)*cc(l,i_ho2)
374     $           + c015(l)*cc(l,i_o3)
375     $           + 2.*c016(l)*cc(l,i_ho2)
376     $           + d003(l)*no(l)             ! ajout 20090401
377     $           + j(l,j_ho2)
378     $           + h001(l)*hetero_ice
379     $           + h003(l)*hetero_dust)
380     $          /( c011(l)*cc(l,i_o2)
381     $           + t001(l)*cc(l,i_h2o)
382     $            /max(cc(l,i_h),dens(l)*1.e-30)   ! ajout 20090401
383     $           )
384
385cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
386c        oh/ho2
387cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
388
389         roh_ho2 = (c001(l)*cc(l,i_o)
390     $            + c003(l)*cc(l,i_o3)*rh_ho2
391     $            + 2.*c004(l)*cc(l,i_h)
392     $            + 2.*c008(l)*cc(l,i_ho2)
393     $            + c015(l)*cc(l,i_o3)
394     $            + d003(l)*no(l)
395     $            + j(l,j_ho2)
396     $            + 2.*b002(l)*cc(l,i_o1d)*cc(l,i_h2o)         ! ajout 20101210
397     $             /max(cc(l,i_ho2),dens(l)*1.e-30)            ! ajout 20101210
398     $            + b003(l)*cc(l,i_o1d)*cc(l,i_h2)             ! ajout 20101210
399     $             /max(cc(l,i_ho2),dens(l)*1.e-30)            ! ajout 20101210
400     $            + j(l,j_h2o)*cc(l,i_h2o)
401     $             /max(cc(l,i_ho2),dens(l)*1.e-30)
402     $            + t001(l)*cc(l,i_h2o)               ! suppression 20090401
403     $             /max(cc(l,i_ho2),dens(l)*1.e-30)   ! suppression 20090401
404     $            )
405     $            /(c002(l)*cc(l,i_o)
406     $            + c007(l)*cc(l,i_ho2)
407     $            + c009(l)*cc(l,i_h2o2)         ! ajout 20090401
408     $            + 2.*c013(l)*cc(l,i_oh)        ! ajout 20090401
409     $            + 2.*c017(l)*cc(l,i_oh)        ! ajout 20090401
410     $            + e001(l)*cc(l,i_co)
411     $            + h002(l)*hetero_ice)
412
413cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
414c        h
415cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
416
417         cc(l,i_h) = cc(l,i_hox)
418     $                 /(1. + (1. + roh_ho2)/rh_ho2)
419
420cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
421c        ho2
422cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
423
424         cc(l,i_ho2) = cc(l,i_h)/rh_ho2
425
426cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
427c        oh
428cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
429
430         cc(l,i_oh) = cc(l,i_ho2)*roh_ho2
431
432      end do
433
434cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
435c        ox species
436cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
437c        day:
438c        - o1d at photochemical equilibrium
439c        - o3 at photochemical equilibrium
440c        - continuity equation for ox
441c        night:
442c        - o1d = 0
443c        - continuity equation for o3
444c        - continuity equation for o
445cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
446       
447      if (sza .le. 95.) then
448
449         do l = 1,lswitch-1
450
451cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
452c        o(1d)
453cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
454
455            cc(l,i_o1d) = (j(l,j_co2_o1d)*cc(l,i_co2)
456     $                   + j(l,j_o2_o1d)*cc(l,i_o2)
457     $                   + j(l,j_o3_o1d)*cc(l,i_o3))
458     $                   /(b001(l)*cc(l,i_co2)
459     $                   + b002(l)*cc(l,i_h2o)
460     $                   + b003(l)*cc(l,i_h2)
461     $                   + b004(l)*cc(l,i_o2)
462     $                   + b005(l)*cc(l,i_o3)
463     $                   + b006(l)*cc(l,i_o3))
464
465cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
466c        o/o3
467cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
468
469            ro_o3 = (j(l,j_o3_o1d) + j(l,j_o3_o)
470     $              + a003(l)*cc(l,i_o)
471     $              + c003(l)*cc(l,i_h)
472     $              + c014(l)*cc(l,i_oh)
473     $              + c015(l)*cc(l,i_ho2)
474     $              )
475     $              /(a001(l)*cc(l,i_o2))
476
477cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
478c        o3
479cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
480
481            cc(l,i_o3) = cc(l,i_ox)/(1. + ro_o3)
482
483cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
484c        o
485cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
486
487            cc(l,i_o) = cc(l,i_o3)*ro_o3
488
489cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
490c        ox = o + o3
491cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
492
493            production(l,i_ox) =
494     $                   + j(l,j_co2_o)*cc(l,i_co2)
495     $                   + j(l,j_co2_o1d)*cc(l,i_co2)
496     $                   + j(l,j_ho2)*cc(l,i_ho2)
497     $                   + 2.*j(l,j_o2_o)*cc(l,i_o2)
498     $                   + 2.*j(l,j_o2_o1d)*cc(l,i_o2)
499     $                   + c006(l)*cc(l,i_h)*cc(l,i_ho2)
500     $                   + c013(l)*cc(l,i_oh)*cc(l,i_oh)
501     $                   + d003(l)*cc(l,i_ho2)*no(l)
502
503            loss(l,i_ox) = 2.*a002(l)*cc(l,i_o)*cc(l,i_o)
504     $                   + 2.*a003(l)*cc(l,i_o)*cc(l,i_o3)
505     $                   + c001(l)*cc(l,i_ho2)*cc(l,i_o)
506     $                   + c002(l)*cc(l,i_oh)*cc(l,i_o)
507     $                   + c003(l)*cc(l,i_h)*cc(l,i_o3)
508     $                   + c012(l)*cc(l,i_o)*cc(l,i_h2o2)
509     $                   + c014(l)*cc(l,i_o3)*cc(l,i_oh)
510     $                   + c015(l)*cc(l,i_o3)*cc(l,i_ho2)
511     $                   + d001(l)*cc(l,i_o)*no2(l)
512     $                   + e002(l)*cc(l,i_o)*cc(l,i_co)
513
514            loss(l,i_ox) = loss(l,i_ox)/cc(l,i_ox)
515
516         end do
517
518      else
519
520         do l = 1,lswitch-1
521
522cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
523c        o(1d)
524cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
525
526            cc(l,i_o1d) = 0.
527
528cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
529c        o3
530cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
531
532            production(l,i_o3) = a001(l)*cc(l,i_o2)*cc(l,i_o)
533
534            loss(l,i_o3) = a003(l)*cc(l,i_o)
535     $                   + c003(l)*cc(l,i_h)
536     $                   + c014(l)*cc(l,i_oh)
537     $                   + c015(l)*cc(l,i_ho2)
538
539cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
540c        o
541cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
542
543            production(l,i_o) = c006(l)*cc(l,i_h)*cc(l,i_ho2)
544     $                        + c013(l)*cc(l,i_oh)*cc(l,i_oh)
545
546            loss(l,i_o)  =  a001(l)*cc(l,i_o2)
547     $                   + 2.*a002(l)*cc(l,i_o)
548     $                   + a003(l)*cc(l,i_o3)
549     $                   + c001(l)*cc(l,i_ho2)
550     $                   + c002(l)*cc(l,i_oh)
551     $                   + c012(l)*cc(l,i_h2o2)
552     $                   + e002(l)*cc(l,i_co)
553
554         end do
555      end if
556
557cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
558c     other species
559cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
560
561      do l = 1,lswitch-1
562
563cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
564c        co2
565cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
566
567         production(l,i_co2) = e001(l)*cc(l,i_oh)*cc(l,i_co)
568     $                       + e002(l)*cc(l,i_o)*cc(l,i_co)
569     $                       + t002(l)*cc(l,i_ch4)*16./44. ! ajout 20090401
570     $                       + t003(l)*cc(l,i_co2)*8./44.  ! ajout 20090401
571
572         loss(l,i_co2) = j(l,j_co2_o)
573     $                 + j(l,j_co2_o1d)
574     $                 + t003(l)
575
576cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
577c        co
578cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
579
580         production(l,i_co) = j(l,j_co2_o)*cc(l,i_co2)
581     $                      + j(l,j_co2_o1d)*cc(l,i_co2)
582     $                      + t003(l)*cc(l,i_co2)
583
584         loss(l,i_co) = e001(l)*cc(l,i_oh)
585     $                + e002(l)*cc(l,i_o)
586
587cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
588c        o2
589cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
590
591         production(l,i_o2) =
592     $                  j(l,j_o3_o)*cc(l,i_o3)
593     $                + j(l,j_o3_o1d)*cc(l,i_o3)
594     $                + a002(l)*cc(l,i_o)*cc(l,i_o)
595     $                + 2.*a003(l)*cc(l,i_o)*cc(l,i_o3)
596     $                + 2.*b005(l)*cc(l,i_o1d)*cc(l,i_o3)
597     $                + b006(l)*cc(l,i_o1d)*cc(l,i_o3)
598     $                + c001(l)*cc(l,i_o)*cc(l,i_ho2)
599     $                + c002(l)*cc(l,i_o)*cc(l,i_oh)
600     $                + c003(l)*cc(l,i_h)*cc(l,i_o3)
601     $                + c005(l)*cc(l,i_h)*cc(l,i_ho2)
602     $                + c007(l)*cc(l,i_oh)*cc(l,i_ho2)
603     $                + c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
604     $                + c014(l)*cc(l,i_o3)*cc(l,i_oh)
605     $                + 2.*c015(l)*cc(l,i_o3)*cc(l,i_ho2)
606     $                + c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
607     $                + d001(l)*cc(l,i_o)*no2(l)
608
609         loss(l,i_o2) = j(l,j_o2_o)
610     $                + j(l,j_o2_o1d)
611     $                + a001(l)*cc(l,i_o)
612     $                + c011(l)*cc(l,i_h)
613
614cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
615c        h2
616cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
617
618         production(l,i_h2) = c005(l)*cc(l,i_h)*cc(l,i_ho2)
619     $                      + c018(l)*cc(l,i_h)*cc(l,i_h)
620
621         loss(l,i_h2) = b003(l)*cc(l,i_o1d)
622     $                + c010(l)*cc(l,i_oh)
623
624cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
625c        h2o
626cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
627
628         production(l,i_h2o) =
629     $                   c006(l)*cc(l,i_h)*cc(l,i_ho2)
630     $                 + c007(l)*cc(l,i_oh)*cc(l,i_ho2)
631     $                 + c009(l)*cc(l,i_oh)*cc(l,i_h2o2)
632     $                 + c010(l)*cc(l,i_oh)*cc(l,i_h2)
633     $                 + c013(l)*cc(l,i_oh)*cc(l,i_oh)
634     $                 + h004(l)*cc(l,i_h2o2)*hetero_ice
635
636         loss(l,i_h2o) = j(l,j_h2o)
637     $                 + b002(l)*cc(l,i_o1d)
638     $                 + t001(l)
639
640cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
641c        h2o2
642cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
643
644         production(l,i_h2o2) =
645     $                    c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
646     $                  + c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
647     $                  + c017(l)*cc(l,i_oh)*cc(l,i_oh)
648c    $                  + 0.5*h001(l)*cc(l,i_ho2)*hetero_ice
649c    $                  + 0.5*h002(l)*cc(l,i_oh)*hetero_ice
650
651         loss(l,i_h2o2) = j(l,j_h2o2)
652     $                  + c009(l)*cc(l,i_oh)
653     $                  + c012(l)*cc(l,i_o)
654     $                  + h004(l)*hetero_ice
655     $                  + h005(l)*hetero_dust
656
657cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
658c        hox = h + oh + ho2
659cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
660
661         production(l,i_hox) =
662     $                   2.*j(l,j_h2o)*cc(l,i_h2o)
663     $                 + 2.*j(l,j_h2o2)*cc(l,i_h2o2)
664     $                 + 2.*b002(l)*cc(l,i_o1d)*cc(l,i_h2o)
665     $                 + 2.*b003(l)*cc(l,i_o1d)*cc(l,i_h2)
666     $                 + 2.*c012(l)*cc(l,i_o)*cc(l,i_h2o2)
667     $                 + 2.*t001(l)*cc(l,i_h2o)
668
669         loss(l,i_hox) = 2.*c005(l)*cc(l,i_h)*cc(l,i_ho2)
670     $                 + 2.*c006(l)*cc(l,i_h)*cc(l,i_ho2)
671     $                 + 2.*c007(l)*cc(l,i_oh)*cc(l,i_ho2)
672     $                 + 2.*c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
673     $                 + 2.*c013(l)*cc(l,i_oh)*cc(l,i_oh)
674     $                 + 2.*c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
675     $                 + 2.*c017(l)*cc(l,i_oh)*cc(l,i_oh)
676     $                 + 2.*c018(l)*cc(l,i_h)*cc(l,i_h)
677     $                 + h001(l)*cc(l,i_ho2)*hetero_ice
678     $                 + h002(l)*cc(l,i_oh)*hetero_ice
679     $                 + h003(l)*cc(l,i_ho2)*hetero_dust
680
681         loss(l,i_hox) = loss(l,i_hox)/cc(l,i_hox)
682
683cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
684c        ch4
685cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
686
687         production(l,i_ch4) = 0.
688
689         if (igcm_ch4/=0) then
690           loss(l,i_ch4) = j(l,j_ch4_ch3_h)
691     $                   + j(l,j_ch4_1ch2_h2)
692     $                   + j(l,j_ch4_3ch2_h_h)
693     $                   + j(l,j_ch4_ch_h2_h)
694     $                   + b007(l)*cc(l,i_o1d)
695     $                   + b008(l)*cc(l,i_o1d)
696     $                   + b009(l)*cc(l,i_o1d)
697     $                   + e003(l)*cc(l,i_oh)
698         else
699           loss(l,i_ch4) = 0.
700         endif
701
702      end do
703
704cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
705c     update number densities
706cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
707
708c     long-lived species
709
710      do l = 1,lswitch-1
711         cc(l,i_co2) = (cc0(l,i_co2) + production(l,i_co2)*dt)
712     $              /(1. + loss(l,i_co2)*dt)
713         cc(l,i_co)  = (cc0(l,i_co) + production(l,i_co)*dt)
714     $              /(1. + loss(l,i_co)*dt)
715         cc(l,i_o2)  = (cc0(l,i_o2) + production(l,i_o2)*dt)
716     $              /(1. + loss(l,i_o2)*dt)
717         cc(l,i_h2)  = (cc0(l,i_h2) + production(l,i_h2)*dt)
718     $              /(1. + loss(l,i_h2)*dt)
719         cc(l,i_h2o2)= (cc0(l,i_h2o2) + production(l,i_h2o2)*dt)
720     $              /(1. + loss(l,i_h2o2)*dt)
721         cc(l,i_h2o) = (cc0(l,i_h2o) + production(l,i_h2o)*dt)
722     $              /(1. + loss(l,i_h2o)*dt)
723         cc(l,i_hox) = (cc0(l,i_hox) + production(l,i_hox)*dt)
724     $              /(1. + loss(l,i_hox)*dt)
725         cc(l,i_ch4) = (cc0(l,i_ch4) + production(l,i_ch4)*dt)
726     $              /(1. + loss(l,i_ch4)*dt)
727      end do
728
729c     ox species
730
731      if (sza .le. 95.) then
732         do l = 1,lswitch-1
733            cc(l,i_ox) = (cc0(l,i_ox) + production(l,i_ox)*dt)
734     $                   /(1. + loss(l,i_ox)*dt)
735         end do
736      else
737         do l = 1,lswitch-1
738            cc(l,i_o)  = (cc0(l,i_o) + production(l,i_o)*dt)
739     $                   /(1. + loss(l,i_o)*dt)
740            cc(l,i_o3) = (cc0(l,i_o3) + production(l,i_o3)*dt)
741     $                   /(1. + loss(l,i_o3)*dt)
742            cc(l,i_ox) = cc(l,i_o) + cc(l,i_o3)
743         end do
744      end if
745
746cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
747c     end of loop over iterations
748cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
749
750      end do
751
752cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
753c     density -> volume mixing ratio conversion
754cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
755
756      do iesp = 1,nesp
757         do l = 1,lswitch-1
758            rm(l,iesp) = max(cc(l,iesp)/dens(l), 1.e-30)
759         end do
760      end do
761
762      end subroutine chimie
763
764c*****************************************************************
765
766      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
767
768c*****************************************************************
769
770      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
771     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
772     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
773     &                      igcm_ch4
774
775      implicit none
776
777#include "callkeys.h"
778
779cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
780c     inputs:
781cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
782     
783      integer, intent(in) :: nlayer ! number of atmospheric layers
784      integer, intent(in) :: nq     ! number of tracers
785      real :: zycol(nlayer,nq)      ! species volume mixing ratio in the gcm
786
787      integer :: nesp               ! number of species in the chemistry
788      integer :: lswitch            ! interface level between chemistries
789
790cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
791c     outputs:
792cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
793
794      real :: rm(nlayer,nesp)   ! species volume mixing ratio
795     
796cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
797c     local:
798cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
799
800      integer :: l, iq
801     
802c     tracer indexes in the chemistry:
803
804      integer,parameter :: i_co2  =  1
805      integer,parameter :: i_co   =  2
806      integer,parameter :: i_o    =  3
807      integer,parameter :: i_o1d  =  4
808      integer,parameter :: i_o2   =  5
809      integer,parameter :: i_o3   =  6
810      integer,parameter :: i_h    =  7
811      integer,parameter :: i_h2   =  8
812      integer,parameter :: i_oh   =  9
813      integer,parameter :: i_ho2  = 10
814      integer,parameter :: i_h2o2 = 11
815      integer,parameter :: i_ch4  = 12
816      integer,parameter :: i_h2o  = 13
817      integer,parameter :: i_n2   = 14
818      integer,parameter :: i_hox  = 15
819      integer,parameter :: i_ox   = 16
820
821cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
822c     initialise chemical species
823cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
824
825      do l = 1,lswitch-1
826         rm(l,i_co2)  = max(zycol(l, igcm_co2),     1.e-30)
827         rm(l,i_co)   = max(zycol(l, igcm_co),      1.e-30)
828         rm(l,i_o)    = max(zycol(l, igcm_o),       1.e-30)
829         rm(l,i_o1d)  = max(zycol(l, igcm_o1d),     1.e-30)
830         rm(l,i_o2)   = max(zycol(l, igcm_o2),      1.e-30)
831         rm(l,i_o3)   = max(zycol(l, igcm_o3),      1.e-30)
832         rm(l,i_h)    = max(zycol(l, igcm_h),       1.e-30)
833         rm(l,i_h2)   = max(zycol(l, igcm_h2),      1.e-30)
834         rm(l,i_oh)   = max(zycol(l, igcm_oh),      1.e-30)
835         rm(l,i_ho2)  = max(zycol(l, igcm_ho2),     1.e-30)
836         rm(l,i_h2o2) = max(zycol(l, igcm_h2o2),    1.e-30)
837         rm(l,i_n2)   = max(zycol(l, igcm_n2),      1.e-30)
838         rm(l,i_h2o)  = max(zycol(l, igcm_h2o_vap), 1.e-30)
839      end do
840
841      if (igcm_ch4 .eq. 0) then
842         do l = 1,lswitch-1
843            rm(l,i_ch4) = 0.
844         end do
845      else
846         do l = 1,lswitch-1
847            rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30)
848         end do
849      end if
850
851cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
852c     initialise chemical families                     c
853cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
854
855      do l = 1,lswitch-1
856         rm(l,i_hox) = rm(l,i_h)
857     $               + rm(l,i_oh)
858     $               + rm(l,i_ho2)
859         rm(l,i_ox)  = rm(l,i_o)
860     $               + rm(l,i_o3)
861      end do
862
863      end subroutine gcmtochim
864
865c*****************************************************************
866c
867      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
868c
869c*****************************************************************
870c
871      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
872     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
873     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
874     &                      igcm_ch4
875
876      implicit none
877
878#include "callkeys.h"
879
880cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
881c     inputs:
882cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
883
884      integer, intent(in) :: nlayer ! number of atmospheric layers
885      integer, intent(in) :: nq     ! number of tracers
886      integer :: nesp               ! number of species in the chemistry
887      integer :: lswitch            ! interface level between chemistries
888     
889      real :: rm(nlayer,nesp)       ! species volume mixing ratio
890
891cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
892c     output:
893cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
894
895      real :: zycol(nlayer,nq)  ! species volume mixing ratio in the gcm
896
897cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
898c     local:
899cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
900
901      integer l, iq
902     
903c     tracer indexes in the chemistry:
904
905      integer,parameter :: i_co2  =  1
906      integer,parameter :: i_co   =  2
907      integer,parameter :: i_o    =  3
908      integer,parameter :: i_o1d  =  4
909      integer,parameter :: i_o2   =  5
910      integer,parameter :: i_o3   =  6
911      integer,parameter :: i_h    =  7
912      integer,parameter :: i_h2   =  8
913      integer,parameter :: i_oh   =  9
914      integer,parameter :: i_ho2  = 10
915      integer,parameter :: i_h2o2 = 11
916      integer,parameter :: i_ch4  = 12
917      integer,parameter :: i_h2o  = 13
918      integer,parameter :: i_n2   = 14
919      integer,parameter :: i_hox  = 15
920      integer,parameter :: i_ox   = 16
921
922cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
923c     save mixing ratios for the gcm
924cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
925
926      do l = 1,lswitch-1
927         zycol(l, igcm_co2)     = rm(l,i_co2)
928         zycol(l, igcm_co)      = rm(l,i_co)
929         zycol(l, igcm_o)       = rm(l,i_o)
930         zycol(l, igcm_o1d)     = rm(l,i_o1d)
931         zycol(l, igcm_o2)      = rm(l,i_o2)
932         zycol(l, igcm_o3)      = rm(l,i_o3)
933         zycol(l, igcm_h)       = rm(l,i_h) 
934         zycol(l, igcm_h2)      = rm(l,i_h2)
935         zycol(l, igcm_oh)      = rm(l,i_oh)
936         zycol(l, igcm_ho2)     = rm(l,i_ho2)
937         zycol(l, igcm_h2o2)    = rm(l,i_h2o2)
938         zycol(l, igcm_n2)      = rm(l,i_n2)
939         zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
940      end do
941
942      if (igcm_ch4 .ne. 0) then
943         do l = 1,lswitch-1
944            zycol(l,igcm_ch4) = rm(l,i_ch4)
945         end do
946      end if
947
948      end subroutine chimtogcm
949
950c*****************************************************************
951
952      subroutine chemrates(nlayer,
953     $                     lswitch, dens, press, t,
954     $                     surfdust1d, surfice1d,
955     $                     a001, a002, a003,
956     $                     b001, b002, b003, b004, b005, b006,
957     $                     b007, b008, b009,
958     $                     c001, c002, c003, c004, c005, c006,
959     $                     c007, c008, c009, c010, c011, c012,
960     $                     c013, c014, c015, c016, c017, c018,
961     $                     d001, d002, d003,
962     $                     e001, e002, e003,
963     $                     h001, h002, h003, h004, h005,
964     $                     t001, t002, t003, tau)
965
966c*****************************************************************
967
968      USE comcstfi_h
969      implicit none
970
971cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
972c     inputs:                                                    c
973cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
974
975      integer, intent(in) :: nlayer   ! number of atmospheric layers
976      integer :: lswitch              ! interface level between chemistries
977
978      real :: dens(nlayer)            ! density (cm-3)
979      real :: press(nlayer)           ! pressure (hpa)
980      real :: t(nlayer)               ! temperature (k)
981      real :: surfdust1d(nlayer)      ! dust surface area (cm^2/cm^3)
982      real :: surfice1d(nlayer)       ! ice surface area (cm^2/cm^3)
983      real :: tau                     ! dust opacity at 7 hpa
984
985      real, parameter :: tribo   = 0. ! switch for triboelectricity
986
987cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
988c     outputs:                                                   c
989cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
990
991      real :: a001(nlayer), a002(nlayer), a003(nlayer)
992      real :: b001(nlayer), b002(nlayer), b003(nlayer),
993     $        b004(nlayer), b005(nlayer), b006(nlayer),
994     $        b007(nlayer), b008(nlayer), b009(nlayer)
995      real :: c001(nlayer), c002(nlayer), c003(nlayer),
996     $        c004(nlayer), c005(nlayer), c006(nlayer),
997     $        c007(nlayer), c008(nlayer), c009(nlayer),
998     $        c010(nlayer), c011(nlayer), c012(nlayer),
999     $        c013(nlayer), c014(nlayer), c015(nlayer),
1000     $        c016(nlayer), c017(nlayer), c018(nlayer)
1001      real :: d001(nlayer), d002(nlayer), d003(nlayer)
1002      real :: e001(nlayer), e002(nlayer), e003(nlayer)
1003      real :: h001(nlayer), h002(nlayer), h003(nlayer),
1004     $        h004(nlayer), h005(nlayer)
1005      real :: t001(nlayer), t002(nlayer), t003(nlayer)
1006
1007cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1008c     local:                                                     c
1009cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1010
1011      real :: ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
1012      real :: ef, efmax, lossh2o, lossch4, lossco2
1013
1014      integer :: l
1015      real :: k1a, k1b, k1a0, k1b0, k1ainf
1016      real :: x, y, fc, fx
1017
1018cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1019c     compute reaction rates
1020cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1021
1022      do l = 1,lswitch-1
1023
1024cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1025c        oxygen compounds
1026cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1027c
1028ccc      a001: o + o2 + co2 -> o3 + co2
1029c
1030c        jpl 2003
1031c
1032c        co2 efficiency as a third body (2.075)
1033c        from sehested et al., j. geophys. res., 100, 1995.
1034c
1035         a001(l) = 2.075
1036     $             *6.0e-34*(t(l)/300.)**(-2.4)*dens(l)
1037c
1038c        mulcahy and williams, 1968
1039c
1040c        a001(l) = 2.68e-33*(t(l)/298.)**(-2.4)*dens(l)
1041c
1042c        nair et al., 1994
1043c
1044c        a001(l) = 1.3e-34*exp(724./t(l))*dens(l)
1045c
1046ccc      a002: o + o + co2 -> o2 + co2
1047c
1048c        Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
1049c
1050c        a002(l) = 2.5*5.2e-35*exp(900./t(l))*dens(l)
1051c
1052c        Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
1053c
1054c        a002(l) = 1.2e-32*(300./t(l))**(2.0)*dens(l)  ! yung expression
1055c
1056         a002(l) = 2.5*9.46e-34*exp(485./t(l))*dens(l) ! nist expression
1057c
1058c        baulch et al., 1976 confirmed by smith and robertson, 2008
1059c
1060c        a002(l) = 2.5*2.76e-34*exp(720./t(l))*dens(l)
1061c
1062ccc      a003: o + o3 -> o2 + o2
1063c
1064c        jpl 2003
1065c
1066         a003(l) = 8.0e-12*exp(-2060./t(l))
1067c
1068cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1069c        reactions with o(1d)
1070cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1071c
1072ccc      b001: o(1d) + co2  -> o + co2
1073c
1074c        jpl 2003
1075c
1076c        b001(l) = 7.4e-11*exp(120./t(l))
1077c
1078c        jpl 2006
1079c
1080         b001(l) = 7.5e-11*exp(115./t(l))
1081c
1082ccc      b002: o(1d) + h2o  -> oh + oh
1083c
1084c        jpl 2003
1085c
1086c        b002(l) = 2.2e-10
1087c
1088c        jpl 2006
1089c
1090         b002(l) = 1.63e-10*exp(60./t(l))
1091c   
1092ccc      b003: o(1d) + h2  -> oh + h
1093c
1094c        jpl 2011
1095c
1096         b003(l) = 1.2e-10
1097c   
1098ccc      b004: o(1d) + o2  -> o + o2
1099c
1100c        jpl 2003
1101c
1102c        b004(l) = 3.2e-11*exp(70./t(l))
1103c
1104c        jpl 2006
1105c
1106         b004(l) = 3.3e-11*exp(55./t(l))
1107c   
1108ccc      b005: o(1d) + o3  -> o2 + o2
1109c
1110c        jpl 2003
1111c
1112         b005(l) = 1.2e-10
1113c   
1114ccc      b006: o(1d) + o3  -> o2 + o + o
1115c
1116c        jpl 2003
1117c
1118         b006(l) = 1.2e-10
1119c   
1120ccc      b007: o(1d) + ch4 -> ch3 + oh
1121c
1122c        jpl 2003
1123c
1124         b007(l) = 1.5e-10*0.75
1125c
1126ccc      b008: o(1d) + ch4 -> ch3o + h
1127c
1128c        jpl 2003
1129c
1130         b008(l) = 1.5e-10*0.20
1131c
1132ccc      b009: o(1d) + ch4 -> ch2o + h2
1133c
1134c        jpl 2003
1135c
1136         b009(l) = 1.5e-10*0.05
1137c
1138cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1139c        hydrogen compounds   
1140cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1141c
1142ccc      c001: o + ho2 -> oh + o2
1143c
1144c        jpl 2003
1145c
1146         c001(l) = 3.0e-11*exp(200./t(l))
1147c
1148ccc      c002: o + oh -> o2 + h
1149c
1150c        jpl 2011
1151c
1152         c002(l) = 1.8e-11*exp(180./t(l))
1153c
1154c        robertson and smith, j. chem. phys. a 110, 6673, 2006
1155c
1156c        c002(l) = 11.2e-11*t(l)**(-0.32)*exp(177./t(l))
1157c
1158ccc      c003: h + o3 -> oh + o2
1159c
1160c        jpl 2003
1161c
1162         c003(l) = 1.4e-10*exp(-470./t(l))
1163c
1164ccc      c004: h + ho2 -> oh + oh
1165c
1166c        jpl 2003
1167c
1168c        c004(l) = 8.1e-11*0.90
1169c
1170c        jpl 2006
1171c
1172         c004(l) = 7.2e-11
1173c
1174ccc      c005: h + ho2 -> h2 + o2
1175c
1176c        jpl 2003
1177c
1178c        c005(l) = 8.1e-11*0.08
1179c
1180c        jpl 2006
1181c
1182         c005(l) = 6.9e-12
1183c
1184ccc      c006: h + ho2 -> h2o + o
1185c
1186c        jpl 2003
1187c
1188c        c006(l) = 8.1e-11*0.02
1189c
1190c        jpl 2006
1191c
1192         c006(l) = 1.6e-12
1193c
1194ccc      c007: oh + ho2 -> h2o + o2
1195c
1196c        jpl 2003
1197c
1198         c007(l) = 4.8e-11*exp(250./t(l))
1199c
1200c        jpl 2003 +20% d'apres canty et al., grl, 2006
1201c
1202c        c007(l) = 4.8e-11*exp(250./t(l))*1.2
1203c
1204ccc      c008: ho2 + ho2 -> h2o2 + o2
1205c
1206c        jpl 2003
1207c
1208c        c008(l) = 2.3e-13*exp(600./t(l))
1209c
1210c        christensen et al., grl, 13, 2002
1211c
1212         c008(l) = 1.5e-12*exp(19./t(l))
1213c
1214ccc      c009: oh + h2o2 -> h2o + ho2
1215c
1216c        jpl 2003
1217c
1218c        c009(l) = 2.9e-12*exp(-160./t(l))
1219c
1220c        jpl 2006
1221c
1222         c009(l) = 1.8e-12
1223c
1224ccc      c010: oh + h2 -> h2o + h
1225c
1226c        jpl 2003
1227c
1228c        c010(l) = 5.5e-12*exp(-2000./t(l))
1229c
1230c        jpl 2006
1231c
1232         c010(l) = 2.8e-12*exp(-1800./t(l))
1233c
1234ccc      c011: h + o2 + co2 -> ho2 + co2
1235c
1236c        jpl 2011
1237c
1238         ak0 = 2.5*4.4e-32*(t(l)/300.)**(-1.3)
1239         ak1 = 7.5e-11*(t(l)/300.)**(0.2)
1240c
1241         rate = (ak0*dens(l))/(1. + ak0*dens(l)/ak1)
1242         xpo = 1./(1. + alog10((ak0*dens(l))/ak1)**2)
1243         c011(l) = rate*0.6**xpo
1244c
1245ccc      c012: o + h2o2 -> oh + ho2
1246c
1247c        jpl 2003
1248c
1249         c012(l) = 1.4e-12*exp(-2000./t(l))
1250c
1251ccc      c013: oh + oh -> h2o + o
1252c
1253c        jpl 2003
1254c
1255c        c013(l) = 4.2e-12*exp(-240./t(l))
1256c
1257c        jpl 2006
1258c
1259         c013(l) = 1.8e-12
1260c
1261ccc      c014: oh + o3 -> ho2 + o2
1262c
1263c        jpl 2003
1264c
1265         c014(l) = 1.7e-12*exp(-940./t(l))
1266c
1267c        jpl 2000
1268c
1269c        c014(l) = 1.5e-12*exp(-880./t(l))
1270c
1271c        nair et al., 1994 (jpl 1997)
1272c
1273c        c014(l) = 1.6e-12*exp(-940./t(l))
1274c
1275ccc      c015: ho2 + o3 -> oh + o2 + o2
1276c
1277c        jpl 2003
1278c
1279         c015(l) = 1.0e-14*exp(-490./t(l))
1280c
1281c        jpl 2000
1282c
1283c        c015(l) = 2.0e-14*exp(-680./t(l))
1284c
1285c        nair et al., 1994 (jpl 1997)
1286c
1287c        c015(l) = 1.1e-14*exp(-500./t(l))
1288c
1289ccc      c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
1290c
1291c        jpl 2011
1292c
1293         c016(l) = 2.5*2.1e-33
1294     $              *exp(920./t(l))*dens(l)
1295c
1296ccc      c017: oh + oh + co2 -> h2o2 + co2
1297c
1298c        jpl 2003
1299c
1300         ak0 = 2.5*6.9e-31*(t(l)/300.)**(-1.0)
1301         ak1 = 2.6e-11*(t(l)/300.)**(0.0)
1302c
1303c        jpl 1997
1304c
1305c        ak0 = 2.5*6.2e-31*(t(l)/300.)**(-1.0)
1306c        ak1 = 2.6e-11*(t(l)/300.)**(0.0)
1307c
1308c        nair et al., 1994
1309c
1310c        ak0 = 2.5*7.1e-31*(t(l)/300.)**(-0.8)
1311c        ak1 = 1.5e-11*(t(l)/300.)**(0.0)
1312c
1313         rate = (ak0*dens(l))/(1. + ak0*dens(l)/ak1)
1314         xpo = 1./(1. + alog10((ak0*dens(l))/ak1)**2)
1315         c017(l) = rate*0.6**xpo
1316c
1317ccc      c018: h + h + co2 -> h2 + co2
1318c
1319c        baulch et al., 1992
1320c
1321c        c018(l) = 2.5*8.85e-33*(t(l)/298.)**(-0.6)*dens(l)
1322c
1323c        baulch et al., 2005
1324c
1325         c018(l) = 2.5*1.8e-30*(t(l)**(-1.0))*dens(l)
1326c
1327cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1328c        nitrogen compounds
1329cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1330c
1331ccc      d001: no2 + o -> no + o2
1332c
1333c        jpl 2003
1334c
1335c        d001(l) = 5.6e-12*exp(180./t(l))
1336c
1337ccc      jpl 2006
1338c
1339         d001(l) = 5.1e-12*exp(210./t(l))
1340c
1341ccc      d002: no + o3 -> no2 + o2
1342c
1343c        jpl 2003
1344c
1345         d002(l) = 3.0e-12*exp(-1500./t(l))
1346c
1347ccc      d003: no + ho2 -> no2 + oh
1348c
1349c        jpl 2011
1350c
1351         d003(l) = 3.3e-12*exp(270./t(l))
1352c
1353cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1354c        carbon compounds
1355cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1356c
1357ccc      e001: oh + co -> co2 + h
1358c
1359c        jpl 2003
1360c
1361c        e001(l) = 1.5e-13*(1 + 0.6*press(l)/1013.)
1362c
1363c        mccabe et al., grl, 28, 3135, 2001
1364c
1365c        e001(l) = 1.57e-13 + 3.54e-33*dens(l)
1366c
1367c        atkinson et al. 2006
1368c
1369c        e001(l) = 1.44e-13 + 3.43e-33*dens(l)
1370c
1371c        joshi et al., 2006
1372c
1373         k1a0 = 1.34*2.5*dens(l)
1374     $         *1/(1/(3.62e-26*t(l)**(-2.739)*exp(-20./t(l)))
1375     $         + 1/(6.48e-33*t(l)**(0.14)*exp(-57./t(l))))   ! corrige de l'erreur publi
1376         k1b0 = 1.17e-19*t(l)**(2.053)*exp(139./t(l))
1377     $        + 9.56e-12*t(l)**(-0.664)*exp(-167./t(l))
1378         k1ainf = 1.52e-17*t(l)**(1.858)*exp(28.8/t(l))
1379     $          + 4.78e-8*t(l)**(-1.851)*exp(-318./t(l))
1380         x = k1a0/(k1ainf - k1b0)
1381         y = k1b0/(k1ainf - k1b0)
1382         fc = 0.628*exp(-1223./t(l)) + (1. - 0.628)*exp(-39./t(l))
1383     $      + exp(-t(l)/255.)
1384         fx = fc**(1./(1. + (alog(x))**2))                   ! corrige de l'erreur publi
1385         k1a = k1a0*((1. + y)/(1. + x))*fx
1386         k1b = k1b0*(1./(1.+x))*fx
1387c
1388         e001(l) = k1a + k1b
1389c
1390ccc      e002: o + co + m -> co2 + m
1391c
1392c        tsang and hampson, 1986.
1393c
1394         e002(l) = 2.5*6.5e-33*exp(-2184./t(l))*dens(l)
1395c
1396c        baulch et al., butterworths, 1976.
1397c
1398c        e002(l) = 1.6e-32*exp(-2184./t(l))*dens(l)
1399c
1400ccc      e003: ch4 + oh -> ch3 + h2o
1401c
1402c        jpl 2003
1403c
1404c        e003(l) = 2.45e-12*exp(-1775./t(l))
1405c
1406c        jpl 2003, three-parameter expression
1407c
1408         e003(l) = 2.80e-14*(t(l)**0.667)*exp(-1575./t(l))
1409c
1410cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1411c        heterogenous chemistry
1412cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1413c
1414c        k = (surface*v*gamma)/4 (s-1)
1415c        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
1416c
1417ccc      h001: ho2 + ice -> products
1418c
1419c        cooper and abbatt, 1996: gamma = 0.025
1420c
1421         h001(l) = surfice1d(l)
1422     $            *100.*sqrt(8.*8.31*t(l)/(33.e-3*pi))*0.025/4.
1423c
1424c        h002: oh + ice -> products
1425c
1426c        cooper and abbatt, 1996: gamma = 0.03
1427c
1428         h002(l) = surfice1d(l)
1429     $            *100.*sqrt(8.*8.31*t(l)/(17.e-3*pi))*0.03/4.
1430c
1431c        h003: ho2 + dust -> products
1432c
1433c        jacob, 2000: gamma = 0.2
1434c        see dereus et al., atm. chem. phys., 2005
1435c
1436c        h003(l) = surfdust1d(l)
1437c    $            *100.*sqrt(8.*8.31*t(l)/(33.e-3*pi))*0.2/4.
1438         h003(l) = 0.     ! advised
1439c
1440ccc      h004: h2o2 + ice -> products
1441c
1442c        gamma = 1.e-3    test value
1443c
1444c        h004(l) = surfice1d(l)
1445c    $            *100.*sqrt(8.*8.31*t(l)/(34.e-3*pi))*0.001/4.
1446         h004(l) = 0.     ! advised
1447c
1448c        h005: h2o2 + dust -> products
1449c
1450c        gamma = 5.e-4
1451c        see dereus et al., atm. chem. phys., 2005
1452c
1453         h005(l) = surfdust1d(l)
1454     $            *100.*sqrt(8.*8.31*t(l)/(34.e-3*pi))*5.e-4/4.
1455         h005(l) = 0.     ! advised
1456c
1457      end do
1458c
1459      if (tribo .eq. 1.) then
1460c
1461c        electrochemical reactions
1462c
1463c        efmax: maximum electric field (kv.m-1)
1464c
1465         efmax = 23.3
1466c
1467c        ef: actual electric field, scaled by tau.
1468c
1469c        if (tau .ge. 1.) then
1470c           ef = efmax
1471c        else
1472c           ef = 0.
1473c        end if
1474c        ef = min(efmax,efmax*tau/1.0)
1475c
1476         ef = (efmax/0.5)*tau - (efmax/0.5)*0.5
1477c
1478         ef = max(ef, 0.)
1479         ef = min(ef, efmax)
1480c
1481ccc      t001: h2o + e -> oh + h-
1482c
1483c        lossh2o: fit of oh/h- production rates
1484c        given by delory et al., astrobiology, 6, 451, 2006
1485c
1486         if (ef .eq. 0.) then
1487            lossh2o = 0.
1488         else if (ef .lt. 10.) then
1489            lossh2o = 0.054136*exp(1.0978*ef)
1490         else if (ef .lt. 16.) then
1491            lossh2o = 64.85*exp(0.38894*ef)
1492         else if (ef .le. 20.) then
1493            lossh2o = 0.2466*exp(0.73719*ef)
1494         else
1495            lossh2o = 2.3269e-8*exp(1.546*ef)
1496         end if
1497c
1498c        production rates are given for h2o = 20 prec. microns.
1499c        t001 is converted to first-order reaction rate
1500c        assuming h2o number density at the surface =  5e13 mol cm-3
1501c
1502         do l = 1,21                     ! 70 km
1503            t001(l) = lossh2o/5.e13      ! s-1
1504         end do
1505         do l = 22,lswitch-1
1506            t001(l) = 0.
1507         end do
1508c
1509ccc      t002: ch4 + e -> products
1510c
1511c        lossch4: fit of ch4 loss rates
1512c        given by farrell et al., grl, 33, 2006
1513c
1514         if (ef .eq. 0.) then
1515            lossch4 = 0.
1516         else if (ef .gt. 20.) then
1517            lossch4 = 1.113e-21*exp(1.6065*ef)
1518         else if (ef .gt. 17.5) then
1519            lossch4 = 1.e-15*exp(0.92103*ef)
1520         else if (ef .gt. 14.) then
1521            lossch4 = 1.e-13*exp(0.65788*ef)
1522         else
1523            lossch4 = 8.9238e-15*exp(0.835*ef)
1524         end if
1525c
1526         do l = 1,21               ! 70 km
1527            t002(l) = lossch4      ! s-1
1528         end do
1529         do l = 22,lswitch-1
1530            t002(l) = 0.
1531         end do
1532c
1533ccc      t003: co2 + e -> co + o-
1534c
1535c        lossco2: fit of co/o- production rates
1536c        given by delory et al., astrobiology, 6, 451, 2006
1537c
1538         if (ef .eq. 0.) then
1539            lossco2 = 0.
1540         else if (ef .lt. 10.) then
1541            lossco2 = 22.437*exp(1.045*ef)
1542         else if (ef .lt. 16.) then
1543            lossco2 = 17518.*exp(0.37896*ef)
1544         else if (ef .lt. 20.) then
1545            lossco2 = 54.765*exp(0.73946*ef)
1546         else
1547            lossco2 = 4.911e-6*exp(1.5508*ef)
1548         end if
1549c
1550c        production rates are assumed to be given for p = 6 hPa
1551c        lossco2 is converted to first-order reaction rate
1552c        assuming co2 number density at the surface =  2e17 mol cm-3
1553c
1554         do l = 1,21                     ! 70 km
1555            t003(l) = lossco2/2.e17      ! s-1
1556         end do
1557         do l = 22,lswitch-1
1558            t003(l) = 0.
1559         end do
1560      else
1561         do l = 1,lswitch-1
1562            t001(l) = 0.
1563            t002(l) = 0.
1564            t003(l) = 0.
1565         end do
1566      end if
1567c
1568      end subroutine chemrates
1569
1570      end subroutine photochemistry
Note: See TracBrowser for help on using the repository browser.