Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 657)
+++ /trunk/LMDZ.MARS/README	(revision 658)
@@ -1665,2 +1665,15 @@
 - aeronomars/conc.h : cleanup of obsolete variables
 - aeronomars/chimiedata.h : cleanup of obsolete variables
+
+== 12/05/12 == FGG+JYC+EM
+- updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h,
+  paramfoto_compact.F, param_read.F , thermosphere.F).
+- minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function
+  which seems to be a problem for g95) .
+- minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if
+  none is computed.
+- in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP"
+  and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines
+  defined in "moldiff_red.F". Added use of automatic-sized array in "tridag"
+  and "tridag_sp" local array "gam" (to avoid using an a priori oversized
+  local array).
Index: /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 658)
@@ -95,5 +95,5 @@
       integer,save :: nbq                   ! number of tracers used in the chemistry
       integer,save :: niq(nqmx)             ! array storing the indexes of the tracers
-      integer :: major(nlayermx)            ! index of major species
+      integer :: iloc(1)            ! index of major species
       integer :: ig,l,i,iq,iqmax
       integer :: foundswitch, lswitch
@@ -171,5 +171,4 @@
             call read_phototable
          end if
-
          ! find index of chemical tracers to use
          ! Listed here are all tracers that can go into photochemistry
@@ -660,5 +659,4 @@
                             zu, zv, zt, zycol, ptimestep, co2ice)
          end if
-
 !=======================================================================
 !     tendencies
@@ -667,10 +665,10 @@
 !     index of the most abundant species at each level
 
-         major(:) = maxloc(zycol, dim = 2)
+!         major(:) = maxloc(zycol, dim = 2)
 
 !     tendency for the most abundant species = - sum of others
-
          do l = 1,nlayermx
-            iqmax = major(l)
+            iloc=maxloc(zycol(l,:))
+            iqmax=iloc(1)
             do i = 1,nbq
                iq = niq(i) ! get tracer index
@@ -683,5 +681,4 @@
             end do
          end do ! of do l = 1,nlayermx
-!         if(ig.eq.800)write(*,*)'calchim/686',dqchim(ig,27,23)
 
 !=======================================================================
Index: /trunk/LMDZ.MARS/libf/aeronomars/iono.h
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/iono.h	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/iono.h	(revision 658)
@@ -19,29 +19,29 @@
      $ taunoplus,taun2plus,tauhplus,tauhco2plus
 
-      real*8    tauco2(nlayermx,nreact)
-      real*8    tauo2(nlayermx,nreact)
-      real*8    tauo3p(nlayermx,nreact)
-      real*8    tauco(nlayermx,nreact)
-      real*8    tauh(nlayermx,nreact)
-      real*8    tauoh(nlayermx,nreact)
-      real*8    tauho2(nlayermx,nreact)
-      real*8    tauh2(nlayermx,nreact)
-      real*8    tauh2o(nlayermx,nreact)
-      real*8    tauo1d(nlayermx,nreact)
-      real*8    tauh2o2(nlayermx,nreact)
-      real*8    tauo3(nlayermx,nreact)
-      real*8    taun(nlayermx,nreact)
-      real*8    tauno(nlayermx,nreact)
-      real*8    taun2(nlayermx,nreact)
-      real*8    taun2d(nlayermx,nreact)
-      real*8    tauno2(nlayermx,nreact)
-      real*8    tauco2plus(nlayermx,nreact)
-      real*8    tauoplus(nlayermx,nreact)
-      real*8    tauo2plus(nlayermx, nreact)
-      real*8    taucoplus(nlayermx, nreact)
-      real*8    taucplus(nlayermx, nreact)
-      real*8    taunplus(nlayermx, nreact) 
-      real*8    taunoplus(nlayermx, nreact)
-      real*8    taun2plus(nlayermx, nreact)
-      real*8    tauhplus(nlayermx, nreact)
-      real*8    tauhco2plus(nlayermx, nreact)
+      real*8    tauco2(nreact,nlayermx)
+      real*8    tauo2(nreact,nlayermx)
+      real*8    tauo3p(nreact,nlayermx)
+      real*8    tauco(nreact,nlayermx)
+      real*8    tauh(nreact,nlayermx)
+      real*8    tauoh(nreact,nlayermx)
+      real*8    tauho2(nreact,nlayermx)
+      real*8    tauh2(nreact,nlayermx)
+      real*8    tauh2o(nreact,nlayermx)
+      real*8    tauo1d(nreact,nlayermx)
+      real*8    tauh2o2(nreact,nlayermx)
+      real*8    tauo3(nreact,nlayermx)
+      real*8    taun(nreact,nlayermx)
+      real*8    tauno(nreact,nlayermx)
+      real*8    taun2(nreact,nlayermx)
+      real*8    taun2d(nreact,nlayermx)
+      real*8    tauno2(nreact,nlayermx)
+      real*8    tauco2plus(nreact,nlayermx)
+      real*8    tauoplus(nreact,nlayermx)
+      real*8    tauo2plus(nreact,nlayermx)
+      real*8    taucoplus(nreact,nlayermx)
+      real*8    taucplus(nreact,nlayermx)
+      real*8    taunplus(nreact,nlayermx) 
+      real*8    taunoplus(nreact,nlayermx)
+      real*8    taun2plus(nreact,nlayermx)
+      real*8    tauhplus(nreact,nlayermx)
+      real*8    tauhco2plus(nreact,nlayermx)
Index: /trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F	(revision 658)
@@ -32,6 +32,4 @@
 c    local parameters and variables
 
-      real       aux1(nlayermx)	                !auxiliar variable
-      real       aux2(nlayermx)                 !auxiliar variable
       real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
       real       o2colx(nlayermx)               !column density of O2(cm^-2)
@@ -57,15 +55,34 @@
 
 
-
 c     variables used in interpolation
 
-      real aux3(nz2), aux4(nz2)
-                                              !auxiliar variables
-      real       limdown                      !limits for interpolation
-      real       limup                        !        ""
+      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*******************************
@@ -117,30 +134,20 @@
 ! in each spectral interval
 
-!AQUI SE PODRIAN AGRUPAR CALCULOS PARA AHORRAR TIEMPO DE CPU
-!P.EJ. LOS CALCULOS DE AUX2 y AUX4 NO ES NECESARIO REPETIRLOS
-!PARA CADA ESPECIE EN UN INTERVALO.
-!REVISAR
-
-!TAMBIEN SE PODRIA PONER, EN LUGAR DE CONDICIONES SOBRE CHEMTHERMOD
-!PARA VER QUE ESPECIES SE CONSIDERAN, PONER CONDICION SOBRE LA EXISTENCIA DE 
-!CADA ESPECIE EN TRACEUR.DEF. ASI SE CALCULARIA LA FOTOABSORCION DE TODAS LAS
-!ESPECIES INCLUIDAS, AUNQUE LUEGO NO SE TENGA EN CUENTA EN LA QUIMICA (P.EJ.,
-!SI HAY NO2 PERO NO NO, NO SE CALCULARIA QUIMICA DEL NITROGENO PERO SE PODRIA
-!TENER EN CUENTA PARA EL CALENTAMIENTO
-!ESTUDIAR EN EL FUTURO
-
-!OTRA POSIBILIDAD ES SUSTITUIR LA ESCRITURA EN DURO DE LOS INDICES
-!DE LAS ESPECIES EN JABSIFOTS, CRSCABSI, ETC. POR INDICES DEL TIPO I_*
-!ESTUDIAR EN EL FUTURO
-c************************************************
-c     co2,0.1,6.0
-c************************************************
-
+
+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
-      limdown=1e-20
-      limup=1e26
-      do i=1,nlayermx
-         aux1(i)=0.
-         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
+      do i=1,nlayermx
+         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
      $        o2colx(i)*crscabsi2(2,indexint) + 
      $        o3pcolx(i)*crscabsi2(3,indexint) + 
@@ -148,132 +155,75 @@
      $        ncolx(i)*crscabsi2(9,indexint)
       end do
+      limdown=1.e-20
+      limup=1.e26
+
+
+c     Interpolations
+
       do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
+         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
-
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
+      !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) = wm(i)*auxjco2(ind+1) +
+     $        wp(i)*auxjco2(ind)
+         !O2 interpolated coefficient
+         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
+     $        wp(i)*auxjo2(ind)
+         !O3p interpolated coefficient
+          jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
+     $        wp(i)*auxjo3p(ind)
+          !H2 interpolated coefficient
+          jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
+     $         wp(i)*auxjh2(ind)
       enddo
-
-c************************************************
-c     o2,0.1,6.0
-c************************************************
-
-      indexint=1
-      limdown=1e-20
-      limup=1e26
-      do i=1,nlayermx
-         aux1(i)=0.
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-      enddo
-
-c**************************************************         
-c     o3p,0.1,6.0
-c**************************************************
-
-      indexint=1
-      limdown=1e-20
-      limup=1e26
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,3,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,3,i) = aux1(nlayermx-i+1)
-      enddo
-
-      
-c**************************************************
-c     h2,0.1,6.0
-c**************************************************
-
-      indexint=1
-      limdown=1e-20
-      limup=1e26
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,5,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,5,i) = aux1(nlayermx-i+1)
-      enddo
-
-      !Only if Nitrogen or ionospheric chemistry requested
+      !Only if chemthermod.ge.2
+      !N interpolated coefficient
       if(chemthermod.ge.2) then
-
-c**************************************************
-c     n,0.1,6.0
-c**************************************************
-
-
-         indexint=1
-         limdown=1e-20
-         limup=1e26
          do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,9,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
