Changeset 2170 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Oct 14, 2019, 11:00:51 PM (5 years ago)
Author:
flefevre
Message:
  • la chimie des ions est désormais optionnelle (FGG)
  • regroupement des options de la photochimie dans calchim_mod (FL)
  • mise à jour cinétique CO + OH, inactive pour le moment (FL)
Location:
trunk/LMDZ.MARS/libf/aeronomars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90

    r2164 r2170  
    2626      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
    2727      use comcstfi_h, only: pi
    28       use photolysis_mod, only: jonline, init_photolysis
     28      use photolysis_mod, only: init_photolysis, nphot
    2929      use iono_h, only: temp_elect
    3030
    … …  
    3838!  Prepare the call for the photochemical module, and send back the
    3939!  tendencies from photochemistry in the chemical species mass mixing ratios
    40 !
    41 !   Author:   Sebastien Lebonnois (08/11/2002)
    42 !   -------
    43 !    update 12/06/2003 for water ice clouds and compatibility with dust
    44 !    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
    45 !    update 03/05/2005 cosmetic changes (Franck Lefevre)
    46 !    update sept. 2008 identify tracers by their names (Ehouarn Millour)
    47 !    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
    48 !    update 16/03/2012 optimization (Franck Lefevre)
    49 !    update 11/12/2013 optimization (Franck Lefevre)
    5040!
    5141!   Arguments:
    … …  
    155145      integer :: ig_vl1
    156146
     147      integer :: nb_reaction_3_max     ! number of quadratic reactions
     148      integer :: nb_reaction_4_max     ! number of bimolecular reactions
     149      integer :: nquench               ! number of quenching + heterogeneous reactions
     150      integer :: nphotion              ! number of photoionizations
     151      integer :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
     152
     153
    157154      real    :: latvl1, lonvl1
    158155      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
    … …  
    164161      real    :: kb                    ! boltzmann constant
    165162
    166       logical,save :: firstcall = .true.
    167       logical,save :: depos = .false.  ! switch for dry deposition
     163      logical, save :: firstcall = .true.
     164      logical, save :: depos       ! switch for dry deposition
     165      logical, save :: ionchem     ! switch for ionospheric chemistry
     166      logical, save :: jonline     ! switch for online photodissociation rates or lookup table
     167      logical, save :: unichim     ! only one unified chemical scheme at all
     168                                   ! layers (default), or upper atmospheric
     169                                   ! scheme in the thermosphere
    168170
    169171!     for each column of atmosphere:
    … …  
    181183      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
    182184      real :: jh2o(nlayer)         ! Photodissociation rate H2O->H+OH (s-1)
    183       real :: em_no(nlayer)        !  NO nightglow emission rate
    184       real :: em_o2(nlayer)        !  O2 nightglow emission rate     
    185 
    186       integer :: iter(nlayer)      !  Number of chemical iterations
    187                                    !  within one physical timestep
     185      real :: em_no(nlayer)        ! NO nightglow emission rate
     186      real :: em_o2(nlayer)        ! O2 nightglow emission rate     
     187
     188      integer :: iter(nlayer)      ! Number of chemical iterations
     189                                   ! within one physical timestep
    188190
    189191!     for output:
    190192
    191193      logical :: output             ! to issue calls to writediagfi and stats
    192       parameter (output = .true.)
    193194      real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
    194195      real :: jh2o_3d(ngrid,nlayer)  ! Photodissociation rate H2O->H+OH (s-1)
    … …  
    198199                                    !  within one physical timestep
    199200
    200       logical :: unichim            ! only one, unified chemical scheme at all
    201                                     ! layers (default), or upper atmospheric
    202                                     ! scheme in the thermosphere?
     201!=======================================================================
     202!     main dashboard for the chemistry
     203!=======================================================================
     204
     205      unichim = .true.     ! true : unified chemistry ! false : separate models in lower and upper atmosphere
     206      jonline = .true.     ! true : on-line calculation of photodissociation rates ! false : lookup table
     207      ionchem = .false.    ! switch for ionospheric chemistry
     208      depos   = .false.    ! switch for dry deposition
     209      output  = .false.    ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc)
    203210
    204211!=======================================================================
    205212!     initialization of the chemistry (first call only)
    206213!=======================================================================
    207 
    208       !call only unified chemistry (default), or also upper atmospheric model
    209       !unichim = .true.   unified chemistry
    210       !unichim = .false.  2 different models
    211       unichim = .true.
    212214
    213215      if (firstcall) then
    … …  
    388390         i_o2plus = igcm_o2plus
    389391         if (i_o2plus == 0) then
    390             write(*,*) "calchim: Error, no O2+ tracer !!!"
    391             stop
     392            write(*,*) "calchim: no O2+ tracer "
     393            write(*,*) "Only neutral chemistry"
    392394         else
    393395            nbq = nbq + 1
    394396            niq(nbq) = i_o2plus
     397            ionchem = .true.
     398            write(*,*) "calchim: O2+ tracer found in traceur.def"
     399            write(*,*) "Ion chemistry included"
    395400         end if
    396401         i_co2plus = igcm_co2plus
    397          if (i_co2plus == 0) then
    398             write(*,*) "calchim: Error, no CO2+ tracer !!!"
    399             stop
    400          else
    401             nbq = nbq + 1
    402             niq(nbq) = i_co2plus
    403          end if
     402         if(ionchem) then
     403            if (i_co2plus == 0) then
     404               write(*,*) "calchim: Error, no CO2+ tracer !!!"
     405               write(*,*) "CO2+ is needed if O2+ is in traceur.def"
     406               stop
     407            else
     408               nbq = nbq + 1
     409               niq(nbq) = i_co2plus
     410            end if
     411         else
     412            if (i_co2plus /= 0) then
     413               write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
     414               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     415               stop
     416            endif
     417         endif
    404418         i_oplus=igcm_oplus
    405          if (i_oplus == 0) then
    406             write(*,*) "calchim: Error, no O+ tracer !!!"
    407             stop
    408          else
    409             nbq = nbq + 1
    410             niq(nbq) = i_oplus
    411          end if
     419         if(ionchem) then
     420            if (i_oplus == 0) then
     421               write(*,*) "calchim: Error, no O+ tracer !!!"
     422               write(*,*) "O+ is needed if O2+ is in traceur.def"
     423               stop
     424            else
     425               nbq = nbq + 1
     426               niq(nbq) = i_oplus
     427            end if
     428         else
     429            if (i_oplus /= 0) then
     430               write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
     431               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     432               stop
     433            endif
     434         endif
    412435         i_noplus=igcm_noplus
    413          if (i_noplus == 0) then
    414             write(*,*) "calchim: Error, no NO+ tracer !!!"
    415             stop
    416          else
    417             nbq = nbq + 1
    418             niq(nbq) = i_noplus
    419          end if
     436         if(ionchem) then
     437            if (i_noplus == 0) then
     438               write(*,*) "calchim: Error, no NO+ tracer !!!"
     439               write(*,*) "NO+ is needed if O2+ is in traceur.def"
     440               stop
     441            else
     442               nbq = nbq + 1
     443               niq(nbq) = i_noplus
     444            end if
     445         else
     446            if (i_noplus /= 0) then
     447               write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
     448               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     449            endif
     450         endif
    420451         i_coplus=igcm_coplus
    421          if (i_coplus == 0) then
    422             write(*,*) "calchim: Error, no CO+ tracer !!!"
    423             stop
    424          else
    425             nbq = nbq + 1
    426             niq(nbq) = i_coplus
    427          end if
     452         if(ionchem) then
     453            if (i_coplus == 0) then
     454               write(*,*) "calchim: Error, no CO+ tracer !!!"
     455               write(*,*) "CO+ is needed if O2+ is in traceur.def"
     456               stop               
     457            else
     458               nbq = nbq + 1
     459               niq(nbq) = i_coplus
     460            end if
     461         else
     462            if (i_coplus /= 0) then
     463               write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
     464               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     465            endif
     466         endif
    428467         i_cplus=igcm_cplus
    429          if (i_cplus == 0) then
    430             write(*,*) "calchim: Error, no C+ tracer !!!"
    431             stop
    432          else
    433             nbq = nbq + 1
    434             niq(nbq) = i_cplus
    435          end if
     468         if(ionchem) then
     469            if (i_cplus == 0) then
     470               write(*,*) "calchim: Error, no C+ tracer !!!"
     471               write(*,*) "C+ is needed if O2+ is in traceur.def"
     472               stop
     473            else
     474               nbq = nbq + 1
     475               niq(nbq) = i_cplus
     476            end if
     477         else
     478            if (i_cplus /= 0) then
     479               write(*,*) "calchim: Error: C+ is present, but O2+ is not!!!"
     480               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     481            endif
     482         endif
    436483         i_n2plus=igcm_n2plus
    437          if (i_n2plus == 0) then
    438             write(*,*) "calchim: Error, no N2+ tracer !!!"
    439             stop
    440          else
    441             nbq = nbq + 1
    442             niq(nbq) = i_n2plus
    443          end if
     484         if(ionchem) then
     485            if (i_n2plus == 0) then
     486               write(*,*) "calchim: Error, no N2+ tracer !!!"
     487               write(*,*) "N2+ is needed if O2+ is in traceur.def"
     488               stop
     489            else
     490               nbq = nbq + 1
     491               niq(nbq) = i_n2plus
     492            end if
     493         else
     494            if (i_n2plus /= 0) then
     495               write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
     496               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     497            endif
     498         endif
    444499         i_nplus=igcm_nplus
    445          if (i_nplus == 0) then
    446             write(*,*) "calchim: Error, no N+ tracer !!!"
    447             stop
    448          else
    449             nbq = nbq + 1
    450             niq(nbq) = i_nplus
    451          end if
     500         if(ionchem) then
     501            if (i_nplus == 0) then
     502               write(*,*) "calchim: Error, no N+ tracer !!!"
     503               write(*,*) "N+ is needed if O2+ is in traceur.def"
     504               stop
     505            else
     506               nbq = nbq + 1
     507               niq(nbq) = i_nplus
     508            end if
     509         else
     510            if (i_nplus /= 0) then
     511               write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
     512               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     513            endif
     514         endif
    452515         i_hplus=igcm_hplus
    453          if (i_hplus == 0) then
    454             write(*,*) "calchim: Error, no H+ tracer !!!"
    455             stop
    456          else
    457             nbq = nbq + 1
    458             niq(nbq) = i_hplus
    459          end if
     516         if(ionchem) then
     517            if (i_hplus == 0) then
     518               write(*,*) "calchim: Error, no H+ tracer !!!"
     519               write(*,*) "H+ is needed if O2+ is in traceur.def"
     520               stop
     521            else
     522               nbq = nbq + 1
     523               niq(nbq) = i_hplus
     524            end if
     525         else
     526            if (i_hplus /= 0) then
     527               write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
     528               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     529            endif
     530         endif
    460531         i_hco2plus=igcm_hco2plus
    461          if (i_hco2plus == 0) then
    462             write(*,*) "calchim: Error, no HCO2+ tracer !!!"
    463             stop
    464          else
    465             nbq = nbq + 1
    466             niq(nbq) = i_hco2plus
    467          end if
     532         if(ionchem) then
     533            if (i_hco2plus == 0) then
     534               write(*,*) "calchim: Error, no HCO2+ tracer !!!"
     535               write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
     536               stop
     537            else
     538               nbq = nbq + 1
     539               niq(nbq) = i_hco2plus
     540            end if
     541         else
     542            if (i_hco2plus /= 0) then
     543               write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
     544               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     545            endif
     546         endif
    468547         i_elec = igcm_elec
    469          if (i_elec == 0) then
    470             write(*,*) "calchim: Error, no e- tracer !!!"
    471             stop
    472          else
    473             nbq = nbq + 1
    474             niq(nbq) = i_elec
    475          end if
    476          
     548         if(ionchem) then
     549            if (i_elec == 0) then
     550               write(*,*) "calchim: Error, no e- tracer !!!"
     551               write(*,*) "e- is needed if O2+ is in traceur.def"
     552               stop
     553            else
     554               nbq = nbq + 1
     555               niq(nbq) = i_elec
     556            end if
     557         else
     558            if (i_elec /= 0) then
     559               write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
     560               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     561            endif
     562         endif
    477563         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
    478564               
    … …  
    523609
    524610!           search for switch index between regions 
    525             if(unichim) then
    526                lswitch=nlayer+1
     611
     612            if (unichim) then
     613               lswitch = nlayer + 1
    527614            else
    528615               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
    … …  
    540627!=======================================================================
    541628
    542 !        chemistry: only one scheme at all layers
    543          
    544629         if (photochem) then
    545             call photochemistry(nlayer,nq,                            &
    546                            ig,lswitch,zycol,szacol,ptimestep,         &
    547                            zpress,zlocal,ztemp,ztelec,zdens,zmmean,   &
    548                            dist_sol,zday,surfdust1d,surfice1d,        &
     630            ! set number of reactions, depending on ion chemistry or not
     631            if (ionchem) then
     632               nb_reaction_4_max = 60   ! set number of bimolecular reactions
     633               nb_reaction_3_max = 6    ! set number of quadratic reactions
     634               nquench           = 9    ! set number of quenching + heterogeneous reactions
     635               nphotion          = 18   ! set number of photoionizations
     636            else
     637               nb_reaction_4_max = 31   ! set number of bimolecular reactions
     638               nb_reaction_3_max = 6    ! set number of quadratic reactions
     639               nquench           = 9    ! set number of quenching + heterogeneous reactions
     640               nphotion          = 0    ! set number of photoionizations
     641            end if
     642
     643!        nb_phot_max is the total number of processes that are treated
     644!        numerically as a photolysis:
     645
     646            nb_phot_max = nphot + nphotion + nquench
     647
     648!        call main photochemistry routine
     649
     650            call photochemistry(nlayer,nq,nbq,ionchem,nb_reaction_3_max,  &
     651                           nb_reaction_4_max, nb_phot_max, nphotion,      &
     652                           jonline, ig,lswitch,zycol,szacol,ptimestep,    &
     653                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,       &
     654                           dist_sol,zday,surfdust1d,surfice1d,            &
    549655                           jo3,jh2o,taucol,iter)
    550656
    … …  
    556662               iter_3d(ig,l) = iter(l)
    557663            end do
     664
    558665!        condensation of h2o2
    559666
    … …  
    562669                         ztemp,zycol,dqcloud,dqscloud)
    563670
    564             if(.not.unichim) then
    565                chemthermod=3   !C/O/H/N/ions
     671!        case of separate photochemical model in the thermosphere
     672
     673            if (.not.unichim) then
     674               chemthermod = 3   !C/O/H/N/ions
    566675               call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
    567676                    zdens,zpress,zlocal,szacol,ptimestep,zday,&
    568677                    em_no,em_o2)
    569                do l=1,nlayer
    570                   emission_no(ig,l)=em_no(l)
    571                   emission_o2(ig,l)=em_o2(l)
    572                enddo
    573             end if
    574          end if
     678               do l = 1,nlayer
     679                  emission_no(ig,l) = em_no(l)
     680                  emission_o2(ig,l) = em_o2(l)
     681               end do
     682            end if
     683
     684         end if  ! photochem
    575685
    576686!        dry deposition
    … …  
    627737            call writediagfi(ngrid,'iter','iterations',  &
    628738                             ' ',3,iter_3d(1,1))
    629             if(.not.unichim) then
     739
     740            if (.not. unichim) then
    630741               call writediagfi(ngrid,'emission_no',        &
    631742                    'NO nightglow emission rate','cm-3 s-1',3,emission_no)

  • trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90

    r2158 r2170  
    1313!*****************************************************************
    1414
    15 subroutine photochemistry(nlayer, nq,                                   &
    16                           ig, lswitch, zycol, sza, ptimestep, press,    &
    17                           alt, temp, temp_elect,dens, zmmean,           &
    18                           dist_sol, zday,                               &
     15subroutine photochemistry(nlayer, nq, nesp, ionchem, nb_reaction_3_max,        &
     16                          nb_reaction_4_max, nb_phot_max, nphotion,            &
     17                          jonline, ig, lswitch, zycol, sza, ptimestep, press,  &
     18                          alt, temp, temp_elect, dens, zmmean,                 &
     19                          dist_sol, zday,                                      &
    1920                          surfdust1d, surfice1d, jo3, jh2o,tau, iter)
    20 
    21 use photolysis_mod, only : nb_phot_max,       &
    22                            nb_reaction_3_max, &
    23                            nb_reaction_4_max, &
    24                            jonline
    2521
    2622use param_v4_h, only: jion
    … …  
    3531
    3632integer, intent(in) :: nlayer ! number of atmospheric layers
    37 integer, intent(in) :: nq     ! number of tracers
     33integer, intent(in) :: nq     ! number of tracers in traceur.def
     34integer, intent(in) :: nesp   ! number of traceurs in chemistry
     35integer, intent(in) :: nb_reaction_3_max   
     36                              ! number of quadratic reactions
     37integer, intent(in) :: nb_reaction_4_max
     38                              ! number of bimolecular reactions
     39integer, intent(in) :: nb_phot_max
     40                              ! number of reactions treated numerically as photodissociations
     41integer, intent(in) :: nphotion
     42                              ! number of photoionizations
     43logical, intent(in) :: ionchem! switch for ion chemistry
     44logical, intent(in) :: jonline! switch for on-line calculation of photolysis rates
    3845integer :: ig                 ! grid point index
    3946     
    … …  
    6976!     local:
    7077!===================================================================
    71 
    72 integer, parameter :: nesp = 28  ! number of species in the chemical code
    7378
    7479integer :: phychemrat            ! (physical timestep)/(nominal chemical timestep)
    … …  
    109114integer,parameter :: i_elec    = 28
    110115
    111 
    112116integer :: ilay
    113117
    … …  
    145149if (firstcall) then
    146150   print*,'photochemistry: initialize indexes'
    147    call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    148                i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    149                i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
    150                i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
     151   call indice(nb_reaction_3_max,nb_reaction_4_max,                 &
     152               nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
     153               i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
     154               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,            &
     155               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
    151156               i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
    152157   firstcall = .false.
    … …  
    157162!===================================================================
    158163
    159 call gcmtochim(nlayer, nq, zycol, lswitch, nesp,               &
     164call gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,      &
    160165               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    161166               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    … …  
    169174!===================================================================
    170175
    171 jionos= .true.
     176jionos = .true.
    172177
    173178if (jonline) then
    174179   if (sza <= 113.) then ! day at 300 km
    175       call photolysis_online(nlayer, alt, press, temp, zmmean,          &
     180      call photolysis_online(nlayer, nb_phot_max, alt, press, temp, zmmean, &
    176181                             i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
    177182                             i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
    … …  
    179184                             tau, sza, dist_sol, v_phot)
    180185
    181       if(jionos) then
     186      if (jionos .and. ionchem) then
    182187         call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
    183188         do ilay=1,lswitch-1
    … …  
    246251else
    247252   tau = tau*7./press(1) ! dust in the lookup table is at 7 hpa
    248    call photolysis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, &
    249                    rm(:,i_co2), rm(:,i_o3), v_phot)
     253   call photolysis(nlayer, nb_phot_max, lswitch, press, temp, sza, tau,  &
     254                   zmmean, dist_sol,rm(:,i_co2), rm(:,i_o3), v_phot)
    250255end if
    251256! save o3 and h2o photolysis for output
    … …  
    267272hetero_ice  = .true.
    268273
    269 call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2),   &
    270                    c(:,i_o), c(:,i_n2), press, temp, temp_elect,      &
    271                    hetero_dust, hetero_ice,                           &
     274call reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
     275                   nb_phot_max, nphotion, lswitch, dens, c(:,i_co2),      &
     276                   c(:,i_o2), c(:,i_o), c(:,i_n2), press, temp,           &
     277                   temp_elect, hetero_dust, hetero_ice,                   &
    272278                   surfdust1d, surfice1d, v_phot, v_3, v_4)
    273279
    … …  
    307313!  first-guess: fill matrix
    308314
    309    call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
     315   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer,              &
     316                    nb_reaction_3_max, nb_reaction_4_max, nb_phot_max,    &
     317                    v_phot, v_3, v_4)
    310318
    311319!  adaptative evaluation of the sub time step
    … …  
    349357   cold(:)   = c(ilev,:)
    350358   c(ilev,:) = cnew(:)
    351    !Mod FGG, July 2019
    352    !Force charge neutrality
    353    if(c(ilev,i_elec).ne.c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
    354         c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+ &
    355         c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)) then
    356       c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
    357            c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+ &
    358            c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)
    359 !      write(*,*)'photochemistry/359'
    360 !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
    361    endif
     359
     360!  force charge neutrality (mod fgg, july 2019)
     361
     362   if (ionchem) then
     363      if(c(ilev,i_elec).ne.c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+&
     364           c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+&
     365           c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)) then
     366         c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
     367              c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+              &
     368              c(ilev,i_n2plus)+c(ilev,i_nplus)+c(ilev,i_hplus)+               &
     369              c(ilev,i_hco2plus)
     370         !      write(*,*)'photochemistry/359'
     371         !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
     372      end if
     373   end if
    362374   cnew(:)   = 0.
    363375
    … …  
    375387!===================================================================
    376388
    377 call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
    378                i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
    379                i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
    380                i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,      &
    381                i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,   &
    382                i_n2plus, i_nplus, i_hplus, i_hco2plus,        &
     389call chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp,      &
     390               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
     391               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
     392               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
     393               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
     394               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
    383395               i_elec, dens, c)
    384396contains
    … …  
    463475
    464476#ifdef LAPACK
    465       call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
     477   call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
    466478#else
    467479   write(*,*) "photochemistry error, missing LAPACK routine dgesv"
    … …  
    505517!======================================================================
    506518
    507  subroutine reactionrates(nlayer,                                  &
    508                           lswitch, dens, co2, o2, o, n2, press, &
    509                           t, t_elect,                              &
    510                           hetero_dust, hetero_ice,                 &
    511                           surfdust1d, surfice1d,                   &
    512                           v_phot, v_3, v_4)
     519 subroutine reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
     520                          nb_phot_max, nphotion, lswitch, dens, co2, o2, o,      &
     521                          n2, press, t, t_elect, hetero_dust, hetero_ice,        &
     522                          surfdust1d, surfice1d, v_phot, v_3, v_4)
    513523 
    514524!================================================================
    … …  
    524534
    525535use comcstfi_h
    526 use photolysis_mod, only : nphot, nphotion, nb_phot_max, &
    527                            nb_reaction_3_max,  &
    528                            nb_reaction_4_max
     536use photolysis_mod, only : nphot
    529537
    530538implicit none
    … …  
    535543
    536544integer, intent(in)     :: nlayer            ! number of atmospheric layers
     545integer, intent(in)     :: nb_reaction_3_max ! number of quadratic reactions
     546integer, intent(in)     :: nb_reaction_4_max ! number of bimolecular reactions
     547integer, intent(in)     :: nb_phot_max       ! number of reactions treated numerically as photodissociations
     548integer, intent(in)     :: nphotion          ! number of photoionizations
     549logical, intent(in)     :: ionchem
    537550integer                 :: lswitch           ! interface level between lower
    538551                                             ! atmosphere and thermosphere chemistries
    … …  
    563576integer          :: ilev
    564577integer          :: nb_phot, nb_reaction_3, nb_reaction_4
    565 real :: ak0, ak1, xpo, rate
     578real :: ak0, ak1, xpo, rate, rate1, rate2
    566579real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
    567580real, dimension(nlayer) :: deq
    … …  
    586599!----------------------------------------------------------------------
    587600
    588       nb_phot       = nphot + nphotion ! initialised to the number of photolysis and photoionization rates
     601      nb_phot       = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates
    589602      nb_reaction_3 = 0
    590603      nb_reaction_4 = 0
    … …  
    10221035!     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
    10231036
    1024 !     jpl 2006
    1025 
    1026 !     ak0 = 1.5e-13*(t(:)/300.)**(0.6)
    1027 !     ak1 = 2.1e-9*(t(:)/300.)**(6.1)
    1028 !     rate1 = ak0/(1. + ak0/(ak1/dens(:)))
    1029 !     xpo1 = 1./(1. + alog10(ak0/(ak1/dens(:)))**2)
    1030 
    1031 !     ak0 = 5.9e-33*(t(:)/300.)**(-1.4)
    1032 !     ak1 = 1.1e-12*(t(:)/300.)**(1.3)
    1033 !     rate2 = (ak0*dens(:))/(1. + ak0*dens(:)/ak1)
    1034 !     xpo2 = 1./(1. + alog10((ak0*dens(:))/ak1)**2)
    1035 
    1036 !     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
     1037!     jpl 2015
     1038
     1039!     do ilev = 1,lswitch-1
     1040
     1041!        branch 1 : oh + co -> h + co2
     1042
     1043!        rate1 = 1.5e-13*(t(ilev)/300.)**(0.0)
     1044
     1045!        branch 2 : oh + co + m -> hoco + m
     1046
     1047!        ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0)
     1048!        ak1 = 1.1e-12*(t(ilev)/300.)**(1.3)
     1049!        rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
     1050!        xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
     1051
     1052!        e001(ilev) = rate1 + rate2*0.6**xpo
     1053!     end do
    10371054
    10381055!     joshi et al., 2006
    … …  
    10531070         k1a = k1a0*((1. + y)/(1. + x))*fx
    10541071         k1b = k1b0*(1./(1.+x))*fx
    1055 
    10561072         e001(ilev) = k1a + k1b
    10571073      end do
    … …  
    10691085      v_4(:,nb_reaction_4) = e002(:)
    10701086
    1071 !---  i001: co2+ + o2 -> o2+ + co2
    1072 
    1073 !     aninich, j. phys. chem. ref. data 1993
    1074 
    1075       i001(:) = 5.5e-11*(300./t_elect(:))**0.5
    1076 
    1077       nb_reaction_4 = nb_reaction_4 + 1
    1078       v_4(:,nb_reaction_4) = i001(:)
    1079 
    1080 !---  i002: co2+ + o -> o+ + co2
    1081 
    1082 !     UMIST database
    1083 
    1084       i002(:) = 9.6e-11
     1087!----------------------------------------------------------------------
     1088!     ionospheric reactions
     1089!     only if ionchem=true
     1090!----------------------------------------------------------------------
     1091
     1092      if (ionchem) then
     1093
     1094!---     i001: co2+ + o2 -> o2+ + co2
     1095
     1096!        aninich, j. phys. chem. ref. data 1993
     1097
     1098         i001(:) = 5.5e-11*(300./t_elect(:))**0.5
     1099
     1100         nb_reaction_4 = nb_reaction_4 + 1
     1101         v_4(:,nb_reaction_4) = i001(:)
     1102
     1103!---     i002: co2+ + o -> o+ + co2
     1104
     1105!        UMIST database
     1106
     1107         i002(:) = 9.6e-11
    10851108     
    1086       nb_reaction_4 = nb_reaction_4 + 1
    1087       v_4(:,nb_reaction_4) = i002(:)
    1088 
    1089 !---  i003: co2+ + o -> o2+ + co
    1090 
    1091 !     UMIST database
    1092 
    1093       i003(:) = 1.64e-10
    1094 
    1095       nb_reaction_4 = nb_reaction_4 + 1
    1096       v_4(:,nb_reaction_4) = i003(:)
    1097 
    1098 !---  i004: o2+ + e- -> o + o
    1099 
    1100 !     Alge et al., J. Phys. B, At. Mol. Phys. 1983
    1101 
    1102       i004(:) = 2.0e-7*(300./t_elect(:))**0.7
    1103 
    1104       nb_reaction_4 = nb_reaction_4 + 1
    1105       v_4(:,nb_reaction_4) = i004(:)
    1106 
    1107 !---  i005: o+ + co2 -> o2+ + co
    1108 
    1109 !     UMIST database
    1110 
    1111       i005(:) = 9.4e-10
    1112 
    1113       nb_reaction_4 = nb_reaction_4 + 1
    1114       v_4(:,nb_reaction_4) = i005(:)
    1115 
    1116 
    1117 !---  i006: co2+ + e- -> co + o
    1118 
    1119 !     UMIST database
    1120 
    1121       i006(:) = 3.8e-7*(300./t_elect(:))**0.5
    1122 
    1123       nb_reaction_4 = nb_reaction_4 + 1
    1124       v_4(:,nb_reaction_4) = i006(:)
    1125 
    1126 
    1127 !---  i007: co2+ + no -> no+ + co2
    1128 
    1129 !     UMIST database
    1130 
    1131       i007(:) = 1.2e-10
    1132 
    1133       nb_reaction_4 = nb_reaction_4 + 1
    1134       v_4(:,nb_reaction_4) = i007(:)
    1135 
    1136 !---  i008: o2+ + no -> no+ + o2
    1137 
    1138 !     UMIST database
    1139 
    1140       i008(:) = 4.6e-10
    1141 
    1142       nb_reaction_4 = nb_reaction_4 + 1
    1143       v_4(:,nb_reaction_4) = i008(:)
    1144 
    1145 !---  i009: o2+ + n2 -> no+ + no
     1109         nb_reaction_4 = nb_reaction_4 + 1
     1110         v_4(:,nb_reaction_4) = i002(:)
     1111
     1112!---     i003: co2+ + o -> o2+ + co
     1113
     1114!        UMIST database
     1115
     1116         i003(:) = 1.64e-10
     1117
     1118         nb_reaction_4 = nb_reaction_4 + 1
     1119         v_4(:,nb_reaction_4) = i003(:)
     1120
     1121!---     i004: o2+ + e- -> o + o
     1122
     1123!        Alge et al., J. Phys. B, At. Mol. Phys. 1983
     1124
     1125         i004(:) = 2.0e-7*(300./t_elect(:))**0.7
     1126
     1127         nb_reaction_4 = nb_reaction_4 + 1
     1128         v_4(:,nb_reaction_4) = i004(:)
     1129
     1130!---     i005: o+ + co2 -> o2+ + co
     1131
     1132!        UMIST database
     1133
     1134         i005(:) = 9.4e-10
     1135
     1136         nb_reaction_4 = nb_reaction_4 + 1
     1137         v_4(:,nb_reaction_4) = i005(:)
     1138
     1139
     1140!---     i006: co2+ + e- -> co + o
     1141
     1142!        UMIST database
     1143
     1144         i006(:) = 3.8e-7*(300./t_elect(:))**0.5
     1145
     1146         nb_reaction_4 = nb_reaction_4 + 1
     1147         v_4(:,nb_reaction_4) = i006(:)
     1148
     1149
     1150!---     i007: co2+ + no -> no+ + co2
     1151
     1152!        UMIST database
     1153
     1154         i007(:) = 1.2e-10
     1155
     1156         nb_reaction_4 = nb_reaction_4 + 1
     1157         v_4(:,nb_reaction_4) = i007(:)
     1158
     1159!---     i008: o2+ + no -> no+ + o2
     1160
     1161!        UMIST database
     1162
     1163         i008(:) = 4.6e-10
     1164
     1165         nb_reaction_4 = nb_reaction_4 + 1
     1166         v_4(:,nb_reaction_4) = i008(:)
     1167
     1168!---     i009: o2+ + n2 -> no+ + no
    11461169     
    1147 !     Fox & Sung 2001
    1148 
    1149       i009(:) = 1.0e-15
     1170!        Fox & Sung 2001
     1171
     1172         i009(:) = 1.0e-15
    11501173     
    1151       nb_reaction_4 = nb_reaction_4 + 1
    1152       v_4(:,nb_reaction_4) = i009(:)
    1153 
    1154 !---  i010: o2+ + n -> no+ + o
    1155 
    1156 !     Fox & Sung 2001
    1157 
    1158       i010(:) = 1.0e-10
    1159 
    1160       nb_reaction_4 = nb_reaction_4 + 1
    1161       v_4(:,nb_reaction_4) = i010(:)
    1162 
    1163 !---  i011: o+ + n2 -> no+ + n
    1164 
    1165 !     Fox & Sung 2001
    1166 
    1167       i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
    1168 
    1169       nb_reaction_4 = nb_reaction_4 + 1
    1170       v_4(:,nb_reaction_4) = i011(:)
    1171 
    1172 !---  i012: no+ + e -> n + o
    1173 
    1174 !     UMIST database
    1175 
    1176       i012(:) = 4.3e-7*(300./t_elect(:))**0.37
    1177 
    1178       nb_reaction_4 = nb_reaction_4 + 1
    1179       v_4(:,nb_reaction_4) = i012(:)
    1180 
    1181 
    1182 !---  i013: co+ + co2 -> co2+ + co
    1183 
    1184 !     UMIST database
    1185 
    1186       i013(:) = 1.0e-9
    1187 
    1188       nb_reaction_4 = nb_reaction_4 + 1
    1189       v_4(:,nb_reaction_4) = i013(:)
    1190 
    1191 
    1192 !---  i014: co+ + o -> o+ + co
    1193 
    1194 !     UMIST database
    1195 
    1196       i014(:) = 1.4e-10
    1197 
    1198       nb_reaction_4 = nb_reaction_4 + 1
    1199       v_4(:,nb_reaction_4) = i014(:)
    1200 
    1201 !---  i015: c+ + co2 -> co+ + co
    1202 
    1203 !     UMIST database
    1204 
    1205       i015(:) = 1.1e-9
    1206 
    1207       nb_reaction_4 = nb_reaction_4 + 1
    1208       v_4(:,nb_reaction_4) = i015(:)
    1209 
    1210 
    1211 !---  i016: N2+ + co2 -> co2+ + N2
    1212 
    1213 !     Fox & Song 2001
    1214 
    1215       i016(:) = 9.0e-10*(300./t_elect(:))**0.23
    1216 
    1217       nb_reaction_4 = nb_reaction_4 + 1
    1218       v_4(:,nb_reaction_4) = i016(:)
    1219 
    1220 
    1221 !---  i017: N2+ + o -> no+ + N
    1222 
    1223 !     Fox & Song 2001
    1224 
    1225       i017(:) = 1.33e-10*(300./t_elect(:))**0.44
    1226 
    1227       nb_reaction_4 = nb_reaction_4 + 1
    1228       v_4(:,nb_reaction_4) = i017(:)
    1229 
    1230 !---  i018: N2+ + co -> co+ + N2
    1231 
    1232 !     UMIST
    1233 
    1234       i018(:) = 7.4e-11
    1235 
    1236       nb_reaction_4 = nb_reaction_4 + 1
    1237       v_4(:,nb_reaction_4) = i018(:)
    1238 
    1239 !---  i019: N2+ + e -> N + N
    1240 
    1241 !     UMIST
    1242 
    1243       i019(:) = 7.7e-7*(300./t_elect(:))**0.3
    1244 
    1245       nb_reaction_4 = nb_reaction_4 + 1
    1246       v_4(:,nb_reaction_4) = i016(:)
    1247 
    1248 !---  i020: N2+ + o -> o+ + N2
    1249 
    1250 !     Fox & Song 2001
    1251 
    1252       i020(:) = 7.0e-12*(300./t_elect(:))**0.23
    1253 
    1254       nb_reaction_4 = nb_reaction_4 + 1
    1255       v_4(:,nb_reaction_4) = i020(:)
    1256 
    1257 !---  i021: N+ + co2 -> co2+ + N
    1258 
    1259 !     UMIST
    1260 
    1261       i021(:) = 7.5e-10
    1262 
    1263       nb_reaction_4 = nb_reaction_4 + 1
    1264       v_4(:,nb_reaction_4) = i021(:)
    1265 
    1266 !---  i022: CO+ + H -> H+ + CO
    1267 
    1268 !     Fox & Sung 2001
    1269 
    1270       i022(:) = 4.0e-10
    1271 
    1272       nb_reaction_4 = nb_reaction_4 + 1
    1273       v_4(:,nb_reaction_4) = i022(:)
    1274 
    1275 !---  i023: O+ + H -> H+ + O
    1276 
    1277 !     UMIST
    1278 
    1279       i023(:) = 5.66e-10*((t_elect(:)/300.)**0.36)*exp(8.6/t_elect(:))
    1280 
    1281       nb_reaction_4 = nb_reaction_4 + 1
    1282       v_4(:,nb_reaction_4) = i023(:)
    1283 
    1284 !---  i024: H+ + O -> O+ + H
    1285 
    1286 !     UMIST
    1287 
    1288       i024(:) = 6.86e-10*((t_elect(:)/300.)**0.26)*exp(-224.3/t_elect(:))
    1289 
    1290       nb_reaction_4 = nb_reaction_4 + 1
    1291       v_4(:,nb_reaction_4) = i024(:)
    1292 
    1293 !---  i025: CO+ + H2 -> HCO2+ + H
    1294 
    1295 !     UMIST
    1296 
    1297       i025(:) = 9.5e-10
    1298 
    1299       nb_reaction_4 = nb_reaction_4 + 1
    1300       v_4(:,nb_reaction_4) = i025(:)
    1301 
    1302 !---  i026: HCO2+ + e -> H + CO2
    1303 
    1304 !     UMIST
    1305 
    1306       i026(:) = 1.75e-8*((300./t_elect(:))**0.5)
    1307 
    1308       nb_reaction_4 = nb_reaction_4 + 1
    1309       v_4(:,nb_reaction_4) = i026(:)
    1310 
    1311 !---  i027+i028: HCO2+ + e -> H + O + CO
    1312 
    1313 !     UMIST
    1314       !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
    1315       !i028: 0.5 (HCO2+ + e-) -> O + CO
    1316 
    1317       i027(:) = 2.38e-7*((300./t_elect(:))**0.5)
    1318 
    1319       nb_reaction_4 = nb_reaction_4 + 1
    1320       v_4(:,nb_reaction_4) = i027(:)
    1321 
    1322       nb_reaction_4 = nb_reaction_4 + 1
    1323       v_4(:,nb_reaction_4) = i027(:)
    1324 
    1325 !---  i029: HCO2+ + e -> H + CO2
    1326 
    1327 !     UMIST
    1328 
    1329       i029(:) = 9.45e-8*((300./t_elect(:))**0.5)
    1330 
    1331       nb_reaction_4 = nb_reaction_4 + 1
    1332       v_4(:,nb_reaction_4) = i029(:)
     1174         nb_reaction_4 = nb_reaction_4 + 1
     1175         v_4(:,nb_reaction_4) = i009(:)
     1176
     1177!---     i010: o2+ + n -> no+ + o
     1178
     1179!        Fox & Sung 2001
     1180
     1181         i010(:) = 1.0e-10
     1182
     1183         nb_reaction_4 = nb_reaction_4 + 1
     1184         v_4(:,nb_reaction_4) = i010(:)
     1185
     1186!---     i011: o+ + n2 -> no+ + n
     1187
     1188!        Fox & Sung 2001
     1189
     1190         i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
     1191
     1192         nb_reaction_4 = nb_reaction_4 + 1
     1193         v_4(:,nb_reaction_4) = i011(:)
     1194
     1195!---     i012: no+ + e -> n + o
     1196
     1197!        UMIST database
     1198
     1199         i012(:) = 4.3e-7*(300./t_elect(:))**0.37
     1200
     1201         nb_reaction_4 = nb_reaction_4 + 1
     1202         v_4(:,nb_reaction_4) = i012(:)
     1203
     1204
     1205!---     i013: co+ + co2 -> co2+ + co
     1206
     1207!        UMIST database
     1208
     1209         i013(:) = 1.0e-9
     1210
     1211         nb_reaction_4 = nb_reaction_4 + 1
     1212         v_4(:,nb_reaction_4) = i013(:)
     1213
     1214
     1215!---     i014: co+ + o -> o+ + co
     1216
     1217!        UMIST database
     1218
     1219         i014(:) = 1.4e-10
     1220
     1221         nb_reaction_4 = nb_reaction_4 + 1
     1222         v_4(:,nb_reaction_4) = i014(:)
     1223
     1224!---     i015: c+ + co2 -> co+ + co
     1225
     1226!        UMIST database
     1227
     1228         i015(:) = 1.1e-9
     1229
     1230         nb_reaction_4 = nb_reaction_4 + 1
     1231         v_4(:,nb_reaction_4) = i015(:)
     1232
     1233
     1234!---     i016: N2+ + co2 -> co2+ + N2
     1235
     1236!        Fox & Song 2001
     1237
     1238         i016(:) = 9.0e-10*(300./t_elect(:))**0.23
     1239
     1240         nb_reaction_4 = nb_reaction_4 + 1
     1241         v_4(:,nb_reaction_4) = i016(:)
     1242
     1243
     1244!---     i017: N2+ + o -> no+ + N
     1245
     1246!        Fox & Song 2001
     1247
     1248         i017(:) = 1.33e-10*(300./t_elect(:))**0.44
     1249
     1250         nb_reaction_4 = nb_reaction_4 + 1
     1251         v_4(:,nb_reaction_4) = i017(:)
     1252
     1253!---     i018: N2+ + co -> co+ + N2
     1254
     1255!        UMIST
     1256
     1257         i018(:) = 7.4e-11
     1258
     1259         nb_reaction_4 = nb_reaction_4 + 1
     1260         v_4(:,nb_reaction_4) = i018(:)
     1261
     1262!---     i019: N2+ + e -> N + N
     1263
     1264!        UMIST
     1265
     1266         i019(:) = 7.7e-7*(300./t_elect(:))**0.3
     1267
     1268         nb_reaction_4 = nb_reaction_4 + 1
     1269         v_4(:,nb_reaction_4) = i016(:)
     1270
     1271!---     i020: N2+ + o -> o+ + N2
     1272
     1273!        Fox & Song 2001
     1274
     1275         i020(:) = 7.0e-12*(300./t_elect(:))**0.23
     1276
     1277         nb_reaction_4 = nb_reaction_4 + 1
     1278         v_4(:,nb_reaction_4) = i020(:)
     1279
     1280!---     i021: N+ + co2 -> co2+ + N
     1281
     1282!        UMIST
     1283
     1284         i021(:) = 7.5e-10
     1285
     1286         nb_reaction_4 = nb_reaction_4 + 1
     1287         v_4(:,nb_reaction_4) = i021(:)
     1288
     1289!---     i022: CO+ + H -> H+ + CO
     1290
     1291!        Fox & Sung 2001
     1292
     1293         i022(:) = 4.0e-10
     1294
     1295         nb_reaction_4 = nb_reaction_4 + 1
     1296         v_4(:,nb_reaction_4) = i022(:)
     1297
     1298!---     i023: O+ + H -> H+ + O
     1299
     1300!        UMIST
     1301
     1302         i023(:) = 5.66e-10*((t_elect(:)/300.)**0.36)*exp(8.6/t_elect(:))
     1303
     1304         nb_reaction_4 = nb_reaction_4 + 1
     1305         v_4(:,nb_reaction_4) = i023(:)
     1306
     1307!---     i024: H+ + O -> O+ + H
     1308
     1309!        UMIST
     1310
     1311         i024(:) = 6.86e-10*((t_elect(:)/300.)**0.26)*exp(-224.3/t_elect(:))
     1312
     1313         nb_reaction_4 = nb_reaction_4 + 1
     1314         v_4(:,nb_reaction_4) = i024(:)
     1315
     1316!---     i025: CO+ + H2 -> HCO2+ + H
     1317
     1318!        UMIST
     1319
     1320         i025(:) = 9.5e-10
     1321
     1322         nb_reaction_4 = nb_reaction_4 + 1
     1323         v_4(:,nb_reaction_4) = i025(:)
     1324
     1325!---     i026: HCO2+ + e -> H + CO2
     1326
     1327!        UMIST
     1328
     1329         i026(:) = 1.75e-8*((300./t_elect(:))**0.5)
     1330
     1331         nb_reaction_4 = nb_reaction_4 + 1
     1332         v_4(:,nb_reaction_4) = i026(:)
     1333
     1334!---     i027+i028: HCO2+ + e -> H + O + CO
     1335
     1336!        UMIST
     1337         !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
     1338         !i028: 0.5 (HCO2+ + e-) -> O + CO
     1339
     1340         i027(:) = 2.38e-7*((300./t_elect(:))**0.5)
     1341
     1342         nb_reaction_4 = nb_reaction_4 + 1
     1343         v_4(:,nb_reaction_4) = i027(:)
     1344
     1345         nb_reaction_4 = nb_reaction_4 + 1
     1346         v_4(:,nb_reaction_4) = i027(:)
     1347
     1348!---     i029: HCO2+ + e -> H + CO2
     1349
     1350!        UMIST
     1351
     1352         i029(:) = 9.45e-8*((300./t_elect(:))**0.5)
     1353
     1354         nb_reaction_4 = nb_reaction_4 + 1
     1355         v_4(:,nb_reaction_4) = i029(:)
     1356
     1357      end if   !ionchem
    13331358
    13341359!----------------------------------------------------------------------
    … …  
    14131438!======================================================================
    14141439
    1415  subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
     1440 subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer,            &
     1441                        nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
     1442                        v_phot, v_3, v_4)
    14161443
    14171444!======================================================================
    … …  
    14201447
    14211448use types_asis
    1422 use photolysis_mod, only : nb_phot_max,       &
    1423                            nb_reaction_3_max, &
    1424                            nb_reaction_4_max
    14251449
    14261450implicit none
    … …  
    14311455integer             :: nesp    ! number of species in the chemistry
    14321456integer, intent(in) :: nlayer  ! number of atmospheric layers
     1457integer, intent(in) :: nb_reaction_3_max
     1458                               ! number of quadratic reactions
     1459integer, intent(in) :: nb_reaction_4_max
     1460                               ! number of bimolecular reactions
     1461integer, intent(in) :: nb_phot_max
     1462                               ! number of processes treated numerically as photodissociations
    14331463
    14341464real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
    … …  
    15311561!================================================================
    15321562
    1533  subroutine indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    1534                    i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    1535                    i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
    1536                    i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
     1563 subroutine indice(nb_reaction_3_max, nb_reaction_4_max,                &
     1564                   nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2, &
     1565                   i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
     1566                   i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,            &
     1567                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
    15371568                   i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
    15381569
    … …  
    15501581
    15511582use types_asis
    1552 use photolysis_mod, only : nb_phot_max,       &
    1553                            nb_reaction_3_max, &
    1554                            nb_reaction_4_max
    15551583
    15561584implicit none
    … …  
    15631591           i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
    15641592           i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
     1593integer, intent(in) :: nb_reaction_3_max
     1594                       ! number of quadratic reactions
     1595integer, intent(in) :: nb_reaction_4_max
     1596                       ! number of bimolecular reactions
     1597integer, intent(in) :: nb_phot_max
     1598                       ! number of processes treated numerically as photodissociations
     1599logical, intent(in) :: ionchem
    15651600
    15661601! local
    … …  
    16831718indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n)
    16841719
     1720!Only if ion chemistry included
     1721if (ionchem) then
     1722
    16851723!===========================================================
    16861724!      CO2 + hv -> CO2+ + e-
    16871725!===========================================================
    16881726
     1727   nb_phot = nb_phot + 1
     1728
     1729   indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
     1730
     1731!===========================================================
     1732!      CO2 + hv -> O+ + CO + e-
     1733!===========================================================
     1734!We divide this reaction in two
     1735
     1736!0.5 CO2 + hv -> CO
     1737   nb_phot = nb_phot + 1
     1738
     1739   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
     1740
     1741!0.5 CO2 + hv -> O+ + e-
     1742   nb_phot = nb_phot + 1
     1743
     1744   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
     1745
     1746!===========================================================
     1747!      CO2 + hv -> CO+ + O + e-
     1748!===========================================================
     1749!We divide this reaction in two
     1750
     1751!0.5 CO2 + hv -> O
     1752   nb_phot = nb_phot + 1
     1753
     1754   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
     1755
     1756!0.5 CO2 + hv -> CO+ + e-
     1757   nb_phot = nb_phot + 1
     1758
     1759   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
     1760
     1761!===========================================================
     1762!      CO2 + hv -> C+ + O2 + e-
     1763!===========================================================
     1764!We divide this reaction in two
     1765
     1766!0.5 CO2 + hv -> O2
     1767   nb_phot = nb_phot + 1
     1768
     1769   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
     1770
     1771!0.5 CO2 + hv -> C+ + e-
     1772   nb_phot = nb_phot + 1
     1773
     1774   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
     1775
     1776!===========================================================
     1777!      O2 + hv -> O2+ + e-
     1778!===========================================================
     1779
     1780   nb_phot = nb_phot + 1
     1781
     1782   indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
     1783
     1784!===========================================================
     1785!      O + hv -> O+ + e-
     1786!===========================================================
     1787
     1788   nb_phot = nb_phot + 1
     1789
     1790   indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
     1791
     1792!===========================================================
     1793!      NO + hv -> NO+ + e-
     1794!===========================================================
     1795
     1796   nb_phot = nb_phot + 1
     1797
     1798   indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
     1799
     1800!===========================================================
     1801!      CO + hv -> CO+ + e-
     1802!===========================================================
     1803
     1804   nb_phot = nb_phot + 1
     1805
     1806   indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
     1807
     1808!===========================================================
     1809!      CO + hv -> C+ + O + e-
     1810!===========================================================
     1811!We divide this reaction in two
     1812
     1813!0.5 CO + hv -> O
     1814   nb_phot = nb_phot + 1
     1815
     1816   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
     1817
     1818!0.5 CO + hv -> C+ + e-
     1819   nb_phot = nb_phot + 1
     1820
     1821   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
     1822
     1823!===========================================================
     1824!      N2 + hv -> N2+ + e-
     1825!===========================================================
     1826
     1827   nb_phot = nb_phot + 1
     1828
     1829   indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
     1830
     1831!===========================================================
     1832!      N2 + hv -> N+ + N + e-
     1833!===========================================================
     1834!We divide this reaction in two
     1835
     1836!0.5 N2 + hv -> N
     1837   nb_phot = nb_phot + 1
     1838
     1839   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
     1840
     1841!0.5 N2 + hv -> N+ + e-
     1842   nb_phot = nb_phot + 1
     1843
     1844   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
     1845
     1846!===========================================================
     1847!      N + hv -> N+ + e-
     1848!===========================================================
     1849
     1850   nb_phot = nb_phot + 1
     1851
     1852   indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
     1853
     1854!===========================================================
     1855!      H + hv -> H+ + e-
     1856!===========================================================
     1857
     1858   nb_phot = nb_phot + 1
     1859
     1860   indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
     1861
     1862end if   !ionchem
     1863
     1864!===========================================================
     1865!      a001 : O + O2 + CO2 -> O3 + CO2
     1866!===========================================================
     1867
     1868nb_reaction_4 = nb_reaction_4 + 1
     1869
     1870indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
     1871
     1872!===========================================================
     1873!      a002 : O + O + CO2 -> O2 + CO2
     1874!===========================================================
     1875
     1876nb_reaction_3 = nb_reaction_3 + 1
     1877
     1878indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
     1879
     1880!===========================================================
     1881!      a003 : O + O3 -> O2 + O2
     1882!===========================================================
     1883
     1884nb_reaction_4 = nb_reaction_4 + 1
     1885
     1886indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
     1887
     1888!===========================================================
     1889!      b001 : O(1D) + CO2 -> O + CO2
     1890!===========================================================
     1891
    16891892nb_phot = nb_phot + 1
    16901893
    1691 indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
    1692 
    1693 !===========================================================
    1694 !      CO2 + hv -> O+ + CO + e-
     1894indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
     1895
     1896!===========================================================
     1897!      b002 : O(1D) + H2O -> OH + OH
     1898!===========================================================
     1899
     1900nb_reaction_4 = nb_reaction_4 + 1
     1901
     1902indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
     1903
     1904!===========================================================
     1905!      b003 : O(1D) + H2 -> OH + H
     1906!===========================================================
     1907
     1908nb_reaction_4 = nb_reaction_4 + 1
     1909
     1910indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
     1911
     1912!===========================================================
     1913!      b004 : O(1D) + O2 -> O + O2
     1914!===========================================================
     1915
     1916nb_phot = nb_phot + 1
     1917
     1918indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
     1919
     1920!===========================================================
     1921!      b005 : O(1D) + O3 -> O2 + O2
     1922!===========================================================
     1923
     1924nb_reaction_4 = nb_reaction_4 + 1
     1925
     1926indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
     1927
     1928!===========================================================
     1929!      b006 : O(1D) + O3 -> O2 + O + O
     1930!===========================================================
     1931
     1932nb_reaction_4 = nb_reaction_4 + 1
     1933
     1934indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
     1935
     1936!===========================================================
     1937!      c001 : O + HO2 -> OH + O2
     1938!===========================================================
     1939
     1940nb_reaction_4 = nb_reaction_4 + 1
     1941
     1942indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
     1943
     1944!===========================================================
     1945!      c002 : O + OH -> O2 + H
     1946!===========================================================
     1947
     1948nb_reaction_4 = nb_reaction_4 + 1
     1949
     1950indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
     1951
     1952!===========================================================
     1953!      c003 : H + O3 -> OH + O2
     1954!===========================================================
     1955
     1956nb_reaction_4 = nb_reaction_4 + 1
     1957
     1958indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
     1959
     1960!===========================================================
     1961!      c004 : H + HO2 -> OH + OH
     1962!===========================================================
     1963
     1964nb_reaction_4 = nb_reaction_4 + 1
     1965
     1966indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
     1967
     1968!===========================================================
     1969!      c005 : H + HO2 -> H2 + O2
     1970!===========================================================
     1971
     1972nb_reaction_4 = nb_reaction_4 + 1
     1973
     1974indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
     1975
     1976!===========================================================
     1977!      c006 : H + HO2 -> H2O + O
     1978!===========================================================
     1979
     1980nb_reaction_4 = nb_reaction_4 + 1
     1981
     1982indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
     1983
     1984!===========================================================
     1985!      c007 : OH + HO2 -> H2O + O2
     1986!===========================================================
     1987
     1988nb_reaction_4 = nb_reaction_4 + 1
     1989
     1990indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
     1991
     1992!===========================================================
     1993!      c008 : HO2 + HO2 -> H2O2 + O2
     1994!===========================================================
     1995
     1996nb_reaction_3 = nb_reaction_3 + 1
     1997
     1998indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
     1999
     2000!===========================================================
     2001!      c009 : OH + H2O2 -> H2O + HO2
     2002!===========================================================
     2003
     2004nb_reaction_4 = nb_reaction_4 + 1
     2005
     2006indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
     2007
     2008!===========================================================
     2009!      c010 : OH + H2 -> H2O + H
     2010!===========================================================
     2011
     2012nb_reaction_4 = nb_reaction_4 + 1
     2013
     2014indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
     2015
     2016!===========================================================
     2017!      c011 : H + O2 + CO2 -> HO2 + CO2
     2018!===========================================================
     2019
     2020nb_reaction_4 = nb_reaction_4 + 1
     2021
     2022indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
     2023
     2024!===========================================================
     2025!      c012 : O + H2O2 -> OH + HO2
     2026!===========================================================
     2027
     2028nb_reaction_4 = nb_reaction_4 + 1
     2029
     2030indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
     2031
     2032!===========================================================
     2033!      c013 : OH + OH -> H2O + O
     2034!===========================================================
     2035
     2036nb_reaction_3 = nb_reaction_3 + 1
     2037
     2038indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
     2039
     2040!===========================================================
     2041!      c014 : OH + O3 -> HO2 + O2
     2042!===========================================================
     2043
     2044nb_reaction_4 = nb_reaction_4 + 1
     2045
     2046indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
     2047
     2048!===========================================================
     2049!      c015 : HO2 + O3 -> OH + O2 + O2
     2050!===========================================================
     2051
     2052nb_reaction_4 = nb_reaction_4 + 1
     2053
     2054indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
     2055
     2056!===========================================================
     2057!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
     2058!===========================================================
     2059
     2060nb_reaction_3 = nb_reaction_3 + 1
     2061
     2062indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
     2063
     2064!===========================================================
     2065!      c017 : OH + OH + CO2 -> H2O2 + CO2
     2066!===========================================================
     2067
     2068nb_reaction_3 = nb_reaction_3 + 1
     2069
     2070indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
     2071
     2072!===========================================================
     2073!      c018 : H + H + CO2 -> H2 + CO2
     2074!===========================================================
     2075
     2076nb_reaction_3 = nb_reaction_3 + 1
     2077
     2078indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
     2079
     2080!===========================================================
     2081!      d001 : NO2 + O -> NO + O2
     2082!===========================================================
     2083
     2084nb_reaction_4 = nb_reaction_4 + 1
     2085
     2086indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
     2087
     2088!===========================================================
     2089!      d002 : NO + O3 -> NO2 + O2
     2090!===========================================================
     2091
     2092nb_reaction_4 = nb_reaction_4 + 1
     2093
     2094indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
     2095
     2096!===========================================================
     2097!      d003 : NO + HO2 -> NO2 + OH
     2098!===========================================================
     2099
     2100nb_reaction_4 = nb_reaction_4 + 1
     2101
     2102indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
     2103
     2104!===========================================================
     2105!      d004 : N + NO -> N2 + O
     2106!===========================================================
     2107
     2108nb_reaction_4 = nb_reaction_4 + 1
     2109
     2110indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
     2111
     2112!===========================================================
     2113!      d005 : N + O2 -> NO + O
     2114!===========================================================
     2115
     2116nb_reaction_4 = nb_reaction_4 + 1
     2117
     2118indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
     2119
     2120!===========================================================
     2121!      d006 : NO2 + H -> NO + OH
     2122!===========================================================
     2123
     2124nb_reaction_4 = nb_reaction_4 + 1
     2125
     2126indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
     2127
     2128!===========================================================
     2129!      d007 : N + O -> NO
     2130!===========================================================
     2131
     2132nb_reaction_4 = nb_reaction_4 + 1
     2133
     2134indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
     2135
     2136!===========================================================
     2137!      d008 : N + HO2 -> NO + OH
     2138!===========================================================
     2139
     2140nb_reaction_4 = nb_reaction_4 + 1
     2141
     2142indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
     2143
     2144!===========================================================
     2145!      d009 : N + OH -> NO + H
     2146!===========================================================
     2147
     2148nb_reaction_4 = nb_reaction_4 + 1
     2149
     2150indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
     2151
     2152!===========================================================
     2153!      d010 : N(2D) + O -> N + O
     2154!===========================================================
     2155
     2156nb_phot = nb_phot + 1
     2157
     2158indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
     2159
     2160!===========================================================
     2161!      d011 : N(2D) + N2 -> N + N2
     2162!===========================================================
     2163
     2164nb_phot = nb_phot + 1
     2165
     2166indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
     2167
     2168!===========================================================
     2169!      d012 : N(2D) + CO2 -> NO + CO
     2170!===========================================================
     2171
     2172nb_reaction_4 = nb_reaction_4 + 1
     2173
     2174indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
     2175
     2176!===========================================================
     2177!      e001 : CO + OH -> CO2 + H
     2178!===========================================================
     2179
     2180nb_reaction_4 = nb_reaction_4 + 1
     2181
     2182indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
     2183
     2184!===========================================================
     2185!      e002 : CO + O + M -> CO2 + M
     2186!===========================================================
     2187
     2188nb_reaction_4 = nb_reaction_4 + 1
     2189
     2190indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
     2191
     2192!Only if ion chemistry
     2193if (ionchem) then
     2194
     2195!===========================================================
     2196!      i001 : CO2+ + O2 -> O2+ + CO2
     2197!===========================================================
     2198
     2199   nb_reaction_4 = nb_reaction_4 + 1
     2200
     2201   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
     2202
     2203!===========================================================
     2204!      i002 : CO2+ + O -> O+ + CO2
     2205!===========================================================
     2206
     2207   nb_reaction_4 = nb_reaction_4 + 1
     2208
     2209   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
     2210
     2211!===========================================================
     2212!      i003 : CO2+ + O -> O2+ + CO
     2213!===========================================================
     2214
     2215   nb_reaction_4 = nb_reaction_4 + 1
     2216
     2217   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
     2218
     2219!===========================================================
     2220!      i004 : O2+ + e- -> O + O
     2221!===========================================================
     2222
     2223   nb_reaction_4 = nb_reaction_4 + 1
     2224
     2225   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
     2226
     2227!===========================================================
     2228!      i005 : O+ + CO2 -> O2+ + CO
     2229!===========================================================
     2230
     2231   nb_reaction_4 = nb_reaction_4 + 1
     2232
     2233   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
     2234
     2235!===========================================================
     2236!      i006 : CO2+ + e -> CO + O
     2237!===========================================================
     2238
     2239   nb_reaction_4 = nb_reaction_4 + 1
     2240
     2241   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
     2242
     2243!===========================================================
     2244!      i007 : CO2+ + NO -> NO+ + CO2
     2245!===========================================================
     2246
     2247   nb_reaction_4 = nb_reaction_4 + 1
     2248
     2249   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
     2250
     2251!===========================================================
     2252!      i008 : O2+ + NO -> NO+ + O2
     2253!===========================================================
     2254
     2255   nb_reaction_4 = nb_reaction_4 + 1
     2256
     2257   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
     2258
     2259!===========================================================
     2260!      i009 : O2+ + N2 -> NO+ + NO
     2261!===========================================================
     2262
     2263   nb_reaction_4 = nb_reaction_4 + 1
     2264
     2265   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
     2266
     2267!===========================================================
     2268!      i010 : O2+ + N -> NO+ + O
     2269!===========================================================
     2270
     2271   nb_reaction_4 = nb_reaction_4 + 1
     2272
     2273   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
     2274
     2275!===========================================================
     2276!      i011 : O+ + N2 -> NO+ + N
     2277!===========================================================
     2278
     2279   nb_reaction_4 = nb_reaction_4 + 1
     2280
     2281   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
     2282
     2283!===========================================================
     2284!      i012 : NO+ + e -> N + O
     2285!===========================================================
     2286
     2287   nb_reaction_4 = nb_reaction_4 + 1
     2288
     2289   indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
     2290
     2291!===========================================================
     2292!      i013 : CO+ + CO2 -> CO2+ + CO
     2293!===========================================================
     2294
     2295   nb_reaction_4 = nb_reaction_4 + 1
     2296
     2297   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
     2298
     2299!===========================================================
     2300!      i014 : CO+ + O -> O+ + CO
     2301!===========================================================
     2302
     2303   nb_reaction_4 = nb_reaction_4 + 1
     2304
     2305   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
     2306
     2307!===========================================================
     2308!      i015 : C+ + CO2 -> CO+ + CO
     2309!===========================================================
     2310
     2311   nb_reaction_4 = nb_reaction_4 + 1
     2312
     2313   indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
     2314
     2315!===========================================================
     2316!      i016 : N2+ + CO2 -> CO2+ + N2
     2317!===========================================================
     2318
     2319   nb_reaction_4 = nb_reaction_4 + 1
     2320
     2321   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
     2322
     2323!===========================================================
     2324!      i017 : N2+ + O -> NO+ + N
     2325!===========================================================
     2326
     2327   nb_reaction_4 = nb_reaction_4 + 1
     2328
     2329   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
     2330
     2331!===========================================================
     2332!      i018 : N2+ + CO -> CO+ + N2
     2333!===========================================================
     2334
     2335   nb_reaction_4 = nb_reaction_4 + 1
     2336
     2337   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
     2338
     2339!===========================================================
     2340!      i019 : N2+ + e -> N + N
     2341!===========================================================
     2342
     2343   nb_reaction_4 = nb_reaction_4 + 1
     2344
     2345   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
     2346
     2347!===========================================================
     2348!      i020 : N2+ + O -> O+ + N2
     2349!===========================================================
     2350
     2351   nb_reaction_4 = nb_reaction_4 + 1
     2352
     2353   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
     2354
     2355!===========================================================
     2356!      i021 : N+ + CO2 -> CO2+ + N
     2357!===========================================================
     2358
     2359   nb_reaction_4 = nb_reaction_4 + 1
     2360
     2361   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
     2362
     2363!===========================================================
     2364!      i022 : CO+ + H -> H+ + CO
     2365!===========================================================
     2366
     2367   nb_reaction_4 = nb_reaction_4 + 1
     2368
     2369   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
     2370
     2371!===========================================================
     2372!      i023 : O+ + H -> H+ + O
     2373!===========================================================
     2374
     2375   nb_reaction_4 = nb_reaction_4 + 1
     2376
     2377   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
     2378
     2379!===========================================================
     2380!      i024 : H+ + O -> O+ + H
     2381!===========================================================
     2382
     2383   nb_reaction_4 = nb_reaction_4 + 1
     2384
     2385   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
     2386
     2387!===========================================================
     2388!      i025 : CO2+ + H2 -> HCO2+ + H
     2389!===========================================================
     2390
     2391   nb_reaction_4 = nb_reaction_4 + 1
     2392
     2393   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
     2394
     2395!===========================================================
     2396!      i026 : HCO2+ + e -> H + CO2
     2397!===========================================================
     2398
     2399   nb_reaction_4 = nb_reaction_4 + 1
     2400
     2401   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
     2402
     2403!===========================================================
     2404!      i027 : HCO2+ + e -> H + O + CO
    16952405!===========================================================
    16962406!We divide this reaction in two
    16972407
    1698 !0.5 CO2 + hv -> CO
    1699 nb_phot = nb_phot + 1
    1700 
    1701 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
    1702 
    1703 !0.5 CO2 + hv -> O+ + e-
    1704 nb_phot = nb_phot + 1
    1705 
    1706 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
    1707 
    1708 
    1709 !===========================================================
    1710 !      CO2 + hv -> CO+ + O + e-
    1711 !===========================================================
    1712 !We divide this reaction in two
    1713 
    1714 !0.5 CO2 + hv -> O
    1715 nb_phot = nb_phot + 1
    1716 
    1717 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
    1718 
    1719 !0.5 CO2 + hv -> CO+ + e-
    1720 nb_phot = nb_phot + 1
    1721 
    1722 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
    1723 
    1724 
    1725 !===========================================================
    1726 !      CO2 + hv -> C+ + O2 + e-
    1727 !===========================================================
    1728 !We divide this reaction in two
    1729 
    1730 !0.5 CO2 + hv -> O2
    1731 nb_phot = nb_phot + 1
    1732 
    1733 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
    1734 
    1735 !0.5 CO2 + hv -> C+ + e-
    1736 nb_phot = nb_phot + 1
    1737 
    1738 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
    1739 
    1740 !===========================================================
    1741 !      O2 + hv -> O2+ + e-
    1742 !===========================================================
    1743 
    1744 nb_phot = nb_phot + 1
    1745 
    1746 indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
    1747 
    1748 !===========================================================
    1749 !      O + hv -> O+ + e-
    1750 !===========================================================
    1751 
    1752 nb_phot = nb_phot + 1
    1753 
    1754 indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
    1755 
    1756 !===========================================================
    1757 !      NO + hv -> NO+ + e-
    1758 !===========================================================
    1759 
    1760 nb_phot = nb_phot + 1
    1761 
    1762 indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
    1763 
    1764 !===========================================================
    1765 !      CO + hv -> CO+ + e-
    1766 !===========================================================
    1767 
    1768 nb_phot = nb_phot + 1
    1769 
    1770 indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
    1771 
    1772 !===========================================================
    1773 !      CO + hv -> C+ + O + e-
    1774 !===========================================================
    1775 !We divide this reaction in two
    1776 
    1777 !0.5 CO + hv -> O
    1778 nb_phot = nb_phot + 1
    1779 
    1780 indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
    1781 
    1782 !0.5 CO + hv -> C+ + e-
    1783 nb_phot = nb_phot + 1
    1784 
    1785 indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
    1786 
    1787 
    1788 !===========================================================
    1789 !      N2 + hv -> N2+ + e-
    1790 !===========================================================
    1791 
    1792 nb_phot = nb_phot + 1
    1793 
    1794 indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
    1795 
    1796 !===========================================================
    1797 !      N2 + hv -> N+ + N + e-
    1798 !===========================================================
    1799 !We divide this reaction in two
    1800 
    1801 !0.5 N2 + hv -> N
    1802 nb_phot = nb_phot + 1
    1803 
    1804 indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
    1805 
    1806 !0.5 N2 + hv -> N+ + e-
    1807 nb_phot = nb_phot + 1
    1808 
    1809 indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
    1810 
    1811 !===========================================================
    1812 !      N + hv -> N+ + e-
    1813 !===========================================================
    1814 
    1815 nb_phot = nb_phot + 1
    1816 
    1817 indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
    1818 
    1819 !===========================================================
    1820 !      H + hv -> H+ + e-
    1821 !===========================================================
    1822 
    1823 nb_phot = nb_phot + 1
    1824 
    1825 indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
    1826 
    1827 !===========================================================
    1828 !      a001 : O + O2 + CO2 -> O3 + CO2
    1829 !===========================================================
    1830 
    1831 nb_reaction_4 = nb_reaction_4 + 1
    1832 
    1833 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy)
    1834 
    1835 !===========================================================
    1836 !      a002 : O + O + CO2 -> O2 + CO2
    1837 !===========================================================
    1838 
    1839 nb_reaction_3 = nb_reaction_3 + 1
    1840 
    1841 indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy)
    1842 
    1843 !===========================================================
    1844 !      a003 : O + O3 -> O2 + O2
    1845 !===========================================================
    1846 
    1847 nb_reaction_4 = nb_reaction_4 + 1
    1848 
    1849 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
    1850 
    1851 !===========================================================
    1852 !      b001 : O(1D) + CO2 -> O + CO2
    1853 !===========================================================
    1854 
    1855 nb_phot = nb_phot + 1
    1856 
    1857 indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
    1858 
    1859 !===========================================================
    1860 !      b002 : O(1D) + H2O -> OH + OH
    1861 !===========================================================
    1862 
    1863 nb_reaction_4 = nb_reaction_4 + 1
    1864 
    1865 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy)
    1866 
    1867 !===========================================================
    1868 !      b003 : O(1D) + H2 -> OH + H
    1869 !===========================================================
    1870 
    1871 nb_reaction_4 = nb_reaction_4 + 1
    1872 
    1873 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h)
    1874 
    1875 !===========================================================
    1876 !      b004 : O(1D) + O2 -> O + O2
    1877 !===========================================================
    1878 
    1879 nb_phot = nb_phot + 1
    1880 
    1881 indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy)
    1882 
    1883 !===========================================================
    1884 !      b005 : O(1D) + O3 -> O2 + O2
    1885 !===========================================================
    1886 
    1887 nb_reaction_4 = nb_reaction_4 + 1
    1888 
    1889 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy)
    1890 
    1891 !===========================================================
    1892 !      b006 : O(1D) + O3 -> O2 + O + O
    1893 !===========================================================
    1894 
    1895 nb_reaction_4 = nb_reaction_4 + 1
    1896 
    1897 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o)
    1898 
    1899 !===========================================================
    1900 !      c001 : O + HO2 -> OH + O2
    1901 !===========================================================
    1902 
    1903 nb_reaction_4 = nb_reaction_4 + 1
    1904 
    1905 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2)
    1906 
    1907 !===========================================================
    1908 !      c002 : O + OH -> O2 + H
    1909 !===========================================================
    1910 
    1911 nb_reaction_4 = nb_reaction_4 + 1
    1912 
    1913 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h)
    1914 
    1915 !===========================================================
    1916 !      c003 : H + O3 -> OH + O2
    1917 !===========================================================
    1918 
    1919 nb_reaction_4 = nb_reaction_4 + 1
    1920 
    1921 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2)
    1922 
    1923 !===========================================================
    1924 !      c004 : H + HO2 -> OH + OH
    1925 !===========================================================
    1926 
    1927 nb_reaction_4 = nb_reaction_4 + 1
    1928 
    1929 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy)
    1930 
    1931 !===========================================================
    1932 !      c005 : H + HO2 -> H2 + O2
    1933 !===========================================================
    1934 
    1935 nb_reaction_4 = nb_reaction_4 + 1
    1936 
    1937 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2)
    1938 
    1939 !===========================================================
    1940 !      c006 : H + HO2 -> H2O + O
    1941 !===========================================================
    1942 
    1943 nb_reaction_4 = nb_reaction_4 + 1
    1944 
    1945 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o)
    1946 
    1947 !===========================================================
    1948 !      c007 : OH + HO2 -> H2O + O2
    1949 !===========================================================
    1950 
    1951 nb_reaction_4 = nb_reaction_4 + 1
    1952 
    1953 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2)
    1954 
    1955 !===========================================================
    1956 !      c008 : HO2 + HO2 -> H2O2 + O2
    1957 !===========================================================
    1958 
    1959 nb_reaction_3 = nb_reaction_3 + 1
    1960 
    1961 indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
    1962 
    1963 !===========================================================
    1964 !      c009 : OH + H2O2 -> H2O + HO2
    1965 !===========================================================
    1966 
    1967 nb_reaction_4 = nb_reaction_4 + 1
    1968 
    1969 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2)
    1970 
    1971 !===========================================================
    1972 !      c010 : OH + H2 -> H2O + H
    1973 !===========================================================
    1974 
    1975 nb_reaction_4 = nb_reaction_4 + 1
    1976 
    1977 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h)
    1978 
    1979 !===========================================================
    1980 !      c011 : H + O2 + CO2 -> HO2 + CO2
    1981 !===========================================================
    1982 
    1983 nb_reaction_4 = nb_reaction_4 + 1
    1984 
    1985 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy)
    1986 
    1987 !===========================================================
    1988 !      c012 : O + H2O2 -> OH + HO2
    1989 !===========================================================
    1990 
    1991 nb_reaction_4 = nb_reaction_4 + 1
    1992 
    1993 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2)
    1994 
    1995 !===========================================================
    1996 !      c013 : OH + OH -> H2O + O
    1997 !===========================================================
    1998 
    1999 nb_reaction_3 = nb_reaction_3 + 1
    2000 
    2001 indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o)
    2002 
    2003 !===========================================================
    2004 !      c014 : OH + O3 -> HO2 + O2
    2005 !===========================================================
    2006 
    2007 nb_reaction_4 = nb_reaction_4 + 1
    2008 
    2009 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2)
    2010 
    2011 !===========================================================
    2012 !      c015 : HO2 + O3 -> OH + O2 + O2
    2013 !===========================================================
    2014 
    2015 nb_reaction_4 = nb_reaction_4 + 1
    2016 
    2017 indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2)
    2018 
    2019 !===========================================================
    2020 !      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2
    2021 !===========================================================
    2022 
    2023 nb_reaction_3 = nb_reaction_3 + 1
    2024 
    2025 indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2)
    2026 
    2027 !===========================================================
    2028 !      c017 : OH + OH + CO2 -> H2O2 + CO2
    2029 !===========================================================
    2030 
    2031 nb_reaction_3 = nb_reaction_3 + 1
    2032 
    2033 indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy)
    2034 
    2035 !===========================================================
    2036 !      c018 : H + H + CO2 -> H2 + CO2
    2037 !===========================================================
    2038 
    2039 nb_reaction_3 = nb_reaction_3 + 1
    2040 
    2041 indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
    2042 
    2043 !===========================================================
    2044 !      d001 : NO2 + O -> NO + O2
    2045 !===========================================================
    2046 
    2047 nb_reaction_4 = nb_reaction_4 + 1
    2048 
    2049 indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
    2050 
    2051 !===========================================================
    2052 !      d002 : NO + O3 -> NO2 + O2
    2053 !===========================================================
    2054 
    2055 nb_reaction_4 = nb_reaction_4 + 1
    2056 
    2057 indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
    2058 
    2059 !===========================================================
    2060 !      d003 : NO + HO2 -> NO2 + OH
    2061 !===========================================================
    2062 
    2063 nb_reaction_4 = nb_reaction_4 + 1
    2064 
    2065 indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
    2066 
    2067 !===========================================================
    2068 !      d004 : N + NO -> N2 + O
    2069 !===========================================================
    2070 
    2071 nb_reaction_4 = nb_reaction_4 + 1
    2072 
    2073 indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
    2074 
    2075 !===========================================================
    2076 !      d005 : N + O2 -> NO + O
    2077 !===========================================================
    2078 
    2079 nb_reaction_4 = nb_reaction_4 + 1
    2080 
    2081 indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
    2082 
    2083 !===========================================================
    2084 !      d006 : NO2 + H -> NO + OH
    2085 !===========================================================
    2086 
    2087 nb_reaction_4 = nb_reaction_4 + 1
    2088 
    2089 indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
    2090 
    2091 !===========================================================
    2092 !      d007 : N + O -> NO
    2093 !===========================================================
    2094 
    2095 nb_reaction_4 = nb_reaction_4 + 1
    2096 
    2097 indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
    2098 
    2099 !===========================================================
    2100 !      d008 : N + HO2 -> NO + OH
    2101 !===========================================================
    2102 
    2103 nb_reaction_4 = nb_reaction_4 + 1
    2104 
    2105 indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
    2106 
    2107 !===========================================================
    2108 !      d009 : N + OH -> NO + H
    2109 !===========================================================
    2110 
    2111 nb_reaction_4 = nb_reaction_4 + 1
    2112 
    2113 indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
    2114 
    2115 !===========================================================
    2116 !      d010 : N(2D) + O -> N + O
    2117 !===========================================================
    2118 
    2119 nb_phot = nb_phot + 1
    2120 
    2121 indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
    2122 
    2123 !===========================================================
    2124 !      d011 : N(2D) + N2 -> N + N2
    2125 !===========================================================
    2126 
    2127 nb_phot = nb_phot + 1
    2128 
    2129 indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
    2130 
    2131 !===========================================================
    2132 !      d012 : N(2D) + CO2 -> NO + CO
    2133 !===========================================================
    2134 
    2135 nb_reaction_4 = nb_reaction_4 + 1
    2136 
    2137 indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
    2138 
    2139 !===========================================================
    2140 !      e001 : CO + OH -> CO2 + H
    2141 !===========================================================
    2142 
    2143 nb_reaction_4 = nb_reaction_4 + 1
    2144 
    2145 indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h)
    2146 
    2147 !===========================================================
    2148 !      e002 : CO + O + M -> CO2 + M
    2149 !===========================================================
    2150 
    2151 nb_reaction_4 = nb_reaction_4 + 1
    2152 
    2153 indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy)
    2154 
    2155 !===========================================================
    2156 !      i001 : CO2+ + O2 -> O2+ + CO2
    2157 !===========================================================
    2158 
    2159 nb_reaction_4 = nb_reaction_4 + 1
    2160 
    2161 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
    2162 
    2163 !===========================================================
    2164 !      i002 : CO2+ + O -> O+ + CO2
    2165 !===========================================================
    2166 
    2167 nb_reaction_4 = nb_reaction_4 + 1
    2168 
    2169 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
    2170 
    2171 !===========================================================
    2172 !      i003 : CO2+ + O -> O2+ + CO
    2173 !===========================================================
    2174 
    2175 nb_reaction_4 = nb_reaction_4 + 1
    2176 
    2177 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
    2178 
    2179 !===========================================================
    2180 !      i004 : O2+ + e- -> O + O
    2181 !===========================================================
    2182 
    2183 nb_reaction_4 = nb_reaction_4 + 1
    2184 
    2185 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
    2186 
    2187 !===========================================================
    2188 !      i005 : O+ + CO2 -> O2+ + CO
    2189 !===========================================================
    2190 
    2191 nb_reaction_4 = nb_reaction_4 + 1
    2192 
    2193 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
    2194 
    2195 !===========================================================
    2196 !      i006 : CO2+ + e -> CO + O
    2197 !===========================================================
    2198 
    2199 nb_reaction_4 = nb_reaction_4 + 1
    2200 
    2201 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
    2202 
    2203 !===========================================================
    2204 !      i007 : CO2+ + NO -> NO+ + CO2
    2205 !===========================================================
    2206 
    2207 nb_reaction_4 = nb_reaction_4 + 1
    2208 
    2209 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
    2210 
    2211 !===========================================================
    2212 !      i008 : O2+ + NO -> NO+ + O2
    2213 !===========================================================
    2214 
    2215 nb_reaction_4 = nb_reaction_4 + 1
    2216 
    2217 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
    2218 
    2219 !===========================================================
    2220 !      i009 : O2+ + N2 -> NO+ + NO
    2221 !===========================================================
    2222 
    2223 nb_reaction_4 = nb_reaction_4 + 1
    2224 
    2225 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
    2226 
    2227 !===========================================================
    2228 !      i010 : O2+ + N -> NO+ + O
    2229 !===========================================================
    2230 
    2231 nb_reaction_4 = nb_reaction_4 + 1
    2232 
    2233 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
    2234 
    2235 !===========================================================
    2236 !      i011 : O+ + N2 -> NO+ + N
    2237 !===========================================================
    2238 
    2239 nb_reaction_4 = nb_reaction_4 + 1
    2240 
    2241 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
    2242 
    2243 !===========================================================
    2244 !      i012 : NO+ + e -> N + O
    2245 !===========================================================
    2246 
    2247 nb_reaction_4 = nb_reaction_4 + 1
    2248 
    2249 indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
    2250 
    2251 !===========================================================
    2252 !      i013 : CO+ + CO2 -> CO2+ + CO
    2253 !===========================================================
    2254 
    2255 nb_reaction_4 = nb_reaction_4 + 1
    2256 
    2257 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
    2258 
    2259 !===========================================================
    2260 !      i014 : CO+ + O -> O+ + CO
    2261 !===========================================================
    2262 
    2263 nb_reaction_4 = nb_reaction_4 + 1
    2264 
    2265 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
    2266 
    2267 !===========================================================
    2268 !      i015 : C+ + CO2 -> CO+ + CO
    2269 !===========================================================
    2270 
    2271 nb_reaction_4 = nb_reaction_4 + 1
    2272 
    2273 indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
    2274 
    2275 !===========================================================
    2276 !      i016 : N2+ + CO2 -> CO2+ + N2
    2277 !===========================================================
    2278 
    2279 nb_reaction_4 = nb_reaction_4 + 1
    2280 
    2281 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
    2282 
    2283 !===========================================================
    2284 !      i017 : N2+ + O -> NO+ + N
    2285 !===========================================================
    2286 
    2287 nb_reaction_4 = nb_reaction_4 + 1
    2288 
    2289 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
    2290 
    2291 !===========================================================
    2292 !      i018 : N2+ + CO -> CO+ + N2
    2293 !===========================================================
    2294 
    2295 nb_reaction_4 = nb_reaction_4 + 1
    2296 
    2297 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
    2298 
    2299 !===========================================================
    2300 !      i019 : N2+ + e -> N + N
    2301 !===========================================================
    2302 
    2303 nb_reaction_4 = nb_reaction_4 + 1
    2304 
    2305 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
    2306 
    2307 !===========================================================
    2308 !      i020 : N2+ + O -> O+ + N2
    2309 !===========================================================
    2310 
    2311 nb_reaction_4 = nb_reaction_4 + 1
    2312 
    2313 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
    2314 
    2315 !===========================================================
    2316 !      i021 : N+ + CO2 -> CO2+ + N
    2317 !===========================================================
    2318 
    2319 nb_reaction_4 = nb_reaction_4 + 1
    2320 
    2321 indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
    2322 
    2323 !===========================================================
    2324 !      i022 : CO+ + H -> H+ + CO
    2325 !===========================================================
    2326 
    2327 nb_reaction_4 = nb_reaction_4 + 1
    2328 
    2329 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
    2330 
    2331 !===========================================================
    2332 !      i023 : O+ + H -> H+ + O
    2333 !===========================================================
    2334 
    2335 nb_reaction_4 = nb_reaction_4 + 1
    2336 
    2337 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
    2338 
    2339 !===========================================================
    2340 !      i024 : H+ + O -> O+ + H
    2341 !===========================================================
    2342 
    2343 nb_reaction_4 = nb_reaction_4 + 1
    2344 
    2345 indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
    2346 
    2347 !===========================================================
    2348 !      i025 : CO2+ + H2 -> HCO2+ + H
    2349 !===========================================================
    2350 
    2351 nb_reaction_4 = nb_reaction_4 + 1
    2352 
    2353 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
    2354 
    2355 !===========================================================
    2356 !      i026 : HCO2+ + e -> H + CO2
    2357 !===========================================================
    2358 
    2359 nb_reaction_4 = nb_reaction_4 + 1
    2360 
    2361 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
    2362 
    2363 !===========================================================
    2364 !      i027 : HCO2+ + e -> H + O + CO
    2365 !===========================================================
    2366 !We divide this reaction in two
    2367 
    23682408!0.5HCO2+ + 0.5e -> H
    23692409
    2370 nb_reaction_4 = nb_reaction_4 + 1
    2371 
    2372 indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
     2410   nb_reaction_4 = nb_reaction_4 + 1
     2411
     2412   indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
    23732413
    23742414!0.5 HCO2+ + 0.5 e -> O + CO
    23752415
    2376 nb_reaction_4 = nb_reaction_4 + 1
    2377 
    2378 indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
     2416   nb_reaction_4 = nb_reaction_4 + 1
     2417
     2418   indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
    23792419
    23802420!===========================================================
    … …  
    23822422!===========================================================
    23832423
    2384 nb_reaction_4 = nb_reaction_4 + 1
    2385 
    2386 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
     2424   nb_reaction_4 = nb_reaction_4 + 1
     2425
     2426   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
     2427
     2428end if    !ionchem
    23872429
    23882430!===========================================================
    … …  
    24552497!*****************************************************************
    24562498
    2457       subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp,         &
     2499      subroutine gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,&
    24582500                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
    24592501                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
    … …  
    24632505                           i_hplus, i_hco2plus, i_elec, dens, rm, c)
    24642506       
    2465 
    24662507!*****************************************************************
    24672508
    … …  
    24852526      integer, intent(in) :: nlayer ! number of atmospheric layers
    24862527      integer, intent(in) :: nq     ! number of tracers in the gcm
     2528      logical, intent(in) :: ionchem
    24872529      integer :: nesp               ! number of species in the chemistry
    24882530      integer :: lswitch            ! interface level between chemistries
    … …  
    25102552      integer      :: l, iesp
    25112553      logical,save :: firstcall = .true.
    2512      
    25132554     
    25142555!     first call initializations
    … …  
    25862627            stop
    25872628         end if
    2588          if (igcm_co2plus == 0) then
    2589             write(*,*) "gcmtochim: Error; no CO2+ tracer !!!"
    2590             stop
    2591          end if
    2592          if (igcm_oplus == 0) then
    2593             write(*,*) "gcmtochim: Error; no O+ tracer !!!"
    2594             stop
    2595          end if
    2596          if (igcm_o2plus == 0) then
    2597             write(*,*) "gcmtochim: Error; no O2+ tracer !!!"
    2598             stop
    2599          end if
    2600          if (igcm_noplus == 0) then
    2601             write(*,*) "gcmtochim: Error; no NO+ tracer !!!"
    2602             stop
    2603          endif
    2604          if (igcm_coplus == 0) then
    2605             write(*,*) "gcmtochim: Error; no CO+ tracer !!!"
    2606             stop
    2607          endif
    2608          if (igcm_cplus == 0) then
    2609             write(*,*) "gcmtochim: Error; no C+ tracer !!!"
    2610             stop
    2611          endif
    2612          if (igcm_n2plus == 0) then
    2613             write(*,*) "gcmtochim: Error; no N2+ tracer !!!"
    2614             stop
    2615          endif
    2616          if (igcm_nplus == 0) then
    2617             write(*,*) "gcmtochim: Error; no N+ tracer !!!"
    2618             stop
    2619          endif
    2620          if (igcm_hplus == 0) then
    2621             write(*,*) "gcmtochim: Error; no H+ tracer !!!"
    2622             stop
    2623          endif
    2624          if (igcm_hco2plus == 0) then
    2625             write(*,*) "gcmtochim: Error; no HCO2+ tracer !!!"
    2626             stop
    2627          endif
    2628          if (igcm_elec == 0) then
    2629             write(*,*) "gcmtochim: Error; no e- tracer !!!"
    2630             stop
    2631          end if
     2629         if (ionchem) then
     2630            if (igcm_co2plus == 0) then
     2631               write(*,*) "gcmtochim: Error; no CO2+ tracer !!!"
     2632               stop
     2633            end if
     2634            if (igcm_oplus == 0) then
     2635               write(*,*) "gcmtochim: Error; no O+ tracer !!!"
     2636               stop
     2637            end if
     2638            if (igcm_o2plus == 0) then
     2639               write(*,*) "gcmtochim: Error; no O2+ tracer !!!"
     2640               stop
     2641            end if
     2642            if (igcm_noplus == 0) then
     2643               write(*,*) "gcmtochim: Error; no NO+ tracer !!!"
     2644               stop
     2645            endif
     2646            if (igcm_coplus == 0) then
     2647               write(*,*) "gcmtochim: Error; no CO+ tracer !!!"
     2648               stop
     2649            endif
     2650            if (igcm_cplus == 0) then
     2651               write(*,*) "gcmtochim: Error; no C+ tracer !!!"
     2652               stop
     2653            endif
     2654            if (igcm_n2plus == 0) then
     2655               write(*,*) "gcmtochim: Error; no N2+ tracer !!!"
     2656               stop
     2657            endif
     2658            if (igcm_nplus == 0) then
     2659               write(*,*) "gcmtochim: Error; no N+ tracer !!!"
     2660               stop
     2661            endif
     2662            if (igcm_hplus == 0) then
     2663               write(*,*) "gcmtochim: Error; no H+ tracer !!!"
     2664               stop
     2665            endif
     2666            if (igcm_hco2plus == 0) then
     2667               write(*,*) "gcmtochim: Error; no HCO2+ tracer !!!"
     2668               stop
     2669            endif
     2670            if (igcm_elec == 0) then
     2671               write(*,*) "gcmtochim: Error; no e- tracer !!!"
     2672               stop
     2673            end if
     2674         end if  ! ionchem
    26322675         firstcall = .false.
    26332676      end if ! of if (firstcall)
    … …  
    26552698         rm(l,i_no2)  = zycol(l, igcm_no2)
    26562699         rm(l,i_n2)   = zycol(l, igcm_n2)
    2657          rm(l,i_co2plus) = zycol(l, igcm_co2plus)
    2658          rm(l,i_oplus) = zycol(l, igcm_oplus)
    2659          rm(l,i_o2plus) = zycol(l, igcm_o2plus)
    2660          rm(l,i_noplus) = zycol(l,igcm_noplus)
    2661          rm(l,i_coplus) = zycol(l,igcm_coplus)
    2662          rm(l,i_cplus)  = zycol(l,igcm_cplus)
    2663          rm(l,i_n2plus) = zycol(l,igcm_n2plus)
    2664          rm(l,i_nplus)  = zycol(l,igcm_nplus)
    2665          rm(l,i_hplus)  = zycol(l,igcm_hplus)
    2666          rm(l,i_hco2plus)=zycol(l,igcm_hco2plus)
    2667          rm(l,i_elec) = zycol(l, igcm_elec)
    26682700      end do
     2701
     2702      if (ionchem) then
     2703         do l = 1,nlayer
     2704            rm(l,i_co2plus)  = zycol(l, igcm_co2plus)
     2705            rm(l,i_oplus)    = zycol(l, igcm_oplus)
     2706            rm(l,i_o2plus)   = zycol(l, igcm_o2plus)
     2707            rm(l,i_noplus)   = zycol(l, igcm_noplus)
     2708            rm(l,i_coplus)   = zycol(l, igcm_coplus)
     2709            rm(l,i_cplus)    = zycol(l, igcm_cplus)
     2710            rm(l,i_n2plus)   = zycol(l, igcm_n2plus)
     2711            rm(l,i_nplus)    = zycol(l, igcm_nplus)
     2712            rm(l,i_hplus)    = zycol(l, igcm_hplus)
     2713            rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
     2714            rm(l,i_elec)     = zycol(l, igcm_elec)
     2715         end do
     2716      end if
    26692717
    26702718      where (rm(:,:) < 1.e-30)
    … …  
    26862734!*****************************************************************
    26872735 
    2688       subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp,        &
    2689                            i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
    2690                            i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
    2691                            i_n, i_n2d, i_no, i_no2, i_n2,            &
    2692                            i_co2plus, i_oplus, i_o2plus, i_noplus,   &
    2693                            i_coplus, i_cplus, i_n2plus, i_nplus,     &
     2736      subroutine chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp, &
     2737                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
     2738                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
     2739                           i_n, i_n2d, i_no, i_no2, i_n2,             &
     2740                           i_co2plus, i_oplus, i_o2plus, i_noplus,    &
     2741                           i_coplus, i_cplus, i_n2plus, i_nplus,      &
    26942742                           i_hplus, i_hco2plus, i_elec, dens, c)
    2695        
    26962743 
    26972744!*****************************************************************
    … …  
    27162763      integer, intent(in) :: nlayer  ! number of atmospheric layers
    27172764      integer, intent(in) :: nq      ! number of tracers in the gcm
     2765      logical, intent(in) :: ionchem
    27182766      integer :: nesp                ! number of species in the chemistry
    27192767      integer :: lswitch             ! interface level between chemistries
    … …  
    27622810         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
    27632811         zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
    2764          zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
    2765          zycol(l, igcm_oplus)   = c(l,i_oplus)/dens(l)
    2766          zycol(l, igcm_o2plus)  = c(l,i_o2plus)/dens(l)
    2767          zycol(l, igcm_noplus)  = c(l,i_noplus)/dens(l)
    2768          zycol(l, igcm_coplus)  = c(l,i_coplus)/dens(l)
    2769          zycol(l, igcm_cplus)   = c(l,i_cplus)/dens(l)
    2770          zycol(l, igcm_n2plus)  = c(l,i_n2plus)/dens(l)
    2771          zycol(l, igcm_nplus)   = c(l,i_nplus)/dens(l)
    2772          zycol(l, igcm_hplus)   = c(l,i_hplus)/dens(l)
    2773          zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
    2774          zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
    27752812      end do
    27762813
     2814      if (ionchem) then
     2815         do l = 1,lswitch-1
     2816            zycol(l, igcm_co2plus) = c(l,i_co2plus)/dens(l)
     2817            zycol(l, igcm_oplus)   = c(l,i_oplus)/dens(l)
     2818            zycol(l, igcm_o2plus)  = c(l,i_o2plus)/dens(l)
     2819            zycol(l, igcm_noplus)  = c(l,i_noplus)/dens(l)
     2820            zycol(l, igcm_coplus)  = c(l,i_coplus)/dens(l)
     2821            zycol(l, igcm_cplus)   = c(l,i_cplus)/dens(l)
     2822            zycol(l, igcm_n2plus)  = c(l,i_n2plus)/dens(l)
     2823            zycol(l, igcm_nplus)   = c(l,i_nplus)/dens(l)
     2824            zycol(l, igcm_hplus)   = c(l,i_hplus)/dens(l)
     2825            zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
     2826            zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
     2827         end do
     2828      end if
     2829
    27772830      end subroutine chimtogcm
    27782831
    27792832end subroutine photochemistry
    2780 

  • trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90

    r2030 r2170  
    11!==========================================================================
    22
    3       subroutine photolysis(nlayer,                              &
     3      subroutine photolysis(nlayer, nb_phot_max,                  &
    44                            lswitch, press, temp, sza, tauref,   &
    55                            zmmean, dist_sol, rmco2, rmo3, v_phot)
    … …  
    88
    99      use comcstfi_h
    10       use photolysis_mod, only : nb_phot_max
    1110
    1211      implicit none
    … …  
    1817!==========================================================================
    1918       
    20       integer, intent(in) :: nlayer ! number of atmospheric layers
    21       integer :: lswitch            ! interface level between chemistries
    22       real :: press(nlayer)         ! pressure (hPa)
    23       real :: temp(nlayer)          ! temperature (K)
    24       real :: sza                   ! solar zenith angle (deg)
    25       real :: tauref                ! optical depth at 7 hpa
    26       real :: zmmean(nlayer)        ! mean molecular mass (g)
    27       real :: dist_sol              ! sun distance (AU)
    28       real :: rmco2(nlayer)         ! co2 mixing ratio
    29       real :: rmo3(nlayer)          ! ozone mixing ratio
     19      integer, intent(in) :: nlayer      ! number of atmospheric layers
     20      integer, intent(in) :: nb_phot_max ! number of processes treated numerically as photodissociations
     21      integer :: lswitch                 ! interface level between chemistries
     22      real :: press(nlayer)              ! pressure (hPa)
     23      real :: temp(nlayer)               ! temperature (K)
     24      real :: sza                        ! solar zenith angle (deg)
     25      real :: tauref                     ! optical depth at 7 hpa
     26      real :: zmmean(nlayer)             ! mean molecular mass (g)
     27      real :: dist_sol                   ! sun distance (AU)
     28      real :: rmco2(nlayer)              ! co2 mixing ratio
     29      real :: rmo3(nlayer)               ! ozone mixing ratio
    3030
    3131!==========================================================================

  • trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90

    r2158 r2170  
    44
    55! photolysis
    6   logical, parameter :: jonline = .true.  ! true: on-line ! false: lookup table
     6
    77  integer, parameter :: nphot = 13        ! number of photolysis
    8   integer, parameter :: nphotion = 18     ! number of photoionizations
    98  integer, parameter :: nabs  = 10        ! number of absorbing gases
    10 
    11 ! number of reactions in chemical solver
    12 
    13   integer, parameter :: nb_phot_max       = nphot + nphotion + 9 ! photolysis + photoionization + quenching/heterogeneous
    14   integer, parameter :: nb_reaction_3_max = 6         ! quadratic
    15   integer, parameter :: nb_reaction_4_max = 60        ! bimolecular
    169
    1710! spectral grid

  • trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F

    r2041 r2170  
    11!==============================================================================
    22
    3       subroutine photolysis_online(nlayer, alt, press, temp,
    4      $           zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
     3      subroutine photolysis_online(nlayer, nb_phot_max, alt, press,
     4     $           temp, zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
    55     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               
    66     $           i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,
    … …  
    1515      integer, intent(in) :: nesp                                    ! total number of chemical species
    1616      integer, intent(in) :: nlayer
     17      integer, intent(in) :: nb_phot_max
    1718      integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
    1819     $                       i_h2, i_oh, i_ho2, i_h2o2, i_h2o,
Note: See TracChangeset for help on using the changeset viewer.