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

Last change on this file since 1430 was 1430, checked in by flefevre, 10 years ago

aeronomars : preparation de l'implementation du solveur chimique ASIS

1) extension de la table de photodissociation "jmars" a la thermosphere

et restriction aux seules especes necessaires.

2) modernisation de la subroutine d'interpolation dans cette

table de photodissociation.

3) chimiedata.h compatible avec solveur ASIS

ATTENTION!! la nouvelle table de photodissociation

jmars.20140930

doit etre presente dans /datadir sinon crash modele...
File size: 50.0 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      return
170      end
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      return
758      end
759
760c*****************************************************************
761
762      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
763
764c*****************************************************************
765
766      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
767     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
768     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
769     &                      igcm_ch4
770
771      implicit none
772
773#include "callkeys.h"
774
775cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
776c     inputs:
777cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
778     
779      integer, intent(in) :: nlayer ! number of atmospheric layers
780      integer, intent(in) :: nq     ! number of tracers
781      real :: zycol(nlayer,nq)      ! species volume mixing ratio in the gcm
782
783      integer :: nesp               ! number of species in the chemistry
784      integer :: lswitch            ! interface level between chemistries
785
786cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
787c     outputs:
788cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
789
790      real :: rm(nlayer,nesp)   ! species volume mixing ratio
791     
792cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
793c     local:
794cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
795
796      integer :: l, iq
797     
798c     tracer indexes in the chemistry:
799
800      integer,parameter :: i_co2  =  1
801      integer,parameter :: i_co   =  2
802      integer,parameter :: i_o    =  3
803      integer,parameter :: i_o1d  =  4
804      integer,parameter :: i_o2   =  5
805      integer,parameter :: i_o3   =  6
806      integer,parameter :: i_h    =  7
807      integer,parameter :: i_h2   =  8
808      integer,parameter :: i_oh   =  9
809      integer,parameter :: i_ho2  = 10
810      integer,parameter :: i_h2o2 = 11
811      integer,parameter :: i_ch4  = 12
812      integer,parameter :: i_h2o  = 13
813      integer,parameter :: i_n2   = 14
814      integer,parameter :: i_hox  = 15
815      integer,parameter :: i_ox   = 16
816
817cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
818c     initialise chemical species
819cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
820
821      do l = 1,lswitch-1
822         rm(l,i_co2)  = max(zycol(l, igcm_co2),     1.e-30)
823         rm(l,i_co)   = max(zycol(l, igcm_co),      1.e-30)
824         rm(l,i_o)    = max(zycol(l, igcm_o),       1.e-30)
825         rm(l,i_o1d)  = max(zycol(l, igcm_o1d),     1.e-30)
826         rm(l,i_o2)   = max(zycol(l, igcm_o2),      1.e-30)
827         rm(l,i_o3)   = max(zycol(l, igcm_o3),      1.e-30)
828         rm(l,i_h)    = max(zycol(l, igcm_h),       1.e-30)
829         rm(l,i_h2)   = max(zycol(l, igcm_h2),      1.e-30)
830         rm(l,i_oh)   = max(zycol(l, igcm_oh),      1.e-30)
831         rm(l,i_ho2)  = max(zycol(l, igcm_ho2),     1.e-30)
832         rm(l,i_h2o2) = max(zycol(l, igcm_h2o2),    1.e-30)
833         rm(l,i_n2)   = max(zycol(l, igcm_n2),      1.e-30)
834         rm(l,i_h2o)  = max(zycol(l, igcm_h2o_vap), 1.e-30)
835      end do
836
837      if (igcm_ch4 .eq. 0) then
838         do l = 1,lswitch-1
839            rm(l,i_ch4) = 0.
840         end do
841      else
842         do l = 1,lswitch-1
843            rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30)
844         end do
845      end if
846
847cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
848c     initialise chemical families                     c
849cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
850
851      do l = 1,lswitch-1
852         rm(l,i_hox) = rm(l,i_h)
853     $               + rm(l,i_oh)
854     $               + rm(l,i_ho2)
855         rm(l,i_ox)  = rm(l,i_o)
856     $               + rm(l,i_o3)
857      end do
858
859      return
860      end
861
862c*****************************************************************
863c
864      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
865c
866c*****************************************************************
867c
868      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
869     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
870     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
871     &                      igcm_ch4
872
873      implicit none
874
875#include "callkeys.h"
876
877cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
878c     inputs:
879cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
880
881      integer, intent(in) :: nlayer ! number of atmospheric layers
882      integer, intent(in) :: nq     ! number of tracers
883      integer :: nesp               ! number of species in the chemistry
884      integer :: lswitch            ! interface level between chemistries
885     
886      real :: rm(nlayer,nesp)       ! species volume mixing ratio
887
888cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
889c     output:
890cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
891
892      real :: zycol(nlayer,nq)  ! species volume mixing ratio in the gcm
893
894cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
895c     local:
896cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
897
898      integer l, iq
899     
900c     tracer indexes in the chemistry:
901
902      integer,parameter :: i_co2  =  1
903      integer,parameter :: i_co   =  2
904      integer,parameter :: i_o    =  3
905      integer,parameter :: i_o1d  =  4
906      integer,parameter :: i_o2   =  5
907      integer,parameter :: i_o3   =  6
908      integer,parameter :: i_h    =  7
909      integer,parameter :: i_h2   =  8
910      integer,parameter :: i_oh   =  9
911      integer,parameter :: i_ho2  = 10
912      integer,parameter :: i_h2o2 = 11
913      integer,parameter :: i_ch4  = 12
914      integer,parameter :: i_h2o  = 13
915      integer,parameter :: i_n2   = 14
916      integer,parameter :: i_hox  = 15
917      integer,parameter :: i_ox   = 16
918
919cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
920c     save mixing ratios for the gcm
921cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
922
923      do l = 1,lswitch-1
924         zycol(l, igcm_co2)     = rm(l,i_co2)
925         zycol(l, igcm_co)      = rm(l,i_co)
926         zycol(l, igcm_o)       = rm(l,i_o)
927         zycol(l, igcm_o1d)     = rm(l,i_o1d)
928         zycol(l, igcm_o2)      = rm(l,i_o2)
929         zycol(l, igcm_o3)      = rm(l,i_o3)
930         zycol(l, igcm_h)       = rm(l,i_h) 
931         zycol(l, igcm_h2)      = rm(l,i_h2)
932         zycol(l, igcm_oh)      = rm(l,i_oh)
933         zycol(l, igcm_ho2)     = rm(l,i_ho2)
934         zycol(l, igcm_h2o2)    = rm(l,i_h2o2)
935         zycol(l, igcm_n2)      = rm(l,i_n2)
936         zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
937      end do
938
939      if (igcm_ch4 .ne. 0) then
940         do l = 1,lswitch-1
941            zycol(l,igcm_ch4) = rm(l,i_ch4)
942         end do
943      end if
944
945      return
946      end
947
948c*****************************************************************
949
950      subroutine chemrates(nlayer,
951     $                     lswitch, dens, press, t,
952     $                     surfdust1d, surfice1d,
953     $                     a001, a002, a003,
954     $                     b001, b002, b003, b004, b005, b006,
955     $                     b007, b008, b009,
956     $                     c001, c002, c003, c004, c005, c006,
957     $                     c007, c008, c009, c010, c011, c012,
958     $                     c013, c014, c015, c016, c017, c018,
959     $                     d001, d002, d003,
960     $                     e001, e002, e003,
961     $                     h001, h002, h003, h004, h005,
962     $                     t001, t002, t003, tau)
963
964c*****************************************************************
965
966      USE comcstfi_h
967      implicit none
968
969cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
970c     inputs:                                                    c
971cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
972
973      integer, intent(in) :: nlayer   ! number of atmospheric layers
974      integer :: lswitch              ! interface level between chemistries
975
976      real :: dens(nlayer)            ! density (cm-3)
977      real :: press(nlayer)           ! pressure (hpa)
978      real :: t(nlayer)               ! temperature (k)
979      real :: surfdust1d(nlayer)      ! dust surface area (cm^2/cm^3)
980      real :: surfice1d(nlayer)       ! ice surface area (cm^2/cm^3)
981      real :: tau                     ! dust opacity at 7 hpa
982
983      real, parameter :: tribo   = 0. ! switch for triboelectricity
984
985cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
986c     outputs:                                                   c
987cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
988
989      real :: a001(nlayer), a002(nlayer), a003(nlayer)
990      real :: b001(nlayer), b002(nlayer), b003(nlayer),
991     $        b004(nlayer), b005(nlayer), b006(nlayer),
992     $        b007(nlayer), b008(nlayer), b009(nlayer)
993      real :: c001(nlayer), c002(nlayer), c003(nlayer),
994     $        c004(nlayer), c005(nlayer), c006(nlayer),
995     $        c007(nlayer), c008(nlayer), c009(nlayer),
996     $        c010(nlayer), c011(nlayer), c012(nlayer),
997     $        c013(nlayer), c014(nlayer), c015(nlayer),
998     $        c016(nlayer), c017(nlayer), c018(nlayer)
999      real :: d001(nlayer), d002(nlayer), d003(nlayer)
1000      real :: e001(nlayer), e002(nlayer), e003(nlayer)
1001      real :: h001(nlayer), h002(nlayer), h003(nlayer),
1002     $        h004(nlayer), h005(nlayer)
1003      real :: t001(nlayer), t002(nlayer), t003(nlayer)
1004
1005cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1006c     local:                                                     c
1007cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1008
1009      real :: ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
1010      real :: ef, efmax, lossh2o, lossch4, lossco2
1011
1012      integer :: l
1013      real :: k1a, k1b, k1a0, k1b0, k1ainf
1014      real :: x, y, fc, fx
1015
1016cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1017c     compute reaction rates
1018cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1019
1020      do l = 1,lswitch-1
1021
1022cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1023c        oxygen compounds
1024cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1025c
1026ccc      a001: o + o2 + co2 -> o3 + co2
1027c
1028c        jpl 2003
1029c
1030c        co2 efficiency as a third body (2.075)
1031c        from sehested et al., j. geophys. res., 100, 1995.
1032c
1033         a001(l) = 2.075
1034     $             *6.0e-34*(t(l)/300.)**(-2.4)*dens(l)
1035c
1036c        mulcahy and williams, 1968
1037c
1038c        a001(l) = 2.68e-33*(t(l)/298.)**(-2.4)*dens(l)
1039c
1040c        nair et al., 1994
1041c
1042c        a001(l) = 1.3e-34*exp(724./t(l))*dens(l)
1043c
1044ccc      a002: o + o + co2 -> o2 + co2
1045c
1046c        Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
1047c
1048c        a002(l) = 2.5*5.2e-35*exp(900./t(l))*dens(l)
1049c
1050c        Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
1051c
1052c        a002(l) = 1.2e-32*(300./t(l))**(2.0)*dens(l)  ! yung expression
1053c
1054         a002(l) = 2.5*9.46e-34*exp(485./t(l))*dens(l) ! nist expression
1055c
1056c        baulch et al., 1976 confirmed by smith and robertson, 2008
1057c
1058c        a002(l) = 2.5*2.76e-34*exp(720./t(l))*dens(l)
1059c
1060ccc      a003: o + o3 -> o2 + o2
1061c
1062c        jpl 2003
1063c
1064         a003(l) = 8.0e-12*exp(-2060./t(l))
1065c
1066cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1067c        reactions with o(1d)
1068cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1069c
1070ccc      b001: o(1d) + co2  -> o + co2
1071c
1072c        jpl 2003
1073c
1074c        b001(l) = 7.4e-11*exp(120./t(l))
1075c
1076c        jpl 2006
1077c
1078         b001(l) = 7.5e-11*exp(115./t(l))
1079c
1080ccc      b002: o(1d) + h2o  -> oh + oh
1081c
1082c        jpl 2003
1083c
1084c        b002(l) = 2.2e-10
1085c
1086c        jpl 2006
1087c
1088         b002(l) = 1.63e-10*exp(60./t(l))
1089c   
1090ccc      b003: o(1d) + h2  -> oh + h
1091c
1092c        jpl 2011
1093c
1094         b003(l) = 1.2e-10
1095c   
1096ccc      b004: o(1d) + o2  -> o + o2
1097c
1098c        jpl 2003
1099c
1100c        b004(l) = 3.2e-11*exp(70./t(l))
1101c
1102c        jpl 2006
1103c
1104         b004(l) = 3.3e-11*exp(55./t(l))
1105c   
1106ccc      b005: o(1d) + o3  -> o2 + o2
1107c
1108c        jpl 2003
1109c
1110         b005(l) = 1.2e-10
1111c   
1112ccc      b006: o(1d) + o3  -> o2 + o + o
1113c
1114c        jpl 2003
1115c
1116         b006(l) = 1.2e-10
1117c   
1118ccc      b007: o(1d) + ch4 -> ch3 + oh
1119c
1120c        jpl 2003
1121c
1122         b007(l) = 1.5e-10*0.75
1123c
1124ccc      b008: o(1d) + ch4 -> ch3o + h
1125c
1126c        jpl 2003
1127c
1128         b008(l) = 1.5e-10*0.20
1129c
1130ccc      b009: o(1d) + ch4 -> ch2o + h2
1131c
1132c        jpl 2003
1133c
1134         b009(l) = 1.5e-10*0.05
1135c
1136cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1137c        hydrogen compounds   
1138cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1139c
1140ccc      c001: o + ho2 -> oh + o2
1141c
1142c        jpl 2003
1143c
1144         c001(l) = 3.0e-11*exp(200./t(l))
1145c
1146ccc      c002: o + oh -> o2 + h
1147c
1148c        jpl 2011
1149c
1150         c002(l) = 1.8e-11*exp(180./t(l))
1151c
1152c        robertson and smith, j. chem. phys. a 110, 6673, 2006
1153c
1154c        c002(l) = 11.2e-11*t(l)**(-0.32)*exp(177./t(l))
1155c
1156ccc      c003: h + o3 -> oh + o2
1157c
1158c        jpl 2003
1159c
1160         c003(l) = 1.4e-10*exp(-470./t(l))
1161c
1162ccc      c004: h + ho2 -> oh + oh
1163c
1164c        jpl 2003
1165c
1166c        c004(l) = 8.1e-11*0.90
1167c
1168c        jpl 2006
1169c
1170         c004(l) = 7.2e-11
1171c
1172ccc      c005: h + ho2 -> h2 + o2
1173c
1174c        jpl 2003
1175c
1176c        c005(l) = 8.1e-11*0.08
1177c
1178c        jpl 2006
1179c
1180         c005(l) = 6.9e-12
1181c
1182ccc      c006: h + ho2 -> h2o + o
1183c
1184c        jpl 2003
1185c
1186c        c006(l) = 8.1e-11*0.02
1187c
1188c        jpl 2006
1189c
1190         c006(l) = 1.6e-12
1191c
1192ccc      c007: oh + ho2 -> h2o + o2
1193c
1194c        jpl 2003
1195c
1196         c007(l) = 4.8e-11*exp(250./t(l))
1197c
1198c        jpl 2003 +20% d'apres canty et al., grl, 2006
1199c
1200c        c007(l) = 4.8e-11*exp(250./t(l))*1.2
1201c
1202ccc      c008: ho2 + ho2 -> h2o2 + o2
1203c
1204c        jpl 2003
1205c
1206c        c008(l) = 2.3e-13*exp(600./t(l))
1207c
1208c        christensen et al., grl, 13, 2002
1209c
1210         c008(l) = 1.5e-12*exp(19./t(l))
1211c
1212ccc      c009: oh + h2o2 -> h2o + ho2
1213c
1214c        jpl 2003
1215c
1216c        c009(l) = 2.9e-12*exp(-160./t(l))
1217c
1218c        jpl 2006
1219c
1220         c009(l) = 1.8e-12
1221c
1222ccc      c010: oh + h2 -> h2o + h
1223c
1224c        jpl 2003
1225c
1226c        c010(l) = 5.5e-12*exp(-2000./t(l))
1227c
1228c        jpl 2006
1229c
1230         c010(l) = 2.8e-12*exp(-1800./t(l))
1231c
1232ccc      c011: h + o2 + co2 -> ho2 + co2
1233c
1234c        jpl 2011
1235c
1236         ak0 = 2.5*4.4e-32*(t(l)/300.)**(-1.3)
1237         ak1 = 7.5e-11*(t(l)/300.)**(0.2)
1238c
1239         rate = (ak0*dens(l))/(1. + ak0*dens(l)/ak1)
1240         xpo = 1./(1. + alog10((ak0*dens(l))/ak1)**2)
1241         c011(l) = rate*0.6**xpo
1242c
1243ccc      c012: o + h2o2 -> oh + ho2
1244c
1245c        jpl 2003
1246c
1247         c012(l) = 1.4e-12*exp(-2000./t(l))
1248c
1249ccc      c013: oh + oh -> h2o + o
1250c
1251c        jpl 2003
1252c
1253c        c013(l) = 4.2e-12*exp(-240./t(l))
1254c
1255c        jpl 2006
1256c
1257         c013(l) = 1.8e-12
1258c
1259ccc      c014: oh + o3 -> ho2 + o2
1260c
1261c        jpl 2003
1262c
1263         c014(l) = 1.7e-12*exp(-940./t(l))
1264c
1265c        jpl 2000
1266c
1267c        c014(l) = 1.5e-12*exp(-880./t(l))
1268c
1269c        nair et al., 1994 (jpl 1997)
1270c
1271c        c014(l) = 1.6e-12*exp(-940./t(l))
1272c
1273ccc      c015: ho2 + o3 -> oh + o2 + o2
1274c
1275c        jpl 2003
1276c
1277         c015(l) = 1.0e-14*exp(-490./t(l))
1278c
1279c        jpl 2000
1280c
1281c        c015(l) = 2.0e-14*exp(-680./t(l))
1282c
1283c        nair et al., 1994 (jpl 1997)
1284c
1285c        c015(l) = 1.1e-14*exp(-500./t(l))
1286c
1287ccc      c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
1288c
1289c        jpl 2011
1290c
1291         c016(l) = 2.5*2.1e-33
1292     $              *exp(920./t(l))*dens(l)
1293c
1294ccc      c017: oh + oh + co2 -> h2o2 + co2
1295c
1296c        jpl 2003
1297c
1298         ak0 = 2.5*6.9e-31*(t(l)/300.)**(-1.0)
1299         ak1 = 2.6e-11*(t(l)/300.)**(0.0)
1300c
1301c        jpl 1997
1302c
1303c        ak0 = 2.5*6.2e-31*(t(l)/300.)**(-1.0)
1304c        ak1 = 2.6e-11*(t(l)/300.)**(0.0)
1305c
1306c        nair et al., 1994
1307c
1308c        ak0 = 2.5*7.1e-31*(t(l)/300.)**(-0.8)
1309c        ak1 = 1.5e-11*(t(l)/300.)**(0.0)
1310c
1311         rate = (ak0*dens(l))/(1. + ak0*dens(l)/ak1)
1312         xpo = 1./(1. + alog10((ak0*dens(l))/ak1)**2)
1313         c017(l) = rate*0.6**xpo
1314c
1315ccc      c018: h + h + co2 -> h2 + co2
1316c
1317c        baulch et al., 1992
1318c
1319c        c018(l) = 2.5*8.85e-33*(t(l)/298.)**(-0.6)*dens(l)
1320c
1321c        baulch et al., 2005
1322c
1323         c018(l) = 2.5*1.8e-30*(t(l)**(-1.0))*dens(l)
1324c
1325cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1326c        nitrogen compounds
1327cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1328c
1329ccc      d001: no2 + o -> no + o2
1330c
1331c        jpl 2003
1332c
1333c        d001(l) = 5.6e-12*exp(180./t(l))
1334c
1335ccc      jpl 2006
1336c
1337         d001(l) = 5.1e-12*exp(210./t(l))
1338c
1339ccc      d002: no + o3 -> no2 + o2
1340c
1341c        jpl 2003
1342c
1343         d002(l) = 3.0e-12*exp(-1500./t(l))
1344c
1345ccc      d003: no + ho2 -> no2 + oh
1346c
1347c        jpl 2011
1348c
1349         d003(l) = 3.3e-12*exp(270./t(l))
1350c
1351cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1352c        carbon compounds
1353cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1354c
1355ccc      e001: oh + co -> co2 + h
1356c
1357c        jpl 2003
1358c
1359c        e001(l) = 1.5e-13*(1 + 0.6*press(l)/1013.)
1360c
1361c        mccabe et al., grl, 28, 3135, 2001
1362c
1363c        e001(l) = 1.57e-13 + 3.54e-33*dens(l)
1364c
1365c        atkinson et al. 2006
1366c
1367c        e001(l) = 1.44e-13 + 3.43e-33*dens(l)
1368c
1369c        joshi et al., 2006
1370c
1371         k1a0 = 1.34*2.5*dens(l)
1372     $         *1/(1/(3.62e-26*t(l)**(-2.739)*exp(-20./t(l)))
1373     $         + 1/(6.48e-33*t(l)**(0.14)*exp(-57./t(l))))   ! corrige de l'erreur publi
1374         k1b0 = 1.17e-19*t(l)**(2.053)*exp(139./t(l))
1375     $        + 9.56e-12*t(l)**(-0.664)*exp(-167./t(l))
1376         k1ainf = 1.52e-17*t(l)**(1.858)*exp(28.8/t(l))
1377     $          + 4.78e-8*t(l)**(-1.851)*exp(-318./t(l))
1378         x = k1a0/(k1ainf - k1b0)
1379         y = k1b0/(k1ainf - k1b0)
1380         fc = 0.628*exp(-1223./t(l)) + (1. - 0.628)*exp(-39./t(l))
1381     $      + exp(-t(l)/255.)
1382         fx = fc**(1./(1. + (alog(x))**2))                   ! corrige de l'erreur publi
1383         k1a = k1a0*((1. + y)/(1. + x))*fx
1384         k1b = k1b0*(1./(1.+x))*fx
1385c
1386         e001(l) = k1a + k1b
1387c
1388ccc      e002: o + co + m -> co2 + m
1389c
1390c        tsang and hampson, 1986.
1391c
1392         e002(l) = 2.5*6.5e-33*exp(-2184./t(l))*dens(l)
1393c
1394c        baulch et al., butterworths, 1976.
1395c
1396c        e002(l) = 1.6e-32*exp(-2184./t(l))*dens(l)
1397c
1398ccc      e003: ch4 + oh -> ch3 + h2o
1399c
1400c        jpl 2003
1401c
1402c        e003(l) = 2.45e-12*exp(-1775./t(l))
1403c
1404c        jpl 2003, three-parameter expression
1405c
1406         e003(l) = 2.80e-14*(t(l)**0.667)*exp(-1575./t(l))
1407c
1408cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1409c        heterogenous chemistry
1410cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1411c
1412c        k = (surface*v*gamma)/4 (s-1)
1413c        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
1414c
1415ccc      h001: ho2 + ice -> products
1416c
1417c        cooper and abbatt, 1996: gamma = 0.025
1418c
1419         h001(l) = surfice1d(l)
1420     $            *100.*sqrt(8.*8.31*t(l)/(33.e-3*pi))*0.025/4.
1421c
1422c        h002: oh + ice -> products
1423c
1424c        cooper and abbatt, 1996: gamma = 0.03
1425c
1426         h002(l) = surfice1d(l)
1427     $            *100.*sqrt(8.*8.31*t(l)/(17.e-3*pi))*0.03/4.
1428c
1429c        h003: ho2 + dust -> products
1430c
1431c        jacob, 2000: gamma = 0.2
1432c        see dereus et al., atm. chem. phys., 2005
1433c
1434c        h003(l) = surfdust1d(l)
1435c    $            *100.*sqrt(8.*8.31*t(l)/(33.e-3*pi))*0.2/4.
1436         h003(l) = 0.     ! advised
1437c
1438ccc      h004: h2o2 + ice -> products
1439c
1440c        gamma = 1.e-3    test value
1441c
1442c        h004(l) = surfice1d(l)
1443c    $            *100.*sqrt(8.*8.31*t(l)/(34.e-3*pi))*0.001/4.
1444         h004(l) = 0.     ! advised
1445c
1446c        h005: h2o2 + dust -> products
1447c
1448c        gamma = 5.e-4
1449c        see dereus et al., atm. chem. phys., 2005
1450c
1451         h005(l) = surfdust1d(l)
1452     $            *100.*sqrt(8.*8.31*t(l)/(34.e-3*pi))*5.e-4/4.
1453         h005(l) = 0.     ! advised
1454c
1455      end do
1456c
1457      if (tribo .eq. 1.) then
1458c
1459c        electrochemical reactions
1460c
1461c        efmax: maximum electric field (kv.m-1)
1462c
1463         efmax = 23.3
1464c
1465c        ef: actual electric field, scaled by tau.
1466c
1467c        if (tau .ge. 1.) then
1468c           ef = efmax
1469c        else
1470c           ef = 0.
1471c        end if
1472c        ef = min(efmax,efmax*tau/1.0)
1473c
1474         ef = (efmax/0.5)*tau - (efmax/0.5)*0.5
1475c
1476         ef = max(ef, 0.)
1477         ef = min(ef, efmax)
1478c
1479ccc      t001: h2o + e -> oh + h-
1480c
1481c        lossh2o: fit of oh/h- production rates
1482c        given by delory et al., astrobiology, 6, 451, 2006
1483c
1484         if (ef .eq. 0.) then
1485            lossh2o = 0.
1486         else if (ef .lt. 10.) then
1487            lossh2o = 0.054136*exp(1.0978*ef)
1488         else if (ef .lt. 16.) then
1489            lossh2o = 64.85*exp(0.38894*ef)
1490         else if (ef .le. 20.) then
1491            lossh2o = 0.2466*exp(0.73719*ef)
1492         else
1493            lossh2o = 2.3269e-8*exp(1.546*ef)
1494         end if
1495c
1496c        production rates are given for h2o = 20 prec. microns.
1497c        t001 is converted to first-order reaction rate
1498c        assuming h2o number density at the surface =  5e13 mol cm-3
1499c
1500         do l = 1,21                     ! 70 km
1501            t001(l) = lossh2o/5.e13      ! s-1
1502         end do
1503         do l = 22,lswitch-1
1504            t001(l) = 0.
1505         end do
1506c
1507ccc      t002: ch4 + e -> products
1508c
1509c        lossch4: fit of ch4 loss rates
1510c        given by farrell et al., grl, 33, 2006
1511c
1512         if (ef .eq. 0.) then
1513            lossch4 = 0.
1514         else if (ef .gt. 20.) then
1515            lossch4 = 1.113e-21*exp(1.6065*ef)
1516         else if (ef .gt. 17.5) then
1517            lossch4 = 1.e-15*exp(0.92103*ef)
1518         else if (ef .gt. 14.) then
1519            lossch4 = 1.e-13*exp(0.65788*ef)
1520         else
1521            lossch4 = 8.9238e-15*exp(0.835*ef)
1522         end if
1523c
1524         do l = 1,21               ! 70 km
1525            t002(l) = lossch4      ! s-1
1526         end do
1527         do l = 22,lswitch-1
1528            t002(l) = 0.
1529         end do
1530c
1531ccc      t003: co2 + e -> co + o-
1532c
1533c        lossco2: fit of co/o- production rates
1534c        given by delory et al., astrobiology, 6, 451, 2006
1535c
1536         if (ef .eq. 0.) then
1537            lossco2 = 0.
1538         else if (ef .lt. 10.) then
1539            lossco2 = 22.437*exp(1.045*ef)
1540         else if (ef .lt. 16.) then
1541            lossco2 = 17518.*exp(0.37896*ef)
1542         else if (ef .lt. 20.) then
1543            lossco2 = 54.765*exp(0.73946*ef)
1544         else
1545            lossco2 = 4.911e-6*exp(1.5508*ef)
1546         end if
1547c
1548c        production rates are assumed to be given for p = 6 hPa
1549c        lossco2 is converted to first-order reaction rate
1550c        assuming co2 number density at the surface =  2e17 mol cm-3
1551c
1552         do l = 1,21                     ! 70 km
1553            t003(l) = lossco2/2.e17      ! s-1
1554         end do
1555         do l = 22,lswitch-1
1556            t003(l) = 0.
1557         end do
1558      else
1559         do l = 1,lswitch-1
1560            t001(l) = 0.
1561            t002(l) = 0.
1562            t003(l) = 0.
1563         end do
1564      end if
1565c
1566      return
1567      end
Note: See TracBrowser for help on using the repository browser.