+            ind=auxind(i)
+            jfotsout(indexint,9,nlayermx-i+1) =  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
-            jfotsout(indexint,9,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif  !Of chemthermod >= 2
-      
-
-c**************************************************
-c     o2, 5-80.5nm
-c**************************************************
-
-      limdown=1e-20
-      limup=1e26
-      do indexint=2,15
-         do i=nlayermx-1,1,-1
-         end do
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
+            auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
      $           o2colx(i)*crscabsi2(2,indexint)+
      $           o3pcolx(i)*crscabsi2(3,indexint)+
@@ -286,286 +236,99 @@
      $           no2colx(i)*crscabsi2(13,indexint)
          end do
+
+c     Interpolations
+
          do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-          call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
+            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
-            jfotsout(indexint,2,i) = aux1(nlayermx-i+1) 
+             ind=auxind(i)
+             auxi = nlayermx-i+1
+             !O2 interpolated coefficient
+             jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
+     $            wp(i)*auxjo2(ind)
+             !O3p interpolated coefficient
+             jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
+     $            wp(i)*auxjo3p(ind)
+             !CO2 interpolated coefficient
+             jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
+     $            wp(i)*auxjco2(ind)
+             !H2 interpolated coefficient
+             jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
+     $            wp(i)*auxjh2(ind)
+             !N2 interpolated coefficient
+             jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
+     $            wp(i)*auxjn2(ind)             
+             !CO interpolated coefficient
+             jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
+     $            wp(i)*auxjco(ind)
+             !H interpolated coefficient
+             jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
+     $            wp(i)*auxjh(ind)             
           enddo
-      end do
-      
-
-c**************************************************
-c     o3p, 5-80.5nm
-c**************************************************
-
-      do indexint=2,15
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,3,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,3,i) = aux1(nlayermx-i+1)
-         enddo
-      end do
-
-c**************************************************
-c     co2, 5-80.5nm
-c**************************************************
-
-      do indexint=2,15
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
-         enddo
-      end do
-
-
-c**************************************************
-c     h2, 5-80.5nm
-c**************************************************
-
-      do indexint=2,15
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,5,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,5,i) = aux1(nlayermx-i+1)
-         enddo
-      end do
-
-
-c**************************************************
-c     n2, 5-80.5nm
-c**************************************************
-
-      do indexint=2,15
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,8,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
-         enddo
-      end do
-
-
-c**************************************************
-c     co, 5-80.5nm
-c**************************************************
-
-      do indexint=2,15
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
-         enddo
-      end do
-
-
-c**************************************************
-c     h, 5-80.5nm
-c**************************************************
-
-      do indexint=2,15
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,12,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,12,i) = aux1(nlayermx-i+1)
-         enddo
-      end do
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-c**************************************************
-c     n, 5-80.5nm
-c**************************************************
-
-         do indexint=2,15
-            do i=1,nlayermx
-               aux2(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
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,9,nz2-i+1)
-               aux4(i) = c1_16(nz2-i+1,indexint)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               jfotsout(indexint,9,i) = aux1(nlayermx-i+1)
-            enddo
-         end do
-
-
-c**************************************************
-c     no, 5-80.5nm
-c**************************************************
-
-         do indexint=2,15
-            do i=1,nlayermx
-               aux2(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
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
-               aux4(i) = c1_16(nz2-i+1,indexint)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
-            enddo
-         end do
-
-c**************************************************
-c     no2, 5-80.5nm
-c**************************************************
-
-         do indexint=2,15
-            do i=1,nlayermx
-               aux2(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
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-               aux4(i) = c1_16(nz2-i+1,indexint)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-            enddo
-         end do
-      
-         endif  !Of chemthermod >= 2
-
-         
-c**************************************************
-c     o2, 80.6-90.8nm
-c**************************************************
-
+          !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) = wm(i)*auxjn(ind+1) +
+     $               wp(i)*auxjn(ind)
+                !NO interpolated coefficient
+                jfotsout(indexint,10,auxi)=wm(i)*auxjno(ind+1) +
+     $               wp(i)*auxjno(ind)
+                !NO2 interpolated coefficient
+                jfotsout(indexint,13,auxi)=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
-         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
+         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
      $        o2colx(i)*crscabsi2(2,indexint)+
      $        o3pcolx(i)*crscabsi2(3,indexint)+
@@ -577,1192 +340,673 @@
      $        no2colx(i)*crscabsi2(13,indexint)
       end do
+
+c     Interpolations
+
       do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
+         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
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,2,i) = aux1(nlayermx-i+1) 
+      !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) = wm(i)*auxjo2(ind+1) +
+     $            wp(i)*auxjo2(ind)
+         !CO2 interpolated coefficient
+         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
+     $            wp(i)*auxjco2(ind) 
+         !O3p interpolated coefficient
+         jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
+     $            wp(i)*auxjo3p(ind) 
+         !N2 interpolated coefficient
+         jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
+     $            wp(i)*auxjn2(ind) 
+         !CO interpolated coefficient
+         jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
+     $            wp(i)*auxjco(ind)  
+         !H interpolated coefficient
+         jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
+     $            wp(i)*auxjh(ind)         
       enddo
-
-
-c**************************************************
-c     co2, 80.6-90.8nm
-c**************************************************
-
-      indexint=16
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,1,i) = aux1(nlayermx-i+1) 
-      enddo
-
-
-c**************************************************
-c     o3p, 80.6-90.8nm
-c**************************************************
-
-      indexint=16
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,3,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,3,i) = aux1(nlayermx-i+1) 
-      enddo
-      
-
-c**************************************************
-c     n2, 80.6-90.8nm
-c**************************************************
-
-      indexint=16
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,8,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,8,i) = aux1(nlayermx-i+1) 
-      enddo
-
-
-c**************************************************
-c     co, 80.6-90.8nm
-c************************************************************
-
-      indexint=16
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,11,i) = aux1(nlayermx-i+1) 
-      enddo
-
-
-c**************************************************
-c     h, 80.6-90.8nm
-c**************************************************
-
-      indexint=16
-      do i=1,nlayermx
-         aux2(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
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,12,nz2-i+1)
-         aux4(i) = c1_16(nz2-i+1,indexint)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,12,i) = aux1(nlayermx-i+1) 
-      enddo
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
+      !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-
-c**************************************************
-c     n, 80.6-90.8nm
-c**************************************************
-
-         indexint=16
          do i=1,nlayermx
-            aux2(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
+            ind=auxind(i)
+            auxi = nlayermx-i+1
+            !N interpolated coefficient
+            jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
+     $           wp(i)*auxjn(ind) 
+            !NO interpolated coefficient
+            jfotsout(indexint,10,auxi) = wm(i)*auxjno(ind+1) +
+     $           wp(i)*auxjno(ind)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,auxi) = 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
-            aux3(i) = jabsifotsintpar(indexint,9,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,9,i) = aux1(nlayermx-i+1) 
-         enddo
-
-
-c**************************************************
-c     no, 80.6-90.8nm
-c**************************************************
-
-         indexint=16
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,10,i) = aux1(nlayermx-i+1) 
-         enddo
-
-
-c***********************************************************
-c     no2, 80.6-90.8nm
-c**************************************************
-
-         indexint=16
-         do i=1,nlayermx
-            aux2(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
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-            aux4(i) = c1_16(nz2-i+1,indexint)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,13,i) = aux1(nlayermx-i+1) 
-         enddo
-      
-         endif !Of chemthermod >= 2
-
-c**************************************************
-c      co2, 90.9-119.5nm
-c**************************************************
-
-      limdown=1e-20
-      limup=1e26
-      do indexint=17,24
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
-     $           nocolx(i) + cocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-            aux4(i) = c17_24(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
-            if(indexint.eq.24) then
-          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-               jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))*
-     $           (1+alfa(indexint,i)*(t2(i)-t0(i)))
+            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 
-                  jfotsout(indexint,1,i)=0.
-               end if
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c     o2, 90.9-119.5nm
-c**************************************************
-      
-      do indexint=17,24
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
-     $           nocolx(i) + cocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-            aux4(i) = c17_24(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-            if(indexint.eq.24) then
-          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-               jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,2,i)=0.
-               end if
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c     n2, 90.9-119.5nm
-c**************************************************
-     
-      do indexint=17,24
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
-     $           nocolx(i) + cocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,8,nz2-i+1)
-            aux4(i) = c17_24(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
-            if(indexint.eq.24) then
-          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-               jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,8,i)=0.
-               end if
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c     co, 90.9-119.5nm
-c**************************************************
-     
-      do indexint=17,24
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
-     $           nocolx(i) + cocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
-            aux4(i) = c17_24(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
-            if(indexint.eq.24) then
-          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-               jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,11,i)=0.
-               end if
-            end if
-         enddo
-      end do
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-
-c**************************************************
-c     no, 90.9-119.5nm
-c**************************************************
-     
-         do indexint=17,24
-            do i=1,nlayermx
-               aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
-     $              nocolx(i) + cocolx(i) + no2colx(i)
-            end do
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
-               aux4(i) = c17_24(nz2-i+1)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
-               if(indexint.eq.24) then
-                  if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $                 .lt.60.) then
-                     jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
-     $             *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-                  else
-                     jfotsout(indexint,10,i)=0.
-                  end if
+                  cortemp(i)=0.
                end if
             enddo
-         end do
-
-
-c**************************************************
-c     no2, 90.9-119.5nm
-c**************************************************
-     
-         do indexint=17,24
+         else
             do i=1,nlayermx
-               aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
-     $              nocolx(i) + cocolx(i) + no2colx(i)
-            end do
+               cortemp(i)=1.
+            enddo
+         end if
+         do i=1,nlayermx           
+            ind=auxind(i)
+            auxi = nlayermx-i+1
+            !O2 interpolated coefficient
+            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
+     $           wp(i)*auxjo2(ind)) * cortemp(i)
+            !CO2 interpolated coefficient
+            jfotsout(indexint,1,auxi) = (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) = (wm(i)*auxjn2(ind+1) +
+     $            wp(i)*auxjn2(ind)) * cortemp(i)            
+            !CO interpolated coefficient
+            jfotsout(indexint,11,auxi) = (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)=(wm(i)*auxjno(ind+1) +
+     $              wp(i)*auxjno(ind)) * cortemp(i)
+               !NO2 interpolated coefficient
+               jfotsout(indexint,13,auxi)=(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
-               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-               aux4(i) = c17_24(nz2-i+1)
+               auxi = nz2-i+1
+               !NO tabulated coefficient
+               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
+               !NO2 tabulated coefficient
+               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
             enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
+         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) = (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) = (wm(i)*auxjo2(ind+1) +
+     $           wp(i)*auxjo2(ind)) * cortemp(i)
+            !H2O interpolated coefficient
+            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
+     $           wp(i)*auxjh2o(ind)) * cortemp(i)
+            !H2O2 interpolated coefficient
+            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
+     $           wp(i)*auxjh2o2(ind)) * cortemp(i)            
+            !CO interpolated coefficient
+            jfotsout(indexint,11,auxi) = (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
-               jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-               if(indexint.eq.24) then
-                  if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $                 .lt.60.) then
-                     jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-                  else
-                     jfotsout(indexint,13,i)=0.
-                  end if
-               end if
+               ind=auxind(i)
+               auxi = nlayermx-i+1
+               !NO interpolated coefficient
+               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
+     $              wp(i)*auxjno(ind)) * cortemp(i)
+               !NO2 interpolated coefficient
+               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
+     $              wp(i)*auxjno2(ind)) * cortemp(i)
             enddo
-         end do
-
-         endif !Of chemthermod >= 2
-
-
-c**************************************************
-c      co2, 119.6-167.0nm
-c**************************************************
-
-      do indexint=25,29
+         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
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
-     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-         end do
+            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) = (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) = (wm(i)*auxjo2(ind+1) +
+     $            wp(i)*auxjo2(ind)) * cortemp(i)
+            !H2O interpolated coefficient
+            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
+     $            wp(i)*auxjh2o(ind)) * cortemp(i)
+            !H2O2 interpolated coefficient
+            jfotsout(indexint,6,auxi) = (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)=(wm(i)*auxjno(ind+1) +
+     $              wp(i)*auxjno(ind)) * cortemp(i)
+               !NO2 interpolated coefficient
+               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
-            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-            aux4(i) = c25_29(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
+            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) = (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) = (wm(i)*auxjo2(ind+1) +
+     $        wp(i)*auxjo2(ind)) * cortemp(i)
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = (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
-          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-               jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-     $           *(1+alfa(indexint,i)*(t2(i)-t0(i)))               
-            else
-               jfotsout(indexint,1,i) = 0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      o2, 119.6-167.0nm
-c**************************************************
-
-      do indexint=25,29
+            auxi = nlayermx-i+1
+            ind=auxind(i)
+            !NO interpolated coefficient
+            jfotsout(indexint,10,auxi) = (wm(i)*auxjno(ind+1) +
+     $           wp(i)*auxjno(ind)) * cortemp(i)
+           !NO2 interpolated coefficient
+            jfotsout(indexint,13,auxi) = (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) = wm(i)*auxjo2(ind+1) +
+     $        wp(i)*auxjo2(ind)
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
+     $        wp(i)*auxjh2o2(ind)         
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
          do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-         end do
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 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
-            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-            aux4(i) = c25_29(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
+            !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) = wm(i)*auxjo2(ind+1) +
+     $        wp(i)*auxjo2(ind)
+         !H2O2 interpolated coefficient
+         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
+     $        wp(i)*auxjh2o2(ind)
+         !O3 interpolated coefficient
+         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
+     $        wp(i)*auxjo3(ind)         
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
          do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,2,i)=0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      h2o, 119.6-167.0nm
