source: trunk/LMDZ.MARS/libf/aeronomars/photochemistry_asis.F90 @ 1499

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

Mars GCM:

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

EM
File size: 55.9 KB
Line 
1!*****************************************************************
2!
3!     Photochemical routine
4!
5!     Author: Franck Lefevre
6!     ------
7!
8!     Version: 10/11/2015
9!
10!*****************************************************************
11
12subroutine photochemistry_asis(nlayer, nq,                            &
13                          ig, lswitch, zycol, sza, ptimestep, press,  &
14                          temp, dens, zmmean, dist_sol, surfdust1d,   &
15                          surfice1d, jo3, tau, iter)
16
17implicit none
18
19#include "chimiedata.h"
20#include "callkeys.h"
21
22!===================================================================
23!     inputs:
24!===================================================================
25
26integer, intent(in) :: nlayer ! number of atmospheric layers
27integer, intent(in) :: nq     ! number of tracers
28integer :: ig                 ! grid point index
29     
30real :: sza                   ! solar zenith angle (deg)
31real :: ptimestep             ! physics timestep (s)
32real :: press(nlayer)         ! pressure (hpa)
33real :: temp(nlayer)          ! temperature (k)
34real :: dens(nlayer)          ! density (cm-3)
35real :: zmmean(nlayer)        ! mean molar mass (g/mole)
36real :: dist_sol              ! sun distance (au)
37real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
38real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
39real :: tau                   ! optical depth at 7 hpa
40
41!===================================================================
42!     input/output:
43!===================================================================
44     
45real :: zycol(nlayer,nq)       ! chemical species volume mixing ratio
46
47!===================================================================
48!     output:
49!===================================================================
50     
51integer :: iter(nlayer)        ! iteration counter
52real    :: jo3(nlayer)         ! photodissociation rate o3 -> o1d
53
54!===================================================================
55!     local:
56!===================================================================
57
58integer :: phychemrat         ! (physical timestep)/(nominal chemical timestep)
59integer :: j_o3_o1d, ilev
60integer :: iesp, nesp
61integer :: lswitch
62
63logical, save :: firstcall = .true.
64
65parameter (nesp = 17)         ! number of species in the chemistry
66
67! tracer indexes in the chemistry:
68
69integer,parameter :: i_co2  =  1
70integer,parameter :: i_co   =  2
71integer,parameter :: i_o    =  3
72integer,parameter :: i_o1d  =  4
73integer,parameter :: i_o2   =  5
74integer,parameter :: i_o3   =  6
75integer,parameter :: i_h    =  7
76integer,parameter :: i_h2   =  8
77integer,parameter :: i_oh   =  9
78integer,parameter :: i_ho2  = 10
79integer,parameter :: i_h2o2 = 11
80integer,parameter :: i_h2o  = 12
81integer,parameter :: i_n    = 13
82integer,parameter :: i_n2d  = 14
83integer,parameter :: i_no   = 15
84integer,parameter :: i_no2  = 16
85integer,parameter :: i_n2   = 17
86
87real :: stimestep           ! standard timestep for the chemistry (s)
88real :: ctimestep           ! real timestep for the chemistry (s)
89real :: dt_guess            ! first-guess timestep (s)
90real :: dt_corrected        ! corrected timestep (s)
91real :: dt_min = 1.         ! minimum allowed timestep (s)
92real :: dtg                 ! correction factor for the timestep (s)
93real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
94real :: time                ! internal time (between 0 and ptimestep, in s)
95
96real, dimension(nlayer,nesp)            :: rm   ! mixing ratios
97real (kind = 8), dimension(nesp)        :: cold ! number densities at previous timestep (molecule.cm-3)
98real (kind = 8), dimension(nlayer,nesp) :: c    ! number densities at current timestep (molecule.cm-3)
99real (kind = 8), dimension(nesp)        :: cnew ! number densities at next timestep (molecule.cm-3)
100 
101! reaction rates
102 
103real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
104real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
105real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
106logical :: hetero_dust, hetero_ice
107
108! matrix
109
110real (kind = 8), dimension(nesp,nesp) :: mat, mat1
111integer, dimension(nesp)              :: indx
112integer                               :: code
113
114! production and loss terms (for first-guess solution only)
115
116real (kind = 8), dimension(nesp) :: prod, loss
117
118! curvatures
119
120real :: ratio, curv, e, e1, e2, e3
121
122!===================================================================
123!     initialisation of the reaction indexes
124!===================================================================
125
126if (firstcall) then
127   print*,'photochemistry: initialize indexes'
128   call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
129               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
130               i_n, i_n2d, i_no, i_no2, i_n2)
131   firstcall = .false.
132end if
133
134!===================================================================
135!     initialisation of mixing ratios and densities       
136!===================================================================
137
138call gcmtochim(nlayer, nq, zycol, lswitch, nesp,               &
139               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
140               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
141               i_n, i_n2d, i_no, i_no2, i_n2, dens, rm, c)
142
143!===================================================================
144!     interpolation of photolysis rates in the lookup table     
145!===================================================================
146
147call photolysis_asis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, &
148                     rm(:,i_co2), rm(:,i_o3), v_phot)
149
150! save o3 photolysis for output
151
152j_o3_o1d = 5
153jo3(:) = v_phot(:,j_o3_o1d)
154
155!===================================================================
156!     reaction rates                                     
157!===================================================================
158!     switches for heterogeneous chemistry
159!     hetero_ice  : reactions on ice clouds
160!     hetero_dust : reactions on dust   
161!===================================================================
162
163hetero_dust = .false.
164hetero_ice  = .true.
165
166call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2), &
167                   press, temp, hetero_dust, hetero_ice,         &
168                   surfdust1d, surfice1d, v_phot, v_3, v_4)
169
170!===================================================================
171!     stimestep : standard chemical timestep (s)                           
172!     ctimestep : real chemical timestep (s),
173!                 taking into account the physical timestep                           
174!===================================================================
175
176stimestep = 600. ! standard value : 10 mn
177
178phychemrat = nint(ptimestep/stimestep)
179phychemrat = 1 ! test 2014_004
180
181ctimestep = ptimestep/real(phychemrat)
182
183!print*, "stimestep  = ", stimestep
184!print*, "ptimestep  = ", ptimestep
185!print*, "phychemrat = ", phychemrat
186!print*, "ctimestep  = ", ctimestep
187!stop
188
189!===================================================================
190!     loop over levels         
191!===================================================================
192
193do ilev = 1,lswitch - 1
194
195!  initializations
196
197   time = 0.
198   iter(ilev) = 0
199   dt_guess = ctimestep
200   cold(:) = c(ilev,:)
201
202!  internal loop for the chemistry
203
204   do while (time < ptimestep)
205
206   iter(ilev) = iter(ilev) + 1    ! iteration counter
207 
208!  first-guess: fill matrix
209
210   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
211
212! first guess : form the matrix identity + mat*dt_guess
213
214   mat(:,:) = mat1(:,:)*dt_guess
215   do iesp = 1,nesp
216      mat(iesp,iesp) = 1. + mat(iesp,iesp)
217   end do
218
219!  first-guess: solve linear system
220     
221   cnew(:) = c(ilev,:)
222
223#ifdef LAPACK
224   call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
225#else
226   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
227   stop
228#endif
229!  eliminate small values
230
231   where (cnew(:)/dens(ilev) < 1.e-30)
232      cnew(:) = 0.
233   end where
234
235!  curvature of the solution and weighted error norm for O
236!
237!  cold : concentrations at time t - dt
238!  c    : concentrations at time t
239!  cnew : concentrations guessed at time t + dt
240
241!  ratio = dtold/dtnew
242!  here ratio = 1 since the timestep for the guess is equal to the previous timestep
243
244   ratio = 1.
245
246   curv = 2.*(ratio*cnew(i_o) - (1. + ratio)*c(ilev,i_o) + cold(i_o)) &
247          /(ratio*(1. + ratio))
248   e1 = (curv/(cnew(i_o) + cnew(i_o3)))*100.
249   e1 = abs(e1)
250
251!  curvature of the solution and weighted error norm for HO2
252
253   curv = 2.*(ratio*cnew(i_ho2) - (1. + ratio)*c(ilev,i_ho2) + cold(i_ho2)) &
254          /(ratio*(1. + ratio))
255   e2 = (curv/(cnew(i_h) + cnew(i_oh) + cnew(i_ho2) + 2.*cnew(i_h2o2)))*100.
256   e2 = abs(e2)
257
258!  curvature of the solution and weighted error norm for NO2
259
260!  curv = 2.*(ratio*cnew(i_no2) - (1. + ratio)*c(ilev,i_no2) + cold(i_no2)) &
261!         /(ratio*(1. + ratio))
262!  e3 = (curv/(cnew(i_no) + cnew(i_no2)))*100.
263!  e3 = abs(e3)
264   e3 = 0.
265
266!  timestep correction
267
268   e = max(e1, e2)
269   e = max(e,  e3)
270   e = max(e, 0.1)
271   dtg = max(0.001, min(2.5,0.8/sqrt(e)))
272!  dtg = 1.                                     ! uncomment this line to turn off the variable timestep scheme
273   dt_corrected = max(dt_min,dtg*dt_guess)      ! minimal timestep: dt_min
274   dt_corrected = min(dt_corrected,ctimestep)   ! maximal timestep: ctimestep
275
276   if (time + dt_corrected > ptimestep) then
277      dt_corrected = ptimestep - time
278   end if
279
280   if (dt_corrected /= dt_guess) then  ! the timestep has been modified
281
282!  form the matrix identity + mat*dt_guess
283
284      mat(:,:) = mat1(:,:)*dt_corrected
285      do iesp = 1,nesp
286         mat(iesp,iesp) = 1. + mat(iesp,iesp)
287      end do
288
289!  solve linear system
290     
291      cnew(:) = c(ilev,:)
292
293#ifdef LAPACK
294      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
295#else
296   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
297   stop
298#endif
299   end if
300
301!  eliminate small values
302
303   where (cnew(:)/dens(ilev) < 1.e-30)
304      cnew(:) = 0.
305   end where
306
307   cold(:)   = c(ilev,:)
308   c(ilev,:) = cnew(:)
309   cnew(:)   = 0.
310
311!  increment internal time
312
313   time = time + dt_corrected
314   dt_guess = dt_corrected     ! first-guess timestep for next iteration
315
316!  if (dt_corrected == dt_min) then
317!     print *, "ilev = ", ilev, " sza = ", sza, " iter = ", iter(ilev)
318!  end if
319
320!  do iesp = 1,nesp
321!     if (c(ilev,iesp) < 0.) then
322!        print *, "sza = ", sza, " species ", iesp, " = ", c(ilev,iesp)/dens(ilev)
323!        print *, "ilev = ", ilev, " iter = ", iter(ilev), " dt_corrected = ", dt_corrected
324!        print *, "e_o = ", e1, " e_ho2 = ", e2, " e_no2 = ", e3
325!     end if
326!  end do
327
328   end do ! while (time < ptimestep)
329
330!  if (iter(ilev) >= 80) then
331!     print *, "*** iter = ", iter(ilev), " ilev = ", ilev, " sza = ", sza
332!     print *, " dt_corrected = ", dt_corrected
333!     stop
334!  end if
335
336!  do iesp = 1,nesp
337!     if (c(ilev,iesp) < 0.) then
338!        print *, "final value is negative!!"
339!        print *, "sza = ", sza, " species ", iesp, " = ", c(ilev,iesp)/dens(ilev)
340!        print *, "ilev = ", ilev, " iter = ", iter(ilev), " dt_corrected = ", dt_corrected
341!        print *, "e_o = ", e1, " e_ho2 = ", e2, " e_no2 = ", e3
342!        stop
343!     end if
344!  end do
345
346!  save cold and dt for next timestep
347
348!  cold3d(ig,ilev,:) = cold(:)
349
350end do ! ilev
351
352!===================================================================
353!     save chemical species for the gcm       
354!===================================================================
355
356call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
357               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
358               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
359               i_n, i_n2d, i_no, i_no2, i_n2, dens, c)
360
361contains
362
363!======================================================================
364
365 subroutine reactionrates(nlayer,                               &
366                          lswitch, dens, co2, o2, press, t,     &
367                          hetero_dust, hetero_ice,              &
368                          surfdust1d, surfice1d,                &
369                          v_phot, v_3, v_4)
370 
371!================================================================
372! compute reaction rates                                        !
373!----------------------------------------------------------------
374! reaction               type                array              !
375!----------------------------------------------------------------
376! A + B    --> C + D     bimolecular         v_4                !
377! A + A    --> B + C     quadratic           v_3                !
378! A + C    --> B + C     quenching           v_phot             !
379! A + ice  --> B + C     heterogeneous       v_phot             !
380!================================================================
381
382use comcstfi_h
383
384implicit none
385
386#include "chimiedata.h"
387
388!----------------------------------------------------------------------
389!     input
390!----------------------------------------------------------------------
391
392integer, intent(in)     :: nlayer            ! number of atmospheric layers
393integer                 :: lswitch           ! interface level between lower
394                                             ! atmosphere and thermosphere chemistries
395real, dimension(nlayer) :: dens              ! total number density (molecule.cm-3)
396real, dimension(nlayer) :: press             ! pressure (hPa)
397real, dimension(nlayer) :: t                 ! temperature (K)
398real, dimension(nlayer) :: surfdust1d        ! dust surface area (cm2.cm-3)
399real, dimension(nlayer) :: surfice1d         ! ice surface area (cm2.cm-3)
400real (kind = 8), dimension(nlayer) :: co2    ! co2 number density (molecule.cm-3)
401real (kind = 8), dimension(nlayer) :: o2     ! o2 number density (molecule.cm-3)
402logical :: hetero_dust, hetero_ice           ! switches for heterogeneous chemistry
403
404!----------------------------------------------------------------------
405!     output
406!----------------------------------------------------------------------
407
408real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
409real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
410real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
411
412!----------------------------------------------------------------------
413!     local
414!----------------------------------------------------------------------
415
416integer          :: ilev
417integer          :: nb_phot, nb_reaction_3, nb_reaction_4
418real :: ak0, ak1, xpo, rate
419real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
420real, dimension(nlayer) :: deq
421real, dimension(nlayer) :: a001, a002, a003,                           &
422                             b001, b002, b003, b004, b005, b006, b007,   &
423                             b008, b009,                                 &
424                             c001, c002, c003, c004, c005, c006, c007,   &
425                             c008, c009, c010, c011, c012, c013, c014,   &
426                             c015, c016, c017, c018,                     &
427                             d001, d002, d003, d004, d005, d006, d007,   &
428                             d008, d009,                                 &
429                             e001, e002,                                 &
430                             h001, h002, h003, h004, h005
431
432!----------------------------------------------------------------------
433!     initialisation
434!----------------------------------------------------------------------
435
436      nb_phot       = 11       ! jmars.20140930 reduit de 13 a 11
437      nb_reaction_3 = 0
438      nb_reaction_4 = 0
439
440!----------------------------------------------------------------------
441!     oxygen reactions
442!----------------------------------------------------------------------
443
444!---  a001: o + o2 + co2 -> o3 + co2
445
446!     jpl 2003
447!
448!     co2 efficiency as a third body (2.075)
449!     from sehested et al., j. geophys. res., 100, 1995.
450
451      a001(:) = 2.075*6.0e-34*(t(:)/300.)**(-2.4)*dens(:)
452
453      nb_reaction_4 = nb_reaction_4 + 1
454      v_4(:,nb_reaction_4) = a001(:)
455
456!---  a002: o + o + co2 -> o2 + co2
457
458!     Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
459
460!     a002(:) = 2.5*5.2e-35*exp(900./t(:))*dens(:)
461
462!     Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
463
464!     a002(:) = 1.2e-32*(300./t(:))**(2.0)*dens(:)  ! yung expression
465      a002(:) = 2.5*9.46e-34*exp(485./t(:))*dens(:) ! nist expression
466
467      nb_reaction_3 = nb_reaction_3 + 1
468      v_3(:,nb_reaction_3) = a002(:)
469
470!---  a003: o + o3 -> o2 + o2
471
472!     jpl 2003
473
474      a003(:) = 8.0e-12*exp(-2060./t(:))
475
476      nb_reaction_4 = nb_reaction_4 + 1
477      v_4(:,nb_reaction_4) = a003(:)
478
479!----------------------------------------------------------------------
480!     o(1d) reactions
481!----------------------------------------------------------------------
482
483!---  b001: o(1d) + co2  -> o + co2
484
485!     jpl 2006
486
487      b001(:) = 7.5e-11*exp(115./t(:))
488   
489      nb_phot = nb_phot + 1
490      v_phot(:,nb_phot) = b001(:)*co2(:)
491
492!---  b002: o(1d) + h2o  -> oh + oh
493
494!     jpl 2006
495 
496      b002(:) = 1.63e-10*exp(60./t(:))
497
498      nb_reaction_4 = nb_reaction_4 + 1
499      v_4(:,nb_reaction_4) = b002(:)
500
501!---  b003: o(1d) + h2  -> oh + h
502
503!     jpl 2011
504
505      b003(:) = 1.2e-10
506
507      nb_reaction_4 = nb_reaction_4 + 1
508      v_4(:,nb_reaction_4) = b003(:)
509
510!---  b004: o(1d) + o2  -> o + o2
511
512!     jpl 2006
513
514      b004(:) = 3.3e-11*exp(55./t(:))
515
516      nb_phot = nb_phot + 1
517      v_phot(:,nb_phot) = b004(:)*o2(:)
518   
519!---  b005: o(1d) + o3  -> o2 + o2
520
521!     jpl 2003
522
523      b005(:) = 1.2e-10
524
525      nb_reaction_4 = nb_reaction_4 + 1
526      v_4(:,nb_reaction_4) = b005(:)
527   
528!---  b006: o(1d) + o3  -> o2 + o + o
529
530!     jpl 2003
531
532      b006(:) = 1.2e-10
533
534      nb_reaction_4 = nb_reaction_4 + 1
535      v_4(:,nb_reaction_4) = b006(:)
536   
537!---  b007: o(1d) + ch4 -> ch3 + oh
538
539!     jpl 2003
540
541      b007(:) = 1.5e-10*0.75
542
543!---  b008: o(1d) + ch4 -> ch3o + h
544
545!     jpl 2003
546
547      b008(:) = 1.5e-10*0.20
548!
549!---  b009: o(1d) + ch4 -> ch2o + h2
550
551!     jpl 2003
552
553      b009(:) = 1.5e-10*0.05
554
555!----------------------------------------------------------------------
556!     hydrogen reactions
557!----------------------------------------------------------------------
558
559!---  c001: o + ho2 -> oh + o2
560
561!     jpl 2003
562
563      c001(:) = 3.0e-11*exp(200./t(:))
564
565      nb_reaction_4 = nb_reaction_4 + 1
566      v_4(:,nb_reaction_4) = c001(:)
567
568!---  c002: o + oh -> o2 + h
569
570!     jpl 2011
571
572      c002(:) = 1.8e-11*exp(180./t(:))
573
574!     robertson and smith, j. chem. phys. a 110, 6673, 2006
575
576!     c002(:) = 11.2e-11*t(:)**(-0.32)*exp(177./t(:))
577
578      nb_reaction_4 = nb_reaction_4 + 1
579      v_4(:,nb_reaction_4) = c002(:)
580
581!---  c003: h + o3 -> oh + o2
582
583!     jpl 2003
584
585      c003(:) = 1.4e-10*exp(-470./t(:))
586
587      nb_reaction_4 = nb_reaction_4 + 1
588      v_4(:,nb_reaction_4) = c003(:)
589
590!---  c004: h + ho2 -> oh + oh
591
592!     jpl 2006
593
594      c004(:) = 7.2e-11
595
596      nb_reaction_4 = nb_reaction_4 + 1
597      v_4(:,nb_reaction_4) = c004(:)
598
599!---  c005: h + ho2 -> h2 + o2
600
601!     jpl 2006
602
603      c005(:) = 6.9e-12
604
605      nb_reaction_4 = nb_reaction_4 + 1
606      v_4(:,nb_reaction_4) = c005(:)
607
608!---  c006: h + ho2 -> h2o + o
609
610!     jpl 2006
611
612      c006(:) = 1.6e-12
613
614      nb_reaction_4 = nb_reaction_4 + 1
615      v_4(:,nb_reaction_4) = c006(:)
616
617!---  c007: oh + ho2 -> h2o + o2
618
619!     jpl 2003
620
621!     canty et al., grl, 2006 suggest to increase this rate
622!     by 20%. not done here.
623
624      c007(:) = 4.8e-11*exp(250./t(:))
625
626      nb_reaction_4 = nb_reaction_4 + 1
627      v_4(:,nb_reaction_4) = c007(:)
628
629!---  c008: ho2 + ho2 -> h2o2 + o2
630
631!     jpl 2006
632
633!     c008(:) = 3.5e-13*exp(430./t(:))
634
635!     christensen et al., grl, 13, 2002
636
637      c008(:) = 1.5e-12*exp(19./t(:))
638
639      nb_reaction_3 = nb_reaction_3 + 1
640      v_3(:,nb_reaction_3) = c008(:)
641
642!---  c009: oh + h2o2 -> h2o + ho2
643
644!     jpl 2006
645
646      c009(:) = 1.8e-12
647
648      nb_reaction_4 = nb_reaction_4 + 1
649      v_4(:,nb_reaction_4) = c009(:)
650
651!---  c010: oh + h2 -> h2o + h
652
653!     jpl 2006
654
655      c010(:) = 2.8e-12*exp(-1800./t(:))
656
657      nb_reaction_4 = nb_reaction_4 + 1
658      v_4(:,nb_reaction_4) = c010(:)
659
660!---  c011: h + o2 + co2 -> ho2 + co2
661
662!     jpl 2011
663
664      do ilev = 1,lswitch-1
665         ak0 = 2.5*4.4e-32*(t(ilev)/300.)**(-1.3)
666         ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
667
668         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
669         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
670         c011(ilev) = rate*0.6**xpo
671      end do
672
673      nb_reaction_4 = nb_reaction_4 + 1
674      v_4(:,nb_reaction_4) = c011(:)
675
676!---  c012: o + h2o2 -> oh + ho2
677
678!     jpl 2003
679
680      c012(:) = 1.4e-12*exp(-2000./t(:))
681
682      nb_reaction_4 = nb_reaction_4 + 1
683      v_4(:,nb_reaction_4) = c012(:)
684
685!---  c013: oh + oh -> h2o + o
686
687!     jpl 2006
688
689      c013(:) = 1.8e-12
690
691      nb_reaction_3 = nb_reaction_3 + 1
692      v_3(:,nb_reaction_3) = c013(:)
693
694!---  c014: oh + o3 -> ho2 + o2
695
696!     jpl 2003
697
698      c014(:) = 1.7e-12*exp(-940./t(:))
699
700      nb_reaction_4 = nb_reaction_4 + 1
701      v_4(:,nb_reaction_4) = c014(:)
702
703!---  c015: ho2 + o3 -> oh + o2 + o2
704
705!     jpl 2003
706
707      c015(:) = 1.0e-14*exp(-490./t(:))
708
709      nb_reaction_4 = nb_reaction_4 + 1
710      v_4(:,nb_reaction_4) = c015(:)
711
712!---  c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
713
714!     jpl 2011
715
716      c016(:) = 2.5*2.1e-33*exp(920./t(:))*dens(:)
717
718      nb_reaction_3 = nb_reaction_3 + 1
719      v_3(:,nb_reaction_3) = c016(:)
720
721!---  c017: oh + oh + co2 -> h2o2 + co2
722
723!     jpl 2003
724
725      do ilev = 1,lswitch-1
726         ak0 = 2.5*6.9e-31*(t(ilev)/300.)**(-1.0)
727         ak1 = 2.6e-11*(t(ilev)/300.)**(0.0)
728
729         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
730         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
731         c017(ilev) = rate*0.6**xpo
732      end do
733
734      nb_reaction_3 = nb_reaction_3 + 1
735      v_3(:,nb_reaction_3) = c017(:)
736
737!---  c018: h + h + co2 -> h2 + co2
738
739!     baulch et al., 2005
740
741      c018(:) = 2.5*1.8e-30*(t(:)**(-1.0))*dens(:)
742
743      nb_reaction_3 = nb_reaction_3 + 1
744      v_3(:,nb_reaction_3) = c018(:)
745
746!----------------------------------------------------------------------
747!     nitrogen reactions
748!----------------------------------------------------------------------
749
750!---  d001: no2 + o -> no + o2
751
752!     jpl 2006
753
754      d001(:) = 5.1e-12*exp(210./t(:))
755
756      nb_reaction_4 = nb_reaction_4 + 1
757      v_4(:,nb_reaction_4) = d001(:)
758
759!---  d002: no + o3 -> no2 + o2
760
761!     jpl 2006
762
763      d002(:) = 3.0e-12*exp(-1500./t(:))
764
765      nb_reaction_4 = nb_reaction_4 + 1
766      v_4(:,nb_reaction_4) = d002(:)
767
768!---  d003: no + ho2 -> no2 + oh
769
770!     jpl 2011
771
772      d003(:) = 3.3e-12*exp(270./t(:))
773
774      nb_reaction_4 = nb_reaction_4 + 1
775      v_4(:,nb_reaction_4) = d003(:)
776
777!---  d004: n + no -> n2 + o
778
779!     jpl 2011
780
781      d004(:) = 2.1e-11*exp(100./t(:))
782
783      nb_reaction_4 = nb_reaction_4 + 1
784      v_4(:,nb_reaction_4) = d004(:)
785
786!---  d005: n + o2 -> no + o
787
788!     jpl 2011
789
790      d005(:) = 1.5e-11*exp(-3600./t(:))
791
792      nb_reaction_4 = nb_reaction_4 + 1
793      v_4(:,nb_reaction_4) = d005(:)
794
795!---  d006: no2 + h -> no + oh
796
797!     jpl 2011
798
799      d006(:) = 4.0e-10*exp(-340./t(:))
800
801      nb_reaction_4 = nb_reaction_4 + 1
802      v_4(:,nb_reaction_4) = d006(:)
803
804!---  d007: n + o -> no
805 
806      d007(:) = 2.8e-17*(300./t(:))**0.5
807
808      nb_reaction_4 = nb_reaction_4 + 1
809      v_4(:,nb_reaction_4) = d007(:)
810
811!---  d008: n + ho2 -> no + oh
812
813!     brune et al., j. chem. phys., 87, 1983
814
815      d008(:) = 2.19e-11
816
817      nb_reaction_4 = nb_reaction_4 + 1
818      v_4(:,nb_reaction_4) = d008(:)
819
820!---  d009: n + oh -> no + h
821
822!     atkinson et al., j. phys. chem. ref. data, 18, 881, 1989
823
824      d009(:) = 3.8e-11*exp(85./t(:))
825
826      nb_reaction_4 = nb_reaction_4 + 1
827      v_4(:,nb_reaction_4) = d009(:)
828
829!----------------------------------------------------------------------
830!     carbon reactions
831!----------------------------------------------------------------------
832
833!---  e001: oh + co -> co2 + h
834
835!     jpl 2003
836
837!     e001(:) = 1.5e-13*(1 + 0.6*press(:)/1013.)
838
839!     mccabe et al., grl, 28, 3135, 2001
840
841!     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
842
843!     jpl 2006
844
845!     ak0 = 1.5e-13*(t(:)/300.)**(0.6)
846!     ak1 = 2.1e-9*(t(:)/300.)**(6.1)
847!     rate1 = ak0/(1. + ak0/(ak1/dens(:)))
848!     xpo1 = 1./(1. + alog10(ak0/(ak1/dens(:)))**2)
849
850!     ak0 = 5.9e-33*(t(:)/300.)**(-1.4)
851!     ak1 = 1.1e-12*(t(:)/300.)**(1.3)
852!     rate2 = (ak0*dens(:))/(1. + ak0*dens(:)/ak1)
853!     xpo2 = 1./(1. + alog10((ak0*dens(:))/ak1)**2)
854
855!     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
856
857!     joshi et al., 2006
858
859      do ilev = 1,lswitch-1
860         k1a0 = 1.34*2.5*dens(ilev)                                  &
861               *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
862               + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
863         k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
864              + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
865         k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
866                + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
867         x = k1a0/(k1ainf - k1b0)
868         y = k1b0/(k1ainf - k1b0)
869         fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
870            + exp(-t(ilev)/255.)
871         fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
872         k1a = k1a0*((1. + y)/(1. + x))*fx
873         k1b = k1b0*(1./(1.+x))*fx
874
875         e001(ilev) = k1a + k1b
876      end do
877
878      nb_reaction_4 = nb_reaction_4 + 1
879      v_4(:,nb_reaction_4) = e001(:)
880
881!---  e002: o + co + m -> co2 + m
882
883!     tsang and hampson, 1986.
884
885      e002(:) = 2.5*6.5e-33*exp(-2184./t(:))*dens(:)
886
887      nb_reaction_4 = nb_reaction_4 + 1
888      v_4(:,nb_reaction_4) = e002(:)
889
890!----------------------------------------------------------------------
891!     heterogeneous chemistry
892!----------------------------------------------------------------------
893
894      if (hetero_ice) then
895
896!        k = (surface*v*gamma)/4 (s-1)
897!        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
898 
899!---     h001: ho2 + ice -> products
900 
901!        cooper and abbatt, 1996: gamma = 0.025
902     
903         gam = 0.025
904         h001(:) = surfice1d(:)       &
905                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
906 
907!        h002: oh + ice -> products
908 
909!        cooper and abbatt, 1996: gamma = 0.03
910 
911         gam = 0.03
912         h002(:) = surfice1d(:)       &
913                   *100.*sqrt(8.*8.31*t(:)/(17.e-3*pi))*gam/4.
914
915!---     h003: h2o2 + ice -> products
916 
917!        gamma = 0.    test value
918 
919         gam = 0.
920         h003(:) = surfice1d(:)       &
921                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
922      else
923         h001(:) = 0.
924         h002(:) = 0.
925         h003(:) = 0.
926      end if
927
928      nb_phot = nb_phot + 1
929      v_phot(:,nb_phot) = h001(:)
930
931      nb_phot = nb_phot + 1
932      v_phot(:,nb_phot) = h002(:)
933
934      nb_phot = nb_phot + 1
935      v_phot(:,nb_phot) = h003(:)
936
937      if (hetero_dust) then
938 
939!---     h004: ho2 + dust -> products
940 
941!        jacob, 2000: gamma = 0.2
942!        see dereus et al., atm. chem. phys., 2005
943 
944         gam = 0.2
945         h004(:) = surfdust1d(:)  &
946                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
947 
948!---     h005: h2o2 + dust -> products
949 
950!        gamma = 5.e-4
951!        see dereus et al., atm. chem. phys., 2005
952 
953         gam = 5.e-4
954         h005(:) = surfdust1d(:)  &
955                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
956      else
957         h004(:) = 0.
958         h005(:) = 0.
959      end if
960 
961      nb_phot = nb_phot + 1
962      v_phot(:,nb_phot) = h004(:)
963
964      nb_phot = nb_phot + 1
965      v_phot(:,nb_phot) = h005(:)
966
967end subroutine reactionrates
968
969!======================================================================
970
971 subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
972
973!======================================================================
974! filling of the jacobian matrix
975!======================================================================
976
977use types_asis
978
979implicit none
980
981#include "chimiedata.h"
982
983! input
984
985integer             :: ilev    ! level index
986integer             :: nesp    ! number of species in the chemistry
987integer, intent(in) :: nlayer  ! number of atmospheric layers
988
989real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
990real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
991real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
992real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
993
994! output
995
996real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
997real (kind = 8), dimension(nesp), intent(out)      :: prod, loss
998
999! local
1000
1001integer :: iesp
1002integer :: ind_phot_2,ind_phot_4,ind_phot_6
1003integer :: ind_3_2,ind_3_4,ind_3_6
1004integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
1005integer :: iphot,i3,i4
1006
1007real(kind = jprb) :: eps, eps_4  ! implicit/explicit coefficient
1008
1009! initialisations
1010
1011mat(:,:) = 0.
1012prod(:)  = 0.
1013loss(:)  = 0.
1014
1015! photodissociations
1016! or reactions a + c -> b + c
1017! or reactions a + ice -> b + c
1018
1019do iphot = 1,nb_phot_max
1020
1021  ind_phot_2 = indice_phot(iphot)%z2
1022  ind_phot_4 = indice_phot(iphot)%z4
1023  ind_phot_6 = indice_phot(iphot)%z6
1024
1025  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
1026  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
1027  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
1028
1029  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
1030  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
1031  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
1032
1033end do
1034
1035! reactions a + a -> b + c
1036
1037do i3 = 1,nb_reaction_3_max
1038
1039  ind_3_2 = indice_3(i3)%z2
1040  ind_3_4 = indice_3(i3)%z4
1041  ind_3_6 = indice_3(i3)%z6
1042
1043  mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
1044  mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)
1045  mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)
1046
1047  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
1048  prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
1049  prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
1050
1051end do
1052
1053! reactions a + b -> c + d
1054
1055eps = 1.d-10
1056
1057do i4 = 1,nb_reaction_4_max
1058
1059  ind_4_2 = indice_4(i4)%z2
1060  ind_4_4 = indice_4(i4)%z4
1061  ind_4_6 = indice_4(i4)%z6
1062  ind_4_8 = indice_4(i4)%z8
1063
1064  eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
1065  eps_4 = min(eps_4,1.0_jprb)
1066
1067  mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
1068  mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
1069  mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
1070  mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)   
1071  mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
1072  mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
1073  mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
1074  mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
1075
1076  loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
1077  loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
1078  prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
1079  prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
1080
1081end do
1082
1083end subroutine fill_matrix
1084
1085!================================================================
1086
1087 subroutine indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
1088                   i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
1089                   i_n, i_n2d, i_no, i_no2, i_n2)
1090
1091!================================================================
1092! set the "indice" arrays used to fill the jacobian matrix      !
1093!----------------------------------------------------------------
1094! reaction                                   type               !
1095!----------------------------------------------------------------
1096! A + hv   --> B + C     photolysis          indice_phot        !
1097! A + B    --> C + D     bimolecular         indice_4           !
1098! A + A    --> B + C     quadratic           indice_3           !
1099! A + C    --> B + C     quenching           indice_phot        !
1100! A + ice  --> B + C     heterogeneous       indice_phot        !
1101!================================================================
1102
1103use types_asis
1104
1105implicit none
1106
1107#include "chimiedata.h"
1108
1109! input
1110
1111integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
1112           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
1113           i_n, i_n2d, i_no, i_no2, i_n2
1114
1115! local
1116
1117integer :: nb_phot, nb_reaction_3, nb_reaction_4
1118integer :: i_dummy
1119
1120allocate (indice_phot(nb_phot_max))
1121allocate (indice_3(nb_reaction_3_max))
1122allocate (indice_4(nb_reaction_4_max))
1123
1124i_dummy = 1
1125
1126nb_phot       = 0
1127nb_reaction_3 = 0
1128nb_reaction_4 = 0
1129
1130!===========================================================
1131!      O2 + hv -> O + O
1132!===========================================================
1133
1134nb_phot = nb_phot + 1
1135
1136indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy)
1137
1138!===========================================================
1139!      O2 + hv -> O + O(1D)
1140!===========================================================
1141
1142nb_phot = nb_phot + 1
1143
1144indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d)
1145
1146!===========================================================
1147!      CO2 + hv -> CO + O
1148!===========================================================
1149
1150nb_phot = nb_phot + 1
1151
1152indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o)
1153
1154!===========================================================
1155!      CO2 + hv -> CO + O(1D)
1156!===========================================================
1157
1158nb_phot = nb_phot + 1
1159
1160indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d)
1161
1162!===========================================================
1163!      O3 + hv -> O2 + O(1D)
1164!===========================================================
1165
1166nb_phot = nb_phot + 1
1167
1168indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d)
1169
1170!===========================================================
1171!      O3 + hv -> O2 + O
1172!===========================================================
1173
1174nb_phot = nb_phot + 1
1175
1176indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o)
1177
1178!===========================================================
1179!      H2O + hv -> H + OH
1180!===========================================================
1181
1182nb_phot = nb_phot + 1
1183
1184indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh)
1185
1186!===========================================================
1187!      H2O2 + hv -> OH + OH
1188!===========================================================
1189
1190nb_phot = nb_phot + 1
1191
1192indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy)
1193
1194!===========================================================
1195!      HO2 + hv -> OH + O
1196!===========================================================
1197
1198nb_phot = nb_phot + 1
1199
1200indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o)
1201
1202!===========================================================
1203!      NO + hv -> N + O
1204!===========================================================
1205
1206nb_phot = nb_phot + 1
1207
1208indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o)
1209
1210!===========================================================
1211!      NO2 + hv -> NO + O
1212!===========================================================
1213
1214nb_phot = nb_phot + 1
1215
1216indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o)
1217
1218!===========================================================
1219!      a001 : O + O2 + CO2 -> O3 + CO2
1220!===========================================================
1221
1222nb_reaction_4 = nb_reaction_4 + 1
1223
1224indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
1225
1226!===========================================================
1227!      a002 : O + O + CO2 -> O2 + CO2
1228!===========================================================
1229
1230nb_reaction_3 = nb_reaction_3 + 1
1231
1232indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
1233
1234!===========================================================
1235!      a003 : O + O3 -> O2 + O2
1236!===========================================================
1237
1238nb_reaction_4 = nb_reaction_4 + 1
1239
1240indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
1241
1242!===========================================================
1243!      b001 : O(1D) + CO2 -> O + CO2
1244!===========================================================
1245
1246nb_phot = nb_phot + 1
1247
1248indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
1249
1250!===========================================================
1251!      b002 : O(1D) + H2O -> OH + OH
1252!===========================================================
1253
1254nb_reaction_4 = nb_reaction_4 + 1
1255
1256indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
1257
1258!===========================================================
1259!      b003 : O(1D) + H2 -> OH + H
1260!===========================================================
1261
1262nb_reaction_4 = nb_reaction_4 + 1
1263
1264indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
1265
1266!===========================================================
1267!      b004 : O(1D) + O2 -> O + O2
1268!===========================================================
1269
1270nb_phot = nb_phot + 1
1271
1272indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
1273
1274!===========================================================
1275!      b005 : O(1D) + O3 -> O2 + O2
1276!===========================================================
1277
1278nb_reaction_4 = nb_reaction_4 + 1
1279
1280indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
1281
1282!===========================================================
1283!      b006 : O(1D) + O3 -> O2 + O + O
1284!===========================================================
1285
1286nb_reaction_4 = nb_reaction_4 + 1
1287
1288indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
1289
1290!===========================================================
1291!      c001 : O + HO2 -> OH + O2
1292!===========================================================
1293
1294nb_reaction_4 = nb_reaction_4 + 1
1295
1296indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
1297
1298!===========================================================
1299!      c002 : O + OH -> O2 + H
1300!===========================================================
1301
1302nb_reaction_4 = nb_reaction_4 + 1
1303
1304indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
1305
1306!===========================================================
1307!      c003 : H + O3 -> OH + O2
1308!===========================================================
1309
1310nb_reaction_4 = nb_reaction_4 + 1
1311
1312indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
1313
1314!===========================================================
1315!      c004 : H + HO2 -> OH + OH
1316!===========================================================
1317
1318nb_reaction_4 = nb_reaction_4 + 1
1319
1320indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
1321
1322!===========================================================
1323!      c005 : H + HO2 -> H2 + O2
1324!===========================================================
1325
1326nb_reaction_4 = nb_reaction_4 + 1
1327
1328indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
1329
1330!===========================================================
1331!      c006 : H + HO2 -> H2O + O
1332!===========================================================
1333
1334nb_reaction_4 = nb_reaction_4 + 1
1335
1336indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
1337
1338!===========================================================
1339!      c007 : OH + HO2 -> H2O + O2
1340!===========================================================
1341
1342nb_reaction_4 = nb_reaction_4 + 1
1343
1344indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
1345
1346!===========================================================
1347!      c008 : HO2 + HO2 -> H2O2 + O2
1348!===========================================================
1349
1350nb_reaction_3 = nb_reaction_3 + 1
1351
1352indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
1353
1354!===========================================================
1355!      c009 : OH + H2O2 -> H2O + HO2
1356!===========================================================
1357
1358nb_reaction_4 = nb_reaction_4 + 1
1359
1360indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
1361
1362!===========================================================
1363!      c010 : OH + H2 -> H2O + H
1364!===========================================================
1365
1366nb_reaction_4 = nb_reaction_4 + 1
1367
1368indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
1369
1370!===========================================================
1371!      c011 : H + O2 + CO2 -> HO2 + CO2
1372!===========================================================
1373
1374nb_reaction_4 = nb_reaction_4 + 1
1375
1376indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
1377
1378!===========================================================
1379!      c012 : O + H2O2 -> OH + HO2
1380!===========================================================
1381
1382nb_reaction_4 = nb_reaction_4 + 1
1383
1384indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
1385
1386!===========================================================
1387!      c013 : OH + OH -> H2O + O
1388!===========================================================
1389
1390nb_reaction_3 = nb_reaction_3 + 1
1391
1392indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
1393
1394!===========================================================
1395!      c014 : OH + O3 -> HO2 + O2
1396!===========================================================
1397
1398nb_reaction_4 = nb_reaction_4 + 1
1399
1400indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
1401
1402!===========================================================
1403!      c015 : HO2 + O3 -> OH + O2 + O2
1404!===========================================================
1405
1406nb_reaction_4 = nb_reaction_4 + 1
1407
1408indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
1409
1410!===========================================================
1411!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
1412!===========================================================
1413
1414nb_reaction_3 = nb_reaction_3 + 1
1415
1416indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
1417
1418!===========================================================
1419!      c017 : OH + OH + CO2 -> H2O2 + CO2
1420!===========================================================
1421
1422nb_reaction_3 = nb_reaction_3 + 1
1423
1424indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
1425
1426!===========================================================
1427!      c018 : H + H + CO2 -> H2 + CO2
1428!===========================================================
1429
1430nb_reaction_3 = nb_reaction_3 + 1
1431
1432indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
1433
1434!===========================================================
1435!      d001 : NO2 + O -> NO + O2
1436!===========================================================
1437
1438nb_reaction_4 = nb_reaction_4 + 1
1439
1440indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
1441
1442!===========================================================
1443!      d002 : NO + O3 -> NO2 + O2
1444!===========================================================
1445
1446nb_reaction_4 = nb_reaction_4 + 1
1447
1448indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
1449
1450!===========================================================
1451!      d003 : NO + HO2 -> NO2 + OH
1452!===========================================================
1453
1454nb_reaction_4 = nb_reaction_4 + 1
1455
1456indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
1457
1458!===========================================================
1459!      d004 : N + NO -> N2 + O
1460!===========================================================
1461
1462nb_reaction_4 = nb_reaction_4 + 1
1463
1464indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
1465
1466!===========================================================
1467!      d005 : N + O2 -> NO + O
1468!===========================================================
1469
1470nb_reaction_4 = nb_reaction_4 + 1
1471
1472indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
1473
1474!===========================================================
1475!      d006 : NO2 + H -> NO + OH
1476!===========================================================
1477
1478nb_reaction_4 = nb_reaction_4 + 1
1479
1480indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
1481
1482!===========================================================
1483!      d007 : N + O -> NO
1484!===========================================================
1485
1486nb_reaction_4 = nb_reaction_4 + 1
1487
1488indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
1489
1490!===========================================================
1491!      d008 : N + HO2 -> NO + OH
1492!===========================================================
1493
1494nb_reaction_4 = nb_reaction_4 + 1
1495
1496indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
1497
1498!===========================================================
1499!      d009 : N + OH -> NO + H
1500!===========================================================
1501
1502nb_reaction_4 = nb_reaction_4 + 1
1503
1504indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
1505
1506!===========================================================
1507!      e001 : CO + OH -> CO2 + H
1508!===========================================================
1509
1510nb_reaction_4 = nb_reaction_4 + 1
1511
1512indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
1513
1514!===========================================================
1515!      e002 : CO + O + M -> CO2 + M
1516!===========================================================
1517
1518nb_reaction_4 = nb_reaction_4 + 1
1519
1520indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
1521
1522!===========================================================
1523!      h001: HO2 + ice -> products
1524!            treated as
1525!            HO2 -> 0.5 H2O + 0.75 O2
1526!===========================================================
1527
1528nb_phot = nb_phot + 1
1529
1530indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
1531
1532!===========================================================
1533!      h002: OH + ice -> products
1534!            treated as
1535!            OH -> 0.5 H2O + 0.25 O2
1536!===========================================================
1537
1538nb_phot = nb_phot + 1
1539
1540indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2)
1541
1542!===========================================================
1543!      h003: H2O2 + ice -> products
1544!            treated as
1545!            H2O2 -> H2O + 0.5 O2
1546!===========================================================
1547
1548nb_phot = nb_phot + 1
1549
1550indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
1551
1552!===========================================================
1553!      h004: HO2 + dust -> products
1554!            treated as
1555!            HO2 -> 0.5 H2O + 0.75 O2
1556!===========================================================
1557
1558nb_phot = nb_phot + 1
1559
1560indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2)
1561
1562!===========================================================
1563!      h005: H2O2 + dust -> products
1564!            treated as
1565!            H2O2 -> H2O + 0.5 O2
1566!===========================================================
1567
1568nb_phot = nb_phot + 1
1569
1570indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
1571
1572!===========================================================
1573!  check dimensions
1574!===========================================================
1575
1576print*, 'nb_phot       = ', nb_phot
1577print*, 'nb_reaction_4 = ', nb_reaction_4
1578print*, 'nb_reaction_3 = ', nb_reaction_3
1579
1580if ((nb_phot /= nb_phot_max)             .or.  &
1581    (nb_reaction_3 /= nb_reaction_3_max) .or.  &
1582    (nb_reaction_4 /= nb_reaction_4_max)) then
1583   print*, 'wrong dimensions in indice'
1584   stop
1585end if 
1586
1587end subroutine indice
1588
1589!*****************************************************************
1590
1591      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp,         &
1592                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
1593                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
1594                           i_n, i_n2d, i_no, i_no2, i_n2,            &
1595                           dens, rm, c)
1596
1597!*****************************************************************
1598
1599      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
1600     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,  &
1601     &                      igcm_ho2, igcm_h2o2, igcm_h2o_vap,           &
1602     &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
1603
1604      implicit none
1605
1606#include "callkeys.h"
1607
1608!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1609!     input:
1610!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1611     
1612      integer, intent(in) :: nlayer ! number of atmospheric layers
1613      integer, intent(in) :: nq     ! number of tracers in the gcm
1614      integer :: nesp               ! number of species in the chemistry
1615      integer :: lswitch            ! interface level between chemistries
1616
1617      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,           &
1618                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                   &
1619                 i_n, i_n2d, i_no, i_no2, i_n2
1620
1621      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
1622      real :: dens(nlayer)          ! total number density (molecule.cm-3)
1623
1624!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1625!     output:
1626!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1627
1628      real, dimension(nlayer,nesp)            :: rm ! volume mixing ratios
1629      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3)
1630     
1631!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1632!     local:
1633!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1634
1635      integer      :: l, iesp
1636      logical,save :: firstcall = .true.
1637     
1638     
1639!     first call initializations
1640
1641      if (firstcall) then
1642
1643!       identify the indexes of the tracers we need
1644
1645         if (igcm_co2 == 0) then
1646            write(*,*) "gcmtochim: Error; no CO2 tracer !!!"
1647            stop
1648         endif
1649         if (igcm_co == 0) then
1650            write(*,*) "gcmtochim: Error; no CO tracer !!!"
1651            stop
1652         end if
1653         if (igcm_o == 0) then
1654            write(*,*) "gcmtochim: Error; no O tracer !!!"
1655            stop
1656         end if
1657         if (igcm_o1d == 0) then
1658            write(*,*) "gcmtochim: Error; no O1D tracer !!!"
1659            stop
1660         end if
1661         if (igcm_o2 == 0) then
1662            write(*,*) "gcmtochim: Error; no O2 tracer !!!"
1663            stop
1664         end if
1665         if (igcm_o3 == 0) then
1666            write(*,*) "gcmtochim: Error; no O3 tracer !!!"
1667            stop
1668         end if
1669         if (igcm_h == 0) then
1670            write(*,*) "gcmtochim: Error; no H tracer !!!"
1671            stop
1672         end if
1673         if (igcm_h2 == 0) then
1674            write(*,*) "gcmtochim: Error; no H2 tracer !!!"
1675            stop
1676         end if
1677         if (igcm_oh == 0) then
1678            write(*,*) "gcmtochim: Error; no OH tracer !!!"
1679            stop
1680         end if
1681         if (igcm_ho2 == 0) then
1682            write(*,*) "gcmtochim: Error; no HO2 tracer !!!"
1683            stop
1684         end if
1685         if (igcm_h2o2 == 0) then
1686            write(*,*) "gcmtochim: Error; no H2O2 tracer !!!"
1687            stop
1688         end if
1689         if (igcm_n == 0) then
1690            write(*,*) "gcmtochim: Error; no N tracer !!!"
1691            stop
1692         end if
1693         if (igcm_n2d == 0) then
1694            write(*,*) "gcmtochim: Error; no N2D tracer !!!"
1695            stop
1696         end if
1697         if (igcm_no == 0) then
1698            write(*,*) "gcmtochim: Error; no NO tracer !!!"
1699            stop
1700         end if
1701         if (igcm_no2 == 0) then
1702            write(*,*) "gcmtochim: Error; no NO2 tracer !!!"
1703            stop
1704         end if
1705         if (igcm_n2 == 0) then
1706            write(*,*) "gcmtochim: Error; no N2 tracer !!!"
1707            stop
1708         end if
1709         if (igcm_h2o_vap == 0) then
1710            write(*,*) "gcmtochim: Error; no water vapor tracer !!!"
1711            stop
1712         end if
1713         firstcall = .false.
1714      end if ! of if (firstcall)
1715
1716!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1717!     initialise mixing ratios
1718!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1719
1720      do l = 1,lswitch-1
1721         rm(l,i_co2)  = zycol(l, igcm_co2)
1722         rm(l,i_co)   = zycol(l, igcm_co)
1723         rm(l,i_o)    = zycol(l, igcm_o)
1724         rm(l,i_o1d)  = zycol(l, igcm_o1d)
1725         rm(l,i_o2)   = zycol(l, igcm_o2)
1726         rm(l,i_o3)   = zycol(l, igcm_o3)
1727         rm(l,i_h)    = zycol(l, igcm_h)
1728         rm(l,i_h2)   = zycol(l, igcm_h2)
1729         rm(l,i_oh)   = zycol(l, igcm_oh)
1730         rm(l,i_ho2)  = zycol(l, igcm_ho2)
1731         rm(l,i_h2o2) = zycol(l, igcm_h2o2)
1732         rm(l,i_h2o)  = zycol(l, igcm_h2o_vap)
1733         rm(l,i_n)    = zycol(l, igcm_n)
1734         rm(l,i_n2d)  = zycol(l, igcm_n2d)
1735         rm(l,i_no)   = zycol(l, igcm_no)
1736         rm(l,i_no2)  = zycol(l, igcm_no2)
1737         rm(l,i_n2)   = zycol(l, igcm_n2)
1738      end do
1739
1740      where (rm(:,:) < 1.e-30)
1741         rm(:,:) = 0.
1742      end where
1743
1744!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1745!     initialise number densities
1746!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1747
1748      do iesp = 1,nesp
1749         do l = 1,lswitch-1
1750            c(l,iesp) = rm(l,iesp)*dens(l)
1751         end do
1752      end do
1753
1754      end subroutine gcmtochim
1755
1756!*****************************************************************
1757 
1758      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp,         &
1759                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
1760                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
1761                           i_n, i_n2d, i_no, i_no2, i_n2, dens, c)
1762 
1763!*****************************************************************
1764 
1765      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,          &
1766                            igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,   &
1767                            igcm_ho2, igcm_h2o2, igcm_h2o_vap,            &
1768                            igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
1769
1770      implicit none
1771
1772#include "callkeys.h"
1773
1774!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1775!     inputs:
1776!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1777 
1778      integer, intent(in) :: nlayer  ! number of atmospheric layers
1779      integer, intent(in) :: nq      ! number of tracers in the gcm
1780      integer :: nesp                ! number of species in the chemistry
1781      integer :: lswitch             ! interface level between chemistries
1782      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
1783                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
1784                 i_n, i_n2d, i_no, i_no2, i_n2
1785
1786      real :: dens(nlayer)     ! total number density (molecule.cm-3)
1787      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3)
1788
1789!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1790!     output:
1791!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1792       
1793      real zycol(nlayer,nq)  ! volume mixing ratios in the gcm
1794
1795!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1796!     local:
1797!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1798
1799      integer l
1800     
1801!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1802!     save mixing ratios for the gcm
1803!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1804
1805      do l = 1,lswitch-1
1806         zycol(l, igcm_co2)     = c(l,i_co2)/dens(l)
1807         zycol(l, igcm_co)      = c(l,i_co)/dens(l)
1808         zycol(l, igcm_o)       = c(l,i_o)/dens(l)
1809         zycol(l, igcm_o1d)     = c(l,i_o1d)/dens(l)
1810         zycol(l, igcm_o2)      = c(l,i_o2)/dens(l)
1811         zycol(l, igcm_o3)      = c(l,i_o3)/dens(l)
1812         zycol(l, igcm_h)       = c(l,i_h)/dens(l) 
1813         zycol(l, igcm_h2)      = c(l,i_h2)/dens(l)
1814         zycol(l, igcm_oh)      = c(l,i_oh)/dens(l)
1815         zycol(l, igcm_ho2)     = c(l,i_ho2)/dens(l)
1816         zycol(l, igcm_h2o2)    = c(l,i_h2o2)/dens(l)
1817         zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
1818         zycol(l, igcm_n)       = c(l,i_n)/dens(l)
1819         zycol(l, igcm_n2d)     = c(l,i_n2d)/dens(l)
1820         zycol(l, igcm_no)      = c(l,i_no)/dens(l)
1821         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
1822         zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
1823      end do
1824
1825      end subroutine chimtogcm
1826
1827end subroutine photochemistry_asis
1828
Note: See TracBrowser for help on using the repository browser.