Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 618)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 635)
@@ -98,4 +98,5 @@
       integer :: ig,l,i,iq,iqmax
       integer :: foundswitch, lswitch
+      integer,save :: chemthermod
 
       integer,save :: i_co2  = 0
@@ -113,4 +114,19 @@
       integer,save :: i_n2   = 0
       integer,save :: i_h2o  = 0
+      integer,save :: i_n    = 0
+      integer,save :: i_no   = 0
+      integer,save :: i_no2  = 0
+      integer,save :: i_n2d  = 0
+      integer,save :: i_co2plus=0
+      integer,save :: i_oplus=0
+      integer,save :: i_o2plus=0
+      integer,save :: i_coplus=0
+      integer,save :: i_cplus=0
+      integer,save :: i_nplus=0
+      integer,save :: i_noplus=0
+      integer,save :: i_n2plus=0
+      integer,save :: i_hplus=0
+      integer,save :: i_hco2plus=0
+      integer,save :: i_elec=0
 
       integer :: ig_vl1
@@ -157,5 +173,7 @@
 
          ! find index of chemical tracers to use
+         ! Listed here are all tracers that can go into photochemistry
          nbq = 0 ! to count number of tracers
+         ! Species ALWAYS present if photochem=.T. or thermochem=.T.
          i_co2 = igcm_co2
          if (i_co2 == 0) then
@@ -248,5 +266,5 @@
          i_ch4 = igcm_ch4
          if (i_ch4 == 0) then
-            write(*,*) "calchim: no CH4 tracer !!!"
+            write(*,*) "calchim: Error; no CH4 tracer !!!"
             write(*,*) "CH4 will be ignored in the chemistry"
          else
@@ -270,7 +288,280 @@
             niq(nbq) = i_h2o
          end if
