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

Last change on this file since 830 was 690, checked in by acolaitis, 13 years ago

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