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

Last change on this file since 616 was 576, checked in by emillour, 13 years ago

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