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

Last change on this file since 1540 was 1499, checked in by emillour, 10 years ago

Mars GCM:

  • Added option "-cpp" so that users can add the definition of specific precompiling flags in makegcm_* scripts (and added default BLAS and LAPACK flags for ifort on Gnome and Ada).
  • Put calls to dgesv (Lapack) routine in photochemistry_asis.F90 under the LAPACK preprocessing flag. Moved secondary routines in photochemistry.F and photochemistry_asis.F90 into the main (via contains instruction) to avoid multiple definitions of routines with identical names.

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