Changeset 2158


Ignore:
Timestamp:
Sep 11, 2019, 9:13:17 AM (5 years ago)
Author:
emillour
Message:

Mars GCM:

  • Updated chemical core to include ionospheric chemistry

FGG

Location:
trunk/LMDZ.MARS
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2157 r2158  
    27312731== 10/09/2019 == JYC
    27322732- Use Wilke's formulation for molecular diffusion
     2733
     2734== 11/09/2019 == FGG
     2735- Updated chemical core to include ionospheric chemistry

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

    r2044 r2158  
    1       subroutine calchim(ngrid,nlayer,nq,                                   &
     1      subroutine calchim(ngrid,nlayer,nq,                           &
    22                         ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
    33                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
    … …  
    1717      use comcstfi_h
    1818      use photolysis_mod
     19      use iono_h, only: temp_elect
    1920
    2021      implicit none
    … …  
    161162      real :: zdens(nlayer)        ! Density  (cm-3)
    162163      real :: ztemp(nlayer)        ! Temperature (K)
     164      real :: ztelec(nlayer)       ! Electronic temperature (K)
    163165      real :: zlocal(nlayer)       ! Altitude (km)
    164166      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
    … …  
    186188                                    !  within one physical timestep
    187189
     190      logical :: unichim            ! only one, unified chemical scheme at all
     191                                    ! layers (default), or upper atmospheric
     192                                    ! scheme in the thermosphere?
     193
    188194!=======================================================================
    189195!     initialization of the chemistry (first call only)
    190196!=======================================================================
     197
     198      !call only unified chemistry (default), or also upper atmospheric model
     199      !unichim = .true.   unified chemistry
     200      !unichim = .false.  2 different models
     201      unichim = .true.
    191202
    192203      if (firstcall) then
    … …  
    201212            end if
    202213         end if
     214
     215         if(.not.unichim) then
     216            !Read reaction rates from external file if the upper atmospheric
     217            !chemistry is called
     218            call chemthermos_readini
     219         endif
     220
    203221         ! find index of chemical tracers to use
    204222         allocate(niq(nq))
    205223         ! Listed here are all tracers that can go into photochemistry
    206224         nbq = 0 ! to count number of tracers
    207          ! Species ALWAYS present if photochem=.T. or thermochem=.T.
     225         ! Species ALWAYS present if photochem=.T.
    208226         i_co2 = igcm_co2
    209227         if (i_co2 == 0) then
    … …  
    358376            niq(nbq) = i_h2o
    359377         end if
    360          !Check tracers needed for thermospheric chemistry
    361          if(thermochem) then
    362             chemthermod=0  !Default: C/O/H chemistry
    363             !Nitrogen chemistry
    364             !NO is used to determine if N chemistry is wanted
    365             !chemthermod=2 -> N chemistry
    366             if (i_no == 0) then
    367                write(*,*) "calchim: no NO tracer"
    368                write(*,*) "C/O/H themosp chemistry only "
    369             else
    370                chemthermod=2
    371                write(*,*) "calchim: NO in traceur.def"
    372                write(*,*) "Nitrogen chemistry included"
    373             end if
    374             ! N
    375             if(chemthermod == 2) then
    376                if (i_n == 0) then
    377                   write(*,*) "calchim: Error; no N tracer !!!"
    378                   write(*,*) "N is needed if NO is in traceur.def"
    379                   stop
    380                end if
    381             ! NO2
    382                if (i_no2 == 0) then
    383                   write(*,*) "calchim: Error; no NO2 tracer !!!"
    384                   write(*,*) "NO2 is needed if NO is in traceur.def"
    385                   stop
    386                end if
    387             ! N(2D)
    388                if (i_n2d == 0) then
    389                   write(*,*) "calchim: Error; no N2D !!!"
    390                   write(*,*) "N2D is needed if NO is in traceur.def"
    391                   stop
    392                end if
    393             endif    !Of if(chemthermod == 2)
    394             ! Ions
    395             ! O2+ is used to determine if ion chemistry is needed
    396             ! chemthermod=3 -> ion chemistry
    397             i_o2plus = igcm_o2plus
    398             if(chemthermod == 2) then
    399                if (i_o2plus == 0) then
    400                   write(*,*) "calchim: no O2+ tracer; no ion chemistry"
    401                else
    402                   nbq = nbq + 1
    403                   niq(nbq) = i_o2plus
    404                   chemthermod = 3
    405                   write(*,*) "calchim: O2+ in traceur.def"
    406                   write(*,*) "Ion chemistry included"
    407                end if
    408             else
    409                if (i_o2plus /= 0) then
    410                   write(*,*) "calchim: O2+ is present, but NO is not!!!"
    411                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    412                   stop
    413                endif
    414             endif   !Of if(chemthermod == 2)
    415             ! CO2+
    416             i_co2plus = igcm_co2plus
    417             if(chemthermod == 3) then
    418                if (i_co2plus == 0) then
    419                   write(*,*) "calchim: Error; no CO2+ tracer !!!"
    420                   write(*,*) "CO2+ is needed if O2+ is in traceur.def"
    421                   stop
    422                else
    423                   nbq = nbq + 1
    424                   niq(nbq) = i_co2plus
    425                end if
    426             else
    427                if (i_co2plus /= 0) then
    428                   write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
    429                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    430                   stop
    431                endif
    432             endif    !Of if(chemthermod == 3)
    433             ! O+
    434             i_oplus = igcm_oplus
    435             if(chemthermod == 3) then
    436                if (i_oplus == 0) then
    437                   write(*,*) "calchim: Error; no O+ tracer !!!"
    438                   write(*,*) "O+ is needed if O2+ is in traceur.def"
    439                   stop
    440                else
    441                   nbq = nbq + 1
    442                   niq(nbq) = i_oplus
    443                end if
    444             else
    445                if (i_oplus /= 0) then
    446                   write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
    447                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    448                   stop
    449                endif
    450             endif   !Of if (chemthermod == 3)
    451             ! CO+
    452             i_coplus = igcm_coplus
    453             if(chemthermod == 3) then
    454                if (i_coplus == 0) then
    455                   write(*,*) "calchim: Error; no CO+ tracer !!!"
    456                   write(*,*) "CO+ is needed if O2+ is in traceur.def"
    457                   stop
    458                else
    459                   nbq = nbq + 1
    460                   niq(nbq) = i_coplus
    461                end if
    462             else
    463                if (i_coplus /= 0) then
    464                   write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
    465                   write(*,*) " Both must be in traceur.def if ionosphere wanted"
    466                   stop
    467                endif
    468             endif   ! Of if (chemthermod == 3)
    469             ! C+
    470             i_cplus = igcm_cplus
    471             if(chemthermod == 3) then
    472                if (i_cplus == 0) then
    473                   write(*,*) "calchim: Error; no C+ tracer !!!"
    474                   write(*,*) "C+ is needed if O2+ is in traceur.def"
    475                   stop
    476                else
    477                   nbq = nbq + 1
    478                   niq(nbq) = i_cplus
    479                end if
    480             else
    481                if (i_cplus /= 0) then
    482                   write(*,*) "calchim: Error; C+ is present, but O2+ is not!!!"
    483                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    484                   stop
    485                endif
    486             endif   ! Of if (chemthermod == 3)
    487             ! N+
    488             i_nplus = igcm_nplus
    489             if(chemthermod == 3) then
    490                if (i_nplus == 0) then
    491                   write(*,*) "calchim: Error; no N+ tracer !!!"
    492                   write(*,*) "N+ is needed if O2+ is in traceur.def"
    493                   stop
    494                else
    495                   nbq = nbq + 1
    496                   niq(nbq) = i_nplus
    497                end if
    498             else
    499                if (i_nplus /= 0) then
    500                   write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
    501                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    502                   stop
    503                endif
    504             endif   !Of if (chemthermod == 3)
    505             ! NO+
    506             i_noplus = igcm_noplus
    507             if(chemthermod == 3) then
    508                if (i_noplus == 0) then
    509                   write(*,*) "calchim: Error; no NO+ tracer !!!"
    510                   write(*,*) "NO+ is needed if O2+ is in traceur.def"
    511                   stop
    512                else
    513                   nbq = nbq + 1
    514                   niq(nbq) = i_noplus
    515                end if
    516             else
    517                if (i_noplus /= 0) then
    518                   write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
    519                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    520                   stop
    521                endif
    522             endif   !Of if (chemthermod == 3)
    523             ! N2+
    524             i_n2plus = igcm_n2plus
    525             if (chemthermod == 3) then
    526                if (i_n2plus == 0) then
    527                   write(*,*) "calchim: Error; no N2+ tracer !!!"
    528                   write(*,*) "N2+ is needed if O2+ is in traceur.def"
    529                   stop
    530                else
    531                   nbq = nbq + 1
    532                   niq(nbq) = i_n2plus
    533                end if
    534             else
    535                if (i_n2plus /= 0) then
    536                   write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
    537                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    538                   stop
    539                endif
    540             endif   !Of if (chemthermod == 3)
    541             !H+
    542             i_hplus = igcm_hplus
    543             if (chemthermod == 3) then
    544                if (i_hplus == 0) then
    545                   write(*,*) "calchim: Error; no H+ tracer !!!"
    546                   write(*,*) "H+ is needed if O2+ is in traceur.def"
    547                   stop
    548                else
    549                   nbq = nbq + 1
    550                   niq(nbq) = i_hplus
    551                end if
    552             else
    553                if (i_hplus /= 0) then
    554                   write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
    555                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    556                   stop
    557                endif
    558             endif   !Of if (chemthermod == 3)
    559             ! HCO2+
    560             i_hco2plus = igcm_hco2plus
    561             if(chemthermod == 3) then
    562                if (i_hco2plus == 0) then
    563                   write(*,*) "calchim: Error; no HCO2+ tracer !!!"
    564                   write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
    565                   stop
    566                else
    567                   nbq = nbq + 1
    568                   niq(nbq) = i_hco2plus
    569                end if
    570             else
    571                if (i_hco2plus /= 0) then
    572                   write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
    573                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    574                   stop
    575                endif
    576             endif    !Of if(chemthermod == 3)
    577             !e-
    578             i_elec = igcm_elec
    579             if(chemthermod == 3) then
    580                if (i_elec == 0) then
    581                   write(*,*) "calchim: Error; no e- tracer !!!"
    582                   write(*,*) "e- is needed if O2+ is in traceur.def"
    583                   stop
    584                else
    585                   nbq = nbq + 1
    586                   niq(nbq) = i_elec
    587                end if
    588             else
    589                if(i_elec /= 0) then
    590                   write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
    591                   write(*,*) "Both must be in traceur.def if ionosphere wanted"
    592                   stop
    593                endif
    594             endif   !Of if (chemthermod == 3)
    595          endif      !Of thermochem
    596 
     378         i_o2plus = igcm_o2plus
     379         if (i_o2plus == 0) then
     380            write(*,*) "calchim: Error, no O2+ tracer !!!"
     381            stop
     382         else
     383            nbq = nbq + 1
     384            niq(nbq) = i_o2plus
     385         end if
     386         i_co2plus = igcm_co2plus
     387         if (i_co2plus == 0) then
     388            write(*,*) "calchim: Error, no CO2+ tracer !!!"
     389            stop
     390         else
     391            nbq = nbq + 1
     392            niq(nbq) = i_co2plus
     393         end if
     394         i_oplus=igcm_oplus
     395         if (i_oplus == 0) then
     396            write(*,*) "calchim: Error, no O+ tracer !!!"
     397            stop
     398         else
     399            nbq = nbq + 1
     400            niq(nbq) = i_oplus
     401         end if
     402         i_noplus=igcm_noplus
     403         if (i_noplus == 0) then
     404            write(*,*) "calchim: Error, no NO+ tracer !!!"
     405            stop
     406         else
     407            nbq = nbq + 1
     408            niq(nbq) = i_noplus
     409         end if
     410         i_coplus=igcm_coplus
     411         if (i_coplus == 0) then
     412            write(*,*) "calchim: Error, no CO+ tracer !!!"
     413            stop
     414         else
     415            nbq = nbq + 1
     416            niq(nbq) = i_coplus
     417         end if
     418         i_cplus=igcm_cplus
     419         if (i_cplus == 0) then
     420            write(*,*) "calchim: Error, no C+ tracer !!!"
     421            stop
     422         else
     423            nbq = nbq + 1
     424            niq(nbq) = i_cplus
     425         end if
     426         i_n2plus=igcm_n2plus
     427         if (i_n2plus == 0) then
     428            write(*,*) "calchim: Error, no N2+ tracer !!!"
     429            stop
     430         else
     431            nbq = nbq + 1
     432            niq(nbq) = i_n2plus
     433         end if
     434         i_nplus=igcm_nplus
     435         if (i_nplus == 0) then
     436            write(*,*) "calchim: Error, no N+ tracer !!!"
     437            stop
     438         else
     439            nbq = nbq + 1
     440            niq(nbq) = i_nplus
     441         end if
     442         i_hplus=igcm_hplus
     443         if (i_hplus == 0) then
     444            write(*,*) "calchim: Error, no H+ tracer !!!"
     445            stop
     446         else
     447            nbq = nbq + 1
     448            niq(nbq) = i_hplus
     449         end if
     450         i_hco2plus=igcm_hco2plus
     451         if (i_hco2plus == 0) then
     452            write(*,*) "calchim: Error, no HCO2+ tracer !!!"
     453            stop
     454         else
     455            nbq = nbq + 1
     456            niq(nbq) = i_hco2plus
     457         end if
     458         i_elec = igcm_elec
     459         if (i_elec == 0) then
     460            write(*,*) "calchim: Error, no e- tracer !!!"
     461            stop
     462         else
     463            nbq = nbq + 1
     464            niq(nbq) = i_elec
     465         end if
     466         
    597467         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
    598468               
    … …  
    634504            zlocal(l) = zzlay(ig,l)/1000.
    635505            zmmean(l) = mmean(ig,l)
     506            ztelec(l) = temp_elect(zlocal(l),ztemp(l),1)
     507            !Electronic temperature. Index 1 -> Viking; Index 2-> MAVEN
    636508
    637509!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
    … …  
    640512            surfice1d(l)  = surfice(ig,l)*1.e-2
    641513
    642 !           search for switch index between regions
    643 
    644             if (photochem .and. thermochem) then
     514!           search for switch index between regions 
     515            if(unichim) then
     516               lswitch=nlayer+1
     517            else
    645518               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
    646519                  lswitch = l
    647520                  foundswitch = 1
    648521               end if
    649             end if
    650             if (.not. photochem) then
    651                lswitch = 22
    652             end if
    653             if (.not. thermochem) then
    654                lswitch = min(50,nlayer+1)
    655             end if
    656 
     522            endif
    657523         end do ! of do l=1,nlayer
    658524
    … …  
    664530!=======================================================================
    665531
    666 !        chemistry in lower atmosphere
    667 
     532!        chemistry: only one scheme at all layers
     533         
    668534         if (photochem) then
    669535            call photochemistry(nlayer,nq,                            &
    670536                           ig,lswitch,zycol,szacol,ptimestep,         &
    671                            zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
    672                            zday,surfdust1d,surfice1d,jo3,jh2o,taucol,iter)
     537                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,   &
     538                           dist_sol,zday,surfdust1d,surfice1d,        &
     539                           jo3,jh2o,taucol,iter)
    673540
    674541!        ozone photolysis, for output
    … …  
    679546               iter_3d(ig,l) = iter(l)
    680547            end do
    681 
    682548!        condensation of h2o2
    683549
    … …  
    685551                         ig,ptimestep,pplev,pplay,                 &
    686552                         ztemp,zycol,dqcloud,dqscloud)
    687          end if
    688 
    689 !        chemistry in upper atmosphere
    690        
    691          if (thermochem) then
    692             call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
    693                              zdens,zpress,zlocal,szacol,ptimestep,zday,&
    694                              em_no,em_o2)
    695             do l=1,nlayer
    696                emission_no(ig,l)=em_no(l)
    697                emission_o2(ig,l)=em_o2(l)
    698             enddo
     553
     554            if(.not.unichim) then
     555               chemthermod=3   !C/O/H/N/ions
     556               call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
     557                    zdens,zpress,zlocal,szacol,ptimestep,zday,&
     558                    em_no,em_o2)
     559               do l=1,nlayer
     560                  emission_no(ig,l)=em_no(l)
     561                  emission_o2(ig,l)=em_o2(l)
     562               enddo
     563            end if
    699564         end if
    700565
    … …  
    737602      end do ! of do ig=1,ngrid
    738603
    739 !      stop
    740604!=======================================================================
    741605!     write outputs
    … …  
    753617            call writediagfi(ngrid,'iter','iterations',  &
    754618                             ' ',3,iter_3d(1,1))
    755             call writediagfi(ngrid,'emission_no',        &
    756                  'NO nightglow emission rate','cm-3 s-1',3,emission_no)
    757             call writediagfi(ngrid,'emission_o2',        &
    758                  'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
     619            if(.not.unichim) then
     620               call writediagfi(ngrid,'emission_no',        &
     621                    'NO nightglow emission rate','cm-3 s-1',3,emission_no)
     622               call writediagfi(ngrid,'emission_o2',        &
     623                    'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
     624            endif
    759625           
    760626            if (callstats) then

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

    r2030 r2158  
    1010     &                      igcm_nplus, igcm_noplus, igcm_n2plus,
    1111     &                      igcm_hplus, igcm_hco2plus, mmol,
    12      &                      igcm_he
     12     &                      igcm_he, igcm_elec
    1313      use conc_mod, only: mmean, Akknew, rnew, cpnew
    1414      implicit none
    … …  
    239239            nbq = nbq + 1
    240240            niq(nbq) = igcm_hco2plus
     241            aki(nbq) = 0.0
     242            cpi(nbq) = 0.0
     243         endif
     244         if(igcm_elec /= 0) then
     245            nbq = nbq + 1
     246            niq(nbq) = igcm_elec
    241247            aki(nbq) = 0.0
    242248            cpi(nbq) = 0.0

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

    r1266 r2158  
    106106         END SUBROUTINE allocate_param_iono
    107107
     108
     109!***********************************************************************
     110      function temp_elect(zkm,tt,origin)
     111
     112!     Computes the electronic temperature, either from Viking (origin=1)
     113!     or MAVEN (origin=2)
     114
     115!***********************************************************************
     116     
     117!     Arguments         
     118
     119      real            tt        ! Temperature
     120      real            zkm       !  Altitude in km
     121      integer         origin    ! Viking (origin=1) or MAVEN (origin=2)
     122
     123! local variables:
     124      real          temp_elect     ! electronic temperatures
     125      real          zhanson(9),tehanson(9)
     126      real          incremento
     127      integer       ii, i1, i2
     128
     129      zhanson(1:9) = (/ 120.,130.,150.,175.,200.,225.,250.,275.,300. /)
     130      tehanson(2:9) = (/ 200.,300.,500.,1250.,2000.,2200.,2400.,2500. /)
     131      tehanson(1) = tt
     132
     133      if(origin.eq.1) then
     134         if ( zkm .le. 120. ) then
     135            temp_elect = tt
     136         else if(zkm .ge.300.) then
     137            temp_elect=tehanson(9)
     138         else
     139            do ii=9,2,-1
     140               if ( zkm .lt. zhanson(ii) ) then
     141                  i1 = ii - 1
     142                  i2 = ii
     143               endif
     144            enddo
     145            incremento=(tehanson(i2)-tehanson(i1))/(zhanson(i2)-zhanson(i1))
     146            temp_elect = tehanson(i1) + (zkm-zhanson(i1)) * incremento
     147        endif
     148      else if(origin.eq.2) then
     149         !MAVEN measured electronic temperature (Ergun et al., GRL 2015)
     150         !Note that the Langmuir probe is not sensitive below ~500K, so
     151         !electronic temperatures in the lower thermosphere (<150 km) may
     152         !be overestimated by this formula
     153         if(zkm.le.120) then
     154            temp_elect = tt
     155         else
     156            temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
     157         endif
     158      else
     159         write(*,*)'Error in function temp_elect:'
     160         write(*,*)'Origin must be either 1 or 2'
     161         write(*,*)'Using neutral temperature instead'
     162         temp_elect = tt
     163      endif
     164
     165      return
     166
     167      end function temp_elect
     168
    108169END MODULE iono_h

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

    r2150 r2158  
    1 !*****************************************************************
     1!****************************************************************
    22!
    33!     Photochemical routine
    … …  
    1515subroutine photochemistry(nlayer, nq,                                   &
    1616                          ig, lswitch, zycol, sza, ptimestep, press,    &
    17                           alt, temp, dens, zmmean, dist_sol, zday,      &
     17                          alt, temp, temp_elect,dens, zmmean,           &
     18                          dist_sol, zday,                               &
    1819                          surfdust1d, surfice1d, jo3, jh2o,tau, iter)
    1920
    … …  
    2324                           jonline
    2425
    25 use param_v4_h, only: jdistot, jdistot_b
     26use param_v4_h, only: jion
    2627
    2728implicit none
    … …  
    4243real :: alt(nlayer)           ! altitude (km)
    4344real :: temp(nlayer)          ! temperature (k)
     45real :: temp_elect(nlayer)    ! electronic temperature (K)
    4446real :: dens(nlayer)          ! density (cm-3)
    4547real :: zmmean(nlayer)        ! mean molar mass (g/mole)
    … …  
    6870!===================================================================
    6971
    70 integer, parameter :: nesp = 17  ! number of species in the chemical code
     72integer, parameter :: nesp = 28  ! number of species in the chemical code
    7173
    7274integer :: phychemrat            ! (physical timestep)/(nominal chemical timestep)
    … …  
    7476integer :: lswitch
    7577logical, save :: firstcall = .true.
    76 logical :: jparam                ! switch for J parameterization
     78logical :: jionos                ! switch for J parameterization
    7779
    7880! tracer indexes in the chemistry:
    … …  
    9597integer,parameter :: i_no2  = 16
    9698integer,parameter :: i_n2   = 17
     99integer,parameter :: i_co2plus = 18
     100integer,parameter :: i_oplus   = 19
     101integer,parameter :: i_o2plus  = 20
     102integer,parameter :: i_noplus  = 21
     103integer,parameter :: i_coplus  = 22
     104integer,parameter :: i_cplus   = 23
     105integer,parameter :: i_n2plus  = 24
     106integer,parameter :: i_nplus   = 25
     107integer,parameter :: i_hplus   = 26
     108integer,parameter :: i_hco2plus= 27
     109integer,parameter :: i_elec    = 28
    97110
    98111
    … …  
    134147   call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    135148               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    136                i_n, i_n2d, i_no, i_no2, i_n2)
     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               i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
    137152   firstcall = .false.
    138153end if
    … …  
    145160               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    146161               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    147                i_n, i_n2d, i_no, i_no2, i_n2, dens, rm, c)
     162               i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
     163               i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
     164               i_n2plus, i_nplus, i_hplus, i_hco2plus,         &
     165               i_elec, dens, rm, c)
    148166
    149167!===================================================================
    … …  
    151169!===================================================================
    152170
    153 jparam= .false.
     171jionos= .true.
    154172
    155173if (jonline) then
    … …  
    160178                             i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,   &
    161179                             tau, sza, dist_sol, v_phot)
     180
     181      if(jionos) then
     182         call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
     183         do ilay=1,lswitch-1
     184            call phdisrate(ig,nlayer,2,sza,ilay)
     185         enddo
     186         v_phot(:,14)=jion(1,:,1)
     187         v_phot(:,15)=jion(1,:,2)
     188         v_phot(:,16)=jion(1,:,2)
     189         v_phot(:,17)=jion(1,:,3)
     190         v_phot(:,18)=jion(1,:,3)
     191         v_phot(:,19)=jion(1,:,4)
     192         v_phot(:,20)=jion(1,:,4)
     193         v_phot(:,21)=jion(2,:,1)
     194         v_phot(:,22)=jion(3,:,1)
     195         v_phot(:,23)=jion(10,:,1)
     196         v_phot(:,24)=jion(11,:,1)
     197         v_phot(:,25)=jion(11,:,2)
     198         v_phot(:,26)=jion(11,:,2)
     199         v_phot(:,27)=jion(8,:,1)
     200         v_phot(:,28)=jion(8,:,2)
     201         v_phot(:,29)=jion(8,:,2)
     202         v_phot(:,30)=jion(9,:,1)
     203         v_phot(:,31)=jion(12,:,1)
     204      endif
     205!      write(*,*)'photochemistry/205',c(:,i_co2),ig
     206!      write(*,*)'photochemistry/206',v_phot(:,3),ig
     207     
    162208   else ! night
    163209      v_phot(:,:) = 0.
    164210   end if
    165 else if(jparam) then
    166    call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
    167    
    168    do ilay=1,lswitch-1
    169       call phdisrate(ig,nlayer,2,sza,ilay)
    170    enddo
    171    v_phot(:,1)=jdistot(2,:)
    172    v_phot(:,2)=jdistot_b(2,:)
    173    v_phot(:,3)=jdistot(1,:)
    174    v_phot(:,4)=jdistot_b(1,:)
    175    v_phot(:,5)=jdistot(7,:)
    176    v_phot(:,6)=jdistot_b(7,:)
    177    v_phot(:,7)=jdistot(4,:)
    178    v_phot(:,8)=jdistot(6,:)
    179    v_phot(:,10)=jdistot(5,:)
    180    v_phot(:,11)=jdistot(10,:)
    181    v_phot(:,12)=jdistot(13,:)
    182    v_phot(:,13)=jdistot(8,:)
     211!else if(jparam) then
     212!   call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
     213!   do ilay=1,lswitch-1
     214!      call phdisrate(ig,nlayer,2,sza,ilay)
     215!   enddo
     216!   v_phot(:,1)=jdistot(2,:)
     217!   v_phot(:,2)=jdistot_b(2,:)
     218!   v_phot(:,3)=jdistot(1,:)
     219!   v_phot(:,4)=jdistot_b(1,:)
     220!   v_phot(:,5)=jdistot(7,:)
     221!   v_phot(:,6)=jdistot_b(7,:)
     222!   v_phot(:,7)=jdistot(4,:)
     223!   v_phot(:,8)=jdistot(6,:)
     224!   v_phot(:,10)=jdistot(5,:)
     225!   v_phot(:,11)=jdistot(10,:)
     226!   v_phot(:,12)=jdistot(13,:)
     227!   v_phot(:,13)=jdistot(8,:)
     228!   v_phot(:,14)=jion(1,:,1)
     229!   v_phot(:,15)=jion(1,:,2)
     230!   v_phot(:,16)=jion(1,:,2)
     231!   v_phot(:,17)=jion(1,:,3)
     232!   v_phot(:,18)=jion(1,:,3)
     233!   v_phot(:,19)=jion(1,:,4)
     234!   v_phot(:,20)=jion(1,:,4)
     235!   v_phot(:,21)=jion(2,:,1)
     236!   v_phot(:,22)=jion(3,:,1)
     237!   v_phot(:,23)=jion(10,:,1)
     238!   v_phot(:,24)=jion(11,:,1)
     239!   v_phot(:,25)=jion(11,:,2)
     240!   v_phot(:,26)=jion(11,:,2)
     241!   v_phot(:,27)=jion(8,:,1)
     242!   v_phot(:,28)=jion(8,:,2)
     243!   v_phot(:,29)=jion(8,:,2)
     244!   v_phot(:,30)=jion(9,:,1)
     245!   v_phot(:,31)=jion(12,:,1)
    183246else
    184247   tau = tau*7./press(1) ! dust in the lookup table is at 7 hpa
    … …  
    186249                   rm(:,i_co2), rm(:,i_o3), v_phot)
    187250end if
    188 
    189251! save o3 and h2o photolysis for output
    190252
    … …  
    205267hetero_ice  = .true.
    206268
    207 call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2),              &
    208                    c(:,i_o), c(:,i_n2), press, temp, hetero_dust, hetero_ice, &
     269call 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,                           &
    209272                   surfdust1d, surfice1d, v_phot, v_3, v_4)
    210273
    … …  
    265328
    266329!  solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
    267      
    268330   cnew(:) = c(ilev,:)
    269331
    … …  
    287349   cold(:)   = c(ilev,:)
    288350   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
    289362   cnew(:)   = 0.
    290363
    … …  
    305378               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
    306379               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
    307                i_n, i_n2d, i_no, i_no2, i_n2, dens, c)
    308 
     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,        &
     383               i_elec, dens, c)
    309384contains
    310385
    … …  
    431506
    432507 subroutine reactionrates(nlayer,                                  &
    433                           lswitch, dens, co2, o2, o, n2, press, t, &
     508                          lswitch, dens, co2, o2, o, n2, press, &
     509                          t, t_elect,                              &
    434510                          hetero_dust, hetero_ice,                 &
    435511                          surfdust1d, surfice1d,                   &
    … …  
    448524
    449525use comcstfi_h
    450 use photolysis_mod, only : nphot, nb_phot_max, &
     526use photolysis_mod, only : nphot, nphotion, nb_phot_max, &
    451527                           nb_reaction_3_max,  &
    452528                           nb_reaction_4_max
    … …  
    464540real, dimension(nlayer) :: press             ! pressure (hPa)
    465541real, dimension(nlayer) :: t                 ! temperature (K)
     542real, dimension(nlayer) :: t_elect           ! electronic temperature (K)
    466543real, dimension(nlayer) :: surfdust1d        ! dust surface area (cm2.cm-3)
    467544real, dimension(nlayer) :: surfice1d         ! ice surface area (cm2.cm-3)
    … …  
    498575                           d008, d009, d010, d011, d012,               &
    499576                           e001, e002,                                 &
     577                           i001, i002, i003, i004, i005, i006,         &
     578                           i007, i008, i009, i010, i011, i012,         &
     579                           i013, i014, i015, i016, i017, i018, i019,   &
     580                           i020, i021, i022, i023, i024, i025, i026,   &
     581                           i027, i028, i029,                           &
    500582                           h001, h002, h003, h004, h005
    501583
    … …  
    504586!----------------------------------------------------------------------
    505587
    506       nb_phot       = nphot ! initialised to the number of photolysis rates
     588      nb_phot       = nphot + nphotion ! initialised to the number of photolysis and photoionization rates
    507589      nb_reaction_3 = 0
    508590      nb_reaction_4 = 0
    … …  
    9861068      nb_reaction_4 = nb_reaction_4 + 1
    9871069      v_4(:,nb_reaction_4) = e002(:)
     1070
     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
     1085     
     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
     1146     
     1147!     Fox & Sung 2001
     1148
     1149      i009(:) = 1.0e-15
     1150     
     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(:)
    9881333
    9891334!----------------------------------------------------------------------
    … …  
    11741519  mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
    11751520
     1521
    11761522  loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
    11771523  loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
    … …  
    11871533 subroutine indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    11881534                   i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    1189                    i_n, i_n2d, i_no, i_no2, i_n2)
     1535                   i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus,       &
     1536                   i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, &
     1537                   i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec)
    11901538
    11911539!================================================================
    … …  
    12121560integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    12131561           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    1214            i_n, i_n2d, i_no, i_no2, i_n2
     1562           i_n, i_n2d, i_no, i_no2, i_n2,                   &
     1563           i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, &
     1564           i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
    12151565
    12161566! local
    … …  
    13341684
    13351685!===========================================================
     1686!      CO2 + hv -> CO2+ + e-
     1687!===========================================================
     1688
     1689nb_phot = nb_phot + 1
     1690
     1691indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
     1692
     1693!===========================================================
     1694!      CO2 + hv -> O+ + CO + e-
     1695!===========================================================
     1696!We divide this reaction in two
     1697
     1698!0.5 CO2 + hv -> CO
     1699nb_phot = nb_phot + 1
     1700
     1701indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
     1702
     1703!0.5 CO2 + hv -> O+ + e-
     1704nb_phot = nb_phot + 1
     1705
     1706indice_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
     1715nb_phot = nb_phot + 1
     1716
     1717indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
     1718
     1719!0.5 CO2 + hv -> CO+ + e-
     1720nb_phot = nb_phot + 1
     1721
     1722indice_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
     1731nb_phot = nb_phot + 1
     1732
     1733indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
     1734
     1735!0.5 CO2 + hv -> C+ + e-
     1736nb_phot = nb_phot + 1
     1737
     1738indice_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
     1744nb_phot = nb_phot + 1
     1745
     1746indice_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
     1752nb_phot = nb_phot + 1
     1753
     1754indice_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
     1760nb_phot = nb_phot + 1
     1761
     1762indice_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
     1768nb_phot = nb_phot + 1
     1769
     1770indice_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
     1778nb_phot = nb_phot + 1
     1779
     1780indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
     1781
     1782!0.5 CO + hv -> C+ + e-
     1783nb_phot = nb_phot + 1
     1784
     1785indice_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
     1792nb_phot = nb_phot + 1
     1793
     1794indice_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
     1802nb_phot = nb_phot + 1
     1803
     1804indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
     1805
     1806!0.5 N2 + hv -> N+ + e-
     1807nb_phot = nb_phot + 1
     1808
     1809indice_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
     1815nb_phot = nb_phot + 1
     1816
     1817indice_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
     1823nb_phot = nb_phot + 1
     1824
     1825indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
     1826
     1827!===========================================================
    13361828!      a001 : O + O2 + CO2 -> O3 + CO2
    13371829!===========================================================
    … …  
    16602152
    16612153indice_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
     2159nb_reaction_4 = nb_reaction_4 + 1
     2160
     2161indice_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
     2167nb_reaction_4 = nb_reaction_4 + 1
     2168
     2169indice_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
     2175nb_reaction_4 = nb_reaction_4 + 1
     2176
     2177indice_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
     2183nb_reaction_4 = nb_reaction_4 + 1
     2184
     2185indice_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
     2191nb_reaction_4 = nb_reaction_4 + 1
     2192
     2193indice_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
     2199nb_reaction_4 = nb_reaction_4 + 1
     2200
     2201indice_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
     2207nb_reaction_4 = nb_reaction_4 + 1
     2208
     2209indice_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
     2215nb_reaction_4 = nb_reaction_4 + 1
     2216
     2217indice_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
     2223nb_reaction_4 = nb_reaction_4 + 1
     2224
     2225indice_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
     2231nb_reaction_4 = nb_reaction_4 + 1
     2232
     2233indice_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
     2239nb_reaction_4 = nb_reaction_4 + 1
     2240
     2241indice_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
     2247nb_reaction_4 = nb_reaction_4 + 1
     2248
     2249indice_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
     2255nb_reaction_4 = nb_reaction_4 + 1
     2256
     2257indice_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
     2263nb_reaction_4 = nb_reaction_4 + 1
     2264
     2265indice_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
     2271nb_reaction_4 = nb_reaction_4 + 1
     2272
     2273indice_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
     2279nb_reaction_4 = nb_reaction_4 + 1
     2280
     2281indice_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
     2287nb_reaction_4 = nb_reaction_4 + 1
     2288
     2289indice_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
     2295nb_reaction_4 = nb_reaction_4 + 1
     2296
     2297indice_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
     2303nb_reaction_4 = nb_reaction_4 + 1
     2304
     2305indice_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
     2311nb_reaction_4 = nb_reaction_4 + 1
     2312
     2313indice_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
     2319nb_reaction_4 = nb_reaction_4 + 1
     2320
     2321indice_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
     2327nb_reaction_4 = nb_reaction_4 + 1
     2328
     2329indice_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
     2335nb_reaction_4 = nb_reaction_4 + 1
     2336
     2337indice_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
     2343nb_reaction_4 = nb_reaction_4 + 1
     2344
     2345indice_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
     2351nb_reaction_4 = nb_reaction_4 + 1
     2352
     2353indice_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
     2359nb_reaction_4 = nb_reaction_4 + 1
     2360
     2361indice_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
     2368!0.5HCO2+ + 0.5e -> H
     2369
     2370nb_reaction_4 = nb_reaction_4 + 1
     2371
     2372indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
     2373
     2374!0.5 HCO2+ + 0.5 e -> O + CO
     2375
     2376nb_reaction_4 = nb_reaction_4 + 1
     2377
     2378indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
     2379
     2380!===========================================================
     2381!      i029 : HCO2+ + e -> OH + CO
     2382!===========================================================
     2383
     2384nb_reaction_4 = nb_reaction_4 + 1
     2385
     2386indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
    16622387
    16632388!===========================================================
    … …  
    17342459                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
    17352460                           i_n, i_n2d, i_no, i_no2, i_n2,            &
    1736                            dens, rm, c)
     2461                           i_co2plus, i_oplus, i_o2plus, i_noplus,   &
     2462                           i_coplus, i_cplus, i_n2plus, i_nplus,     &
     2463                           i_hplus, i_hco2plus, i_elec, dens, rm, c)
     2464       
    17372465
    17382466!*****************************************************************
    … …  
    17412469     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,  &
    17422470     &                      igcm_ho2, igcm_h2o2, igcm_h2o_vap,           &
    1743      &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
     2471     &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2,&
     2472     &                      igcm_co2plus, igcm_oplus, igcm_o2plus,       &
     2473     &                      igcm_noplus, igcm_coplus, igcm_cplus,        &
     2474     &                      igcm_n2plus, igcm_nplus, igcm_hplus,         &
     2475     &                      igcm_hco2plus, igcm_elec
    17442476
    17452477      implicit none
    … …  
    17582490      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,           &
    17592491                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                   &
    1760                  i_n, i_n2d, i_no, i_no2, i_n2
     2492                 i_n, i_n2d, i_no, i_no2, i_n2,                      &
     2493                 i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus,   &
     2494                 i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, i_elec
    17612495
    17622496      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
    … …  
    18522586            stop
    18532587         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
    18542632         firstcall = .false.
    18552633      end if ! of if (firstcall)
    … …  
    18772655         rm(l,i_no2)  = zycol(l, igcm_no2)
    18782656         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)
    18792668      end do
    18802669
    … …  
    18882677
    18892678      do iesp = 1,nesp
    1890          do l = 1,nlayer!-1!lswitch-1
     2679         do l = 1,nlayer
    18912680            c(l,iesp) = rm(l,iesp)*dens(l)
    18922681         end do
    … …  
    19002689                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
    19012690                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
    1902                            i_n, i_n2d, i_no, i_no2, i_n2, dens, c)
     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,     &
     2694                           i_hplus, i_hco2plus, i_elec, dens, c)
     2695       
    19032696 
    19042697!*****************************************************************
    … …  
    19072700                            igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,   &
    19082701                            igcm_ho2, igcm_h2o2, igcm_h2o_vap,            &
    1909                             igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
     2702                            igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2, &
     2703                            igcm_co2plus, igcm_oplus, igcm_o2plus,        &
     2704                            igcm_noplus, igcm_coplus, igcm_cplus,         &
     2705                            igcm_n2plus, igcm_nplus, igcm_hplus,          &
     2706                            igcm_hco2plus, igcm_elec
    19102707
    19112708      implicit none
    … …  
    19232720      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
    19242721                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
    1925                  i_n, i_n2d, i_no, i_no2, i_n2
     2722                 i_n, i_n2d, i_no, i_no2, i_n2,                  &
     2723                 i_co2plus, i_oplus, i_o2plus, i_noplus,         &
     2724                 i_coplus, i_cplus, i_n2plus, i_nplus,           &
     2725                 i_hplus, i_hco2plus, i_elec
    19262726
    19272727      real :: dens(nlayer)     ! total number density (molecule.cm-3)
    … …  
    19622762         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
    19632763         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)
    19642775      end do
    19652776

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

    r2044 r2158  
    44
    55! photolysis
    6 
    76  logical, parameter :: jonline = .true.  ! true: on-line ! false: lookup table
    87  integer, parameter :: nphot = 13        ! number of photolysis
     8  integer, parameter :: nphotion = 18     ! number of photoionizations
    99  integer, parameter :: nabs  = 10        ! number of absorbing gases
    1010
    1111! number of reactions in chemical solver
    1212
    13   integer, parameter :: nb_phot_max       = nphot + 9 ! photolysis + quenching/heterogeneous
     13  integer, parameter :: nb_phot_max       = nphot + nphotion + 9 ! photolysis + photoionization + quenching/heterogeneous
    1414  integer, parameter :: nb_reaction_3_max = 6         ! quadratic
    15   integer, parameter :: nb_reaction_4_max = 31        ! bimolecular
     15  integer, parameter :: nb_reaction_4_max = 60        ! bimolecular
    1616
    1717! spectral grid

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r2149 r2158  
    3131     &                      nuice_ref, rho_ice, rho_dust, ref_r0,
    3232     &                      igcm_he, igcm_stormdust_mass,
    33      &                      igcm_stormdust_number 
     33     &                      igcm_stormdust_number
    3434      use comsoil_h, only: inertiedat, ! soil thermal inertia
    3535     &                     tsoil, nsoilmx ! number of subsurface layers
    … …  
    576576         if(callnlte.and.nltemodel.eq.2) call nlte_setup
    577577         if(callnirco2.and.nircorr.eq.1) call NIR_leedat
    578          if(thermochem) call chemthermos_readini
     578
    579579
    580580        IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
    … …  
    16141614c        photochemistry :
    16151615c        --------------
    1616          IF (photochem .or. thermochem) then
     1616         IF (photochem) then
    16171617
    16181618!           dust and ice surface area
    … …  
    16651665           endif
    16661666
    1667          END IF  ! of IF (photochem.or.thermochem)
     1667         END IF  ! of IF (photochem)
    16681668#endif
    16691669
    … …  
    23962396           endif ! (dustbin.ne.0)
    23972397
    2398            if (thermochem .or. photochem) then
     2398           if (photochem) then
    23992399              do iq=1,nq
    24002400                 if (noms(iq) .ne. "dust_mass" .and.
    … …  
    24282428     &                           (1000*mmol(iq))
    24292429
    2430 !                   call wstats(ngrid,"rho_"//trim(noms(iq)),
    2431 !    $                     "Number density","cm-3",3,rhopart)
    2432 !                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
    2433 !    $                     "Number density","cm-3",3,rhopart)
     2430                   call wstats(ngrid,"rho_"//trim(noms(iq)),
     2431     $                   "Number density","cm-3",3,rhopart)
     2432                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
     2433     $                  "Number density","cm-3",3,rhopart)
    24342434
    24352435!                   vertical column (molecule.cm-2)
    … …  
    24612461                 end if ! of if (noms(iq) .ne. "dust_mass" ...)
    24622462              end do ! of do iq=1,nq
    2463            end if ! of if (thermochem .or. photochem)
     2463           end if ! of if (photochem)
    24642464
    24652465           end if ! of if (tracer)
    … …  
    29392939            call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
    29402940     $           3,zdtnirco2)
     2941
     2942!            call wstats(ngrid,"q15um","15 um cooling","K/s",
     2943!     $           3,zdtnlte)
     2944!            call wstats(ngrid,"quv","UV heating","K/s",
     2945!     $           3,zdteuv)
     2946!            call wstats(ngrid,"cond","Thermal conduction","K/s",
     2947!     $           3,zdtconduc)
     2948!            call wstats(ngrid,"qnir","NIR heating","K/s",
     2949!     $           3,zdtnirco2)
    29412950
    29422951         endif  !(callthermos)
Note: See TracChangeset for help on using the changeset viewer.