Index: trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90	(revision 703)
+++ trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90	(revision 705)
@@ -44,5 +44,5 @@
 !    Local variables :
 !    -----------------
-      integer :: l,iq
+      integer :: l
       integer,save :: nesptherm
       real,allocatable :: rm(:,:) 		!number density (cm-3)
@@ -462,8 +462,13 @@
 
       !Solar flux calculation
-      call flujo(solarcondate)!+zday/365.)
 
       !Photoabsorption coefficients
-      call jthermcalc(ig,chemthermod,rm,nesptherm,ztemp,zlocal,zenit)
+      if(solvarmod.eq.0) then
+         call flujo(solarcondate)!+zday/365.)
+         call jthermcalc(ig,chemthermod,rm,nesptherm,ztemp,zlocal,zenit)
+      else if(solvarmod.eq.1) then
+         call jthermcalc_e107(ig,chemthermod,rm,nesptherm,ztemp,zlocal,zenit,zday)
+      endif
+         
 
       !Chemistry
Index: trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90	(revision 703)
+++ trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90	(revision 705)
@@ -362,5 +362,5 @@
       
       !Solar flux calculation
-      call flujo(solarcondate)
+      if(solvarmod.eq.0) call flujo(solarcondate)
 
                                          ! Not recommended for long runs 
@@ -407,5 +407,5 @@
          enddo
         !Routine to calculate the UV heating
-         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,jtot)
+         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
 
 !        euveff=0.16    !UV heating efficiency. Following Fox et al. ASR 1996
Index: trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F	(revision 703)
+++ trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F	(revision 705)
@@ -1,5 +1,5 @@
 c**********************************************************************
 
-      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,jtot)
+      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday,jtot)
 
 
@@ -18,4 +18,5 @@
       include 'param.h'
       include 'param_v4.h'
+      include "callkeys.h"
 
 
@@ -38,4 +39,5 @@
       real       zenit
       real       iz(nlayermx)
+      real       zday
 
       ! tracer indexes for the EUV heating:
@@ -99,5 +101,12 @@
 
       !Calculation of photoabsortion coefficient
-      call jthermcalc(ig,euvmod,rm,nespeuv,tx,iz,zenit)
+      if(solvarmod.eq.0) then
+         call jthermcalc(ig,euvmod,rm,nespeuv,tx,iz,zenit)
+      else if (solvarmod.eq.1) then
+         call jthermcalc_e107(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday)
+         do indexint=1,ninter
+            fluxtop(indexint)=1.
+         enddo
+      endif
 
       !Total photoabsorption coefficient
