Ignore:
Timestamp:
May 13, 2020, 9:24:14 AM (5 years ago)
Author:
emillour
Message:

Mars GCM:

  • Finalized water ion chemisty (added H3O+ and OH+ ions). Added an example of relevent traceur.def file in deftank.
  • Added the computation of NO nightglow.

FGG

Location:
trunk/LMDZ.MARS/libf/aeronomars
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r2302 r2321  
    2323                            igcm_noplus, igcm_n2plus, igcm_hplus,         &
    2424                            igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
    25                             igcm_elec, mmol
     25                            igcm_h3oplus, igcm_ohplus, igcm_elec, mmol
    2626
    2727      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
    … …  
    144144      integer,save :: i_hcoplus=0
    145145      integer,save :: i_h2oplus=0
     146      integer,save :: i_h3oplus=0
     147      integer,save :: i_ohplus=0
    146148      integer,save :: i_elec=0
    147149
    … …  
    577579            if (i_h2oplus /= 0) then
    578580               write(*,*) "calchim: Error: H2O+ is present, but O2+ is not!!!"
     581               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     582            endif
     583         endif
     584         i_h3oplus=igcm_h3oplus
     585         if(ionchem) then
     586            if (i_h3oplus == 0) then
     587               write(*,*) "calchim: Error, no H3O+ tracer !!!"
     588               write(*,*) "H3O+ is needed if O2+ is in traceur.def"
     589               call abort_physic("calchim","missing h3oplus tracer",1)
     590            else
     591               nbq = nbq + 1
     592               niq(nbq) = i_h3oplus
     593            end if
     594         else
     595            if (i_h3oplus /= 0) then
     596               write(*,*) "calchim: Error: H3O+ is present, but O2+ is not!!!"
     597               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     598            endif
     599         endif
     600         i_ohplus=igcm_ohplus
     601         if(ionchem) then
     602            if (i_ohplus == 0) then
     603               write(*,*) "calchim: Error, no OH+ tracer !!!"
     604               write(*,*) "OH+ is needed if O2+ is in traceur.def"
     605               call abort_physic("calchim","missing ohplus tracer",1)
     606            else
     607               nbq = nbq + 1
     608               niq(nbq) = i_ohplus
     609            end if
     610         else
     611            if (i_ohplus /= 0) then
     612               write(*,*) "calchim: Error: OH+ is present, but O2+ is not!!!"
    579613               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
    580614            endif
    … …  
    665699            ! set number of reactions, depending on ion chemistry or not
    666700            if (ionchem) then
    667                nb_reaction_4_max = 80   ! set number of bimolecular reactions
     701               nb_reaction_4_max = 95   ! set number of bimolecular reactions
    668702               nb_reaction_3_max = 6    ! set number of quadratic reactions
    669703               nquench           = 9    ! set number of quenching + heterogeneous reactions
    … …  
    688722                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,       &
    689723                           dist_sol,zday,surfdust1d,surfice1d,            &
    690                            jo3,jh2o,taucol,iter)
     724                           jo3,jh2o,em_no,em_o2,taucol,iter)
    691725
    692726!        ozone photolysis, for output
    … …  
    711745                    zdens,zpress,zlocal,szacol,ptimestep,zday,&
    712746                    em_no,em_o2)
    713                do l = 1,nlayer
    714                   emission_no(ig,l) = em_no(l)
    715                   emission_o2(ig,l) = em_o2(l)
    716                end do
    717             end if
    718 
     747            end if
     748
     749            do l = 1,nlayer
     750               emission_no(ig,l) = em_no(l)
     751               emission_o2(ig,l) = em_o2(l)
     752            end do
    719753         end if  ! photochem
    720754
    … …  
    773807                             ' ',3,iter_3d(1,1))
    774808
    775             if (.not. unichim) then
     809!            if (.not. unichim) then
    776810               call writediagfi(ngrid,'emission_no',        &
    777811                    'NO nightglow emission rate','cm-3 s-1',3,emission_no)
    778812               call writediagfi(ngrid,'emission_o2',        &
    779813                    'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
    780             endif
     814!            endif
    781815           
    782816            if (callstats) then

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

    r2312 r2321  
    8787      igcm_h2o_vap=0
    8888      igcm_h2o_ice=0
    89       igcm_hdo_vap=0
    90       igcm_hdo_ice=0
    9189      igcm_co2=0
    9290      igcm_co=0
    … …  
    121119      igcm_hcoplus=0
    122120      igcm_h2oplus=0
     121      igcm_h3oplus=0
     122      igcm_ohplus=0
    123123      igcm_elec=0
    124124
    … …  
    278278           count = count + 1
    279279        end if
    280         if (noms(iq) == "hdo_vap") then
    281            igcm_hdo_vap = iq
    282            mmol(igcm_hdo_vap) = 18.
    283            count = count + 1
    284         end if
    285         if (noms(iq) == "hdo_ice") then
    286            igcm_hdo_ice = iq
    287            mmol(igcm_hdo_ice) = 18.
    288            count = count + 1
    289         end if
    290280        if (noms(iq).eq."he") then
    291281          igcm_he=iq
    … …  
    354344           igcm_h2oplus = iq
    355345           mmol(igcm_h2oplus) = 18.
     346           count = count + 1
     347        end if
     348        if (noms(iq) == "h3oplus") then
     349           igcm_h3oplus = iq
     350           mmol(igcm_h3oplus) = 19.
     351           count = count + 1
     352        end if
     353        if (noms(iq) == "ohplus") then
     354           igcm_ohplus = iq
     355           mmol(igcm_ohplus) = 17.
    356356           count = count + 1
    357357        end if
    … …  
    573573                  do n = 1,nspe
    574574                     iq = niq(n)
    575                      if (iq /= igcm_h2o_vap .or. iq/= igcm_hdo_vap .or. flagh2o == 1) then
     575                     if (iq /= igcm_h2o_vap .or. flagh2o == 1) then
    576576                        pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l)
    577577                     end if
    … …  
    617617              .or. igcm_cplus == 0 .or. igcm_nplus == 0 .or. igcm_noplus == 0    &
    618618              .or. igcm_n2plus == 0 .or. igcm_hplus == 0 .or. igcm_hco2plus == 0 &
    619               .or. igcm_hcoplus == 0 .or. igcm_h2oplus == 0 .or. igcm_elec == 0) then
     619              .or. igcm_hcoplus == 0 .or. igcm_h2oplus == 0                      &
     620              .or. igcm_h3oplus == 0.or. igcm_ohplus == 0 .or.igcm_elec == 0) then
    620621            write(*,*)'inichim_newstart error:'
    621622            write(*,*)'if co2plus is in traceur.def, all other ions must also be'
    622623            write(*,*)'o2plus, oplus, coplus, cplus, nplus, noplus, n2plus'
    623             write(*,*)'hplus, hco2plus, hcoplus, h2oplus, and elec'
     624            write(*,*)'hplus, hco2plus, hcoplus, h2oplus, h3oplus, ohplus, and elec'
    624625            write(*,*)'stop'
    625626            call abort_physic("inichim_newstart","missing ions in tracers",1)
    … …  
    642643                  pq(i,j,l,igcm_hcoplus)  = 0.
    643644                  pq(i,j,l,igcm_h2oplus)  = 0.
     645                  pq(i,j,l,igcm_h3oplus)  = 0.
     646                  pq(i,j,l,igcm_ohplus)   = 0.
    644647                  pq(i,j,l,igcm_elec)     = 0.
    645648               end do
    … …  
    661664         qsurf(1:ngrid,igcm_hcoplus)  = 0.
    662665         qsurf(1:ngrid,igcm_h2oplus)  = 0.
     666         qsurf(1:ngrid,igcm_h3oplus)  = 0.
     667         qsurf(1:ngrid,igcm_ohplus)   = 0.
    663668         qsurf(1:ngrid,igcm_elec)     = 0.
    664669
    … …  
    668673              .or. igcm_cplus /= 0 .or. igcm_nplus /= 0 .or. igcm_noplus /= 0    &
    669674              .or. igcm_n2plus /= 0 .or. igcm_hplus /= 0 .or. igcm_hco2plus /= 0 &
    670               .or. igcm_hcoplus /= 0 .or. igcm_h2oplus /= 0 .or. igcm_elec /= 0) then
     675              .or. igcm_hcoplus /= 0 .or. igcm_h2oplus /= 0                      &
     676              .or. igcm_h3oplus /= 0 .or. igcm_ohplus /= 0 .or. igcm_elec /= 0) then
    671677            write(*,*)'inichim_newstart error:'
    672678            write(*,*)'some ions are in traceur.def, but not co2plus'

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

    r2302 r2321  
    1818                          alt, temp, temp_elect, dens, zmmean,                 &
    1919                          dist_sol, zday,                                      &
    20                           surfdust1d, surfice1d, jo3, jh2o,tau, iter)
     20                          surfdust1d, surfice1d, jo3, jh2o,em_no,em_o2,        &
     21                          tau, iter)
    2122
    2223use param_v4_h, only: jion
    … …  
    7273real    :: jo3(nlayer)         ! photodissociation rate o3 -> o1d
    7374real    :: jh2o(nlayer)        ! photodissociation rate h2o -> h + oh
     75real    :: em_no(nlayer)       ! NO nightglow volume emission rate
     76real    :: em_o2(nlayer)       ! O2 nightglow volume emission rate
    7477
    7578!===================================================================
    … …  
    114117integer,parameter :: i_hcoplus = 28
    115118integer,parameter :: i_h2oplus = 29
    116 integer,parameter :: i_elec    = 30
     119integer,parameter :: i_h3oplus = 30
     120integer,parameter :: i_ohplus  = 31
     121integer,parameter :: i_elec    = 32
    117122
    118123integer :: ilay
    … …  
    157162               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
    158163               i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
    159                i_h2oplus, i_elec)
     164               i_h2oplus, i_h3oplus, i_ohplus, i_elec)
    160165   firstcall = .false.
    161166end if
    … …  
    171176               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
    172177               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
    173                i_hcoplus, i_h2oplus, i_elec, dens, rm, c)
     178               i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus,      &
     179               i_elec, dens, rm, c)
    174180
    175181!===================================================================
    … …  
    211217         v_phot(:,31)=jion(12,:,1)
    212218      endif
    213 !      write(*,*)'photochemistry/205',c(:,i_co2),ig
    214 !      write(*,*)'photochemistry/206',v_phot(:,3),ig
    215219     
    216220   else ! night
    … …  
    367371           c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+&
    368372           c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)+                &
    369            c(ilev,i_hcoplus)+c(ilev,i_h2oplus)) then
     373           c(ilev,i_hcoplus)+c(ilev,i_h2oplus)+c(ilev,i_h3oplus)+             &
     374           c(ilev,i_ohplus)) then
    370375         c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
    371376              c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+              &
    372377              c(ilev,i_n2plus)+c(ilev,i_nplus)+c(ilev,i_hplus)+               &
    373               c(ilev,i_hco2plus)+c(ilev,i_hcoplus)+c(ilev,i_h2oplus)
     378              c(ilev,i_hco2plus)+c(ilev,i_hcoplus)+c(ilev,i_h2oplus)+         &
     379              c(ilev,i_h3oplus)+c(ilev,i_ohplus)
    374380         !      write(*,*)'photochemistry/359'
    375381         !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
    376382      end if
    377383   end if
     384
    378385   cnew(:)   = 0.
    379386
    … …  
    386393
    387394end do ! ilev
     395
     396!add calculation of NO and O2 nightglows
     397em_no(:)=c(:,i_o)*c(:,i_n)*v_4(:,26)   !2.8e-17*(300./temp(:)))**0.5
     398em_o2(:)=0.75*c(:,i_o)*c(:,i_o)*c(:,i_co2)*v_3(:,1)   !2.5*9.46e-34*exp(485./temp(:))*dens(:)
    388399
    389400!===================================================================
    … …  
    397408               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
    398409               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
    399                i_hcoplus, i_h2oplus, i_elec, dens, c)
     410               i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus,      &
     411               i_elec, dens, c)
    400412contains
    401413
    … …  
    599611                           i034, i035, i036, i037, i038, i039, i040,   &
    600612                           i041, i042, i043, i044, i045, i046, i047,   &
    601                            i048, i049,                                 &
     613                           i048, i049, i050, i051, i052, i053, i054,   &
     614                           i055, i056, i057, i058, i059, i060, i061,   &
     615                           i062, i063,                                 &
    602616                           h001, h002, h003, h004, h005
    603617
    … …  
    15211535         nb_reaction_4 = nb_reaction_4 + 1
    15221536         v_4(:,nb_reaction_4) = i049(:)
     1537
     1538!---     i050: H2O+ + H2O -> H3O+ + OH
     1539
     1540         !UMIST, from Huntress, W.T. and Pinizzotto, R.F., J. Chem. Phys., 59, 4742 (1973)
     1541
     1542         i050(:) = 2.1e-9*((300./t_elect(:))**0.5)
     1543         nb_reaction_4 = nb_reaction_4 + 1
     1544         v_4(:,nb_reaction_4) = i050(:)
     1545
     1546
     1547!---     i051: H2O+ + H2 -> H3O+ + H
     1548
     1549         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
     1550
     1551         i051(:) = 6.4e-10
     1552         nb_reaction_4 = nb_reaction_4 + 1
     1553         v_4(:,nb_reaction_4) = i051(:)
     1554
     1555!---     i052: HCO+ + H2O -> H3O+ + CO
     1556
     1557         !UMIST, from Adams, N.G., Smith, D., and Grief, D., Int. J. Mass Spectrom. Ion Phys., 26, 405 (1978)
     1558
     1559         i052(:) = 2.5e-9*((300./t_elect(:))**0.5)
     1560         nb_reaction_4 = nb_reaction_4 + 1
     1561         v_4(:,nb_reaction_4) = i052(:)
     1562
     1563!---     i053: H3O+ + e -> OH + H + H
     1564
     1565         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     1566
     1567         i053(:) = 3.05e-7*((300./t_elect(:))**0.5)
     1568         nb_reaction_4 = nb_reaction_4 + 1
     1569         v_4(:,nb_reaction_4) = i053(:)
     1570
     1571!---     i054: H3O + + e -> H2O + H
     1572
     1573         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     1574         
     1575         i054(:) = 7.09e-8*((300./t_elect(:))**0.5)
     1576         nb_reaction_4 = nb_reaction_4 + 1
     1577         v_4(:,nb_reaction_4) = i054(:)
     1578
     1579!---     i055: H3O+ + e -> OH + H2
     1580
     1581         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     1582
     1583         i055(:) = 5.37e-8*((300./t_elect(:))**0.5)
     1584         nb_reaction_4 = nb_reaction_4 + 1
     1585         v_4(:,nb_reaction_4) = i055(:)
     1586
     1587!---     i056: H3O+ + e -> O + H2 + H
     1588
     1589         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     1590
     1591         i056(:) = 5.6e-9*((300./t_elect(:))**0.5)
     1592         nb_reaction_4 = nb_reaction_4 + 1
     1593         v_4(:,nb_reaction_4) = i056(:)
     1594
     1595         nb_reaction_4 = nb_reaction_4 + 1
     1596         v_4(:,nb_reaction_4) = i056(:)
     1597
     1598!---     i057: O+ + H2 -> OH+ + H
     1599
     1600         !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
     1601
     1602         i057(:) = 1.7e-9
     1603         nb_reaction_4 = nb_reaction_4 + 1
     1604         v_4(:,nb_reaction_4) = i057(:)
     1605
     1606!---     i058: OH+ + O -> O2+ + H
     1607
     1608         !UMIST, from Prasad & Huntress, 1980, ApJS, 43, 1
     1609
     1610         i058(:) = 7.1e-10
     1611         nb_reaction_4 = nb_reaction_4 + 1
     1612         v_4(:,nb_reaction_4) = i058(:)
     1613
     1614!---     i059: OH+ + CO2 -> HCO2+ + O
     1615
     1616         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     1617
     1618         i059(:) = 1.44e-9
     1619         nb_reaction_4 = nb_reaction_4 + 1
     1620         v_4(:,nb_reaction_4) = i059(:)
     1621
     1622!---     i060: OH+ + CO -> HCO+ + O
     1623
     1624         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     1625
     1626         i060(:) = 1.05e-9
     1627         nb_reaction_4 = nb_reaction_4 + 1
     1628         v_4(:,nb_reaction_4) = i060(:)
     1629
     1630!---     i061: OH+ + NO -> NO+ + OH (tasa de reacción UMIST 3.59e-10)
     1631
     1632         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     1633
     1634         i061(:) = 3.59e-10
     1635         nb_reaction_4 = nb_reaction_4 + 1
     1636         v_4(:,nb_reaction_4) = i061(:)
     1637
     1638!---     i062: OH+ + H2 -> H2O+ + H (tasa de reacción UMIST 1.01e-9,
     1639
     1640         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     1641
     1642         i062(:) = 1.01e-9
     1643         nb_reaction_4 = nb_reaction_4 + 1
     1644         v_4(:,nb_reaction_4) = i062(:)
     1645
     1646!---     i063: OH+ + O2 -> O2+ + OH (tasa de reacción UMIST 5.9e-10
     1647
     1648         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     1649
     1650         i063(:) = 5.9e-10
     1651         nb_reaction_4 = nb_reaction_4 + 1
     1652         v_4(:,nb_reaction_4) = i063(:)
    15231653
    15241654      end if   !ionchem
    … …  
    17341864                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
    17351865                   i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
    1736                    i_h2oplus, i_elec)
     1866                   i_h2oplus, i_h3oplus, i_ohplus, i_elec)
    17371867
    17381868!================================================================
    … …  
    17591889           i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
    17601890           i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, &
    1761            i_hcoplus, i_h2oplus, i_elec
     1891           i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec
    17621892integer, intent(in) :: nb_reaction_3_max
    17631893                       ! number of quadratic reactions
    … …  
    27442874   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o)
    27452875
     2876!===========================================================
     2877!      i050 : H2O+ + H2O -> H3O+ + OH
     2878!===========================================================
     2879
     2880   nb_reaction_4 = nb_reaction_4 + 1
     2881   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_oh)
     2882
     2883!===========================================================
     2884!      i051 : H2O+ + H2 -> H3O+ + H
     2885!===========================================================
     2886
     2887   nb_reaction_4 = nb_reaction_4 + 1
     2888   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2, 1.0, i_h3oplus, 1.0, i_h)
     2889
     2890!===========================================================
     2891!      i052 : HCO+ + H2O -> H3O+ + CO
     2892!===========================================================
     2893
     2894   nb_reaction_4 = nb_reaction_4 + 1
     2895   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_co)
     2896
     2897!===========================================================
     2898!      i053: H3O+ + e -> OH + H + H
     2899!===========================================================
     2900
     2901   nb_reaction_4 = nb_reaction_4 + 1
     2902   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 2.0, i_h)
     2903
     2904!===========================================================
     2905!      i054: H3O+ + e -> H2O + H
     2906!===========================================================
     2907
     2908   nb_reaction_4 = nb_reaction_4 + 1
     2909   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_h2o, 1.0, i_h)
     2910
     2911!===========================================================
     2912!      i055: H3O+ + e -> HO + H2
     2913!===========================================================
     2914
     2915   nb_reaction_4 = nb_reaction_4 + 1
     2916   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 1.0, i_h2)
     2917
     2918!===========================================================
     2919!      i056: H3O+ + e -> O + H2 + H
     2920!===========================================================
     2921!We divide this reaction in two
     2922
     2923!0.5H3O+ + 0.5e -> O
     2924
     2925   nb_reaction_4 = nb_reaction_4 + 1
     2926   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_o, 0.0, i_dummy)
     2927
     2928!0.5H3O+ + 0.5e -> H2 + H
     2929
     2930   nb_reaction_4 = nb_reaction_4 + 1
     2931   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_h2, 1.0, i_h)
     2932
     2933!===========================================================
     2934!      i057: O+ + H2 -> OH+ + H
     2935!===========================================================
     2936
     2937   nb_reaction_4 = nb_reaction_4 + 1
     2938   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2, 1.0, i_ohplus, 1.0, i_h)
     2939
     2940!===========================================================
     2941!      i058: OH+ + O -> O2+ + H
     2942!===========================================================
     2943
     2944   nb_reaction_4 = nb_reaction_4 + 1
     2945   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h)
     2946
     2947!===========================================================
     2948!      i059: OH+ + CO2 -> HCO2+ + O
     2949!===========================================================
     2950
     2951   nb_reaction_4 = nb_reaction_4 + 1
     2952   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co2, 1.0, i_hco2plus, 1.0, i_o)
     2953
     2954!===========================================================
     2955!      i060: OH+ + CO -> HCO+ + O
     2956!===========================================================
     2957
     2958   nb_reaction_4 = nb_reaction_4 + 1
     2959   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_o)
     2960
     2961!===========================================================
     2962!      i061: OH+ + NO -> NO+ + OH
     2963!===========================================================
     2964
     2965   nb_reaction_4 = nb_reaction_4 + 1
     2966   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_oh)
     2967
     2968!===========================================================
     2969!      i062: OH+ + H2 -> H2O+ + H
     2970!===========================================================
     2971
     2972   nb_reaction_4 = nb_reaction_4 + 1
     2973   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_h2, 1.0, i_h2oplus, 1.0, i_h)
     2974
     2975!===========================================================
     2976!      i063: OH+ + O2 -> O2+ + OH
     2977!===========================================================
     2978
     2979   nb_reaction_4 = nb_reaction_4 + 1
     2980   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_oh)
     2981
    27462982end if    !ionchem
    27472983
    … …  
    28223058                           i_coplus, i_cplus, i_n2plus, i_nplus,     &
    28233059                           i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,&
    2824                            i_elec, dens, rm, c)
     3060                           i_h3oplus, i_ohplus, i_elec, dens, rm, c)
    28253061       
    28263062!*****************************************************************
    … …  
    28343070     &                      igcm_n2plus, igcm_nplus, igcm_hplus,         &
    28353071     &                      igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,   &
    2836      &                      igcm_elec
     3072     &                      igcm_h3oplus, igcm_ohplus, igcm_elec
    28373073
    28383074      implicit none
    … …  
    28553091                 i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus,   &
    28563092                 i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus,    &
    2857                  i_hcoplus, i_h2oplus, i_elec
     3093                 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec
    28583094
    28593095      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
    … …  
    29973233               call abort_physic("gcmtochim","missing h2oplus tracer",1)
    29983234            endif
     3235            if (igcm_h3oplus == 0) then
     3236               write(*,*) "gcmtochim: Error; no H3O+ tracer !!!"
     3237               call abort_physic("gcmtochim","missing h3oplus tracer",1)
     3238            endif
     3239            if (igcm_ohplus == 0) then
     3240               write(*,*) "gcmtochim: Error; no OH+ tracer !!!"
     3241               call abort_physic("gcmtochim","missing ohplus tracer",1)
     3242            endif
    29993243            if (igcm_elec == 0) then
    30003244               write(*,*) "gcmtochim: Error; no e- tracer !!!"
    … …  
    30433287            rm(l,i_hcoplus)  = zycol(l, igcm_hcoplus)
    30443288            rm(l,i_h2oplus)  = zycol(l, igcm_h2oplus)
     3289            rm(l,i_h3oplus)  = zycol(l, igcm_h3oplus)
     3290            rm(l,i_ohplus)   = zycol(l, igcm_ohplus)
    30453291            rm(l,i_elec)     = zycol(l, igcm_elec)
    30463292         end do
    … …  
    30723318                           i_coplus, i_cplus, i_n2plus, i_nplus,      &
    30733319                           i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, &
    3074                            i_elec, dens, c)
     3320                           i_h3oplus, i_ohplus, i_elec, dens, c)
    30753321 
    30763322!*****************************************************************
    … …  
    30843330                            igcm_n2plus, igcm_nplus, igcm_hplus,          &
    30853331                            igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
    3086                             igcm_elec
     3332                            igcm_h3oplus, igcm_ohplus, igcm_elec
    30873333
    30883334      implicit none
    … …  
    31043350                 i_co2plus, i_oplus, i_o2plus, i_noplus,         &
    31053351                 i_coplus, i_cplus, i_n2plus, i_nplus,           &
    3106                  i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, i_elec
     3352                 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,      &
     3353                 i_h3oplus, i_ohplus, i_elec
    31073354
    31083355      real :: dens(nlayer)     ! total number density (molecule.cm-3)
    … …  
    31593406            zycol(l, igcm_hcoplus) = c(l,i_hcoplus)/dens(l)
    31603407            zycol(l, igcm_h2oplus) = c(l,i_h2oplus)/dens(l)
     3408            zycol(l, igcm_h3oplus) = c(l,i_h3oplus)/dens(l)
     3409            zycol(l, igcm_ohplus)  = c(l,i_ohplus)/dens(l)
    31613410            zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
    31623411         end do
Note: See TracChangeset for help on using the changeset viewer.