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

Last change on this file since 1198 was 1125, checked in by flefevre, 11 years ago

calchim.F90 et photochemistry.F:

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