+         !Check tracers needed for thermospheric chemistry
+         if(thermochem) then
+            chemthermod=0  !Default: C/O/H chemistry
+            !Nitrogen chemistry
+            !NO is used to determine if N chemistry is wanted
+            !chemthermod=2 -> N chemistry
+            i_no = igcm_no
+            if (i_no == 0) then
+               write(*,*) "calchim: no NO tracer"
+               write(*,*) "C/O/H themosp chemistry only "
+            else
+               nbq = nbq + 1
+               niq(nbq) = i_no
+               chemthermod=2
+               write(*,*) "calchim: NO in traceur.def"
+               write(*,*) "Nitrogen chemistry included"
+            end if
+            ! N
+            i_n = igcm_n
+            if(chemthermod == 2) then
+               if (i_n == 0) then
+                  write(*,*) "calchim: Error; no N tracer !!!"
+                  write(*,*) "N is needed if NO is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_n
+               end if
+            else 
+               if (i_n /= 0) then
+                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
+                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
+                  stop
+               endif
+            endif    !Of if(chemthermod == 2) 
+            ! NO2
+            i_no2 = igcm_no2
+            if(chemthermod == 2) then
+               if (i_no2 == 0) then
+                  write(*,*) "calchim: Error; no NO2 tracer !!!"
+                  write(*,*) "NO2 is needed if NO is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_no2
+               end if
+            else
+               if (i_no2 /= 0) then
+                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
+                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
+                  stop
+               endif
+            endif     !Of if(chemthermod == 2)
+            ! N(2D)
+            if(chemthermod == 2) then
+               i_n2d = igcm_n2d
+               if (i_n2d == 0) then
+                  write(*,*) "calchim: Error; no N2D !!!"
+                  write(*,*) "N2D is needed if NO is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_n2d
+               end if
+            else
+               if (i_n2d /= 0) then
+                  write(*,*) "calchim: Error: N2D is present, but NO is not!!!"
+                  write(*,*) "Both must be in traceur.def if N chemistry wanted"
+                  stop
+               endif
+            endif    !Of if(chemthermod == 2)
+            ! Ions
+            ! O2+ is used to determine if ion chemistry is needed
+            ! chemthermod=3 -> ion chemistry
+            i_o2plus = igcm_o2plus
+            if(chemthermod == 2) then
+               if (i_o2plus == 0) then
+                  write(*,*) "calchim: no O2+ tracer; no ion chemistry"
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_o2plus
+                  chemthermod = 3
+                  write(*,*) "calchim: O2+ in traceur.def"
+                  write(*,*) "Ion chemistry included"
+               end if
+            else
+               if (i_o2plus /= 0) then
+                  write(*,*) "calchim: O2+ is present, but NO is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if(chemthermod == 2)
+            ! CO2+
+            i_co2plus = igcm_co2plus
+            if(chemthermod == 3) then
+               if (i_co2plus == 0) then
+                  write(*,*) "calchim: Error; no CO2+ tracer !!!"
+                  write(*,*) "CO2+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_co2plus
+               end if
+            else
+               if (i_co2plus /= 0) then
+                  write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif    !Of if(chemthermod == 3)
+            ! O+
+            i_oplus = igcm_oplus
+            if(chemthermod == 3) then
+               if (i_oplus == 0) then
+                  write(*,*) "calchim: Error; no O+ tracer !!!"
+                  write(*,*) "O+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_oplus
+               end if
+            else
+               if (i_oplus /= 0) then
+                  write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if (chemthermod == 3)
+            ! CO+
+            i_coplus = igcm_coplus
+            if(chemthermod == 3) then
+               if (i_coplus == 0) then
+                  write(*,*) "calchim: Error; no CO+ tracer !!!"
+                  write(*,*) "CO+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_coplus
+               end if
+            else
+               if (i_coplus /= 0) then
+                  write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
+                  write(*,*) " Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   ! Of if (chemthermod == 3)
+            ! C+
+            i_cplus = igcm_cplus
+            if(chemthermod == 3) then
+               if (i_cplus == 0) then
+                  write(*,*) "calchim: Error; no C+ tracer !!!"
+                  write(*,*) "C+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_cplus
+               end if
+            else
+               if (i_cplus /= 0) then
+                  write(*,*) "calchim: Error; C+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   ! Of if (chemthermod == 3)
+            ! N+
+            i_nplus = igcm_nplus
+            if(chemthermod == 3) then
+               if (i_nplus == 0) then
+                  write(*,*) "calchim: Error; no N+ tracer !!!"
+                  write(*,*) "N+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_nplus
+               end if
+            else
+               if (i_nplus /= 0) then
+                  write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if (chemthermod == 3)
+            ! NO+
+            i_noplus = igcm_noplus
+            if(chemthermod == 3) then
+               if (i_noplus == 0) then
+                  write(*,*) "calchim: Error; no NO+ tracer !!!"
+                  write(*,*) "NO+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_noplus
+               end if
+            else
+               if (i_noplus /= 0) then
+                  write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if (chemthermod == 3)
+            ! N2+
+            i_n2plus = igcm_n2plus
+            if (chemthermod == 3) then
+               if (i_n2plus == 0) then
+                  write(*,*) "calchim: Error; no N2+ tracer !!!"
+                  write(*,*) "N2+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_n2plus
+               end if
+            else
+               if (i_n2plus /= 0) then
+                  write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if (chemthermod == 3)
+            !H+
+            i_hplus = igcm_hplus
+            if (chemthermod == 3) then
+               if (i_hplus == 0) then
+                  write(*,*) "calchim: Error; no H+ tracer !!!"
+                  write(*,*) "H+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_hplus
+               end if
+            else
+               if (i_hplus /= 0) then
+                  write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if (chemthermod == 3)
+            ! HCO2+
+            i_hco2plus = igcm_hco2plus
+            if(chemthermod == 3) then
+               if (i_hco2plus == 0) then
+                  write(*,*) "calchim: Error; no HCO2+ tracer !!!"
+                  write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_hco2plus
+               end if
+            else
+               if (i_hco2plus /= 0) then
+                  write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif    !Of if(chemthermod == 3)
+            !e-
+            i_elec = igcm_elec
+            if(chemthermod == 3) then
+               if (i_elec == 0) then
+                  write(*,*) "calchim: Error; no e- tracer !!!"
+                  write(*,*) "e- is needed if O2+ is in traceur.def"
+                  stop
+               else
+                  nbq = nbq + 1
+                  niq(nbq) = i_elec
+               end if
+            else
+               if(i_elec /= 0) then
+                  write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
+                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
+                  stop
+               endif
+            endif   !Of if (chemthermod == 3)
+         endif      !Of thermochem
 
          write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
-         
+               
          firstcall = .false.
       end if ! if (firstcall)
@@ -279,6 +570,6 @@
 
       zycol(:,:)    = 0.
-      dqchim(:,:,:) = 0
-      dqschim(:,:)  = 0
+      dqchim(:,:,:) = 0.
+      dqschim(:,:)  = 0.
 
 !     latvl1= 22.27
@@ -292,5 +583,5 @@
 
       do ig = 1,ngridmx
-
+         
          foundswitch = 0
          do l = 1,nlayermx
@@ -316,5 +607,5 @@
 
             if (photochem .and. thermochem) then
-               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
+               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
                   lswitch = l
                   foundswitch = 1
@@ -357,8 +648,8 @@
 
 !        chemistry in upper atmosphere
-
+        
          if (thermochem) then
-            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,  &
-                             zlocal,szacol,ptimestep,zday)
+            call chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens,  &
+                             zpress,zlocal,szacol,ptimestep,zday)
          end if
 
@@ -392,4 +683,5 @@
             end do
          end do ! of do l = 1,nlayermx
+!         if(ig.eq.800)write(*,*)'calchim/686',dqchim(ig,27,23)
 
 !=======================================================================
@@ -419,2 +711,3 @@
       return
       end
+