Index: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F	(revision 705)
+++ trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F	(revision 705)
@@ -0,0 +1,1117 @@
+c**********************************************************************
+
+      subroutine jthermcalc_e107
+     $     (ig,chemthermod,rm,nesptherm,tx,iz,zenit,zday)
+
+
+c     feb 2002        fgg           first version
+c     nov 2002        fgg           second version
+c
+c modified from paramhr.F
+c MAC July 2003
+c**********************************************************************
+
+      implicit none
+
+c     common variables and constants
+      include "dimensions.h"
+      include "dimphys.h"
+      include 'param.h'
+      include 'param_v4.h'
+
+c     input and output variables
+
+      integer    ig
+      integer    chemthermod
+      integer    nesptherm                      !Number of species considered
+      real       rm(nlayermx,nesptherm)         !Densities (cm-3)
+      real       tx(nlayermx)                   !temperature
+      real       zenit                          !SZA
+      real       iz(nlayermx)                   !Local altitude
+      real       zday                           !Martian day after Ls=0
+
+
+c    local parameters and variables
+
+      real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
+      real       o2colx(nlayermx)               !column density of O2(cm^-2)
+      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2)
+      real       h2colx(nlayermx)               !H2 column density (cm-2)
+      real       h2ocolx(nlayermx)              !H2O column density (cm-2)
+      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2)
+      real       o3colx(nlayermx)               !O3 column density (cm-2)
+      real       n2colx(nlayermx)               !N2 column density (cm-2)
+      real       ncolx(nlayermx)                !N column density (cm-2)
+      real       nocolx(nlayermx)               !NO column density (cm-2)
+      real       cocolx(nlayermx)               !CO column density (cm-2)
+      real       hcolx(nlayermx)                !H column density (cm-2)
+      real       no2colx(nlayermx)              !NO2 column density (cm-2)
+      real       t2(nlayermx)
+      real       coltemp(nlayermx)
+      real       sigma(ninter,nlayermx)
+      real       alfa(ninter,nlayermx)
+      real       realday
+      
+      integer    i,j,k,indexint                 !indexes
+      character  dn
+      integer    tinf,tsup
+
+
+
+c     variables used in interpolation
+
+      real*8      auxcoltab(nz2)
+      real*8      auxjco2(nz2)
+      real*8      auxjo2(nz2)
+      real*8      auxjo3p(nz2)
+      real*8      auxjh2o(nz2)
+      real*8      auxjh2(nz2)
+      real*8      auxjh2o2(nz2)
+      real*8      auxjo3(nz2)
+      real*8      auxjn2(nz2)
+      real*8      auxjn(nz2)
+      real*8      auxjno(nz2)
+      real*8      auxjco(nz2)
+      real*8      auxjh(nz2)
+      real*8      auxjno2(nz2)
+      real*8      wp(nlayermx),wm(nlayermx)
+      real*8      auxcolinp(nlayermx)
+      integer     auxind(nlayermx)
+      integer     auxi
+      integer     ind
+      real*8      cortemp(nlayermx)
+
+      real*8      limdown                      !limits for interpolation
+      real*8      limup                        !        ""
+
+      !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
+      !!!If the value is changed there, if has to be changed also here !!!!
+      integer,parameter :: i_co2=1
+
+
+c*************************PROGRAM STARTS*******************************
+      
+      if(zenit.gt.140.) then
+         dn='n'
+         else
+         dn='d'
+      end if
+      if(dn.eq.'n') then
+        return
+      endif
+      
+      !Initializing the photoabsorption coefficients
+      jfotsout(:,:,:)=0.
+
+      !Auxiliar temperature to take into account the temperature dependence
+      !of CO2 cross section
+      do i=1,nlayermx
+         t2(i)=tx(i)
+         if(t2(i).lt.195.0) t2(i)=195.0
+         if(t2(i).gt.295.0) t2(i)=295.0
+      end do
+
+      !Calculation of column amounts 
+      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
+     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
+     $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
+
+      !Auxiliar column to include the temperature dependence 
+      !of CO2 cross section
+      coltemp(nlayermx)=co2colx(nlayermx)*abs(t2(nlayermx)-t0(nlayermx))
+      do i=nlayermx-1,1,-1
+        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR 
+     $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5 
+     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
+      end do
+      
+      !Calculation of CO2 cross section at temperature t0(i)
+      do i=1,nlayermx
+         do indexint=24,32
+           sigma(indexint,i)=co2crsc195(indexint-23)
+           alfa(indexint,i)=((co2crsc295(indexint-23)
+     $          /sigma(indexint,i))-1.)/(295.-t0(i))
+        end do
+      end do
+
+      !E10.7 for the day: linear interpolation to tabulated values
+      realday=mod(zday,669.)
+      if(realday.lt.date_e107(1)) then
+         e107=e107_tab(1)
+      else if(realday.ge.date_e107(1).and.
+     $        realday.le.date_e107(669)) then
+         do i=1,669
+            if(realday.ge.date_e107(i).and.
+     $           realday.lt.date_e107(i+1)) then
+               tinf=i
+               tsup=i+1
+               e107=e107_tab(tinf)+(zday-date_e107(tinf))*
+     $              (e107_tab(tsup)-e107_tab(tinf))
+            endif
+         enddo
+      endif
+      !Photoabsorption coefficients at TOA as a function of E10.7
+      do j=1,nabs
+         do indexint=1,ninter
+            jfotsout(indexint,j,nlayermx)=coefit0(indexint,j)+
+     $           coefit1(indexint,j)*e107+coefit2(indexint,j)*e107**2+
+     $           coefit3(indexint,j)*e107**3+coefit4(indexint,j)*e107**4
+         enddo
+      enddo
+
+! Interpolation to the tabulated photoabsorption coefficients for each species
+! in each spectral interval
+
+
+c     auxcolinp-> Actual atmospheric column
+c     auxj*-> Tabulated photoabsorption coefficients
+c     auxcoltab-> Tabulated atmospheric columns
+
+ccccccccccccccccccccccccccccccc
+c     0.1,5.0 (int 1)
+c
+c     Absorption by: 
+c     CO2, O2, O, H2, N
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+      indexint=1
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
+     $        o2colx(i)*crscabsi2(2,indexint) + 
+     $        o3pcolx(i)*crscabsi2(3,indexint) + 
+     $        h2colx(i)*crscabsi2(5,indexint) + 
+     $        ncolx(i)*crscabsi2(9,indexint)
+      end do
+      limdown=1.e-20
+      limup=1.e26
+
+
+c     Interpolations
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !CO2 tabulated coefficient
+         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+         !O2 tabulated coefficient
+         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+         !O3p tabulated coefficient
+         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
+         !H2 tabulated coefficient
+         auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
+         !Tabulated column
+         auxcoltab(i) = c1_16(auxi,indexint)
+      enddo
+      !Only if chemthermod.ge.2
+      !N tabulated coefficient
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
+         enddo
+      endif
+
+      call interfast 
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         auxi=nlayermx-i+1
+         !CO2 interpolated coefficient
+         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayermx) *
+     $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
+         !O2 interpolated coefficient
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) *
+     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
+         !O3p interpolated coefficient
+         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayermx) *
+     $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
+         !H2 interpolated coefficient
+         jfotsout(indexint,5,auxi) = jfotsout(indexint,5,auxi) *
+     $        (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
+      enddo
+      !Only if chemthermod.ge.2
+      !N interpolated coefficient
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            ind=auxind(i)
+            jfotsout(indexint,9,nlayermx-i+1) =  
+     $           jfotsout(indexint,9,nlayermx) * 
+     $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
+         enddo
+      endif
+         
+
+c     End interval 1
+
+
+ccccccccccccccccccccccccccccccc
+c     5-80.5nm (int 2-15)
+c
+c     Absorption by:
+c     CO2, O2, O, H2, N2, N, 
+c     NO, CO, H, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+      do indexint=2,15
+         do i=1,nlayermx
+            auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
+     $           o2colx(i)*crscabsi2(2,indexint)+
+     $           o3pcolx(i)*crscabsi2(3,indexint)+
+     $           h2colx(i)*crscabsi2(5,indexint)+
+     $           n2colx(i)*crscabsi2(8,indexint)+
+     $           ncolx(i)*crscabsi2(9,indexint)+
+     $           nocolx(i)*crscabsi2(10,indexint)+
+     $           cocolx(i)*crscabsi2(11,indexint)+
+     $           hcolx(i)*crscabsi2(12,indexint)+
+     $           no2colx(i)*crscabsi2(13,indexint)
+         end do
+
+c     Interpolations
+
+         do i=1,nz2
+            auxi = nz2-i+1
+            !O2 tabulated coefficient
+            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+            !O3p tabulated coefficient
+            auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
+            !CO2 tabulated coefficient
+            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+            !H2 tabulated coefficient
+            auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
+            !N2 tabulated coefficient
+            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
+            !CO tabulated coefficient
+            auxjco(i) = jabsifotsintpar(auxi,11,indexint)
+            !H tabulated coefficient
+            auxjh(i) = jabsifotsintpar(auxi,12,indexint)
+            !tabulated column
+            auxcoltab(i) = c1_16(auxi,indexint)
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nz2
+               auxi = nz2-i+1
+               !N tabulated coefficient
+               auxjn(i) = jabsifotsintpar(auxi,9,indexint)
+               !NO tabulated coefficient
+               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+               !NO2 tabulated coefficient
+               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+            enddo
+         endif
+
+          call interfast(wm,wp,auxind,auxcolinp,nlayermx,
+     $        auxcoltab,nz2,limdown,limup)
+          do i=1,nlayermx
+             ind=auxind(i)
+             auxi = nlayermx-i+1
+             !O2 interpolated coefficient
+             jfotsout(indexint,2,auxi) = 
+     $            jfotsout(indexint,2,nlayermx) *
+     $            (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
+             !O3p interpolated coefficient
+             jfotsout(indexint,3,auxi) = 
+     $            jfotsout(indexint,3,nlayermx) *
+     $            (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
+             !CO2 interpolated coefficient
+             jfotsout(indexint,1,auxi) = 
+     $            jfotsout(indexint,1,nlayermx) *
+     $            (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
+             !H2 interpolated coefficient
+             jfotsout(indexint,5,auxi) = 
+     $            jfotsout(indexint,5,nlayermx) *
+     $            (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
+             !N2 interpolated coefficient
+             jfotsout(indexint,8,auxi) = 
+     $            jfotsout(indexint,8,nlayermx) *
+     $            (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
+             !CO interpolated coefficient
+             jfotsout(indexint,11,auxi) = 
+     $            jfotsout(indexint,11,nlayermx) *
+     $            (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
+             !H interpolated coefficient
+             jfotsout(indexint,12,auxi) = 
+     $            jfotsout(indexint,12,nlayermx) *
+     $            (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
+          enddo
+          !Only if chemthermod.ge.2
+          if(chemthermod.ge.2) then
+             do i=1,nlayermx
+                ind=auxind(i)
+                auxi = nlayermx-i+1
+                !N interpolated coefficient
+                jfotsout(indexint,9,auxi) = 
+     $               jfotsout(indexint,9,nlayermx) *
+     $               (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
+                !NO interpolated coefficient
+                jfotsout(indexint,10,auxi)=
+     $               jfotsout(indexint,10,nlayermx) *
+     $               (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
+                !NO2 interpolated coefficient
+                jfotsout(indexint,13,auxi)=
+     $               jfotsout(indexint,13,nlayermx) *
+     $               (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
+             enddo
+          endif   
+      end do
+
+c     End intervals 2-15
+
+
+ccccccccccccccccccccccccccccccc
+c     80.6-90.8nm (int16)
+c
+c     Absorption by:
+c     CO2, O2, O, N2, N, NO,
+c     CO, H, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+      indexint=16
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
+     $        o2colx(i)*crscabsi2(2,indexint)+
+     $        o3pcolx(i)*crscabsi2(3,indexint)+
+     $        n2colx(i)*crscabsi2(8,indexint)+
+     $        ncolx(i)*crscabsi2(9,indexint)+
+     $        nocolx(i)*crscabsi2(10,indexint)+
+     $        cocolx(i)*crscabsi2(11,indexint)+
+     $        hcolx(i)*crscabsi2(12,indexint)+
+     $        no2colx(i)*crscabsi2(13,indexint)
+      end do
+
+c     Interpolations
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !O2 tabulated coefficient
+         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+         !CO2 tabulated coefficient
+         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+         !O3p tabulated coefficient
+         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
+         !N2 tabulated coefficient
+         auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
+         !CO tabulated coefficient
+         auxjco(i) = jabsifotsintpar(auxi,11,indexint)
+         !H tabulated coefficient
+         auxjh(i) = jabsifotsintpar(auxi,12,indexint)
+         !NO2 tabulated coefficient
+         auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+         !Tabulated column
+         auxcoltab(i) = c1_16(auxi,indexint)
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            auxi = nz2-i+1
+            !N tabulated coefficient
+            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
+            !NO tabulated coefficient
+            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+            !NO2 tabulated coefficient
+            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+         enddo
+      endif
+
+      call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         auxi = nlayermx-i+1
+         !O2 interpolated coefficient
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) *
+     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
+         !CO2 interpolated coefficient
+         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayermx) *
+     $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
+         !O3p interpolated coefficient
+         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayermx) *
+     $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
+         !N2 interpolated coefficient
+         jfotsout(indexint,8,auxi) = jfotsout(indexint,8,nlayermx) *
+     $        (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
+         !CO interpolated coefficient
+         jfotsout(indexint,11,auxi) = 
+     $        jfotsout(indexint,11,nlayermx) *
+     $        (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
+         !H interpolated coefficient
+         jfotsout(indexint,12,auxi) = 
+     $        jfotsout(indexint,12,nlayermx) *
+     $        (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            ind=auxind(i)
+            auxi = nlayermx-i+1
+            !N interpolated coefficient
+            jfotsout(indexint,9,auxi) = 
+     $           jfotsout(indexint,9,nlayermx) *
+     $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
+            !NO interpolated coefficient
+            jfotsout(indexint,10,auxi) = 
+     $           jfotsout(indexint,10,nlayermx) * 
+     $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,auxi) = 
+     $           jfotsout(indexint,13,nlayermx) *
+     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
+         enddo
+      endif
+c     End interval 16
+
+
+ccccccccccccccccccccccccccccccc
+c     90.9-119.5nm (int 17-24)
+c
+c     Absorption by:
+c     CO2, O2, N2, NO, CO, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input column
+
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
+     $        nocolx(i) + cocolx(i) + no2colx(i)
+      end do
+
+      do indexint=17,24
+
+c     Interpolations
+
+         do i=1,nz2
+            auxi = nz2-i+1
+            !CO2 tabulated coefficient
+            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+            !O2 tabulated coefficient
+            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+            !N2 tabulated coefficient
+            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
+            !CO tabulated coefficient
+            auxjco(i) = jabsifotsintpar(auxi,11,indexint)            
+            !Tabulated column
+            auxcoltab(i) = c17_24(auxi)
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nz2
+               auxi = nz2-i+1
+               !NO tabulated coefficient
+               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+               !NO2 tabulated coefficient
+               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+            enddo
+         endif
+
+         call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+         !Correction to include T variation of CO2 cross section
+         if(indexint.eq.24) then
+            do i=1,nlayermx
+               auxi = nlayermx-i+1
+               if(sigma(indexint,auxi)*
+     $              alfa(indexint,auxi)*coltemp(auxi)
+     $              .lt.60.) then
+                  cortemp(i)=exp(-sigma(indexint,auxi)*
+     $                alfa(indexint,auxi)*coltemp(auxi))
+               else 
+                  cortemp(i)=0.
+               end if
+            enddo
+         else
+            do i=1,nlayermx
+               cortemp(i)=1.
+            enddo
+         end if
+         do i=1,nlayermx           
+            ind=auxind(i)
+            auxi = nlayermx-i+1
+            !O2 interpolated coefficient
+            jfotsout(indexint,2,auxi) = 
+     $           jfotsout(indexint,2,nlayermx) *
+     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
+     $           cortemp(i)
+            !CO2 interpolated coefficient
+            jfotsout(indexint,1,auxi) = 
+     $           jfotsout(indexint,1,nlayermx) *
+     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) 
+     $           * cortemp(i)
+            if(indexint.eq.24) jfotsout(indexint,1,auxi)=
+     $           jfotsout(indexint,1,auxi)*
+     $           (1+alfa(indexint,auxi)*
+     $           (t2(auxi)-t0(auxi)))
+            !N2 interpolated coefficient
+            jfotsout(indexint,8,auxi) = 
+     $           jfotsout(indexint,8,nlayermx) *
+     $           (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind)) * 
+     $           cortemp(i)            
+            !CO interpolated coefficient
+            jfotsout(indexint,11,auxi) = 
+     $           jfotsout(indexint,11,nlayermx) *
+     $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) * 
+     $           cortemp(i)            
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nlayermx
+               ind=auxind(i)
+               auxi = nlayermx-i+1
+               !NO interpolated coefficient
+               jfotsout(indexint,10,auxi)=
+     $              jfotsout(indexint,10,nlayermx) *
+     $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) * 
+     $              cortemp(i)
+               !NO2 interpolated coefficient
+               jfotsout(indexint,13,auxi)=
+     $              jfotsout(indexint,13,nlayermx) *
+     $              (wm(i)*auxjno2(ind+1)+ wp(i)*auxjno2(ind)) * 
+     $              cortemp(i)
+            enddo
+         endif               
+      end do
+c     End intervals 17-24
+
+
+ccccccccccccccccccccccccccccccc
+c     119.6-167.0nm (int 25-29)
+c
+c     Absorption by:
+c     CO2, O2, H2O, H2O2, NO,
+c     CO, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
+     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
+      end do
+
+      do indexint=25,29
+
+c     Interpolations
+
+         do i=1,nz2
+            auxi = nz2-i+1
+            !CO2 tabulated coefficient
+            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+            !O2 tabulated coefficient
+            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+            !H2O tabulated coefficient
+            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
+            !H2O2 tabulated coefficient
+            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)            
+            !CO tabulated coefficient
+            auxjco(i) = jabsifotsintpar(auxi,11,indexint)            
+            !Tabulated column
+            auxcoltab(i) = c25_29(auxi)
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nz2
+               auxi = nz2-i+1
+               !NO tabulated coefficient
+               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+               !NO2 tabulated coefficient
+               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+            enddo
+         endif
+         call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayermx
+            ind=auxind(i)
+            auxi = nlayermx-i+1
+            !Correction to include T variation of CO2 cross section
+            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
+     $           coltemp(auxi).lt.60.) then
+               cortemp(i)=exp(-sigma(indexint,auxi)*
+     $              alfa(indexint,auxi)*coltemp(auxi))
+            else 
+               cortemp(i)=0.
+            end if
+            !CO2 interpolated coefficient
+            jfotsout(indexint,1,auxi) = 
+     $           jfotsout(indexint,1,nlayermx) *
+     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) * 
+     $           cortemp(i) *
+     $           (1+alfa(indexint,auxi)*
+     $           (t2(auxi)-t0(auxi)))
+            !O2 interpolated coefficient
+            jfotsout(indexint,2,auxi) = 
+     $           jfotsout(indexint,2,nlayermx) *
+     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
+     $           cortemp(i)
+            !H2O interpolated coefficient
+            jfotsout(indexint,4,auxi) = 
+     $           jfotsout(indexint,4,nlayermx) *
+     $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) * 
+     $           cortemp(i)
+            !H2O2 interpolated coefficient
+            jfotsout(indexint,6,auxi) = 
+     $           jfotsout(indexint,6,nlayermx) *
+     $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
+     $           cortemp(i)            
+            !CO interpolated coefficient
+            jfotsout(indexint,11,auxi) = 
+     $           jfotsout(indexint,11,nlayermx) *
+     $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) * 
+     $           cortemp(i)
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nlayermx
+               ind=auxind(i)
+               auxi = nlayermx-i+1
+               !NO interpolated coefficient
+               jfotsout(indexint,10,auxi)=
+     $              jfotsout(indexint,10,nlayermx) *
+     $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) * 
+     $              cortemp(i)
+               !NO2 interpolated coefficient
+               jfotsout(indexint,13,auxi)=
+     $              jfotsout(indexint,13,nlayermx) *
+     $              (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) * 
+     $              cortemp(i)
+            enddo
+         endif
+
+      end do
+
+c     End intervals 25-29
+
+
+cccccccccccccccccccccccccccccccc
+c     167.1-202.5nm (int 30-31)
+c    
+c     Absorption by:
+c     CO2, O2, H2O, H2O2, NO,
+c     NO2
+cccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
+     $        h2o2colx(i) + nocolx(i) + no2colx(i)
+      end do
+
+c     Interpolation
+
+      do indexint=30,31
+
+         do i=1,nz2
+            auxi = nz2-i+1
+            !CO2 tabulated coefficient
+            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+            !O2 tabulated coefficient
+            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+            !H2O tabulated coefficient
+            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
+            !H2O2 tabulated coefficient
+            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)            
+            !Tabulated column
+            auxcoltab(i) = c30_31(auxi)
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nz2
+               auxi = nz2-i+1
+               !NO tabulated coefficient
+               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+               !NO2 tabulated coefficient
+               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+            enddo
+         endif
+
+         call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayermx
+            ind=auxind(i)
+            auxi = nlayermx-i+1
+            !Correction to include T variation of CO2 cross section
+            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
+     $           coltemp(auxi).lt.60.) then
+               cortemp(i)=exp(-sigma(indexint,auxi)*
+     $              alfa(indexint,auxi)*coltemp(auxi))
+            else 
+               cortemp(i)=0.
+            end if
+            !CO2 interpolated coefficient
+            jfotsout(indexint,1,auxi) = 
+     $           jfotsout(indexint,1,nlayermx) *
+     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) * 
+     $           cortemp(i) *
+     $           (1+alfa(indexint,auxi)*
+     $           (t2(auxi)-t0(auxi)))
+            !O2 interpolated coefficient
+            jfotsout(indexint,2,auxi) = 
+     $           jfotsout(indexint,2,nlayermx) *
+     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
+     $           cortemp(i)
+            !H2O interpolated coefficient
+            jfotsout(indexint,4,auxi) = 
+     $           jfotsout(indexint,4,nlayermx) *
+     $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) * 
+     $           cortemp(i)
+            !H2O2 interpolated coefficient
+            jfotsout(indexint,6,auxi) = 
+     $           jfotsout(indexint,6,nlayermx) *
+     $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
+     $           cortemp(i)            
+         enddo
+         !Only if chemthermod.ge.2
+         if(chemthermod.ge.2) then
+            do i=1,nlayermx 
+               ind=auxind(i)
+               auxi = nlayermx-i+1
+               !NO interpolated coefficient
+               jfotsout(indexint,10,auxi)=
+     $              jfotsout(indexint,10,nlayermx) *
+     $              (wm(i)*auxjno(ind+1) +wp(i)*auxjno(ind)) * 
+     $              cortemp(i)
+               !NO2 interpolated coefficient
+               jfotsout(indexint,13,auxi)=
+     $              jfotsout(indexint,13,auxi) *
+     $              (wm(i)*auxjno2(ind+1)+wp(i)*auxjno2(ind)) * 
+     $              cortemp(i)
+            enddo
+         endif
+
+      end do
+
+c     End intervals 30-31
+
+
+ccccccccccccccccccccccccccccccc
+c     202.6-210.0nm (int 32)
+c
+c     Absorption by:
+c     CO2, O2, H2O2, NO, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      indexint=32
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) + 
+     $        nocolx(i) + no2colx(i)
+      end do
+
+c     Interpolation
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !CO2 tabulated coefficient
+         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
+         !O2 tabulated coefficient
+         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+         !H2O2 tabulated coefficient
+         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)         
+         !Tabulated column
+         auxcoltab(i) = c32(auxi)
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            auxi = nz2-i+1
+            !NO tabulated coefficient
+            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+            !NO2 tabulated coefficient
+            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
+         enddo
+      endif
+      call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         auxi = nlayermx-i+1
+         !Correction to include T variation of CO2 cross section
+         if(sigma(indexint,nlayermx-i+1)*alfa(indexint,auxi)*
+     $        coltemp(auxi).lt.60.) then
+            cortemp(i)=exp(-sigma(indexint,auxi)*
+     $           alfa(indexint,auxi)*coltemp(auxi))
+         else 
+            cortemp(i)=0.
+         end if
+         !CO2 interpolated coefficient
+         jfotsout(indexint,1,auxi) = 
+     $        jfotsout(indexint,1,nlayermx) *
+     $        (wm(i)*auxjco2(ind+1)+wp(i)*auxjco2(ind)) * 
+     $        cortemp(i) *
+     $        (1+alfa(indexint,auxi)*
+     $        (t2(auxi)-t0(auxi)))
+         !O2 interpolated coefficient
+         jfotsout(indexint,2,auxi) = 
+     $        jfotsout(indexint,2,nlayermx) *
+     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
+     $        cortemp(i)
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = 
+     $        jfotsout(indexint,6,nlayermx) *
+     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
+     $        cortemp(i)         
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            auxi = nlayermx-i+1
+            ind=auxind(i)
+            !NO interpolated coefficient
+            jfotsout(indexint,10,auxi) = 
+     $           jfotsout(indexint,10,nlayermx) * 
+     $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) * 
+     $           cortemp(i)
+           !NO2 interpolated coefficient
+            jfotsout(indexint,13,auxi) = 
+     $           jfotsout(indexint,13,nlayermx) *
+     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) * 
+     $           cortemp(i)
+         enddo
+      endif
+
+c     End of interval 32
+
+
+ccccccccccccccccccccccccccccccc
+c     210.1-231.0nm (int 33)
+c     
+c     Absorption by:
+c     O2, H2O2, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      indexint=33
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
+      end do
+
+c     Interpolation
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !O2 tabulated coefficient
+         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+         !H2O2 tabulated coefficient
+         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
+         !Tabulated column
+         auxcoltab(i) = c33(auxi)
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            !NO2 tabulated coefficient
+            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
+         enddo
+      endif
+      call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         auxi = nlayermx-i+1
+         !O2 interpolated coefficient
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) *
+     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayermx) *
+     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 
+     $           jfotsout(indexint,13,nlayermx) *
+     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
+         enddo
+      endif
+
+c     End of interval 33
+
+
+ccccccccccccccccccccccccccccccc
+c     231.1-240.0nm (int 34)
+c
+c     Absorption by:
+c     O2, H2O2, O3, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      indexint=34
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
+     $        no2colx(i)
+      end do
+
+c     Interpolation
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !O2 tabulated coefficient
+         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
+         !H2O2 tabulated coefficient
+         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
+         !O3 tabulated coefficient
+         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
+         !Tabulated column
+         auxcoltab(i) = c34(nz2-i+1)
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            !NO2 tabulated coefficient
+            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
+         enddo
+      endif
+      call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         auxi = nlayermx-i+1
+         !O2 interpolated coefficient
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) * 
+     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayermx) *
+     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
+         !O3 interpolated coefficient
+         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayermx) *
+     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 
+     $           jfotsout(indexint,13,nlayermx) *
+     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
+         enddo
+      endif
+
+c     End of interval 34      
+
+
+ccccccccccccccccccccccccccccccc
+c     240.1-337.7nm (int 35)
+c
+c     Absorption by:
+c     H2O2, O3, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      indexint=35
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
+      end do
+
+c     Interpolation
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !H2O2 tabulated coefficient
+         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
+         !O3 tabulated coefficient
+         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
+         !Tabulated column
+         auxcoltab(i) = c35(auxi)
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            !NO2 tabulated coefficient
+            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
+         enddo
+      endif
+      call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         auxi = nlayermx-i+1
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayermx) *
+     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
+         !O3 interpolated coefficient
+         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayermx) *
+     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
+      enddo
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 
+     $           jfotsout(indexint,13,nlayermx) *
+     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
+         enddo
+      endif
+
+c     End of interval 35
+
+ccccccccccccccccccccccccccccccc
+c     337.8-800.0 nm (int 36)
+c     
+c     Absorption by:
+c     O3, NO2
+ccccccccccccccccccccccccccccccc
+
+c     Input atmospheric column
+
+      indexint=36
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = o3colx(i) + no2colx(i)
+      end do
+
+c     Interpolation
+
+      do i=1,nz2
+         auxi = nz2-i+1
+         !O3 tabulated coefficient
+         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
+         !Tabulated column
+         auxcoltab(i) = c36(auxi)
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nz2
+            !NO2 tabulated coefficient
+            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
+         enddo
+      endif
+      call interfast
+     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayermx
+         ind=auxind(i)
+         !O3 interpolated coefficient
+         jfotsout(indexint,7,nlayermx-i+1) = 
+     $        jfotsout(indexint,7,nlayermx) *
+     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
+         do i=1,nlayermx
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 
+     $           jfotsout(indexint,13,nlayermx) *
+     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
+         enddo
+      endif
+
+c     End of interval 36
+
+c     End of interpolation to obtain photoabsorption rates
+
+
+      return
+
+      end
+
+
+
+
+
+
Index: trunk/LMDZ.MARS/libf/aeronomars/param_read.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_read.F	(revision 703)
+++ trunk/LMDZ.MARS/libf/aeronomars/param_read.F	(revision 705)
@@ -15,5 +15,5 @@
 
       integer    i,j,k,inter                          !indexes
