source: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90 @ 2031

Last change on this file since 2031 was 2031, checked in by flefevre, 6 years ago

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