Ignore:
Timestamp:
Apr 28, 2020, 9:53:07 AM (5 years ago)
Author:
emillour
Message:

Mars GCM:

  • Add H2O+ ion chemistry

FGG

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

Legend:

Unmodified
Added
Removed

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

    r2284 r2302  
    2222                            igcm_coplus, igcm_cplus, igcm_nplus,          &
    2323                            igcm_noplus, igcm_n2plus, igcm_hplus,         &
    24                             igcm_hco2plus, igcm_hcoplus, igcm_elec, mmol
     24                            igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
     25                            igcm_elec, mmol
    2526
    2627      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
    … …  
    142143      integer,save :: i_hco2plus=0
    143144      integer,save :: i_hcoplus=0
     145      integer,save :: i_h2oplus=0
    144146      integer,save :: i_elec=0
    145147
    … …  
    240242         if (i_co2 == 0) then
    241243            write(*,*) "calchim: Error; no CO2 tracer !!!"
    242             stop
     244            call abort_physic("calchim","missing co2 tracer",1)
    243245         else
    244246            nbq = nbq + 1
    … …  
    248250         if (i_co == 0) then
    249251            write(*,*) "calchim: Error; no CO tracer !!!"
    250             stop
     252            call abort_physic("calchim","missing co tracer",1)
    251253         else
    252254            nbq = nbq + 1
    … …  
    256258         if (i_o == 0) then
    257259            write(*,*) "calchim: Error; no O tracer !!!"
    258             stop
     260            call abort_physic("calchim","missing o tracer",1)
    259261         else
    260262            nbq = nbq + 1
    … …  
    264266         if (i_o1d == 0) then
    265267            write(*,*) "calchim: Error; no O1D tracer !!!"
    266             stop
     268            call abort_physic("calchim","missing o1d tracer",1)
    267269         else
    268270            nbq = nbq + 1
    … …  
    272274         if (i_o2 == 0) then
    273275            write(*,*) "calchim: Error; no O2 tracer !!!"
    274             stop
     276            call abort_physic("calchim","missing o2 tracer",1)
    275277         else
    276278            nbq = nbq + 1
    … …  
    280282         if (i_o3 == 0) then
    281283            write(*,*) "calchim: Error; no O3 tracer !!!"
    282             stop
     284            call abort_physic("calchim","missing o3 tracer",1)
    283285         else
    284286            nbq = nbq + 1
    … …  
    288290         if (i_h == 0) then
    289291            write(*,*) "calchim: Error; no H tracer !!!"
    290             stop
     292            call abort_physic("calchim","missing h tracer",1)
    291293         else
    292294            nbq = nbq + 1
    … …  
    296298         if (i_h2 == 0) then
    297299            write(*,*) "calchim: Error; no H2 tracer !!!"
    298             stop
     300            call abort_physic("calchim","missing h2 tracer",1)
    299301         else
    300302            nbq = nbq + 1
    … …  
    304306         if (i_oh == 0) then
    305307            write(*,*) "calchim: Error; no OH tracer !!!"
    306             stop
     308            call abort_physic("calchim","missing oh tracer",1)
    307309         else
    308310            nbq = nbq + 1
    … …  
    312314         if (i_ho2 == 0) then
    313315            write(*,*) "calchim: Error; no HO2 tracer !!!"
    314             stop
     316            call abort_physic("calchim","missing ho2 tracer",1)
    315317         else
    316318            nbq = nbq + 1
    … …  
    320322         if (i_h2o2 == 0) then
    321323            write(*,*) "calchim: Error; no H2O2 tracer !!!"
    322             stop
     324            call abort_physic("calchim","missing h2o2 tracer",1)
    323325         else
    324326            nbq = nbq + 1
    … …  
    336338         if (i_n2 == 0) then
    337339            write(*,*) "calchim: Error; no N2 tracer !!!"
    338             stop
     340            call abort_physic("calchim","missing n2 tracer",1)
    339341         else
    340342            nbq = nbq + 1
    … …  
    345347            if (photochem) then
    346348               write(*,*) "calchim: Error; no N tracer !!!"
    347                stop
     349               call abort_physic("calchim","missing n tracer",1)
    348350            end if
    349351         else
    … …  
    355357            if (photochem) then
    356358               write(*,*) "calchim: Error; no N2D tracer !!!"
    357                stop
     359               call abort_physic("calchim","missing n2d tracer",1)
    358360            end if
    359361         else
    … …  
    365367            if (photochem) then
    366368               write(*,*) "calchim: Error; no NO tracer !!!"
    367                stop
     369               call abort_physic("calchim","missing no tracer",1)
    368370            end if
    369371         else
    … …  
    375377            if (photochem) then
    376378               write(*,*) "calchim: Error; no NO2 tracer !!!"
    377                stop
     379               call abort_physic("calchim","missing no2 tracer",1)
    378380            end if
    379381         else
    … …  
    384386         if (i_h2o == 0) then
    385387            write(*,*) "calchim: Error; no water vapor tracer !!!"
    386             stop
     388            call abort_physic("calchim","missing h2o_vap tracer",1)
    387389         else
    388390            nbq = nbq + 1
    … …  
    405407               write(*,*) "calchim: Error, no CO2+ tracer !!!"
    406408               write(*,*) "CO2+ is needed if O2+ is in traceur.def"
    407                stop
     409               call abort_physic("calchim","missing co2plus tracer",1)
    408410            else
    409411               nbq = nbq + 1
    … …  
    414416               write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
    415417               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
    416                stop
     418               call abort_physic("calchim","missing o2plus tracer",1)
    417419            endif
    418420         endif
    … …  
    422424               write(*,*) "calchim: Error, no O+ tracer !!!"
    423425               write(*,*) "O+ is needed if O2+ is in traceur.def"
    424                stop
     426               call abort_physic("calchim","missing oplus tracer",1)
    425427            else
    426428               nbq = nbq + 1
    … …  
    431433               write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
    432434               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
    433                stop
     435               call abort_physic("calchim","missing o2plus tracer",1)
    434436            endif
    435437         endif
    … …  
    439441               write(*,*) "calchim: Error, no NO+ tracer !!!"
    440442               write(*,*) "NO+ is needed if O2+ is in traceur.def"
    441                stop
     443               call abort_physic("calchim","missing noplus tracer",1)
    442444            else
    443445               nbq = nbq + 1
    … …  
    455457               write(*,*) "calchim: Error, no CO+ tracer !!!"
    456458               write(*,*) "CO+ is needed if O2+ is in traceur.def"
    457                stop               
     459               call abort_physic("calchim","missing coplus tracer",1)               
    458460            else
    459461               nbq = nbq + 1
    … …  
    471473               write(*,*) "calchim: Error, no C+ tracer !!!"
    472474               write(*,*) "C+ is needed if O2+ is in traceur.def"
    473                stop
     475               call abort_physic("calchim","missing cplus tracer",1)
    474476            else
    475477               nbq = nbq + 1
    … …  
    487489               write(*,*) "calchim: Error, no N2+ tracer !!!"
    488490               write(*,*) "N2+ is needed if O2+ is in traceur.def"
    489                stop
     491               call abort_physic("calchim","missing n2plus tracer",1)
    490492            else
    491493               nbq = nbq + 1
    … …  
    503505               write(*,*) "calchim: Error, no N+ tracer !!!"
    504506               write(*,*) "N+ is needed if O2+ is in traceur.def"
    505                stop
     507               call abort_physic("calchim","missing nplus tracer",1)
    506508            else
    507509               nbq = nbq + 1
    … …  
    519521               write(*,*) "calchim: Error, no H+ tracer !!!"
    520522               write(*,*) "H+ is needed if O2+ is in traceur.def"
    521                stop
     523               call abort_physic("calchim","missing hplus tracer",1)
    522524            else
    523525               nbq = nbq + 1
    … …  
    535537               write(*,*) "calchim: Error, no HCO2+ tracer !!!"
    536538               write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
    537                stop
     539               call abort_physic("calchim","missing hco2plus tracer",1)
    538540            else
    539541               nbq = nbq + 1
    … …  
    551553               write(*,*) "calchim: Error, no HCO+ tracer !!!"
    552554               write(*,*) "HCO+ is needed if O2+ is in traceur.def"
    553                stop
     555               call abort_physic("calchim","missing hcoplus tracer",1)
    554556            else
    555557               nbq = nbq + 1
    … …  
    559561            if (i_hcoplus /= 0) then
    560562               write(*,*) "calchim: Error: HCO+ is present, but O2+ is not!!!"
     563               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
     564            endif
     565         endif
     566         i_h2oplus=igcm_h2oplus
     567         if(ionchem) then
     568            if (i_h2oplus == 0) then
     569               write(*,*) "calchim: Error, no H2O+ tracer !!!"
     570               write(*,*) "H2O+ is needed if O2+ is in traceur.def"
     571               call abort_physic("calchim","missing h2oplus tracer",1)
     572            else
     573               nbq = nbq + 1
     574               niq(nbq) = i_h2oplus
     575            end if
     576         else
     577            if (i_h2oplus /= 0) then
     578               write(*,*) "calchim: Error: H2O+ is present, but O2+ is not!!!"
    561579               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
    562580            endif
    … …  
    567585               write(*,*) "calchim: Error, no e- tracer !!!"
    568586               write(*,*) "e- is needed if O2+ is in traceur.def"
    569                stop
     587               call abort_physic("calchim","missing elec tracer",1)
    570588            else
    571589               nbq = nbq + 1
    … …  
    647665            ! set number of reactions, depending on ion chemistry or not
    648666            if (ionchem) then
    649                nb_reaction_4_max = 67   ! set number of bimolecular reactions
     667               nb_reaction_4_max = 80   ! set number of bimolecular reactions
    650668               nb_reaction_3_max = 6    ! set number of quadratic reactions
    651669               nquench           = 9    ! set number of quenching + heterogeneous reactions

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

    r2284 r2302  
    118118      igcm_hco2plus=0
    119119      igcm_hcoplus=0
     120      igcm_h2oplus=0
    120121      igcm_elec=0
    121122
    … …  
    338339           count = count + 1
    339340        end if
     341        if (noms(iq) == "h2oplus") then
     342           igcm_h2oplus = iq
     343           mmol(igcm_h2oplus) = 18.
     344           count = count + 1
     345        end if
    340346        if (noms(iq) == "elec") then
    341347           igcm_elec = iq
    … …  
    362368            write(*,*) '      ', iq, ' ', trim(noms(iq))
    363369         end do
    364          stop
     370         call abort_physic("inichim_newstart","tracer mismatch",1)
    365371      else
    366372         write(*,*) "inichim_newstart: found all expected tracers"
    … …  
    378384            write(*,*)'N, NO2 and/or N2D are in traceur.def, but not NO'
    379385            write(*,*)'stop'
    380             stop
     386            call abort_physic("inichim_newstart","missing no tracer",1)
    381387         endif
    382388         flagnitro = .false.
    … …  
    388394            write(*,*)'if NO is in traceur.def, N, NO2 and N2D must also be'
    389395            write(*,*)'stop'
    390             stop
     396            call abort_physic("inichim_newstart","missing n* tracer",1)
    391397         endif
    392398         flagnitro = .true.
    … …  
    470476         write(*,*)'   can be obtained online on:'
    471477         write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
    472          stop
     478         call abort_physic("inichim_newstart","missing input file",1)
    473479      end if
    474480      open(220, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_min.dat')
    … …  
    482488         write(*,*)'   can be obtained online on:'
    483489         write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
    484          stop
     490         call abort_physic("inichim_newstart","missing input file",1)
    485491      end if
    486492      if(flagnitro) then
    … …  
    495501            write(*,*)'   can be obtained online on:'
    496502            write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
    497             STOP
     503            call abort_physic("inichim_newstart","missing input file",1)
    498504         endif
    499505      endif   ! Of if(flagnitro)
    … …  
    599605              .or. igcm_cplus == 0 .or. igcm_nplus == 0 .or. igcm_noplus == 0    &
    600606              .or. igcm_n2plus == 0 .or. igcm_hplus == 0 .or. igcm_hco2plus == 0 &
    601               .or. igcm_hcoplus == 0 .or. igcm_elec == 0) then
     607              .or. igcm_hcoplus == 0 .or. igcm_h2oplus == 0 .or. igcm_elec == 0) then
    602608            write(*,*)'inichim_newstart error:'
    603609            write(*,*)'if co2plus is in traceur.def, all other ions must also be'
    604610            write(*,*)'o2plus, oplus, coplus, cplus, nplus, noplus, n2plus'
    605             write(*,*)'hplus, hco2plus, hcoplus, and elec'
     611            write(*,*)'hplus, hco2plus, hcoplus, h2oplus, and elec'
    606612            write(*,*)'stop'
    607             stop
     613            call abort_physic("inichim_newstart","missing ions in tracers",1)
    608614         end if
    609615
    … …  
    623629                  pq(i,j,l,igcm_hco2plus) = 0.
    624630                  pq(i,j,l,igcm_hcoplus)  = 0.
     631                  pq(i,j,l,igcm_h2oplus)  = 0.
    625632                  pq(i,j,l,igcm_elec)     = 0.
    626633               end do
    … …  
    641648         qsurf(1:ngrid,igcm_hco2plus) = 0.
    642649         qsurf(1:ngrid,igcm_hcoplus)  = 0.
     650         qsurf(1:ngrid,igcm_h2oplus)  = 0.
    643651         qsurf(1:ngrid,igcm_elec)     = 0.
    644652
    … …  
    648656              .or. igcm_cplus /= 0 .or. igcm_nplus /= 0 .or. igcm_noplus /= 0    &
    649657              .or. igcm_n2plus /= 0 .or. igcm_hplus /= 0 .or. igcm_hco2plus /= 0 &
    650               .or. igcm_hcoplus /= 0 .or. igcm_elec /= 0) then
     658              .or. igcm_hcoplus /= 0 .or. igcm_h2oplus /= 0 .or. igcm_elec /= 0) then
    651659            write(*,*)'inichim_newstart error:'
    652660            write(*,*)'some ions are in traceur.def, but not co2plus'
    653661            write(*,*)'stop'
    654             stop
     662            call abort_physic("inichim_newstart","missing ions in tracers",1)
    655663         end if
    656664      end if    ! of if(igcm_co2 /= 0)

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

    r2284 r2302  
    2424implicit none
    2525
    26 #include "callkeys.h"
     26include "callkeys.h"
    2727
    2828!===================================================================
    … …  
    113113integer,parameter :: i_hco2plus= 27
    114114integer,parameter :: i_hcoplus = 28
    115 integer,parameter :: i_elec    = 29
     115integer,parameter :: i_h2oplus = 29
     116integer,parameter :: i_elec    = 30
    116117
    117118integer :: ilay
    … …  
    156157               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
    157158               i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
    158                i_elec)
     159               i_h2oplus, i_elec)
    159160   firstcall = .false.
    160161end if
    … …  
    170171               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
    171172               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
    172                i_hcoplus, i_elec, dens, rm, c)
     173               i_hcoplus, i_h2oplus, i_elec, dens, rm, c)
    173174
    174175!===================================================================
    … …  
    344345#else
    345346   write(*,*) "photochemistry error, missing LAPACK routine dgesv"
    346    stop
     347   call abort_physic("photochemistry","missing LAPACK routine",1)
    347348#endif
    348349
    … …  
    366367           c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+c(ilev,i_n2plus)+&
    367368           c(ilev,i_nplus)+c(ilev,i_hplus)+c(ilev,i_hco2plus)+                &
    368            c(ilev,i_hcoplus)) then
     369           c(ilev,i_hcoplus)+c(ilev,i_h2oplus)) then
    369370         c(ilev,i_elec) = c(ilev,i_co2plus)+c(ilev,i_oplus)+c(ilev,i_o2plus)+ &
    370371              c(ilev,i_noplus)+c(ilev,i_coplus)+c(ilev,i_cplus)+              &
    371372              c(ilev,i_n2plus)+c(ilev,i_nplus)+c(ilev,i_hplus)+               &
    372               c(ilev,i_hco2plus)+c(ilev,i_hcoplus)
     373              c(ilev,i_hco2plus)+c(ilev,i_hcoplus)+c(ilev,i_h2oplus)
    373374         !      write(*,*)'photochemistry/359'
    374375         !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
    … …  
    396397               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
    397398               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
    398                i_hcoplus, i_elec, dens, c)
     399               i_hcoplus, i_h2oplus, i_elec, dens, c)
    399400contains
    400401
    … …  
    481482#else
    482483   write(*,*) "photochemistry error, missing LAPACK routine dgesv"
    483    stop
     484   call abort_physic("photochemistry","missing LAPACK routine",1)
    484485#endif
    485486
    … …  
    596597                           i020, i021, i022, i023, i024, i025, i026,   &
    597598                           i027, i028, i029, i030, i031, i032, i033,   &
    598                            i034, i035, i036,                           &
     599                           i034, i035, i036, i037, i038, i039, i040,   &
     600                           i041, i042, i043, i044, i045, i046, i047,   &
     601                           i048, i049,                                 &
    599602                           h001, h002, h003, h004, h005
    600603
    … …  
    14131416         nb_reaction_4 = nb_reaction_4 + 1
    14141417         v_4(:,nb_reaction_4) = i036(:)
     1418
     1419!---     i037: CO2+ + H2O -> H2O+ + CO2
     1420
     1421         !UMIST, from Karpas, Z., Anicich, V.G., and Huntress, W.T., Chem. Phys. Lett., 59, 84 (1978)
     1422
     1423         i037(:) = 2.04e-9 *((300./t_elect(:))**0.5)
     1424         nb_reaction_4 = nb_reaction_4 + 1
     1425         v_4(:,nb_reaction_4) = i037(:)
     1426
     1427!---     i038: CO+ + H2O -> H2O+ + CO
     1428
     1429         !UMIST, from Huntress, W.T., McEwan, M.J., Karpas, Z., and Anicich, V.G., Astrophys. J. Supp. Series, 44, 481 (1980)
     1430
     1431         i038(:) = 1.72e-9*((300./t_elect(:))**0.5)
     1432         nb_reaction_4 = nb_reaction_4 + 1
     1433         v_4(:,nb_reaction_4) = i038(:)
     1434
     1435!---     i039: O+ + H2O -> H2O+ + O
     1436
     1437         !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)
     1438
     1439         i039(:) = 3.2e-9*((300./t_elect(:))**0.5)
     1440         nb_reaction_4 = nb_reaction_4 + 1
     1441         v_4(:,nb_reaction_4) = i039(:)
     1442
     1443!---     i040: N2+ + H2O -> H2O+ + N2
     1444
     1445         !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)
     1446
     1447         i040(:) = 2.3e-9*((300./t_elect(:))**0.5)
     1448         nb_reaction_4 = nb_reaction_4 + 1
     1449         v_4(:,nb_reaction_4) = i040(:)
     1450
     1451!---     i041: N+ + H2O -> H2O+ + N
     1452
     1453         !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)
     1454
     1455         i041(:) = 2.8e-9*((300./t_elect(:))**0.5)
     1456         nb_reaction_4 = nb_reaction_4 + 1
     1457         v_4(:,nb_reaction_4) = i041(:)
     1458
     1459
     1460!---     i042: H+ + H2O -> H2O+ + H
     1461
     1462         !UMIST, from D. Smith, P. Spanel and C. A. Mayhew, Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
     1463
     1464         i042(:) = 6.9e-9*((300./t_elect(:))**0.5)
     1465         nb_reaction_4 = nb_reaction_4 + 1
     1466         v_4(:,nb_reaction_4) = i042(:)
     1467
     1468!---     i043: H2O+ + O2 -> O2+ + H2O
     1469
     1470         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
     1471
     1472         i043(:) = 4.6e-10
     1473         nb_reaction_4 = nb_reaction_4 + 1
     1474         v_4(:,nb_reaction_4) = i043(:)
     1475
     1476!---     i044: H2O+ + CO -> HCO+ + OH
     1477
     1478         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     1479
     1480         i044(:) = 5.0e-10
     1481         nb_reaction_4 = nb_reaction_4 + 1
     1482         v_4(:,nb_reaction_4) = i044(:)
     1483
     1484!---     i045: H2O+ + O -> O2+ + H2
     1485
     1486         !UMIST, from Viggiano, A.A, Howarka, F., Albritton, D.L., Fehsenfeld, F.C., Adams, N.G., and Smith, D., Astrophys. J., 236, 492 (1980)
     1487
     1488         i045(:) = 4.0e-11
     1489         nb_reaction_4 = nb_reaction_4 + 1
     1490         v_4(:,nb_reaction_4) = i045(:)
     1491
     1492!---     i046: H2O+ + NO -> NO+ + H2O
     1493
     1494         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
     1495
     1496         i046(:) = 2.7e-10
     1497         nb_reaction_4 = nb_reaction_4 + 1
     1498         v_4(:,nb_reaction_4) = i046(:)
     1499
     1500!---     i047: H2O+ + e- -> H + H + O
     1501
     1502         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
     1503         
     1504         i047(:) = 3.05e-7*((300./t_elect(:))**0.5)
     1505         nb_reaction_4 = nb_reaction_4 + 1
     1506         v_4(:,nb_reaction_4) = i047(:)
     1507
     1508!---     i048: H2O+ + e- -> H + OH
     1509         
     1510         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
     1511
     1512         i048(:) = 8.6e-8*((300./t_elect(:))**0.5)
     1513         nb_reaction_4 = nb_reaction_4 + 1
     1514         v_4(:,nb_reaction_4) = i048(:)
     1515
     1516!---     i049: H2O+ + e- -> O + H2
     1517
     1518         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
     1519
     1520         i049(:) = 3.9e-8*((300./t_elect(:))**0.5)
     1521         nb_reaction_4 = nb_reaction_4 + 1
     1522         v_4(:,nb_reaction_4) = i049(:)
    14151523
    14161524      end if   !ionchem
    … …  
    16261734                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus,      &
    16271735                   i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus,   &
    1628                    i_elec)
     1736                   i_h2oplus, i_elec)
    16291737
    16301738!================================================================
    … …  
    16511759           i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
    16521760           i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, &
    1653            i_hcoplus, i_elec
     1761           i_hcoplus, i_h2oplus, i_elec
    16541762integer, intent(in) :: nb_reaction_3_max
    16551763                       ! number of quadratic reactions
    … …  
    25452653   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_elec, 1.0, i_co, 1.0, i_h)
    25462654
     2655!===========================================================
     2656!      i037 : CO2+ + H2O -> H2O+ + CO2
     2657!===========================================================
     2658
     2659   nb_reaction_4 = nb_reaction_4 + 1
     2660   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co2)
     2661
     2662!===========================================================
     2663!      i038 : CO+ + H2O -> H2O+ + CO
     2664!===========================================================
     2665
     2666   nb_reaction_4 = nb_reaction_4 + 1
     2667   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co)
     2668
     2669!===========================================================
     2670!      i039 : O+ + H2O -> H2O+ + O
     2671!===========================================================
     2672
     2673   nb_reaction_4 = nb_reaction_4 + 1
     2674   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_o)
     2675
     2676!===========================================================
     2677!      i040 : N2+ + H2O -> H2O+ + N2
     2678!===========================================================
     2679
     2680   nb_reaction_4 = nb_reaction_4 + 1
     2681   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n2)
     2682
     2683!===========================================================
     2684!      i041 : N+ + H2O -> H2O+ + N
     2685!===========================================================
     2686
     2687   nb_reaction_4 = nb_reaction_4 + 1
     2688   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n)
     2689
     2690!===========================================================
     2691!      i042 : H+ + H2O -> H2O+ + H
     2692!===========================================================
     2693
     2694   nb_reaction_4 = nb_reaction_4 + 1
     2695   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_h)
     2696
     2697!===========================================================
     2698!      i043 : H2O+ + O2 -> O2+ + H2O
     2699!===========================================================
     2700
     2701   nb_reaction_4 = nb_reaction_4 + 1
     2702   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_h2o)
     2703
     2704!===========================================================
     2705!      i044 : H2O+ + CO -> HCO+ + OH
     2706!===========================================================
     2707
     2708   nb_reaction_4 = nb_reaction_4 + 1
     2709   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_oh)
     2710
     2711!===========================================================
     2712!      i045 : H2O+ + O -> O2+ + H2
     2713!===========================================================
     2714
     2715   nb_reaction_4 = nb_reaction_4 + 1
     2716   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h2)
     2717
     2718!===========================================================
     2719!      i046 : H2O+ + NO -> NO+ + H2O
     2720!===========================================================
     2721
     2722   nb_reaction_4 = nb_reaction_4 + 1
     2723   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_h2o)
     2724
     2725!===========================================================
     2726!      i047 : H2O+ + e- -> H + H + O
     2727!===========================================================
     2728
     2729   nb_reaction_4 = nb_reaction_4 + 1
     2730   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 2.0, i_h, 1.0, i_o)
     2731
     2732!===========================================================
     2733!      i048 : H2O+ + e- -> H + OH
     2734!===========================================================
     2735
     2736   nb_reaction_4 = nb_reaction_4 + 1
     2737   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h, 1.0, i_oh)
     2738
     2739!===========================================================
     2740!      i049 : H2O+ + e- -> H2 + O
     2741!===========================================================
     2742
     2743   nb_reaction_4 = nb_reaction_4 + 1
     2744   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o)
    25472745
    25482746end if    !ionchem
    … …  
    26102808    (nb_reaction_4 /= nb_reaction_4_max)) then
    26112809   print*, 'wrong dimensions in indice'
    2612    stop
     2810   call abort_physic("indice","wrong array dimensions",1)
    26132811end if 
    26142812
    … …  
    26232821                           i_co2plus, i_oplus, i_o2plus, i_noplus,   &
    26242822                           i_coplus, i_cplus, i_n2plus, i_nplus,     &
    2625                            i_hplus, i_hco2plus, i_hcoplus, i_elec,   &
    2626                            dens, rm, c)
     2823                           i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,&
     2824                           i_elec, dens, rm, c)
    26272825       
    26282826!*****************************************************************
    … …  
    26352833     &                      igcm_noplus, igcm_coplus, igcm_cplus,        &
    26362834     &                      igcm_n2plus, igcm_nplus, igcm_hplus,         &
    2637      &                      igcm_hco2plus, igcm_hcoplus, igcm_elec
     2835     &                      igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,   &
     2836     &                      igcm_elec
    26382837
    26392838      implicit none
    26402839
    2641 #include "callkeys.h"
     2840      include "callkeys.h"
    26422841
    26432842!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    26562855                 i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus,   &
    26572856                 i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus,    &
    2658                  i_hcoplus, i_elec
     2857                 i_hcoplus, i_h2oplus, i_elec
    26592858
    26602859      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
    … …  
    26832882         if (igcm_co2 == 0) then
    26842883            write(*,*) "gcmtochim: Error; no CO2 tracer !!!"
    2685             stop
     2884            call abort_physic("gcmtochim","missing co2 tracer",1)
    26862885         endif
    26872886         if (igcm_co == 0) then
    26882887            write(*,*) "gcmtochim: Error; no CO tracer !!!"
    2689             stop
     2888            call abort_physic("gcmtochim","missing co tracer",1)
    26902889         end if
    26912890         if (igcm_o == 0) then
    26922891            write(*,*) "gcmtochim: Error; no O tracer !!!"
    2693             stop
     2892            call abort_physic("gcmtochim","missing o tracer",1)
    26942893         end if
    26952894         if (igcm_o1d == 0) then
    26962895            write(*,*) "gcmtochim: Error; no O1D tracer !!!"
    2697             stop
     2896            call abort_physic("gcmtochim","missing o1d tracer",1)
    26982897         end if
    26992898         if (igcm_o2 == 0) then
    27002899            write(*,*) "gcmtochim: Error; no O2 tracer !!!"
    2701             stop
     2900            call abort_physic("gcmtochim","missing o2 tracer",1)
    27022901         end if
    27032902         if (igcm_o3 == 0) then
    27042903            write(*,*) "gcmtochim: Error; no O3 tracer !!!"
    2705             stop
     2904            call abort_physic("gcmtochim","missing o3 tracer",1)
    27062905         end if
    27072906         if (igcm_h == 0) then
    27082907            write(*,*) "gcmtochim: Error; no H tracer !!!"
    2709             stop
     2908            call abort_physic("gcmtochim","missing h tracer",1)
    27102909         end if
    27112910         if (igcm_h2 == 0) then
    27122911            write(*,*) "gcmtochim: Error; no H2 tracer !!!"
    2713             stop
     2912            call abort_physic("gcmtochim","missing h2 tracer",1)
    27142913         end if
    27152914         if (igcm_oh == 0) then
    27162915            write(*,*) "gcmtochim: Error; no OH tracer !!!"
    2717             stop
     2916            call abort_physic("gcmtochim","missing oh tracer",1)
    27182917         end if
    27192918         if (igcm_ho2 == 0) then
    27202919            write(*,*) "gcmtochim: Error; no HO2 tracer !!!"
    2721             stop
     2920            call abort_physic("gcmtochim","missing ho2 tracer",1)
    27222921         end if
    27232922         if (igcm_h2o2 == 0) then
    27242923            write(*,*) "gcmtochim: Error; no H2O2 tracer !!!"
    2725             stop
     2924            call abort_physic("gcmtochim","missing h2o2 tracer",1)
    27262925         end if
    27272926         if (igcm_n == 0) then
    27282927            write(*,*) "gcmtochim: Error; no N tracer !!!"
    2729             stop
     2928            call abort_physic("gcmtochim","missing n tracer",1)
    27302929         end if
    27312930         if (igcm_n2d == 0) then
    27322931            write(*,*) "gcmtochim: Error; no N2D tracer !!!"
    2733             stop
     2932            call abort_physic("gcmtochim","missing n2d tracer",1)
    27342933         end if
    27352934         if (igcm_no == 0) then
    27362935            write(*,*) "gcmtochim: Error; no NO tracer !!!"
    2737             stop
     2936            call abort_physic("gcmtochim","missing no tracer",1)
    27382937         end if
    27392938         if (igcm_no2 == 0) then
    27402939            write(*,*) "gcmtochim: Error; no NO2 tracer !!!"
    2741             stop
     2940            call abort_physic("gcmtochim","missing no2 tracer",1)
    27422941         end if
    27432942         if (igcm_n2 == 0) then
    27442943            write(*,*) "gcmtochim: Error; no N2 tracer !!!"
    2745             stop
     2944            call abort_physic("gcmtochim","missing n2 tracer",1)
    27462945         end if
    27472946         if (igcm_h2o_vap == 0) then
    27482947            write(*,*) "gcmtochim: Error; no water vapor tracer !!!"
    2749             stop
     2948            call abort_physic("gcmtochim","missing h2o_vap tracer",1)
    27502949         end if
    27512950         if (ionchem) then
    27522951            if (igcm_co2plus == 0) then
    27532952               write(*,*) "gcmtochim: Error; no CO2+ tracer !!!"
    2754                stop
     2953               call abort_physic("gcmtochim","missing co2plus tracer",1)
    27552954            end if
    27562955            if (igcm_oplus == 0) then
    27572956               write(*,*) "gcmtochim: Error; no O+ tracer !!!"
    2758                stop
     2957               call abort_physic("gcmtochim","missing oplus tracer",1)
    27592958            end if
    27602959            if (igcm_o2plus == 0) then
    27612960               write(*,*) "gcmtochim: Error; no O2+ tracer !!!"
    2762                stop
     2961               call abort_physic("gcmtochim","missing o2plus tracer",1)
    27632962            end if
    27642963            if (igcm_noplus == 0) then
    27652964               write(*,*) "gcmtochim: Error; no NO+ tracer !!!"
    2766                stop
     2965               call abort_physic("gcmtochim","missing noplus tracer",1)
    27672966            endif
    27682967            if (igcm_coplus == 0) then
    27692968               write(*,*) "gcmtochim: Error; no CO+ tracer !!!"
    2770                stop
     2969               call abort_physic("gcmtochim","missing coplus tracer",1)
    27712970            endif
    27722971            if (igcm_cplus == 0) then
    27732972               write(*,*) "gcmtochim: Error; no C+ tracer !!!"
    2774                stop
     2973               call abort_physic("gcmtochim","missing cplus tracer",1)
    27752974            endif
    27762975            if (igcm_n2plus == 0) then
    27772976               write(*,*) "gcmtochim: Error; no N2+ tracer !!!"
    2778                stop
     2977               call abort_physic("gcmtochim","missing n2plus tracer",1)
    27792978            endif
    27802979            if (igcm_nplus == 0) then
    27812980               write(*,*) "gcmtochim: Error; no N+ tracer !!!"
    2782                stop
     2981               call abort_physic("gcmtochim","missing nplus tracer",1)
    27832982            endif
    27842983            if (igcm_hplus == 0) then
    27852984               write(*,*) "gcmtochim: Error; no H+ tracer !!!"
    2786                stop
     2985               call abort_physic("gcmtochim","missing hplus tracer",1)
    27872986            endif
    27882987            if (igcm_hco2plus == 0) then
    27892988               write(*,*) "gcmtochim: Error; no HCO2+ tracer !!!"
    2790                stop
     2989               call abort_physic("gcmtochim","missing hco2plus tracer",1)
    27912990            endif
    27922991            if (igcm_hcoplus == 0) then
    27932992               write(*,*) "gcmtochim: Error; no HCO+ tracer !!!"
    2794                stop
     2993               call abort_physic("gcmtochim","missing hcoplus tracer",1)
     2994            endif
     2995            if (igcm_h2oplus == 0) then
     2996               write(*,*) "gcmtochim: Error; no H2O+ tracer !!!"
     2997               call abort_physic("gcmtochim","missing h2oplus tracer",1)
    27952998            endif
    27962999            if (igcm_elec == 0) then
    27973000               write(*,*) "gcmtochim: Error; no e- tracer !!!"
    2798                stop
     3001               call abort_physic("gcmtochim","missing elec tracer",1)
    27993002            end if
    28003003         end if  ! ionchem
    … …  
    28393042            rm(l,i_hco2plus) = zycol(l, igcm_hco2plus)
    28403043            rm(l,i_hcoplus)  = zycol(l, igcm_hcoplus)
     3044            rm(l,i_h2oplus)  = zycol(l, igcm_h2oplus)
    28413045            rm(l,i_elec)     = zycol(l, igcm_elec)
    28423046         end do
    … …  
    28673071                           i_co2plus, i_oplus, i_o2plus, i_noplus,    &
    28683072                           i_coplus, i_cplus, i_n2plus, i_nplus,      &
    2869                            i_hplus, i_hco2plus, i_hcoplus, i_elec,    &
    2870                            dens, c)
     3073                           i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, &
     3074                           i_elec, dens, c)
    28713075 
    28723076!*****************************************************************
    … …  
    28793083                            igcm_noplus, igcm_coplus, igcm_cplus,         &
    28803084                            igcm_n2plus, igcm_nplus, igcm_hplus,          &
    2881                             igcm_hco2plus, igcm_hcoplus, igcm_elec
     3085                            igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
     3086                            igcm_elec
    28823087
    28833088      implicit none
    … …  
    28993104                 i_co2plus, i_oplus, i_o2plus, i_noplus,         &
    29003105                 i_coplus, i_cplus, i_n2plus, i_nplus,           &
    2901                  i_hplus, i_hco2plus, i_hcoplus, i_elec
     3106                 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, i_elec
    29023107
    29033108      real :: dens(nlayer)     ! total number density (molecule.cm-3)
    … …  
    29533158            zycol(l, igcm_hco2plus)= c(l,i_hco2plus)/dens(l)
    29543159            zycol(l, igcm_hcoplus) = c(l,i_hcoplus)/dens(l)
     3160            zycol(l, igcm_h2oplus) = c(l,i_h2oplus)/dens(l)
    29553161            zycol(l, igcm_elec)    = c(l,i_elec)/dens(l)
    29563162         end do
Note: See TracChangeset for help on using the changeset viewer.