-      integer ierr,lnblnk 
+      integer ierr
       real       nada
   
@@ -89,14 +89,13 @@
       !Tabulated column amount
       open(210, status = 'old',
-c    $file=datafile(1:lnblnk(datafile))//'/EUVDAT/coln.dat',iostat=ierr)
-     $file=datafile
-     $   (1:lnblnk(datafile))//'/EUVDAT/param_v5/coln.dat',iostat=ierr)
+c    $file=trim(datafile)//'/EUVDAT/coln.dat',iostat=ierr)
+     $file=trim(datafile)//'/EUVDAT/param_v5/coln.dat',iostat=ierr)
 
       IF (ierr.NE.0) THEN 
-       write(*,*)'cant find directory EUVDAT and content coln.dat'
+       write(*,*)'cant find directory EUVDAT containing param_v5 subdir'
        write(*,*)'(in aeronomars/param_read.F)'
-       write(*,*)'It should be in :', datafile,'/'
+       write(*,*)'It should be in :', trim(datafile),'/'
        write(*,*)'1) You can change this directory address in '
-       write(*,*)'   file phymars/datafile.h'
+       write(*,*)'   callphys.def with datadir=/path/to/dir'
        write(*,*)'2) If necessary, EUVDAT (and other datafiles)'
        write(*,*)'   can be obtained online on:'
@@ -106,56 +105,30 @@
  
       !Tabulated photoabsorption coefficients
-      open(220,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j2_an.dat')
-      open(230,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j3_an.dat')
-      open(240,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j1_an.dat')
-      open(250,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j2_bn.dat')
-      open(260,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j2_cn.dat')
-      open(300,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j2_dn.dat')
-      open(270,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j1_bn.dat')
-      open(280,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j1_cn.dat')
-      open(290,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j1_dn.dat')
-      open(150,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j4n.dat')
-      open(160,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j5n.dat')
-      open(170,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j6n.dat')
-      open(180,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j7n.dat')
-      open(390,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j8_an.dat')
-      open(400,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j8_bn.dat')
-      open(410,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j9n.dat')
-      open(420,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j10_an.dat')
-      open(430,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j10_bn.dat')
-      open(440,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j10_cn.dat')
-      open(450,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j11_an.dat')
-      open(460,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j11_bn.dat')
-      open(470,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j11_cn.dat')
-      open(480,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j12n.dat')
-      open(490,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j13_an.dat')
-      open(500,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j13_bn.dat')
-      open(510,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j13_cn.dat')
+      open(220,file=trim(datafile)//'/EUVDAT/param_v5/j2_an.dat')
+      open(230,file=trim(datafile)//'/EUVDAT/param_v5/j3_an.dat')
+      open(240,file=trim(datafile)//'/EUVDAT/param_v5/j1_an.dat')
+      open(250,file=trim(datafile)//'/EUVDAT/param_v5/j2_bn.dat')
+      open(260,file=trim(datafile)//'/EUVDAT/param_v5/j2_cn.dat')
+      open(300,file=trim(datafile)//'/EUVDAT//param_v5/j2_dn.dat')
+      open(270,file=trim(datafile)//'/EUVDAT//param_v5/j1_bn.dat')
+      open(280,file=trim(datafile)//'/EUVDAT//param_v5/j1_cn.dat')
+      open(290,file=trim(datafile)//'/EUVDAT//param_v5/j1_dn.dat')
+      open(150,file=trim(datafile)//'/EUVDAT//param_v5/j4n.dat')
+      open(160,file=trim(datafile)//'/EUVDAT//param_v5/j5n.dat')
+      open(170,file=trim(datafile)//'/EUVDAT//param_v5/j6n.dat')
+      open(180,file=trim(datafile)//'/EUVDAT//param_v5/j7n.dat')
+      open(390,file=trim(datafile)//'/EUVDAT//param_v5/j8_an.dat')
+      open(400,file=trim(datafile)//'/EUVDAT//param_v5/j8_bn.dat')
+      open(410,file=trim(datafile)//'/EUVDAT//param_v5/j9n.dat')
+      open(420,file=trim(datafile)//'/EUVDAT//param_v5/j10_an.dat')
+      open(430,file=trim(datafile)//'/EUVDAT//param_v5/j10_bn.dat')
+      open(440,file=trim(datafile)//'/EUVDAT//param_v5/j10_cn.dat')
+      open(450,file=trim(datafile)//'/EUVDAT//param_v5/j11_an.dat')
+      open(460,file=trim(datafile)//'/EUVDAT//param_v5/j11_bn.dat')
+      open(470,file=trim(datafile)//'/EUVDAT//param_v5/j11_cn.dat')
+      open(480,file=trim(datafile)//'/EUVDAT//param_v5/j12n.dat')
+      open(490,file=trim(datafile)//'/EUVDAT//param_v5/j13_an.dat')
+      open(500,file=trim(datafile)//'/EUVDAT//param_v5/j13_bn.dat')
+      open(510,file=trim(datafile)//'/EUVDAT//param_v5/j13_cn.dat')
 
       
@@ -310,6 +283,5 @@
 
       !Parameters for the variation of the solar flux with 11 years cycle
-      open(100,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/varflujo.dat')
+      open(100,file=trim(datafile)//'/EUVDAT/param_v5/varflujo.dat')
       read(100,*)
       do i=1,24
@@ -349,6 +321,5 @@
 c     CO2, O2, NO
 
-      open(120,file=datafile
-     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/efdis_inter.dat')
+      open(120,file=trim(datafile)//'/EUVDAT/param_v5/efdis_inter.dat')
       read(120,*)
 !      do i=1,21
Index: trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F	(revision 705)
+++ trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F	(revision 705)
@@ -0,0 +1,419 @@
+      subroutine param_read_e107
+
+      implicit none
+
+
+c     common variables and constants
+      include "dimensions.h"
+      include "dimphys.h"
+      include 'param.h'
+      include 'param_v4.h'
+      include 'datafile.h'
+      include "callkeys.h"
+ 
+ 
+c     local variables
+
+      integer    i,j,k,inter                          !indexes
+      integer ierr,lnblnk 
+      real       nada
+      character*13  filename
+  
+      
+c*************************+PROGRAM STARTS**************************
+
+
+c     Reads tabulated functions
+
+      !Tabulated column amount
+      open(210, status = 'old',
+c    $file=trim(datafile)//'/EUVDAT/coln.dat',iostat=ierr)
+     $file=trim(datafile)//'/EUVDAT/param_v6/coln.dat',iostat=ierr)
+
+      IF (ierr.NE.0) THEN 
+       write(*,*)'cant find directory EUVDAT containing param_v6 subdir'
+       write(*,*)'(in aeronomars/param_read.F)'
+       write(*,*)'It should be in :', trim(datafile),'/'
+       write(*,*)'1) You can change this directory address in '
+       write(*,*)'   callphys.def with datadir=/path/to/dir'
+       write(*,*)'2) If necessary, EUVDAT (and other datafiles)'
+       write(*,*)'   can be obtained online on:'
+       write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
+       STOP
+      ENDIF
+ 
+      !Tabulated photoabsorption coefficients
+      open(220,file=trim(datafile)//'/EUVDAT/param_v6/trans2_an.dat')
+      open(230,file=trim(datafile)//'/EUVDAT/param_v6/trans3_an.dat')
+      open(240,file=trim(datafile)//'/EUVDAT/param_v6/trans1_an.dat')
+      open(250,file=trim(datafile)//'/EUVDAT/param_v6/trans2_bn.dat')
+      open(260,file=trim(datafile)//'/EUVDAT/param_v6/trans2_cn.dat')
+      open(300,file=trim(datafile)//'/EUVDAT/param_v6/trans2_dn.dat')
+      open(270,file=trim(datafile)//'/EUVDAT/param_v6/trans1_bn.dat')
+      open(280,file=trim(datafile)//'/EUVDAT/param_v6/trans1_cn.dat')
+      open(290,file=trim(datafile)//'/EUVDAT/param_v6/trans1_dn.dat')
+      open(150,file=trim(datafile)//'/EUVDAT/param_v6/trans4n.dat')
+      open(160,file=trim(datafile)//'/EUVDAT/param_v6/trans5n.dat')
+      open(170,file=trim(datafile)//'/EUVDAT/param_v6/trans6n.dat')
+      open(180,file=trim(datafile)//'/EUVDAT/param_v6/trans7n.dat')
+      open(390,file=trim(datafile)//'/EUVDAT/param_v6/trans8_an.dat')
+      open(400,file=trim(datafile)//'/EUVDAT/param_v6/trans8_bn.dat')
+      open(410,file=trim(datafile)//'/EUVDAT/param_v6/trans9n.dat')
+      open(420,file=trim(datafile)//'/EUVDAT/param_v6/trans10_an.dat')
+      open(430,file=trim(datafile)//'/EUVDAT/param_v6/trans10_bn.dat')
+      open(440,file=trim(datafile)//'/EUVDAT/param_v6/trans10_cn.dat')
+      open(450,file=trim(datafile)//'/EUVDAT/param_v6/trans11_an.dat')
+      open(460,file=trim(datafile)//'/EUVDAT/param_v6/trans11_bn.dat')
+      open(470,file=trim(datafile)//'/EUVDAT/param_v6/trans11_cn.dat')
+      open(480,file=trim(datafile)//'/EUVDAT/param_v6/trans12n.dat')
+      open(490,file=trim(datafile)//'/EUVDAT/param_v6/trans13_an.dat')
+      open(500,file=trim(datafile)//'/EUVDAT/param_v6/trans13_bn.dat')
+      open(510,file=trim(datafile)//'/EUVDAT/param_v6/trans13_cn.dat')
+
+      
+      do i=210,300,10
+         read(i,*)
+         read(i,*)
+      end do
+
+      do i=150,180,10
+         read(i,*)
+         read(i,*)
+      end do
+
+      do i=390,510,10
+         read(i,*)
+         read(i,*)
+      enddo
+
+      do i=nz2,1,-1
+         read(210,*) (c1_16(i,j),j=1,16),c17_24(i),c25_29(i),c30_31(i),
+     $        c32(i),c33(i),c34(i),c35(i),c36(i)
+      end do
+
+      do i=nz2,1,-1
+         read(220,*) (jabsifotsintpar(i,2,j),j=1,16)
+      end do
+      
+      do i=nz2,1,-1
+         read(230,*) (jabsifotsintpar(i,3,j),j=1,16)
+      end do
+
+      do i=nz2,1,-1
+         read(240,*) (jabsifotsintpar(i,1,j),j=1,16)
+      end do
+
+      do i=nz2,1,-1
+         read(250,*) (jabsifotsintpar(i,2,j),j=17,24)
+      end do
+
+
+      do i=nz2,1,-1
+         read(260,*) (jabsifotsintpar(i,2,j),j=25,31)
+      end do
+
+      do i=nz2,1,-1
+         read(270,*) (jabsifotsintpar(i,1,j),j=17,24)
+      end do
+
+      do i=nz2,1,-1
+         read(280,*) (jabsifotsintpar(i,1,j),j=25,31)
+      end do
+
+      do i=nz2,1,-1
+         read(290,*) jabsifotsintpar(i,1,32)
+      end do
+
+      do i=nz2,1,-1
+         read(300,*) (jabsifotsintpar(i,2,j),j=32,34)
+      end do
+
+      do i=nz2,1,-1
+         read(160,*) (jabsifotsintpar(i,5,j),j=1,15)
+      end do
+
+      do i=nz2,1,-1
+         read(150,*) (jabsifotsintpar(i,4,j),j=25,31)
+      end do
+
+      do i=nz2,1,-1
+         read(170,*) (jabsifotsintpar(i,6,j),j=25,35)
+      end do
+
+      do i=nz2,1,-1
+         read(180,*) (jabsifotsintpar(i,7,j),j=34,36)
+      end do
+
+      do i=nz2,1,-1
+         read(390,*) (jabsifotsintpar(i,8,j),j=2,16)
+      enddo
+
+      do i=nz2,1,-1
+         read(400,*) (jabsifotsintpar(i,8,j),j=17,24)
+      enddo
+
+      do i=nz2,1,-1
+         read(410,*) (jabsifotsintpar(i,9,j),j=1,16)
+      enddo
+
+      do i=nz2,1,-1
+         read(420,*) (jabsifotsintpar(i,10,j),j=2,16)
+      enddo
+
+      do i=nz2,1,-1
+         read(430,*) (jabsifotsintpar(i,10,j),j=17,24)
+      enddo
+
+      do i=nz2,1,-1
+         read(440,*) (jabsifotsintpar(i,10,j),j=25,32)
+      enddo
+
+      do i=nz2,1,-1
+         read(450,*) (jabsifotsintpar(i,11,j),j=2,16)
+      enddo
+
+      do i=nz2,1,-1
+         read(460,*) (jabsifotsintpar(i,11,j),j=17,24)
+      enddo
+
+      do i=nz2,1,-1
+         read(470,*) (jabsifotsintpar(i,11,j),j=25,29)
+      enddo
+      
+      do i=nz2,1,-1
+         read(480,*) (jabsifotsintpar(i,12,j),j=2,16)
+      enddo
+
+      do i=nz2,1,-1
+         read(490,*) (jabsifotsintpar(i,13,j),j=2,16)
+      enddo
+      
+      do i=nz2,1,-1
+         read(500,*) (jabsifotsintpar(i,13,j),j=17,24)
+      enddo
+      
+      do i=nz2,1,-1
+         read(510,*) (jabsifotsintpar(i,13,j),j=25,36)
+      enddo
+
+      do i=210,300,10
+         close(i)
+      end do
+
+      do i=150,180,10
+         close(i)
+      end do
+
+      do i=390,510,10
+         close(i)
+      enddo
+
+
+c     set t0
+
+      do i=1,nz2
+         t0(i)=195.
+      end do
+
+
+      do i=1,ninter
+         fluxtop(i)=1.
+      end do
+
+      !Parameters for the variation of the solar flux with 11 years cycle
+      open(620,file=trim(datafile)//'/EUVDAT/param_v6/fit_js_e107.dat')
+      do i=1,ninter
+         read(620,*)
+         do j=1,nabs
+            read(620,*)coefit0(i,j),coefit1(i,j),coefit2(i,j),
+     $           coefit3(i,j),coefit4(i,j)
+         enddo
+      enddo
+      close(620)
+
+
+      !Tabulated values for E10.7
+      if(solvaryear.eq.23) then
+         filename="e107_MY23.dat"
+      else if(solvaryear.eq.24) then
+         filename="e107_MY24.dat"
+      else if(solvaryear.eq.25) then
+         filename="e107_MY25.dat"
+      else if (solvaryear.eq.26) then
+         filename="e107_MY26.dat"
+      else if (solvaryear.eq.27) then
+         filename="e107_MY27.dat"
+      else if (solvaryear.eq.28) then
+         filename="e107_MY28.dat"
+      else if (solvaryear.eq.29) then
+         filename="e107_MY29.dat"
+      else if(solvaryear.eq.30) then
+         filename="e107_MY30.dat"
+      else
+         write(*,*)"param_read_e107: "
+         write(*,*)"bad value for solvaryear in callphys.def"
+         write(*,*)"solvaryear must be between 24 and 30"
+         stop
+      endif
+      
+      open(640,file=trim(datafile)//'/EUVDAT/param_v6/'//filename)
+      read(640,*)
+      do i=1,669
+         read(640,*)date_e107(i),e107_tab(i)
+         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
+      enddo
+      close(640)
+         
+
+c     dissociation and ionization efficiencies
+
+      do inter=1,ninter
+         efdisco2(inter)=0.
+         efdiso2(inter)=0.
+         efdish2(inter)=0.
+         efdish2o(inter)=0.
+         efdish2o2(inter)=0.
+         efdiso3(inter)=0.
+         efdisco(inter)=0.
+         efdisn2(inter)=0.
+         efdisno(inter)=0.
+         efdisno2(inter)=0.
+         efionco2(inter,1)=0.
+         efionco2(inter,2)=0.
+         efionco2(inter,3)=0.
+         efionco2(inter,4)=0.
+         efiono3p(inter)=0.
+         efionn2(inter,1)=0.
+         efionn2(inter,2)=0.
+         efionco(inter,1)=0.
+         efionco(inter,2)=0.
+         efionco(inter,3)=0.
+         efionn(inter)=0.
+         efionh(inter)=0.
+         efionno(inter)=0.
+      enddo
+
+
+c     CO2, O2, NO
+
+      open(120,file=trim(datafile)//'/EUVDAT/param_v5/efdis_inter.dat')
+      read(120,*)
+!      do i=1,21
+!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
+      do inter=8,28
+         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
+      enddo
+      do inter=29,ninter
+         efdisco2(inter)=1.
+         efdiso2(inter)=1.
+         efdisno(inter)=1.
+      enddo
+
+
+c     N2
+
+      efdisn2(15)=0.1
+      do inter=16,ninter
+         efdisn2(inter)=1.
+      enddo
+
+
+c     CO
+
+      efdisco(16)=0.5
+      do inter=17,ninter
+         efdisco(inter)=1.
+      enddo
+
+      
+c     O, N, H
+
+      do inter=1,ninter
+         efdiso(inter)=0.
+         efdisn(inter)=0.
+         efdish(inter)=0.
+      enddo
+
+
+c     H2O, H2O2, O3, NO2
+
+      do inter=25,31
+         efdish2o(inter)=1.
+      enddo
+      do inter=25,35
+         efdish2o2(inter)=1.
+      enddo
+      do inter=34,36
+         efdiso3(inter)=1.
+      enddo
+      do inter=27,36
+         efdisno2(inter)=1.
+      enddo
+      do inter=1,15
+         efdish2(inter)=1.
+      enddo
+         
+      !4 possible channels for CO2 ionization
+      do inter=14,16
+         efionco2(inter,1)=1.-efdisco2(inter)
+      enddo
+      efionco2(13,1)=0.805*(1.-efdisco2(13))
+      efionco2(13,2)=0.195*(1.-efdisco2(13))
+      do inter=11,12
+         efionco2(inter,3)=1.-efdisco2(inter)
+      enddo
+      efionco2(10,3)=0.9*(1.-efdisco2(10))
+      efionco2(10,4)=0.1*(1.-efdisco2(10))
+      do inter=2,9
+         efionco2(inter,4)=1.-efdisco2(inter)
+      enddo
+
+      !For O(3p) total ionization under 91.1 nm
+      do inter=1,16
+         efiono3p(inter)=1.d0
+      enddo
+
+      !2 channels for N2 ionization
+      do inter=9,15
+         efionn2(inter,1)=1.-efdisn2(inter)
+      enddo
+      do inter=2,8
+         efionn2(inter,2)=1.-efdisn2(inter)
+      enddo
+      
+      !3 channels for CO ionization
+      do inter=11,16
+         efionco(inter,1)=1.-efdisco(inter)
+      enddo
+      efionco(10,1)=0.87*(1.-efdisco(10))
+      efionco(10,2)=0.13*(1.-efdisco(10))
+      do inter=8,9
+         efionco(inter,2)=1.-efdisco(inter)
+      enddo
+      efionco(7,2)=0.1*(1.-efdisco(7))
+      efionco(7,3)=0.9*(1.-efdisco(7))
+      do inter=2,6
+         efionco(inter,3)=1.-efdisco(inter)
+      enddo
+
+      !Total ionization under 85 nm for N
+      do inter=1,16
+         efionn(inter)=1.
+      enddo
+
+      !NO
+      do inter=2,28
+         efionno(inter)=1.-efdisno(inter)
+      enddo
+
+      !Total ionization under 90 nm for H
+      do inter=3,16
+         efionh(inter)=1.
+      enddo
+
+
+      return
+
+
+      end
+
Index: trunk/LMDZ.MARS/libf/aeronomars/param_v4.h
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_v4.h	(revision 703)
+++ trunk/LMDZ.MARS/libf/aeronomars/param_v4.h	(revision 705)
@@ -9,5 +9,6 @@
       common/uv/ crscabsi2,c1_16,c17_24,c25_29,c30_31,c32,                         &
      &   c33,c34,c35,c36,jfotsout,co2crsc195,co2crsc295,                           &
-     &   t0,fluxtop,freccen,jabsifotsintpar
+     &   t0,fluxtop,freccen,jabsifotsintpar,e107,date_e107,e107_tab,               &
+     &   coefit0,coefit1,coefit2,coefit3,coefit4
      
      
@@ -30,4 +31,8 @@
       real freccen(ninter)          !representative wavelenght
       real jabsifotsintpar(nz2,nabs,ninter)
+      real e107,date_e107(669),e107_tab(669)
+      real coefit0(ninter,nabs),coefit1(ninter,nabs)
+      real coefit2(ninter,nabs)
+      real coefit3(ninter,nabs),coefit4(ninter,nabs)
 
 
@@ -172,5 +177,5 @@
         real*8 tminnplus(nlayermx),tminnoplus(nlayermx)
         real*8 tminn2plus(nlayermx)
-        real*8 tminhplus(nlayermx),tminhcoplus(nlayermx)
+        real*8 tminhplus(nlayermx)
         real*8 tminhco2plus(nlayermx)
 
Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 703)
+++ trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 705)
@@ -16,5 +16,5 @@
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
-     &   ,dustbin,nltemodel,nircorr
+     &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear
      
       COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,euveff,     &
@@ -47,4 +47,6 @@
       integer ilwn
       integer ncouche
+      integer solvarmod   ! model for solar EUV variation
+      integer solvaryear  ! mars year for realisticly varying solar EUV 
 
       logical rayleigh
Index: trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 703)
+++ trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 705)
@@ -584,6 +584,14 @@
          write(*,*) " thermochem = ",thermochem
 
+         write(*,*) "Method to include solar variability"
+         write(*,*) "0-> old method (using solarcondate); ",
+     &                  "1-> variability wit E10.7"
+         solvarmod=1
+         call getin("solvarmod",solvarmod)
+         write(*,*) " solvarmod = ",solvarmod
+
          write(*,*) "date for solar flux calculation:",
-     &   " (1985 < date < 2002)"
+     &   " (1985 < date < 2002)",
+     $   " (Only used if solvarmod=0)"
          write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
          solarcondate=1993.4 ! default value
@@ -591,4 +599,10 @@
          write(*,*) " solarcondate = ",solarcondate
          
+         write(*,*) "Solar variability as observed for MY: "
+         write(*,*) "Only if solvarmod=1"
+         solvaryear=24
+         call getin("solvaryear",solvaryear)
+         write(*,*) " solvaryear = ",solvaryear
+
          write(*,*) "UV heating efficiency:",
      &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 703)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 705)
@@ -308,4 +308,5 @@
       real rho(ngridmx,nlayermx)  ! density
       real vmr(ngridmx,nlayermx)  ! volume mixing ratio
+      real rhopart(ngridmx,nlayermx) ! number density of a given species
       real colden(ngridmx,nqmx)   ! vertical column of tracers
       REAL mtot(ngridmx)          ! Total mass of water vapor (kg/m2)
@@ -429,5 +430,8 @@
 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-         if (callthermos) call param_read
+         if (callthermos) then
+            if(solvarmod.eq.0) call param_read
+            if(solvarmod.eq.1) call param_read_e107
+         endif
 #endif
 c        Initialize R and Cp as constant
@@ -561,16 +565,16 @@
                  CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte)
               else if(nltemodel.eq.2) then
-                 do ig=1,ngrid
-                    do l=1,nlayer
-                       co2vmr_gcm(ig,l)=pq(ig,l,igcm_co2)*
-     $                      mmean(ig,l)/mmol(igcm_co2)
-                       n2vmr_gcm(ig,l)=pq(ig,l,igcm_n2)*
-     $                      mmean(ig,l)/mmol(igcm_n2)
-                       covmr_gcm(ig,l)=pq(ig,l,igcm_co)*
-     $                      mmean(ig,l)/mmol(igcm_co)
-                       ovmr_gcm(ig,l)=pq(ig,l,igcm_o)*
-     $                      mmean(ig,l)/mmol(igcm_o)
-                    enddo
-                 enddo
+                co2vmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_co2)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
+                n2vmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_n2)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
+                covmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_co)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
+                ovmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_o)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
                  
                  CALL NLTEdlvr09_TCOOL(ngrid,nlayer,pplay*9.869e-6,
@@ -578,9 +582,6 @@
      $                ovmr_gcm,  zdtnlte )
 
-                 do ig=1,ngrid
-                    do l=1,nlayer
-                       zdtnlte(ig,l)=zdtnlte(ig,l)/86400.
-                    enddo
-                 enddo
+                 zdtnlte(1:ngrid,1:nlayer)=
+     &                             zdtnlte(1:ngrid,1:nlayer)/86400.
               endif
            else
@@ -1549,8 +1550,8 @@
         call wstats(ngrid,"v","Meridional (North-South) wind",
      &                "m.s-1",3,zv)
-c           call wstats(ngrid,"w","Vertical (down-up) wind",
-c     &                "m.s-1",3,pw)
+        call wstats(ngrid,"w","Vertical (down-up) wind",
+     &                "m.s-1",3,pw)
         call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
-c           call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
+        call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
 c          call wstats(ngrid,"q2",
 c    &                "Boundary layer eddy kinetic energy",
@@ -1624,11 +1625,15 @@
      $                 noms(iq) .ne. "ccn_mass" .and.
      $                 noms(iq) .ne. "ccn_number") then
-                   do l=1,nlayer
-                      do ig=1,ngrid
-                         vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
-                      end do
-                   end do
+                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
+     &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
+                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
+     &                          *rho(1:ngrid,1:nlayer)*n_avog/
+     &                           (1000*mmol(iq))
                    call wstats(ngrid,"vmr_"//trim(noms(iq)),
      $                         "Volume mixing ratio","mol/mol",3,vmr)
+!                   call wstats(ngrid,"rho_"//trim(noms(iq)),
+!     $                     "Number density","cm-3",3,rhopart)
+!                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
+!     $                     "Number density","cm-3",3,rhopart)
                    if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
      $                 (noms(iq).eq."o3")) then                     
@@ -1928,4 +1933,22 @@
 c          endif ! (submicron)
          end if  ! (tracer.and.(dustbin.ne.0))
+
+
+c        ----------------------------------------------------------
+c        Thermospheric outputs
+c        ----------------------------------------------------------
+
+         if(callthermos) then 
+
+            call WRITEDIAGFI(ngridmx,"q15um","15 um cooling","K/s",
+     $           3,zdtnlte)
+            call WRITEDIAGFI(ngridmx,"quv","UV heating","K/s",
+     $           3,zdteuv)
+            call WRITEDIAGFI(ngridmx,"cond","Thermal conduction","K/s",
+     $           3,zdtconduc)
+            call WRITEDIAGFI(ngridmx,"qnir","NIR heating","K/s",
+     $           3,zdtnirco2)
+
+         endif  !(callthermos)
 
 c        ----------------------------------------------------------
