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

Last change on this file since 1574 was 1574, checked in by flefevre, 9 years ago

correction de type (real*8 -> real) pour variable densite

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