-c**************************************************
-
-      do indexint=25,29
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 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) = wm(i)*auxjh2o2(ind+1) +
+     $        wp(i)*auxjh2o2(ind)
+         !O3 interpolated coefficient
+         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
+     $        wp(i)*auxjo3(ind)         
+      enddo
+      if(chemthermod.ge.2) then
          do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-         end do
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = 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
-            aux3(i) = jabsifotsintpar(indexint,4,nz2-i+1)
-            aux4(i) = c25_29(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
+            !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) = wm(i)*auxjo3(ind+1) +
+     $        wp(i)*auxjo3(ind)         
+      enddo
+      !Only if chemthermod.ge.2
+      if(chemthermod.ge.2) then
          do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,4,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,4,i) = 0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      h2o2,119.6-167.0nm
-c**************************************************
-
-      do indexint=25,29
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
-            aux4(i) = c25_29(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,6,i) = 0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      co, 119.6-167.0nm
-c**************************************************
-
-      do indexint=25,29
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
-            aux4(i) = c25_29(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,11,i) = 0.
-            end if
-         enddo
-      end do
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-
-c**************************************************
-c      no, 119.6-167.0nm
-c**************************************************
-
-         do indexint=25,29
-            do i=1,nlayermx
-               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $              h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-            end do
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
-               aux4(i) = c25_29(nz2-i+1)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $              .lt.60.) then
-                  jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
-     $             *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,10,i) = 0.
-               end if
-            enddo
-         end do
-
-
-c**************************************************
-c      no2, 119.6-167.0nm
-c**************************************************
-
-         do indexint=25,29
-            do i=1,nlayermx
-               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $              h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-            end do
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-               aux4(i) = c25_29(nz2-i+1)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $              .lt.60.) then
-                  jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-     $             *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,13,i) = 0.
-               end if
-            enddo
-         end do
-
-         endif !Of chemthermod >= 2
-
-
-c**************************************************
-c      co2, 167.0-202.5nm
-c**************************************************
-
-      do indexint=30,31
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
-     $           h2o2colx(i) + nocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-            aux4(i) = c30_31(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-               jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-     $           *(1+alfa(indexint,i)*(t2(i)-t0(i)))               
-            else
-               jfotsout(indexint,1,i) = 0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      o2, 167.0-202.5nm
-c**************************************************
-
-      do indexint=30,31
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-            aux4(i) = c30_31(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,2,i)=0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      h2o, 167.0-202.5nm
-c**************************************************
-
-      do indexint=30,31
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,4,nz2-i+1)
-            aux4(i) = c30_31(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,4,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,4,i) = 0.
-            end if
-         enddo
-      end do
-
-
-c**************************************************
-c      h2o2, 167.0-202.5nm
-c**************************************************
-
-      do indexint=30,31
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $           h2o2colx(i) + nocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
-            aux4(i) = c30_31(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
-     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,6,i) = 0.
-            end if
-         enddo
-      end do
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-      
-c**************************************************
-c      no, 167.0-202.5nm
-c**************************************************
-
-         do indexint=30,31
-            do i=1,nlayermx
-               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $              h2o2colx(i) + nocolx(i) + no2colx(i)
-            end do
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
-               aux4(i) = c30_31(nz2-i+1)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $              .lt.60.) then
-                  jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
-     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,10,i) = 0.
-               end if
-            enddo
-         end do
-
-
-c**************************************************
-c      no2, 167.0-202.5nm
-c**************************************************
-
-         do indexint=30,31
-            do i=1,nlayermx
-               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
-     $              h2o2colx(i) + nocolx(i) + no2colx(i)
-            end do
-            do i=1,nz2
-               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-               aux4(i) = c30_31(nz2-i+1)
-            enddo
-            call interpfast 
-     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-            do i=1,nlayermx
-               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $              .lt.60.) then
-                  jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-               else
-                  jfotsout(indexint,13,i) = 0.
-               end if
-            enddo
-         end do
-         
-      endif !Of chemthermod >= 2
-
-c**************************************************
-c     co2, 202.6-210.0nm
-c**************************************************
-
-      indexint=32
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) + 
-     $        nocolx(i) + no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
-         aux4(i) = c32(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
-     $           *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-     $           *(1+alfa(indexint,i)*(t2(i)-t0(i)))
-         else
-            jfotsout(indexint,1,i)=0.
-         end if
-      enddo
-
-
-c**************************************************
-c     o2, 202.6-210.0nm
-c**************************************************
-
-      indexint=32
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) + 
-     $        nocolx(i) + no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-         aux4(i) = c32(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-     $           *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-         else
-            jfotsout(indexint,2,i)=0.
-         end if
-      enddo
-
-
-c**************************************************
-c     h2o2, 202.6-210.0nm
-c**************************************************
-
-      indexint=32
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) + 
-     $        nocolx(i) + no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
-         aux4(i) = c32(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
-            jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
-     $           *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-         else
-            jfotsout(indexint,6,i)=0.
-         end if
-      enddo
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.eq.2) then
-
-c**************************************************
-c     no, 202.6-210.0nm
-c**************************************************
-
-         indexint=32
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) + 
-     $           nocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
-            aux4(i) = c32(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $           .lt.60.) then
-               jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
-     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,10,i)=0.
-            end if
-         enddo
-
-
-c**************************************************
-c     no2, 202.6-210.0nm
-c**************************************************
-
-         indexint=32
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) + 
-     $           nocolx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-            aux4(i) = c32(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
-     $           .lt.60.) then
-               jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
-            else
-               jfotsout(indexint,13,i)=0.
-            end if
-         enddo
-         
-      endif !Of chemthermod >= 2
-
-
-c**************************************************
-c     o2, 210.0-231.0nm
-c**************************************************
-
-      indexint=33
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-         aux4(i) = c33(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-      enddo
-
-
-c**************************************************
-c     h2o2, 210.1-231.0nm
-c**************************************************
-
-      indexint=33
-      limdown=1.e-20
-      limup=1e25
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
-         aux4(i) = c33(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
-      enddo
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-
-c**************************************************
-c     no2, 210.1-231.0nm
-c**************************************************
-
-         indexint=33
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-            aux4(i) = c33(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif !Of chemthermod >= 2
-
-
-c**************************************************
-c     o2, 231.-240.nm
-c**************************************************
-
-      indexint=34
-      limdown=1.e-20
-      limup=1e25
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
-     $        no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
-         aux4(i) = c34(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
-      enddo
-
-
-c**************************************************
-c     h2o2, 231.0-240.nm
-c**************************************************
-
-      indexint=34
-      limdown=1.e-20
-      limup=1e25
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
-     $        no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
-         aux4(i) = c34(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
-      enddo
-
-
-      !Only if Ozone, Nitrogen or ionospheric chem. requested
-      if(chemthermod.ge.1) then
-
-c**************************************************
-c     o3, 231.0-240.nm
-c**************************************************
-
-         indexint=34
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
-     $           no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,7,nz2-i+1)
-            aux4(i) = c34(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,7,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif !Of chemthermod.ge.1
-
-
-      !Only if nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-
-c**************************************************
-c     no2, 231.0-240.nm
-c**************************************************
-
-         indexint=34
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
-     $           no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-            aux4(i) = c34(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif !Of chemthermod >= 2
-
-
-c**************************************************
-c     h2o2, 240.0-337.7nm
-c**************************************************
-
-      indexint=35
-      limdown=1.e-20
-      limup=1e25
-      do i=1,nlayermx
-         aux2(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
-      end do
-      do i=1,nz2
-         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
-         aux4(i) = c35(nz2-i+1)
-      enddo
-      call interpfast 
-     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-      do i=1,nlayermx
-         jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
-      enddo
-      
-
-      !Only if Ozone, Nitrogen or ionospheric chem. requested
-      if(chemthermod.ge.1) then
-
-c**************************************************
-c     o3, 240.0-337.7nm
-c**************************************************
-
-         indexint=35
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,7,nz2-i+1)
-            aux4(i) = c35(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,7,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif !Of chemthermod.ge.1
-
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-
-c**************************************************
-c     no2, 240.0-337.7nm
-c**************************************************
-
-         indexint=35
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-            aux4(i) = c35(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif !Of chemthermod.ge.2
-
-
-      !Only if Ozone, Nitrogen or ionospheric chem. requested
-      if(chemthermod.ge.1) then
-
-c**************************************************
-c     o3, 337.7-800.0nm
-c**************************************************
-
-         indexint=36
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = o3colx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,7,nz2-i+1)
-            aux4(i) = c36(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,7,i) = aux1(nlayermx-i+1)
-         enddo
-
+            ind=auxind(i)
+            !NO2 interpolated coefficient
+            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
+     $           wp(i)*auxjno2(ind)
+         enddo
       endif
 
-
-      !Only if Nitrogen or ionospheric chemistry requested
-      if(chemthermod.ge.2) then
-
-c**************************************************
-c     no2, 337.7-800.0nm
-c**************************************************
-
-         indexint=36
-         limdown=1.e-20
-         limup=1e25
-         do i=1,nlayermx
-            aux2(nlayermx-i+1) = o3colx(i) + no2colx(i)
-         end do
-         do i=1,nz2
-            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
-            aux4(i) = c36(nz2-i+1)
-         enddo
-         call interpfast 
-     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
-         do i=1,nlayermx
-            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
-         enddo
-
-      endif !Of chemthermod.ge.2
+c     End of interval 36
+
+c     End of interpolation to obtain photoabsorption rates
+
+
       return
 
@@ -1885,5 +1129,4 @@
 
       nz3 = nlayermx*2
-      szadeg = zenit*180./3.141592
       do i=1,nlayermx
          xx = ( radio + iz(i) ) * 1.e5
@@ -1911,4 +1154,5 @@
          Hno2  = xx / mmol(igcm_no2)
       endif
+      ! first loop in altitude : initialisation
       do i=nlayermx,1,-1
          !Column initialisation
@@ -1945,4 +1189,7 @@
             no2x(i)  = rm(i,i_no2)
          endif
+      enddo
+      ! second loop in altitude : column calculations
+      do i=nlayermx,1,-1
          !Routine to calculate the geometrical length of each layer
          call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
@@ -1971,5 +1218,5 @@
                !Top layer
                if(jj.eq.nlayermx) then
-                  if(szadeg.le.60.) then 
+                  if(zenit.le.60.) then 
                      o3pcolx(i)=o3pcolx(i)+o3px(nlayermx)*Ho3p*esp(j)
      $                    *1.e-5
@@ -2003,5 +1250,5 @@
      $                       *1.e-5
                      endif
-                  else if(szadeg.gt.60.) then
+                  else if(zenit.gt.60.) then
                      espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
                      espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
@@ -2040,5 +1287,5 @@
                         no2colx(i) = no2colx(i) + espno2*no2x(nlayermx)
                      endif
-                  endif !Of if szadeg.lt.60
+                  endif !Of if zenit.lt.60
                !Other layers
                else 
@@ -2086,5 +1333,6 @@
 c**********************************************************************
 c**********************************************************************
-      subroutine interpfast(escout,p,nlayer,escin,pin,nl,limdown,limup)
+
+      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
 C
 C subroutine to perform linear interpolation in pressure from 1D profile 
@@ -2092,28 +1340,29 @@
 C escout(nlayer) on pressure grid p(nlayer).
 C
-      real escout(nlayer),p(nlayer)
-      real escin(nl),pin(nl),wm,wp
-      real limup,limdown
-      integer nl,nlayer,n1,n,nm,np
-      nm=1
-      do 5 n1=1,nlayer
+      real*8 wm(nlayer),wp(nlayer),p(nlayer)
+      integer nm(nlayer)
+      real*8 pin(nl)
+      real*8 limup,limdown
+      integer nl,nlayer,n1,n,np,nini
+      nini=1
+      do n1=1,nlayer
          if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
-            escout(n1) = 0.0
+            wm(n1) = 0.d0
+            wp(n1) = 0.d0
          else
-            do n = nm,nl-1
+            do n = nini,nl-1
                if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
-                  nm=n
+                  nm(n1)=n
                   np=n+1
-                  wm=abs(pin(nm)-p(n1))/(pin(np)-pin(nm))
-                  wp=1.0 - wm
-                  goto 33
+                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
+                  wp(n1)=1.d0 - wm(n1)
+                  nini = n
+                  exit
                endif
             enddo
- 33         escout(n1) = escin(np)*wm + escin(nm)*wp
          endif
-    5 continue
+      enddo
       return
       end
-
 
 
Index: /trunk/LMDZ.MARS/libf/aeronomars/moldiff.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/moldiff.F	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/moldiff.F	(revision 658)
@@ -279,9 +279,9 @@
 c Inverting the alfa matrix
 c
-        call ludcmp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
+        call ludcmp_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
 
 c       TEMPORAIRE *****************************
         if (ierr.ne.0) then
-            write(*,*) 'In moldiff: Problem in LUDCMP with matrix alf'
+            write(*,*)'In moldiff: Problem in LUDCMP_SP with matrix alf'
             write(*,*) 'Singular matrix ?'
 c           write(*,*) 'Matrix alf = ', alf
@@ -294,5 +294,5 @@
 c       *******************************************
         do n=1,ncompmoldiff-1
-          call lubksb(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
+          call lubksb_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
           do nn=1,ncompmoldiff-1
             alfinv(l,nn,n)=y(nn,n)/hh
@@ -415,5 +415,5 @@
         btri(1)=btri(1)+atri(1)
 
-        call tridag(atri,btri,ctri,rtri,qtri,nz-2)
+        call tridag_sp(atri,btri,ctri,rtri,qtri,nz-2)
 
         do l=2,nz-1
@@ -464,10 +464,10 @@
 c    ********************************************************************
  
-      subroutine tridag(a,b,c,r,u,n)
-      parameter (nmax=100)
+      subroutine tridag_sp(a,b,c,r,u,n)
+c      parameter (nmax=100)
 c      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
-      real gam(nmax),a(n),b(n),c(n),r(n),u(n)
+      real gam(n),a(n),b(n),c(n),r(n),u(n)
       if(b(1).eq.0.)then
-        stop 'tridag: error: b(1)=0 !!! '
+        stop 'tridag_sp: error: b(1)=0 !!! '
       endif
       bet=b(1)
@@ -477,5 +477,5 @@
         bet=b(j)-a(j)*gam(j)
         if(bet.eq.0.) then
-          stop 'tridag: error: bet=0 !!! '
+          stop 'tridag_sp: error: bet=0 !!! '
         endif
         u(j)=(r(j)-a(j)*u(j-1))/bet
@@ -491,5 +491,5 @@
 c    ********************************************************************
 
-      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
+      SUBROUTINE LUBKSB_SP(A,N,NP,INDX,B)
 
       implicit none
@@ -530,5 +530,5 @@
 c    ********************************************************************
 
-      SUBROUTINE LUDCMP(A,N,NP,INDX,D,ierr)
+      SUBROUTINE LUDCMP_SP(A,N,NP,INDX,D,ierr)
 
       implicit none
@@ -550,5 +550,5 @@
 11      CONTINUE
         IF (AAMAX.EQ.0.) then
-            write(*,*) 'In moldiff: Problem in LUDCMP with matrix A'
+            write(*,*) 'In moldiff: Problem in LUDCMP_SP with matrix A'
             write(*,*) 'Singular matrix ?'
 c           write(*,*) 'Matrix A = ', A
Index: /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 658)
@@ -726,7 +726,7 @@
 
       subroutine tridag(a,b,c,r,u,n)
-      parameter (nmax=4000)
+!      parameter (nmax=4000)
 !      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
-      real*8 gam(nmax),a(n),b(n),c(n),r(n),u(n)
+      real*8 gam(n),a(n),b(n),c(n),r(n),u(n)
       if(b(1).eq.0.)then
         stop 'tridag: error: b(1)=0 !!! '
Index: /trunk/LMDZ.MARS/libf/aeronomars/param_read.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/param_read.F	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/param_read.F	(revision 658)
@@ -95,5 +95,5 @@
       IF (ierr.NE.0) THEN 
        write(*,*)'cant find directory EUVDAT and content coln.dat'
-       write(*,*)'(in aeronomars/param_read_iono.F)'
+       write(*,*)'(in aeronomars/param_read.F)'
        write(*,*)'It should be in :', datafile,'/'
        write(*,*)'1) You can change this directory address in '
@@ -181,106 +181,106 @@
 
       do i=nz2,1,-1
-         read(220,*) (jabsifotsintpar(j,2,i),j=1,16)
-      end do
-      
-      do i=nz2,1,-1
-         read(230,*) (jabsifotsintpar(j,3,i),j=1,16)
-      end do
-
-      do i=nz2,1,-1
-         read(240,*) (jabsifotsintpar(j,1,i),j=1,16)
-      end do
-
-      do i=nz2,1,-1
-         read(250,*) (jabsifotsintpar(j,2,i),j=17,24)
-      end do
-
-
-      do i=nz2,1,-1
-         read(260,*) (jabsifotsintpar(j,2,i),j=25,31)
-      end do
-
-      do i=nz2,1,-1
-         read(270,*) (jabsifotsintpar(j,1,i),j=17,24)
-      end do
-
-      do i=nz2,1,-1
-         read(280,*) (jabsifotsintpar(j,1,i),j=25,31)
-      end do
-
-      do i=nz2,1,-1
-         read(290,*) jabsifotsintpar(32,1,i)
-      end do
-
-      do i=nz2,1,-1
-         read(300,*) (jabsifotsintpar(j,2,i),j=32,34)
-      end do
-
-      do i=nz2,1,-1
-         read(160,*) (jabsifotsintpar(j,5,i),j=1,15)
-      end do
-
-      do i=nz2,1,-1
-         read(150,*) (jabsifotsintpar(j,4,i),j=25,31)
-      end do
-
-      do i=nz2,1,-1
-         read(170,*) (jabsifotsintpar(j,6,i),j=25,35)
-      end do
-
-      do i=nz2,1,-1
-         read(180,*) (jabsifotsintpar(j,7,i),j=34,36)
-      end do
-
-      do i=nz2,1,-1
-         read(390,*) (jabsifotsintpar(j,8,i),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(400,*) (jabsifotsintpar(j,8,i),j=17,24)
-      enddo
-
-      do i=nz2,1,-1
-         read(410,*) (jabsifotsintpar(j,9,i),j=1,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(420,*) (jabsifotsintpar(j,10,i),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(430,*) (jabsifotsintpar(j,10,i),j=17,24)
-      enddo
-
-      do i=nz2,1,-1
-         read(440,*) (jabsifotsintpar(j,10,i),j=25,32)
-      enddo
-
-      do i=nz2,1,-1
-         read(450,*) (jabsifotsintpar(j,11,i),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(460,*) (jabsifotsintpar(j,11,i),j=17,24)
-      enddo
-
-      do i=nz2,1,-1
-         read(470,*) (jabsifotsintpar(j,11,i),j=25,29)
-      enddo
-      
-      do i=nz2,1,-1
-         read(480,*) (jabsifotsintpar(j,12,i),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(490,*) (jabsifotsintpar(j,13,i),j=2,16)
-      enddo
-      
-      do i=nz2,1,-1
-         read(500,*) (jabsifotsintpar(j,13,i),j=17,24)
-      enddo
-      
-      do i=nz2,1,-1
-         read(510,*) (jabsifotsintpar(j,13,i),j=25,36)
+         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
 
Index: /trunk/LMDZ.MARS/libf/aeronomars/param_v4.h
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/param_v4.h	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/param_v4.h	(revision 658)
@@ -29,5 +29,5 @@
       real fluxtop(ninter)
       real freccen(ninter)          !representative wavelenght
-      real jabsifotsintpar(ninter,nabs,nz2)
+      real jabsifotsintpar(nz2,nabs,ninter)
 
 
Index: /trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F	(revision 658)
@@ -1,3 +1,2 @@
-
 c**********************************************************************
 
@@ -1449,5 +1448,5 @@
 	   enddo
 	   incremento = ( tehanson(i2)-tehanson(i1) ) / 
-	1	( zhanson(i2)-zhanson(i1) )
+     $          ( zhanson(i2)-zhanson(i1) )
 	   t_elect = tehanson(i1) + (zkm-zhanson(i1)) * incremento
 	endif
@@ -1814,141 +1813,141 @@
          !Initialization of lifetimes
          do j = 1, nreact
-            tauco2(i,j)  = 1.d30
-            tauo2(i,j)   = 1.d30
-            tauo3p(i,j)  = 1.d30
-            tauco(i,j)   = 1.d30
-            tauh(i,j)    = 1.d30
-            tauoh(i,j)   = 1.d30
-            tauho2(i,j)  = 1.d30
-            tauh2(i,j)   = 1.d30
-            tauh2o(i,j)  = 1.d30
-            tauo1d(i,j)  = 1.d30
-            tauh2o2(i,j) = 1.d30
-            tauo3(i,j)   = 1.d30
-            taun2(i,j)   = 1.d30
-            taun(i,j)    = 1.d30
-            tauno(i,j)   = 1.d30
-            taun2d(i,j)  = 1.d30
-            tauno2(i,j)  = 1.d30
-            tauco2plus(i,j) = 1.d30
-            tauoplus(i,j)   = 1.d30
-	    tauo2plus(i,j)  = 1.d30
-            taucoplus(i,j)  = 1.d30
-	    taucplus(i,j)   = 1.d30
-	    taunplus(i,j)   = 1.d30
-	    taun2plus(i,j)  = 1.d30
-            taunoplus(i,j)  = 1.d30
-	    tauhplus(i,j)   = 1.d30
-            tauhco2plus(i,j)= 1.d30
+            tauco2(j,i)  = 1.d30
+            tauo2(j,i)   = 1.d30
+            tauo3p(j,i)  = 1.d30
+            tauco(j,i)   = 1.d30
+            tauh(j,i)    = 1.d30
+            tauoh(j,i)   = 1.d30
+            tauho2(j,i)  = 1.d30
+            tauh2(j,i)   = 1.d30
+            tauh2o(j,i)  = 1.d30
+            tauo1d(j,i)  = 1.d30
+            tauh2o2(j,i) = 1.d30
+            tauo3(j,i)   = 1.d30
+            taun2(j,i)   = 1.d30
+            taun(j,i)    = 1.d30
+            tauno(j,i)   = 1.d30
+            taun2d(j,i)  = 1.d30
+            tauno2(j,i)  = 1.d30
+            tauco2plus(j,i) = 1.d30
+            tauoplus(j,i)   = 1.d30
+	    tauo2plus(j,i)  = 1.d30
+            taucoplus(j,i)  = 1.d30
+	    taucplus(j,i)   = 1.d30
+	    taunplus(j,i)   = 1.d30
+	    taun2plus(j,i)  = 1.d30
+            taunoplus(j,i)  = 1.d30
+	    tauhplus(j,i)   = 1.d30
+            tauhco2plus(j,i)= 1.d30
          end do
 
          !Lifetime of each species in each reaction
-         if(jdistot8(1,i).gt.1.d-30) tauco2(i,1) = 1.d0 / jdistot8(1,i)
+         if(jdistot8(1,i).gt.1.d-30) tauco2(1,i) = 1.d0 / jdistot8(1,i)
 		 
          if(ch2*o2xini*co2xini.gt.1.d-30) 
-     $      tauh(i,2) = 1.d0 / (ch2 * o2xini * co2xini)  
+     $      tauh(2,i) = 1.d0 / (ch2 * o2xini * co2xini)  
          if(ch2*hxini*co2xini.gt.1.d-30)
-     $      tauo2(i,2) = 1.d0 / (ch2 * hxini * co2xini)
+     $      tauo2(2,i) = 1.d0 / (ch2 * hxini * co2xini)
  
-         if(ch3*o3pxini.gt.1.d-30) tauho2(i,3) = 1.d0 / 
+         if(ch3*o3pxini.gt.1.d-30) tauho2(3,i) = 1.d0 / 
      $        (ch3 * o3pxini)
-         if(ch3*ho2xini.gt.1.d-30) tauo3p(i,3) = 1.d0 / 
+         if(ch3*ho2xini.gt.1.d-30) tauo3p(3,i) = 1.d0 / 
      $        (ch3 * ho2xini)
  
-         if(ch4*coxini.gt.1.d-30) tauoh(i,4) = 1.d0 / 
+         if(ch4*coxini.gt.1.d-30) tauoh(4,i) = 1.d0 / 
      $        (ch4 * coxini)
-         if(ch4*ohxini.gt.1.d-30) tauco(i,4) = 1.d0 / 
+         if(ch4*ohxini.gt.1.d-30) tauco(4,i) = 1.d0 / 
      $        (ch4 * ohxini)
  
-         if(ch5*ho2xini.gt.1.d-30)tauho2(i,5)=1.d0 / 
+         if(ch5*ho2xini.gt.1.d-30)tauho2(5,i)=1.d0 / 
      $        (2.d0*ch5*ho2xini)
 
 		 
-	 if(jdistot8(6,i).gt.1.d-30) tauh2o2(i,6) = 1.d0 / jdistot8(6,i)
-
-         if(ch7*ohxini.gt.1.d-30) tauho2(i,7) = 1.d0 / 
+	 if(jdistot8(6,i).gt.1.d-30) tauh2o2(6,i) = 1.d0 / jdistot8(6,i)
+
+         if(ch7*ohxini.gt.1.d-30) tauho2(7,i) = 1.d0 / 
      $        (ch7 * ohxini)
-         if(ch7*ho2xini.gt.1.d-30) tauoh(i,7) = 1.d0 / 
+         if(ch7*ho2xini.gt.1.d-30) tauoh(7,i) = 1.d0 / 
      $        (ch7 * ho2xini)
  
-         if(jdistot8(4,i).gt.1.d-30) tauh2o(i,8) = 1.d0 / jdistot8(4,i)
-
-	 if(ch9*o1dxini.gt.1.d-30) tauh2o(i,9) = 1.d0 / 
+         if(jdistot8(4,i).gt.1.d-30) tauh2o(8,i) = 1.d0 / jdistot8(4,i)
+
+	 if(ch9*o1dxini.gt.1.d-30) tauh2o(9,i) = 1.d0 / 
      $        (ch9 * o1dxini)
-         if(ch9*h2oxini.gt.1.d-30) tauo1d(i,9) = 1.d0 / 
+         if(ch9*h2oxini.gt.1.d-30) tauo1d(9,i) = 1.d0 / 
      $        (ch9 * h2oxini)
 
          if(ch10*o3pxini*co2xini.gt.1.d-30) 
-     $     tauo3p(i,10) = 1.d0 / 
+     $     tauo3p(10,i) = 1.d0 / 
      $        (2.d0 * ch10 * o3pxini * co2xini)
 
-         if(ch11*o3pxini.gt.1.d-30) tauoh(i,11)=1.d0 / 
+         if(ch11*o3pxini.gt.1.d-30) tauoh(11,i)=1.d0 / 
      $        (ch11 * o3pxini)
-         if(ch11*ohxini.gt.1.d-30) tauo3p(i,11) = 1.d0 / 
+         if(ch11*ohxini.gt.1.d-30) tauo3p(11,i) = 1.d0 / 
      $        (ch11 * ohxini)
 
-         if(jdistot8(2,i).gt.1.d-30) tauo2(i,12) = 1.d0 / jdistot8(2,i)
-
-         if(ch13*hxini.gt.1.d-30) tauho2(i,13) = 1.d0 / 
+         if(jdistot8(2,i).gt.1.d-30) tauo2(12,i) = 1.d0 / jdistot8(2,i)
+
+         if(ch13*hxini.gt.1.d-30) tauho2(13,i) = 1.d0 / 
      $        (ch13 * hxini)
-         if(ch13*ho2xini.gt.1.d-30) tauh(i,13) = 1.d0 / 
+         if(ch13*ho2xini.gt.1.d-30) tauh(13,i) = 1.d0 / 
      $        (ch13 * ho2xini)
 
-         if(ch14*o1dxini.gt.1.d-30) tauh2(i,14) = 1.d0 / 
+         if(ch14*o1dxini.gt.1.d-30) tauh2(14,i) = 1.d0 / 
      $        (ch14 * o1dxini)
-         if(ch14*h2xini.gt.1.d-30) tauo1d(i,14) = 1.d0 / 
+         if(ch14*h2xini.gt.1.d-30) tauo1d(14,i) = 1.d0 / 
      $        (ch14 * h2xini)
 
-         if(ch15*ohxini.gt.1.d-30) tauh2(i,15) = 1.d0 / 
+         if(ch15*ohxini.gt.1.d-30) tauh2(15,i) = 1.d0 / 
      $        (ch15 * ohxini)
-         if(ch15*h2xini.gt.1.d-30) tauoh(i,15) = 1.d0 / 
+         if(ch15*h2xini.gt.1.d-30) tauoh(15,i) = 1.d0 / 
      $        (ch15 * h2xini)
 
-         if(jdistot8_b(1,i).gt.1.d-30) tauco2(i,16)=1.d0/jdistot8_b(1,i)
-
-         if(jdistot8_b(2,i).gt.1.d-30) tauo2(i,17)=1.d0/jdistot8_b(2,i)
-
-         if(ch18*ohxini.gt.1.d-30) tauh2o2(i,18)=1.d0 / 
+         if(jdistot8_b(1,i).gt.1.d-30) tauco2(16,i)=1.d0/jdistot8_b(1,i)
+
+         if(jdistot8_b(2,i).gt.1.d-30) tauo2(17,i)=1.d0/jdistot8_b(2,i)
+
+         if(ch18*ohxini.gt.1.d-30) tauh2o2(18,i)=1.d0 / 
      $        (ch18 * ohxini)
-         if(ch18*h2o2xini.gt.1.d-30) tauoh(i,18)=1.d0 / 
+         if(ch18*h2o2xini.gt.1.d-30) tauoh(18,i)=1.d0 / 
      $        (ch18 * h2o2xini)
 
-         if(ch19*co2xini.gt.1.d-30)tauo1d(i,19)=1.d0 / 
+         if(ch19*co2xini.gt.1.d-30)tauo1d(19,i)=1.d0 / 
      $        (ch19 * co2xini)
 
-         if(ch20*o2xini.gt.1.d-30)tauo1d(i,20)= 1.d0 / 
+         if(ch20*o2xini.gt.1.d-30)tauo1d(20,i)= 1.d0 / 
      $        (ch20 * o2xini)
 
-         if(ch21*o2xini*co2xini.gt.1.d-30)tauo3p(i,21)= 1.d0 / 
+         if(ch21*o2xini*co2xini.gt.1.d-30)tauo3p(21,i)= 1.d0 / 
      $        (ch21 * o2xini * co2xini)
-         if(ch21*o3pxini*co2xini.gt.1.d-30) tauo2(i,21) = 1.d0 / 
+         if(ch21*o3pxini*co2xini.gt.1.d-30) tauo2(21,i) = 1.d0 / 
      $        (ch21 * o3pxini * co2xini)
 
          !Only if O3, N or ion chemistry requested
          if(chemthermod.ge.1) then
-            if(ch22*hxini.gt.1.d-30) tauo3(i,22) = 1.d0 / 
+            if(ch22*hxini.gt.1.d-30) tauo3(22,i) = 1.d0 / 
      $           (ch22 * hxini)
-            if(ch22*o3xini.gt.1.d-30) tauh(i,22) = 1.d0 / 
+            if(ch22*o3xini.gt.1.d-30) tauh(22,i) = 1.d0 / 
      $           (ch22 * o3xini)
             
-            if(ch23*ohxini.gt.1.d-30) tauo3(i,23) = 1.d0 / 
+            if(ch23*ohxini.gt.1.d-30) tauo3(23,i) = 1.d0 / 
      $           (ch23 * ohxini)
-            if(ch23*o3xini.gt.1.d-30) tauoh(i,23) = 1.d0 / 
+            if(ch23*o3xini.gt.1.d-30) tauoh(23,i) = 1.d0 / 
      $           (ch23 * o3xini)
             
-            if(ch24 * ho2xini.gt.1.d-30)tauo3(i,24)= 1.d0 / 
+            if(ch24 * ho2xini.gt.1.d-30)tauo3(24,i)= 1.d0 / 
      $           (ch24 * ho2xini)
-            if(ch24 * o3xini.gt.1.d-30)tauho2(i,24)= 1.d0 / 
+            if(ch24 * o3xini.gt.1.d-30)tauho2(24,i)= 1.d0 / 
      $           (ch24 * o3xini)
             
-            if(jdistot8(7,i).gt.1.d-30) tauo3(i,25)=1.d0 /
+            if(jdistot8(7,i).gt.1.d-30) tauo3(25,i)=1.d0 /
      $           jdistot8(7,i)
 
-            if(jdistot8_b(7,i).gt.1.d-30) tauo3(i,26)= 1.d0 /
+            if(jdistot8_b(7,i).gt.1.d-30) tauo3(26,i)= 1.d0 /
      $           jdistot8_b(7,i)
 
          endif    !Of chemthermod.ge.1
 
-         if(jdistot8(5,i).gt.1.d-30) tauh2(i,27)= 1.d0/jdistot8(5,i)
+         if(jdistot8(5,i).gt.1.d-30) tauh2(27,i)= 1.d0/jdistot8(5,i)
 
 !         if(ch28 * h2oxini.gt.1.d-30) tauo3(i,28) = 1.d0/(ch28*h2oxini)
@@ -1957,86 +1956,86 @@
          !Only if N or ion chemistry requested
          if(chemthermod.ge.2) then
-            if(jdistot8(8,i).gt.1.d-30) taun2(i,28) = 1.d0 / 
+            if(jdistot8(8,i).gt.1.d-30) taun2(28,i) = 1.d0 / 
      $           jdistot8(8,i)
 
-            if(jdistot8(10,i).gt.1.d-30) tauno(i,29) = 1.d0 /
+            if(jdistot8(10,i).gt.1.d-30) tauno(29,i) = 1.d0 /
      $           jdistot8(10,i)
 
-            if(ch30 * noxini.gt.1.d-30) taun(i,30) = 1.d0 / 
+            if(ch30 * noxini.gt.1.d-30) taun(30,i) = 1.d0 / 
      $           (ch30 * noxini)
-            if(ch30 * nxini.gt.1.d-30) tauno(i,30) = 1.d0 / 
+            if(ch30 * nxini.gt.1.d-30) tauno(30,i) = 1.d0 / 
      $           (ch30 * nxini)
 
-            if(ch31 * o1dxini.gt.1.d-30) taun2(i,31) = 1.d0 / 
+            if(ch31 * o1dxini.gt.1.d-30) taun2(31,i) = 1.d0 / 
      $           (ch31 * o1dxini)
-            if(ch31 * n2xini.gt.1.d-30) tauo1d(i,31) = 1.d0 /
+            if(ch31 * n2xini.gt.1.d-30) tauo1d(31,i) = 1.d0 /
      $           (ch31 * n2xini)
 
-            if(ch32 * o2xini.gt.1.d-30) taun(i,32) = 1.d0 /
+            if(ch32 * o2xini.gt.1.d-30) taun(32,i) = 1.d0 /
      $           (ch32 * o2xini)
-            if(ch32 * nxini.gt.1.d-30) tauo2(i,32) = 1.d0 /
+            if(ch32 * nxini.gt.1.d-30) tauo2(32,i) = 1.d0 /
      $           (ch32 * nxini)
             
-            if(ch33 * ohxini.gt.1.d-30) taun(i,33) = 1.d0 /
+            if(ch33 * ohxini.gt.1.d-30) taun(33,i) = 1.d0 /
      $           (ch33 * ohxini)
-            if(ch33 * nxini.gt.1.d-30) tauoh(i,33) = 1.d0 /
+            if(ch33 * nxini.gt.1.d-30) tauoh(33,i) = 1.d0 /
      $           (ch33 * nxini)
 
-            if(ch34 * o3xini.gt.1.d-30) taun(i,34) = 1.d0 /
+            if(ch34 * o3xini.gt.1.d-30) taun(34,i) = 1.d0 /
      $           (ch34 * o3xini)
-            if(ch34 * nxini.gt.1.d-30) tauo3(i,34) = 1.d0 /
+            if(ch34 * nxini.gt.1.d-30) tauo3(34,i) = 1.d0 /
      $           (ch34 * nxini)
             
-            if(ch35 * ho2xini.gt.1.d-30) taun(i,35) = 1.d0 /
+            if(ch35 * ho2xini.gt.1.d-30) taun(35,i) = 1.d0 /
      $           (ch35 * ho2xini)
-            if(ch35 * nxini.gt.1.d-30) tauho2(i,35) = 1.d0 /
+            if(ch35 * nxini.gt.1.d-30) tauho2(35,i) = 1.d0 /
      $           (ch35 * nxini)
             
-            if(ch36 * o3pxini.gt.1.d-30) taun2d(i,36) = 1.d0 /
+            if(ch36 * o3pxini.gt.1.d-30) taun2d(36,i) = 1.d0 /
      $           (ch36 * o3pxini)
-            if(ch36 * n2dxini.gt.1.d-30) tauo3p(i,36) = 1.d0 /
+            if(ch36 * n2dxini.gt.1.d-30) tauo3p(36,i) = 1.d0 /
      $           (ch36 * n2dxini)
             
-            if(ch37 * n2xini.gt.1.d-30) taun2d(i,37) = 1.d0 /
+            if(ch37 * n2xini.gt.1.d-30) taun2d(37,i) = 1.d0 /
      $           (ch37 * n2xini)
-            if(ch37 * n2dxini.gt.1.d-30) taun2(i,37) = 1.d0 /
+            if(ch37 * n2dxini.gt.1.d-30) taun2(37,i) = 1.d0 /
      $           (ch37 * n2dxini)
             
-            if(ch38 * co2xini.gt.1.d-30) taun2d(i,38) = 1.d0 /
+            if(ch38 * co2xini.gt.1.d-30) taun2d(38,i) = 1.d0 /
      $           (ch38 * co2xini)
-            if(ch38 * n2dxini.gt.1.d-30) tauco2(i,38) = 1.d0 /
+            if(ch38 * n2dxini.gt.1.d-30) tauco2(38,i) = 1.d0 /
      $           (ch38 * n2dxini)
             
-            if(ch39 * ho2xini.gt.1.d-30) tauno(i,39) = 1.d0 /
+            if(ch39 * ho2xini.gt.1.d-30) tauno(39,i) = 1.d0 /
      $           (ch39 * ho2xini)
-            if(ch39 * noxini.gt.1.d-30) tauho2(i,39) = 1.d0 /
+            if(ch39 * noxini.gt.1.d-30) tauho2(39,i) = 1.d0 /
      $           (ch39 * noxini)
             
-            if(ch40 * noxini * co2xini.gt.1.d-30) tauo3p(i,40) = 1.d0 /
+            if(ch40 * noxini * co2xini.gt.1.d-30) tauo3p(40,i) = 1.d0 /
      $           (ch40 * noxini * co2xini)
-            if(ch40 * o3pxini * co2xini.gt.1.d-30) tauno(i,40) = 1.d0 /
+            if(ch40 * o3pxini * co2xini.gt.1.d-30) tauno(40,i) = 1.d0 /
      $           (ch40 * o3pxini * co2xini)
             
-            if(ch41 * no2xini.gt.1.d-30) tauo3p(i,41) = 1.d0 /
+            if(ch41 * no2xini.gt.1.d-30) tauo3p(41,i) = 1.d0 /
      $           (ch41 * no2xini)
-            if(ch41 * o3pxini.gt.1.d-30) tauno2(i,41) = 1.d0 /
+            if(ch41 * o3pxini.gt.1.d-30) tauno2(41,i) = 1.d0 /
      $           (ch41 * o3pxini)
             
-            if(ch42 * noxini.gt.1.d-30) tauo3(i,42) = 1.d0 /
+            if(ch42 * noxini.gt.1.d-30) tauo3(42,i) = 1.d0 /
      $           (ch42 * noxini)
-            if(ch42 * o3xini.gt.1.d-30) tauno(i,42) = 1.d0 /
+            if(ch42 * o3xini.gt.1.d-30) tauno(42,i) = 1.d0 /
      $           (ch42 * o3xini)
             
-            if(ch43 * no2xini.gt.1.d-30) tauh(i,43) = 1.d0 /
+            if(ch43 * no2xini.gt.1.d-30) tauh(43,i) = 1.d0 /
      $           (ch43 * no2xini)
-            if(ch43 * hxini.gt.1.d-30) tauno2(i,43) = 1.d0 /
+            if(ch43 * hxini.gt.1.d-30) tauno2(43,i) = 1.d0 /
      $           (ch43 * hxini)
 
-            if(jdistot8(13,i).gt.1.d-30) tauno2(i,44) = 1.d0 / 
+            if(jdistot8(13,i).gt.1.d-30) tauno2(44,i) = 1.d0 / 
      $           jdistot8(13,i)
 
-            if(ch45 * nxini.gt.1.d-30) tauo3p(i,45) = 1.d0 / 
+            if(ch45 * nxini.gt.1.d-30) tauo3p(45,i) = 1.d0 / 
      $           (ch45 * nxini)
-            if(ch45 * o3pxini.gt.1.d-30) taun(i,45) = 1.d0 /
+            if(ch45 * o3pxini.gt.1.d-30) taun(45,i) = 1.d0 /
      $           (ch45 * o3pxini)
 
@@ -2047,198 +2046,198 @@
          !Only if ion chemistry requested
          if(chemthermod.eq.3) then
-            if(ch46 * co2plusxini .gt.1.d-30) tauo2(i,46) = 
+            if(ch46 * co2plusxini .gt.1.d-30) tauo2(46,i) = 
      @           1.d0/(ch46*co2plusxini)
-            if(ch46 * o2xini .gt.1.d-30) tauco2plus(i,46) = 
+            if(ch46 * o2xini .gt.1.d-30) tauco2plus(46,i) = 
      @           1.d0/(ch46*o2xini)
 
-            if ( ch47*o3pxini .gt. 1.d-30 ) tauco2plus(i,47) = 
+            if ( ch47*o3pxini .gt. 1.d-30 ) tauco2plus(47,i) = 
      @           1.d0/( ch47*o3pxini )
-            if ( ch47*co2plusxini .gt. 1.d-30 ) tauo3p(i,47) = 
+            if ( ch47*co2plusxini .gt. 1.d-30 ) tauo3p(47,i) = 
      @           1.d0/( ch47*co2plusxini )
             
-            if ( ch48*o3pxini .gt. 1.d-30 ) tauco2plus(i,48) =
+            if ( ch48*o3pxini .gt. 1.d-30 ) tauco2plus(48,i) =
      @           1.d0/(ch48*o3pxini)	   
-            if ( ch48*co2plusxini .gt. 1.d-30 ) tauo3p(i,48) =
+            if ( ch48*co2plusxini .gt. 1.d-30 ) tauo3p(48,i) =
      @           1.d0/(ch48*co2plusxini)	 
 
-            if ( ch49*electxini .gt. 1.d-30 ) tauo2plus(i,49) =
+            if ( ch49*electxini .gt. 1.d-30 ) tauo2plus(49,i) =
      @           1.d0/(ch49*electxini)	 
 
-            if ( ch50*co2xini  .gt. 1.d-30 ) tauoplus(i,50) =
+            if ( ch50*co2xini  .gt. 1.d-30 ) tauoplus(50,i) =
      @           1.d0/(ch50*co2xini)
-            if ( ch50*oplusxini .gt. 1.d-30 ) tauco2(i,50) =
+            if ( ch50*oplusxini .gt. 1.d-30 ) tauco2(50,i) =
      @           1.d0/(ch50*oplusxini)
 
-            if ( jion8(1,i,1).gt.1.d-30 ) tauco2(i,51) = 
+            if ( jion8(1,i,1).gt.1.d-30 ) tauco2(51,i) = 
      $           1.d0 / jion8(1,i,1)
 
-            if ( jion8(1,i,2).gt.1.d-30 ) tauco2(i,52) = 
+            if ( jion8(1,i,2).gt.1.d-30 ) tauco2(52,i) = 
      $           1.d0 / jion8(1,i,2)
 
-            if ( jion8(1,i,3).gt.1.d-30 ) tauco2(i,53) = 
+            if ( jion8(1,i,3).gt.1.d-30 ) tauco2(53,i) = 
      $           1.d0 / jion8(1,i,3)
 
-            if ( jion8(1,i,4).gt.1.d-30 ) tauco2(i,54) = 
+            if ( jion8(1,i,4).gt.1.d-30 ) tauco2(54,i) = 
      $           1.d0 / jion8(1,i,4) 
 		
-            if ( ch55*electxini .gt. 1.d-30 ) tauco2plus(i,55) =
+            if ( ch55*electxini .gt. 1.d-30 ) tauco2plus(55,i) =
      @           1.d0/(ch55*electxini)	 
 
-            if ( ch56*oplusxini .gt. 1.d-30 ) tauco2(i,56) =
+            if ( ch56*oplusxini .gt. 1.d-30 ) tauco2(56,i) =
      @           1.d0/(ch56*oplusxini)
-            if ( ch56*co2xini .gt. 1.d-30 ) tauoplus(i,56) =
+            if ( ch56*co2xini .gt. 1.d-30 ) tauoplus(56,i) =
      @           1.d0/(ch56*co2xini)	
     
-            if ( ch57*coplusxini .gt. 1.d-30 ) tauco2(i,57) =
+            if ( ch57*coplusxini .gt. 1.d-30 ) tauco2(57,i) =
      @           1.d0/(ch57*coplusxini)	 
-            if ( ch57*co2xini .gt. 1.d-30 ) taucoplus(i,57) =
+            if ( ch57*co2xini .gt. 1.d-30 ) taucoplus(57,i) =
      @           1.d0/(ch57*co2xini)
 
-            if ( ch58*coplusxini .gt. 1.d-30 ) tauo3p(i,58) =
+            if ( ch58*coplusxini .gt. 1.d-30 ) tauo3p(58,i) =
      @           1.d0/(ch58*coplusxini)	
-            if ( ch58*o3pxini .gt. 1.d-30 ) taucoplus(i,58) =
+            if ( ch58*o3pxini .gt. 1.d-30 ) taucoplus(58,i) =
      @           1.d0/(ch58*o3pxini)
 
-            if ( ch59*cplusxini .gt. 1.d-30 ) tauco2(i,59) =
+            if ( ch59*cplusxini .gt. 1.d-30 ) tauco2(59,i) =
      @           1.d0/(ch59*cplusxini)	
-            if ( ch59*co2xini .gt. 1.d-30 ) taucplus(i,59) =
+            if ( ch59*co2xini .gt. 1.d-30 ) taucplus(59,i) =
      @           1.d0/(ch59*co2xini)	   	 	 
          
-            if ( jion8(2,i,1).gt.1.d-30 ) tauo2(i,60) = 
+            if ( jion8(2,i,1).gt.1.d-30 ) tauo2(60,i) = 
      $           1.d0 / jion8(2,i,1)
 
-            if ( jion8(3,i,1).gt.1.d-30 ) tauo3p(i,61) = 
+            if ( jion8(3,i,1).gt.1.d-30 ) tauo3p(61,i) = 
      $           1.d0 / jion8(3,i,1)  
 
-            if ( ch62*co2plusxini .gt. 1.d-30 ) tauno(i,62) =
+            if ( ch62*co2plusxini .gt. 1.d-30 ) tauno(62,i) =
      @           1.d0/(ch62*co2plusxini)	
-            if ( ch62*noxini .gt. 1.d-30 ) tauco2plus(i,62) =
+            if ( ch62*noxini .gt. 1.d-30 ) tauco2plus(62,i) =
      @           1.d0/(ch62*noxini)	   	 	 
  
-            if ( ch63*co2plusxini .gt. 1.d-30 ) taun(i,63) =
+            if ( ch63*co2plusxini .gt. 1.d-30 ) taun(63,i) =
      @           1.d0/(ch63*cplusxini)	
-            if ( ch63*nxini .gt. 1.d-30 ) tauco2plus(i,63) =
+            if ( ch63*nxini .gt. 1.d-30 ) tauco2plus(63,i) =
      @           1.d0/(ch63*nxini)	   	 	 
      
-            if ( ch64*o2plusxini .gt. 1.d-30 ) tauno(i,64) =
+            if ( ch64*o2plusxini .gt. 1.d-30 ) tauno(64,i) =
      @           1.d0/(ch64*o2plusxini)	
-            if ( ch64*noxini .gt. 1.d-30 ) tauo2plus(i,64) =
+            if ( ch64*noxini .gt. 1.d-30 ) tauo2plus(64,i) =
      @           1.d0/(ch64*noxini)	   	 	 
    
-            if ( ch65*o2plusxini .gt. 1.d-30 ) taun2(i,65) =
+            if ( ch65*o2plusxini .gt. 1.d-30 ) taun2(65,i) =
      @           1.d0/(ch65*o2plusxini)	
-            if ( ch65*n2xini .gt. 1.d-30 ) tauo2plus(i,65) =
+            if ( ch65*n2xini .gt. 1.d-30 ) tauo2plus(65,i) =
      @           1.d0/(ch65*n2xini)	   	 	 
    
-            if ( ch66*o2plusxini .gt. 1.d-30 ) taun(i,66) =
+            if ( ch66*o2plusxini .gt. 1.d-30 ) taun(66,i) =
      @           1.d0/(ch66*o2plusxini)	
-            if ( ch66*nxini .gt. 1.d-30 ) tauo2plus(i,66) =
+            if ( ch66*nxini .gt. 1.d-30 ) tauo2plus(66,i) =
      @           1.d0/(ch66*nxini)	   	 	 
    
-            if ( ch67*oplusxini .gt. 1.d-30 ) taun2(i,67) =
+            if ( ch67*oplusxini .gt. 1.d-30 ) taun2(67,i) =
      @           1.d0/(ch67*oplusxini)	
-            if ( ch67*n2xini .gt. 1.d-30 ) tauoplus(i,67) =
+            if ( ch67*n2xini .gt. 1.d-30 ) tauoplus(67,i) =
      @           1.d0/(ch67*n2xini)	   	 	 
 
-            if ( ch68*n2plusxini .gt. 1.d-30 ) tauco2(i,68) =
+            if ( ch68*n2plusxini .gt. 1.d-30 ) tauco2(68,i) =
      @           1.d0/(ch68*n2plusxini)
-            if ( ch68*co2xini .gt. 1.d-30 ) taun2plus(i,68) =
+            if ( ch68*co2xini .gt. 1.d-30 ) taun2plus(68,i) =
      @           1.d0/(ch68*co2xini)	   	 	 
     
-            if ( ch69*cplusxini .gt. 1.d-30 ) tauco2(i,69) =
+            if ( ch69*cplusxini .gt. 1.d-30 ) tauco2(69,i) =
      @           1.d0/(ch69*cplusxini)	
-            if ( ch69*co2xini .gt. 1.d-30 ) taucplus(i,69) =
+            if ( ch69*co2xini .gt. 1.d-30 ) taucplus(69,i) =
      @           1.d0/(ch69*co2xini)	   	 	 
             
-            if ( ch70*n2plusxini .gt. 1.d-30 ) tauco(i,70) =
+            if ( ch70*n2plusxini .gt. 1.d-30 ) tauco(70,i) =
      @           1.d0/(ch70*n2plusxini)	
-            if ( ch70*coxini .gt. 1.d-30 ) taun2plus(i,70) =
+            if ( ch70*coxini .gt. 1.d-30 ) taun2plus(70,i) =
      @           1.d0/(ch70*coxini)	   	 	 
             
-            if ( ch71*electxini .gt. 1.d-30 ) taun2plus(i,71) =
+            if ( ch71*electxini .gt. 1.d-30 ) taun2plus(71,i) =
      @           1.d0/(ch71*electxini)	
             
-            if ( ch72*n2plusxini .gt. 1.d-30 ) tauo3p(i,72) =
+            if ( ch72*n2plusxini .gt. 1.d-30 ) tauo3p(72,i) =
      @           1.d0/(ch72*n2plusxini)	
-            if ( ch72*o3pxini .gt. 1.d-30 ) taun2plus(i,72) =
+            if ( ch72*o3pxini .gt. 1.d-30 ) taun2plus(72,i) =
      @           1.d0/(ch72*o3pxini)	   	 	 
      	 	 
-            if ( ch73*coplusxini .gt. 1.d-30 ) tauh(i,73) =
+            if ( ch73*coplusxini .gt. 1.d-30 ) tauh(73,i) =
      @           1.d0/(ch73*coplusxini)	
-            if ( ch73*hxini .gt. 1.d-30 ) taucoplus(i,73) =
+            if ( ch73*hxini .gt. 1.d-30 ) taucoplus(73,i) =
      @           1.d0/(ch73*hxini)	   	 	 
          
-            if ( ch74*oplusxini .gt. 1.d-30 ) tauh(i,74) =
+            if ( ch74*oplusxini .gt. 1.d-30 ) tauh(74,i) =
      @           1.d0/(ch74*oplusxini)	
-            if ( ch74*hxini .gt. 1.d-30 ) tauoplus(i,74) =
+            if ( ch74*hxini .gt. 1.d-30 ) tauoplus(74,i) =
      @           1.d0/(ch74*hxini)	
  
-            if ( ch75*electxini .gt. 1.d-30 ) taunoplus(i,75) =
+            if ( ch75*electxini .gt. 1.d-30 ) taunoplus(75,i) =
      @           1.d0/(ch75*electxini)	   	 	 
             
-            if ( ch76*hplusxini .gt. 1.d-30 ) tauo3p(i,76) =
+            if ( ch76*hplusxini .gt. 1.d-30 ) tauo3p(76,i) =
      @           1.d0/(ch76*hplusxini)	
-            if ( ch76*o3pxini .gt. 1.d-30 ) tauhplus(i,76) =
+            if ( ch76*o3pxini .gt. 1.d-30 ) tauhplus(76,i) =
      @           1.d0/(ch76*o3pxini)	
 	 
-            if( jion8(11,i,1).gt.1.d-30 )  tauco(i,77) = 
+            if( jion8(11,i,1).gt.1.d-30 )  tauco(77,i) = 
      $           1.d0 / jion8(11,i,1)
 
-            if( jion8(11,i,2).gt.1.d-30 )  tauco(i,78) = 
+            if( jion8(11,i,2).gt.1.d-30 )  tauco(78,i) = 
      $           1.d0 / jion8(11,i,2)  
 
-         !if( jion8(11,i,3).gt.1.d-30 ) tauco(i,79) = 
+         !if( jion8(11,i,3).gt.1.d-30 ) tauco(79,i) = 
 !     $        1.d0 / jion8(11,i,3) 
 
-            if( jion8(10,i,1).gt.1.d-30 )  tauno(i,80) = 
+            if( jion8(10,i,1).gt.1.d-30 )  tauno(80,i) = 
      $           1.d0 / jion8(10,i,1) 
          
-            if( jion8(8,i,1).gt.1.d-30 )   taun2(i,81)  = 
+            if( jion8(8,i,1).gt.1.d-30 )   taun2(81,i)  = 
      $           1.d0 / jion8(8,i,1) 
 
-            if( jion8(8,i,2).gt.1.d-30 )   taun2(i,82)  = 
+            if( jion8(8,i,2).gt.1.d-30 )   taun2(82,i)  = 
      $           1.d0 / jion8(8,i,2)
 		 
-            if( jion8(12,i,1).gt.1.d-30 )  tauh(i,83)  = 
+            if( jion8(12,i,1).gt.1.d-30 )  tauh(83,i)  = 
      $           1.d0 / jion8(12,i,1)
 
-            if( jion8(9,i,1).gt.1.d-30 )   taun(i,84)  = 
+            if( jion8(9,i,1).gt.1.d-30 )   taun(84,i)  = 
      $           1.d0 / jion8(9,i,1)
 
-            if ( ch85*nplusxini .gt. 1.d-30 ) tauco2(i,85) =
+            if ( ch85*nplusxini .gt. 1.d-30 ) tauco2(85,i) =
      @           1.d0/(ch85*nplusxini)	
-            if ( ch85*co2xini .gt. 1.d-30 ) taunplus(i,85) =
+            if ( ch85*co2xini .gt. 1.d-30 ) taunplus(85,i) =
      @           1.d0/(ch85*co2xini)
 
-            if ( ch86*co2plusxini .gt. 1.d-30) tauh2(i,86) =
+            if ( ch86*co2plusxini .gt. 1.d-30) tauh2(86,i) =
      $           1.d0/(ch86*co2plusxini)
-            if ( ch86*h2xini .gt. 1.d-30) tauco2plus(i,86) =
+            if ( ch86*h2xini .gt. 1.d-30) tauco2plus(86,i) =
      $           1.d0/(ch86*h2xini)
          
-            if ( ch87*electxini .gt. 1.d-30) tauhco2plus(i,87) =
+            if ( ch87*electxini .gt. 1.d-30) tauhco2plus(87,i) =
      $                  1.d0/(ch87*electxini)
 
             if ( jion8(9,i,1)*ionsec_nplus(zenit,zx(i)).gt.1.d-30) 
-     $           taun(i,88) = 1.d0 / 
+     $           taun(88,i) = 1.d0 / 
      $           (jion8(9,i,1)*ionsec_nplus(zenit,zx(i)))
 
             if ( jion8(8,i,1)*ionsec_n2plus(zenit,zx(i)).gt.1.d-30) 
-     $           taun2(i,89) = 1.d0 / 
+     $           taun2(89,i) = 1.d0 / 
      $           (jion8(8,i,1)*ionsec_n2plus(zenit,zx(i))) 
             
             if ( jion8(3,i,1)*ionsec_oplus(zenit,zx(i)).gt.1.d-30) 
-     $           tauo3p(i,90) = 1.d0 /
+     $           tauo3p(90,i) = 1.d0 /
      $           (jion8(3,i,1)*ionsec_oplus(zenit,zx(i)))
             
             if (jion8(11,i,1)*ionsec_coplus(zenit,zx(i)).gt.1.d-30) 
-     $           tauco(i,91) = 1.d0 /
+     $           tauco(91,i) = 1.d0 /
      $           (jion8(11,i,1)*ionsec_coplus(zenit,zx(i)))
 
             if (jion8(1,i,1)*ionsec_co2plus(zenit,zx(i)).gt.1.d-30) 
-     $           tauco2(i,92) = 1.d0 /
+     $           tauco2(92,i) = 1.d0 /
      $           (jion8(1,i,1)*ionsec_co2plus(zenit,zx(i)))
 
             if ( jion8(2,i,1)*ionsec_o2plus(zenit,zx(i)).gt.1.d-30) 
-     $           tauo2(i,93) = 1.d0 /
+     $           tauo2(93,i) = 1.d0 /
      $           (jion8(2,i,1)*ionsec_o2plus(zenit,zx(i)))
 
@@ -2276,32 +2275,32 @@
 
          do j=1,nreact
-            tminco2(i)  = min(tminco2(i),tauco2(i,j))
-            tmino2(i)   = min(tmino2(i),tauo2(i,j))
-            tmino3p(i)  = min(tmino3p(i),tauo3p(i,j))
-            tminco(i)   = min(tminco(i),tauco(i,j))
-            tminh(i)    = min(tminh(i),tauh(i,j))
-            tminoh(i)   = min(tminoh(i),tauoh(i,j))
-            tminho2(i)  = min(tminho2(i),tauho2(i,j))
-            tminh2(i)   = min(tminh2(i),tauh2(i,j))
-            tminh2o(i)  = min(tminh2o(i),tauh2o(i,j))
-            tmino1d(i)  = min(tmino1d(i),tauo1d(i,j))
-            tminh2o2(i) = min(tminh2o2(i),tauh2o2(i,j))
-            tmino3(i)   = min(tmino3(i),tauo3(i,j))
-            tminn(i)    = min(tminn(i),taun(i,j))
-            tminno(i)   = min(tminno(i),tauno(i,j))
-	    tminn2(i)   = min(tminn2(i),taun2(i,j))
-	    tminn2d(i)  = min(tminn2d(i),taun2d(i,j))
-            tminno2(i)  = min(tminno2(i),tauno2(i,j))
+            tminco2(i)  = min(tminco2(i),tauco2(j,i))
+            tmino2(i)   = min(tmino2(i),tauo2(j,i))
+            tmino3p(i)  = min(tmino3p(i),tauo3p(j,i))
+            tminco(i)   = min(tminco(i),tauco(j,i))
+            tminh(i)    = min(tminh(i),tauh(j,i))
+            tminoh(i)   = min(tminoh(i),tauoh(j,i))
+            tminho2(i)  = min(tminho2(i),tauho2(j,i))
+            tminh2(i)   = min(tminh2(i),tauh2(j,i))
+            tminh2o(i)  = min(tminh2o(i),tauh2o(j,i))
+            tmino1d(i)  = min(tmino1d(i),tauo1d(j,i))
+            tminh2o2(i) = min(tminh2o2(i),tauh2o2(j,i))
+            tmino3(i)   = min(tmino3(i),tauo3(j,i))
+            tminn(i)    = min(tminn(i),taun(j,i))
+            tminno(i)   = min(tminno(i),tauno(j,i))
+	    tminn2(i)   = min(tminn2(i),taun2(j,i))
+	    tminn2d(i)  = min(tminn2d(i),taun2d(j,i))
+            tminno2(i)  = min(tminno2(i),tauno2(j,i))
             !
-            tminco2plus(i) = min(tminco2plus(i),tauco2plus(i,j))
-	    tminoplus(i)   = min(tminoplus(i),tauoplus(i,j))
-            tmino2plus(i)  = min(tmino2plus(i),tauo2plus(i,j))
-            tmincoplus(i)  = min(tmincoplus(i),taucoplus(i,j))
-            tmincplus(i)   = min(tmincplus(i),taucplus(i,j))
-	    tminnplus(i)   = min(tminnplus(i),taunplus(i,j))
-            tminn2plus(i)  = min(tminn2plus(i),taun2plus(i,j))
-            tminnoplus(i)  = min(tminnoplus(i),taunoplus(i,j))
-            tminhplus(i)   = min(tminhplus(i),tauhplus(i,j))
-            tminhco2plus(i)= min(tminhco2plus(i),tauhco2plus(i,j))
+            tminco2plus(i) = min(tminco2plus(i),tauco2plus(j,i))
+	    tminoplus(i)   = min(tminoplus(i),tauoplus(j,i))
+            tmino2plus(i)  = min(tmino2plus(i),tauo2plus(j,i))
+            tmincoplus(i)  = min(tmincoplus(i),taucoplus(j,i))
+            tmincplus(i)   = min(tmincplus(i),taucplus(j,i))
+	    tminnplus(i)   = min(tminnplus(i),taunplus(j,i))
+            tminn2plus(i)  = min(tminn2plus(i),taun2plus(j,i))
+            tminnoplus(i)  = min(tminnoplus(i),taunoplus(j,i))
+            tminhplus(i)   = min(tminhplus(i),tauhplus(j,i))
+            tminhco2plus(i)= min(tminhco2plus(i),tauhco2plus(j,i))
          end do
 
@@ -2441,5 +2440,5 @@
 
       subroutine timemarching( ig,i,chemthermod,n_comp_en_EQ,
-     $     compmin,tmin,timefrac_sec, deltat ) 
+     $     compmin,tmin,timefrac_sec, deltat,fmargin1 ) 
 
   
@@ -2478,5 +2477,5 @@
 ccccccccccccccc CODE STARTS 
 
-      fmargin1=1.
+!      fmargin1=1.
       fmargin2=5.
 
@@ -3661,5 +3660,4 @@
       Lo1dtot(i) = Lo1d(i,9) + Lo1d(i,14) + Lo1d(i,19) + Lo1d(i,20) +
      $     Lo1d(i,31)
-      
 
 c     H2O2
Index: /trunk/LMDZ.MARS/libf/aeronomars/perosat.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/perosat.F	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/perosat.F	(revision 658)
@@ -127,4 +127,6 @@
 c
            zynew(1)=zysat(1)
+        else
+          pdqscloud(ig,igcm_h2o2)=0
         end if
 c
Index: /trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F	(revision 657)
+++ /trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F	(revision 658)
@@ -79,5 +79,5 @@
       if (callmoldiff) THEN
         call zerophys(ngridmx*nlayermx*nqmx,zdqmoldiff)
-        call moldiff(pplay,pplev,pt,pdt,pq,pdq,ptimestep,
+        call moldiff_red(pplay,pplev,pt,pdt,pq,pdq,ptimestep,
      &                   zzlay,zdteuv,zdtconduc,zdqmoldiff)
       endif
