Index: trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90	(revision 2464)
@@ -13,5 +13,5 @@
                  i_cocl2, i_s, i_so, i_so2, i_so3,       &
                  i_s2o2, i_ocs, i_hso3, i_h2so4, i_s2,   &
-                 i_clso2, i_oscl, i_n2 
+                 i_clso2, i_oscl, i_n2, i_he 
                  
 INTEGER, SAVE :: i_h2oliq, i_h2so4liq
@@ -834,4 +834,7 @@
 		i_n2=i
 		M_tr(i_n2)=28.013
+		CASE('he')
+		i_he=i
+		M_tr(i_he)=4.0026
 ! MICROPHYSICAL TRACERS FOR CL_SCHEME=1
   		CASE('h2oliq')
Index: trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 2464)
@@ -23,5 +23,5 @@
        REAL    z0, lmixmin
        REAL    ksta, inertie
-       REAL    euveff, solarcondate
+       REAL    euveff, fixed_euv_value
 
        COMMON/clesphys_l/ cycle_diurne, soil_model,                     &
@@ -37,4 +37,4 @@
 
        COMMON/clesphys_r/ ecriphy, solaire, z0, lmixmin,                &
-     &     ksta, inertie, euveff, solarcondate
+     &     ksta, inertie, euveff, fixed_euv_value
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F	(revision 2464)
@@ -96,5 +96,5 @@
      $            + d_tr_chem(:,:,i_h2so4liq)*ptimestep
 
-      wgt_SA(:,:) = wh2so4(:,:)
+      wgt_SA(:,1:n_lev) = wh2so4(:,:)
 
 c     Init F_sed
Index: trunk/LMDZ.VENUS/libf/phyvenus/compo_hedin_mod2.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/compo_hedin_mod2.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/compo_hedin_mod2.F90	(revision 2464)
@@ -97,6 +97,6 @@
          enddo
 
-!            print*, pres_hedin(1,:)
-!	    print*, ' '
+!           print*, pres_hedin(1,:)
+!           print*, ' '
 !           print*, T(musize,:)  
 !	    print*, ' '
@@ -112,15 +112,15 @@
 !		print*, " "
 !		print*, ntot(1,:)
-!		print*, co2_hedin(:,10)
+!               print*, co2_hedin(1,:)
 !		print*, mu_hedin(:)
-!		print*, " "
-!		print*, n2_hedin(1,:)
-!		print*, " "
-!		print*, o_hedin(1,:)
-!		print*, " "
-!		print*, co_hedin(1,:)
-!		print*, ' '
-!	 stop
-		
+!               print*, " "
+!               print*, n2_hedin(1,:)
+!               print*, " "
+!               print*, o_hedin(1,:)
+!               print*, " "
+!               print*, co_hedin(1,:)
+!               print*, ' '
+!           stop
+
        end subroutine compo_hedin83_init2
 
@@ -165,7 +165,11 @@
           enddo
 
-        factp(k) = (log10(pression(i,k))-log10(pres_hedin(mu_ok,z_ok-1)))/(log10(pres_hedin(mu_ok,z_ok))-log10(pres_hedin(mu_ok,z_ok-1)))
+        if (pres_hedin(mu_ok,zsize).ge.pression(i,k)) then
+! avoid the extrapolation... Problem when p < 1.E-6 Pa...
+         factp(k) = 1.
+        else
+         factp(k) = (log10(pression(i,k))-log10(pres_hedin(mu_ok,z_ok-1)))/(log10(pres_hedin(mu_ok,z_ok))-log10(pres_hedin(mu_ok,z_ok-1)))
+        endif
 
-      	
         ovmr_gcm(i,k) = (10.**(log10(o_hedin(mu_ok,z_ok))*factp(k)+log10(o_hedin(mu_ok,z_ok-1))*(1.-factp(k))))
         n2vmr_gcm(i,k) = 10.**(log10(n2_hedin(mu_ok,z_ok))*factp(k)+log10(n2_hedin(mu_ok,z_ok-1))*(1.-factp(k)))
Index: trunk/LMDZ.VENUS/libf/phyvenus/conc.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/conc.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/conc.F90	(revision 2464)
@@ -14,52 +14,4 @@
       REAL, ALLOCATABLE, SAVE :: cpnew(:,:)
 !$OMP THREADPRIVATE(cpnew)
-      REAL, ALLOCATABLE, SAVE :: jfotsout(:,:,:)
-!$OMP THREADPRIVATE(jfotsout)
-      REAL, ALLOCATABLE, SAVE :: coefit4(:,:)
-!$OMP THREADPRIVATE(coefit4)
-      REAL, ALLOCATABLE, SAVE :: coefit3(:,:)
-!$OMP THREADPRIVATE(coefit3)
-      REAL, ALLOCATABLE, SAVE :: coefit2(:,:)
-!$OMP THREADPRIVATE(coefit2)
-      REAL, ALLOCATABLE, SAVE :: coefit1(:,:)
-!$OMP THREADPRIVATE(coefit1)
-      REAL, ALLOCATABLE, SAVE :: coefit0(:,:)
-!$OMP THREADPRIVATE(coefit0)
-      REAL, ALLOCATABLE, SAVE :: fluxtop(:)
-!$OMP THREADPRIVATE(fluxtop)
-!      REAL, ALLOCATABLE, SAVE :: freccen(:)
-!$OMP THREADPRIVATE(freccen)
-      REAL, ALLOCATABLE, SAVE :: t0(:)
-!$OMP THREADPRIVATE(t0)
-!      REAL, ALLOCATABLE, SAVE :: crscabsi2(:,:)
-!$OMP THREADPRIVATE(crscabsi2)
-      REAL, ALLOCATABLE, SAVE :: jabsifotsintpar(:,:,:)
-!$OMP THREADPRIVATE(crscabsi2)
-      REAL, ALLOCATABLE, SAVE :: c1_16(:,:)
-!$OMP THREADPRIVATE(c1_16)
-      REAL, ALLOCATABLE, SAVE :: c17_24(:)
-!$OMP THREADPRIVATE(c17_24)
-      REAL, ALLOCATABLE, SAVE :: c25_29(:)
-!$OMP THREADPRIVATE(c25_29)
-      REAL, ALLOCATABLE, SAVE :: c30_31(:)
-!$OMP THREADPRIVATE(c30_31)
-      REAL, ALLOCATABLE, SAVE :: c32(:)
-!$OMP THREADPRIVATE(c32)
-      REAL, ALLOCATABLE, SAVE :: c33(:)
-!$OMP THREADPRIVATE(c33)
-      REAL, ALLOCATABLE, SAVE :: c34(:)
-!$OMP THREADPRIVATE(c34)
-      REAL, ALLOCATABLE, SAVE :: c35(:)
-!$OMP THREADPRIVATE(c35)
-      REAL, ALLOCATABLE, SAVE :: c36(:)
-!$OMP THREADPRIVATE(c36)
-      REAL, ALLOCATABLE, SAVE :: ct1(:)
-!$OMP THREADPRIVATE(ct1)
-      REAL, ALLOCATABLE, SAVE :: ct2(:)
-!$OMP THREADPRIVATE(ct2)
-      REAL, ALLOCATABLE, SAVE :: p1(:)
-!$OMP THREADPRIVATE(p1)
-      REAL, ALLOCATABLE, SAVE :: p2(:)
-!$OMP THREADPRIVATE(p2)
 
 
@@ -69,5 +21,4 @@
      USE dimphy
      implicit none
-#include"param.h"
 
      ALLOCATE(mmean(klon,klev)) 
@@ -76,28 +27,4 @@
      ALLOCATE(rnew(klon,klev))  
      ALLOCATE(cpnew(klon,klev)) 
-     ALLOCATE(jfotsout(ninter,nabs,klev))
-     ALLOCATE(coefit0(ninter,nabs))
-     ALLOCATE(coefit1(ninter,nabs))  
-     ALLOCATE(coefit2(ninter,nabs))  
-     ALLOCATE(coefit3(ninter,nabs))  
-     ALLOCATE(coefit4(ninter,nabs)) 
-     ALLOCATE(fluxtop(ninter)) 
-!     ALLOCATE(freccen(ninter)) 
-     ALLOCATE(t0(nz2))   
-!     ALLOCATE(crscabsi2(nabs,16))
-     ALLOCATE(c1_16(nz2,16))
-     ALLOCATE(c17_24(nz2) )
-     ALLOCATE(c25_29(nz2) )
-     ALLOCATE(c30_31(nz2) )
-     ALLOCATE(c32(nz2) )
-     ALLOCATE(c33(nz2) )
-     ALLOCATE(c34(nz2) )
-     ALLOCATE(c35(nz2) )
-     ALLOCATE(c36(nz2) )
-     ALLOCATE(ct2(ninter) )
-     ALLOCATE(ct1(ninter) )
-     ALLOCATE(p1(ninter) )
-     ALLOCATE(p2(ninter) )
-     ALLOCATE(jabsifotsintpar(nz2,nabs,ninter))
      end subroutine conc_init
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F	(revision 2464)
@@ -1,3 +1,3 @@
-      SUBROUTINE concentrations2(pplay,t_seri,pdt,tr_seri, nqmx)
+      SUBROUTINE concentrations2(pplay,t_seri,tr_seri, nqmx)
 
       use dimphy
@@ -31,12 +31,7 @@
 
       real pplay(klon,klev)
-c      real pt(klon,klev)
       integer,intent(in) :: nqmx    ! number of tracers
       real t_seri(klon, klev)
-      real pdt(klon,klev)
-      real n2vmr_gcm(klon,klev),nvmr_gcm(klon,klev)
       real tr_seri(klon,klev,nqmx)
-c      real pdq(klon,klev,nqmx)
-      real ptimestep
 
 !     local variables
Index: trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 2464)
@@ -454,11 +454,11 @@
 
 !
-!Config Key  = solarcondate
+!Config Key  = fixed_euv_value
 !Config Desc = 
-!Config Def  = 1993.4  ## Average solar cycle condition
-!Config Help =
-!
-  solarcondate = 1993.4
-  call getin('solarcondate',solarcondate)
+!Config Def  = 140.  ## Average solar cycle condition
+!Config Help =
+!
+  fixed_euv_value =140. 
+  call getin('fixed_euv_value',fixed_euv_value)
 
 !
@@ -528,5 +528,5 @@
   write(lunout,*)' callthermos = ',callthermos
   write(lunout,*)' solvarmod = ',solvarmod
-  write(lunout,*)' solarcondate = ',solarcondate
+  write(lunout,*)' fixed_euv_value = ',fixed_euv_value
   write(lunout,*)' euveff = ',euveff
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/euvheat.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/euvheat.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/euvheat.F90	(revision 2464)
@@ -1,8 +1,7 @@
       SUBROUTINE euvheat(nlon, nlev,nqmx, pt,pplev,pplay,zzlay, &
-        	   mu0,ptimestep,ptime,zday,       &
-                   pq, pdq, pdteuv)
+                   mu0,ptimestep,ptime,pq, pdq, pdteuv)
 
         use chemparam_mod
-	use dimphy
+        use dimphy
         use conc, only:  rnew, cpnew
 
@@ -34,7 +33,4 @@
 #include "YOMCST.h"
 #include "clesphys.h"
-#include "param.h"
-#include "param_v4.h"
-!#include "chimiedata.h"
 #include "mmol.h"
 !-----------------------------------------------------------------------
@@ -54,5 +50,4 @@
 
       real :: ptimestep,ptime
-      real :: zday
       real :: pq(nlon,nlev,nqmx)
       real :: pdq(nlon,nlev,nqmx)
@@ -336,11 +331,4 @@
       
       !Solar flux calculation
-      if(solvarmod.eq.0) call flujo(solarcondate) 
-
-                                         ! Not recommended for long runs 
-                                         ! (e.g. to build the MCD), as the
-                                         ! solar conditions at the end will 
-                                         ! be different to the ones initially
-                                         ! set
       
       do ig=1,nlon
@@ -377,5 +365,5 @@
 
         !Routine to calculate the UV heating
-         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
+         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,jtot)
 	 
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F	(revision 2464)
@@ -1,5 +1,5 @@
 c**********************************************************************
 
-      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday,jtot)
+      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,jtot)
 
 
@@ -9,5 +9,6 @@
 c**********************************************************************
       use dimphy
-      use conc
+      use param_v4_h, only: ninter,nabs,jfotsout,fluxtop,freccen
+
       implicit none
 
@@ -15,6 +16,4 @@
 
 
-#include "param.h"
-#include "param_v4.h"
 #include "clesphys.h"
 
@@ -38,5 +37,4 @@
       real       zenit
       real       iz(klev)
-      real       zday
 
       ! tracer indexes for the EUV heating:
@@ -95,12 +93,5 @@
 
       !Calculation of photoabsortion coefficient
-      if(solvarmod.eq.0) then
-         call jthermcalc(ig,euvmod,rm,nespeuv,tx,iz,zenit)
-      else if (solvarmod.eq.1) then
-         call jthermcalc_e107(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday)
-         do indexint=1,ninter
-            fluxtop(indexint)=1.
-         enddo
-      endif
+      call jthermcalc_e107(ig,klev,euvmod,rm,nespeuv,tx,iz,zenit)
 
       !Total photoabsorption coefficient    !  erg/(s*cm3) 
Index: trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc.F	(revision 2463)
+++ 	(revision )
@@ -1,1743 +1,0 @@
-c**********************************************************************
-
-      subroutine jthermcalc(ig,chemthermod,rm,nesptherm,tx,iz,zenit)
-
-
-c     feb 2002        fgg           first version
-c     nov 2002        fgg           second version
-c     
-c     mar 2014        gg            update for Venus GCM 
-c
-c modified from paramhr.F
-c MAC July 2003
-c**********************************************************************
-      use dimphy
-      use conc
-c      use chemparam_mod
-      implicit none
-
-c     common variables and constants
-      include "param.h"
-      include "param_v4.h"
-
-c     input and output variables
-
-      integer    ig
-      integer    chemthermod 
-      integer    nesptherm                  !Number of species considered
-      real       rm(klev,nesptherm)         !Densities (cm-3)
-      real       tx(klev)                   !temperature
-      real       zenit                      !SZA
-      real       iz(klev)                   !Local altitude
-
-
-c    local parameters and variables
-
-      real       co2colx(klev)              !column density of CO2 (cm^-2)
-      real       o3pcolx(klev)              !column density of O(3P)(cm^-2)
-      real       n2colx(klev)               !N2 column density (cm-2)
-      real       cocolx(klev)               !CO column density (cm-2)
-c      real       o2colx(klev)               !column density of O2(cm^-2)
-c      real       h2colx(klev)               !H2 column density (cm-2)
-c      real       h2ocolx(klev)              !H2O column density (cm-2)
-c      real       h2o2colx(klev)             !column density of H2O2(cm^-2)
-c      real       o3colx(klev)               !O3 column density (cm-2)
-c     real       hcolx(klev)                !H column density (cm-2)
-c     real       no2colx(klev)              !NO2 column density (cm-2)
-c     real       nocolx(klev)               !NO column density (cm-2)
-      real       t2(klev)
-      real       coltemp(klev)
-      real       sigma(ninter,klev)
-      real       alfa(ninter,klev)
-    
-      integer    i,j,k,indexint                 !indexes
-      character  dn
-
-
-
-c     variables used in interpolation
-
-      real*8      auxcoltab(nz2)
-      real*8      auxjco2(nz2)
-c     real*8      auxjo2(nz2)
-      real*8      auxjo3p(nz2)
-c     real*8      auxjh2o(nz2)
-c     real*8      auxjh2(nz2)
-c     real*8      auxjh2o2(nz2)
-c     real*8      auxjo3(nz2)
-      real*8      auxjn2(nz2)
-c      real*8      auxjn(nz2)
-c     real*8      auxjno(nz2)
-      real*8      auxjco(nz2)
-c     real*8      auxjh(nz2)
-c     real*8      auxjno2(nz2)
-      real*8      wp(klev),wm(klev)
-      real*8      auxcolinp(klev)
-      integer     auxind(klev)
-      integer     auxi
-      integer     ind
-      real*8      cortemp(klev)
-
-      real*8      limdown                      !limits for interpolation
-      real*8      limup                        !        ""
-
- 
-! Tracer indexes in the thermospheric chemistry:
-!!! ATTENTION. These values have to be identical to those in euvheat.F90
-!!! If the values are changed there, the same has to be done here  !!!
-
-      integer,parameter :: ix_co2=1
-      integer,parameter :: ix_n2=13
-c      integer,parameter :: i_n=14
-      integer,parameter :: ix_o=3
-      integer,parameter :: ix_co=4
-
-
-c*************************PROGRAM STARTS*******************************
-      
-      if(zenit.gt.140.) then
-         dn='n'
-         else
-         dn='d'
-      end if
-      if(dn.eq.'n') then
-        return
-      endif
-      
-      !Initializing the photoabsorption coefficients
-      jfotsout(:,:,:)=0.
-
-      !Auxiliar temperature to take into account the temperature dependence
-      !of CO2 cross section
-      do i=1,klev
-         t2(i)=tx(i)
-         if(t2(i).lt.195.0) t2(i)=195.0
-         if(t2(i).gt.295.0) t2(i)=295.0
-      end do
-
-      !Calculation of column amounts 
-      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
-     $     co2colx,o3pcolx, n2colx,cocolx)
-
-      !Auxiliar column to include the temperature dependence 
-      !of CO2 cross section
-      coltemp(klev)=co2colx(klev)*abs(t2(klev)-t0(klev))
-      do i=klev-1,1,-1
-        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR 
-     $         ( rm(i,ix_co2) + rm(i+1,ix_co2) ) * 0.5 
-     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
-      end do
-      
-      !Calculation of CO2 cross section at temperature t0(i)
-      do i=1,klev
-         do indexint=24,32
-           sigma(indexint,i)=co2crsc195(indexint-23)
-           alfa(indexint,i)=((co2crsc295(indexint-23)
-     $          /sigma(indexint,i))-1.)/(295.-t0(i))
-        end do
-      end do
-
-! Interpolation to the tabulated photoabsorption coefficients for each species
-! in each spectral interval
-
-
-c     auxcolinp-> Actual atmospheric column
-c     auxj*-> Tabulated photoabsorption coefficients
-c     auxcoltab-> Tabulated atmospheric columns
-
-ccccccccccccccccccccccccccccccc
-c     0.1,5.0 (int 1)
-c
-c     Absorption by: 
-c     CO2, O2, O, H2, N
-ccccccccccccccccccccccccccccccc
-
-c     Input atmospheric column
-      indexint=1
-      do i=1,klev
-         auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint) +
-c     $        o2colx(i)*crscabsi2(2,indexint) + 
-     $        o3pcolx(i)*crscabsi2(3,indexint) 
-c     $        h2colx(i)*crscabsi2(5,indexint) + 
-      end do
-      limdown=1.e-20
-      limup=1.e26
-
-
-c     Interpolations
-
-      do i=1,nz2
-         auxi = nz2-i+1
-         !CO2 tabulated coefficient
-         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
-         !O2 tabulated coefficient
-c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
-         !O3p tabulated coefficient
-         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
-         !H2 tabulated coefficient
-c         auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
-         !N tabulated coefficient
-c         auxjn(i) = jabsifotsintpar(auxi,9,indexint)
-         !Tabulated column
-         auxcoltab(i) = c1_16(auxi,indexint)
-      enddo
-
-      !Only if chemthermod.ge.2
-      !N tabulated coefficient
-c      if(chemthermod.ge.2) then
-c         do i=1,nz2
-c            auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
-c         enddo
-c      endif
-
-      call interfast 
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
-         ind=auxind(i)
-         auxi=klev-i+1
-         !CO2 interpolated coefficient
-         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
-     $        wp(i)*auxjco2(ind)
-         !O2 interpolated coefficient
-c         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
-c     $        wp(i)*auxjo2(ind)
-         !O3p interpolated coefficient
-          jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
-     $        wp(i)*auxjo3p(ind)
-          !H2 interpolated coefficient
-c          jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
-c     $         wp(i)*auxjh2(ind)
-c      !N interpolated coefficient
-c            jfotsout(indexint,9,auxi) =  wm(i)*auxjn(ind+1) +
-c     $         wp(i)*auxjn(ind)
-
-C      print*, '--- L214 jthermcal.F ---'
-C      print*, jfotsout(indexint,1,auxi)
-c      STOP         
-
-      enddo
-
-      !Only if chemthermod.ge.2
-      !N interpolated coefficient
-c      if(chemthermod.ge.2) then
-c         do i=1,klev
-c            ind=auxind(i)
-c            jfotsout(indexint,9,klev-i+1) =  wm(i)*auxjn(ind+1) +
-c     $         wp(i)*auxjn(ind)
-c         enddo
-c      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,klev
-            auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
-     $           o3pcolx(i)*crscabsi2(3,indexint)+
-     $           n2colx(i)*crscabsi2(8,indexint)+
-     $           cocolx(i)*crscabsi2(11,indexint)
-
-c     $           o2colx(i)*crscabsi2(2,indexint)+
-c     $           h2colx(i)*crscabsi2(5,indexint)+
-c     $           nocolx(i)*crscabsi2(10,indexint)+
-c     $           hcolx(i)*crscabsi2(12,indexint)+
-c     $           no2colx(i)*crscabsi2(13,indexint)
-         end do
-
-c     Interpolations
-
-         do i=1,nz2
-            auxi = nz2-i+1
-            !O2 tabulated coefficient
-c            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
-c            auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
-            !N2 tabulated coefficient
-            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
-            !N tabulated coefficient
-c            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
-            !CO tabulated coefficient
-            auxjco(i) = jabsifotsintpar(auxi,11,indexint)
-            !H tabulated coefficient
-c            auxjh(i) = jabsifotsintpar(auxi,12,indexint)
-            !tabulated column
-            auxcoltab(i) = c1_16(auxi,indexint)
-         enddo
-c         !Only if chemthermod.ge.2
-c         if(chemthermod.ge.2) then
-c            do i=1,nz2
-c               auxi = nz2-i+1
-c               !N tabulated coefficient
-c               auxjn(i) = jabsifotsintpar(auxi,9,indexint)
-c               !NO tabulated coefficient
-c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-c               !NO2 tabulated coefficient
-c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-c            enddo
-c         endif
-
-          call interfast(wm,wp,auxind,auxcolinp,klev,
-     $        auxcoltab,nz2,limdown,limup)
-          do i=1,klev
-             ind=auxind(i)
-             auxi = klev-i+1
-             !O2 interpolated coefficient
-c             jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
-c     $            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
-c             jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
-c     $            wp(i)*auxjh2(ind)
-             !N2 interpolated coefficient
-             jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
-     $            wp(i)*auxjn2(ind)             
-             !N interpolated coefficient
-c             jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
-c     $            wp(i)*auxjn(ind)
-             !CO interpolated coefficient
-             jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
-     $            wp(i)*auxjco(ind)
-             !H interpolated coefficient
-c             jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
-c     $            wp(i)*auxjh(ind)             
-          enddo
-          !Only if chemthermod.ge.2
-c          if(chemthermod.ge.2) then
-c             do i=1,klev
-c                ind=auxind(i)
-c                auxi = klev-i+1
-c                !N interpolated coefficient
-c                jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
-c     $               wp(i)*auxjn(ind)
-c                !NO interpolated coefficient
-c                jfotsout(indexint,10,auxi)=wm(i)*auxjno(ind+1) +
-c     $               wp(i)*auxjno(ind)
-c                !NO2 interpolated coefficient
-c                jfotsout(indexint,13,auxi)=wm(i)*auxjno2(ind+1)+
-c     $               wp(i)*auxjno2(ind)
-c             enddo
-c          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,klev
-         auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
-c     $        o2colx(i)*crscabsi2(2,indexint)+
-     $        o3pcolx(i)*crscabsi2(3,indexint)+
-     $        n2colx(i)*crscabsi2(8,indexint)+
-     $        cocolx(i)*crscabsi2(11,indexint)
-c     $        hcolx(i)*crscabsi2(12,indexint)+
-c     $        no2colx(i)*crscabsi2(13,indexint)
-      end do
-
-c     Interpolations
-
-      do i=1,nz2
-         auxi = nz2-i+1
-         !O2 tabulated coefficient
-c         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)
-c        !N tabulated coefficient
-c         auxjn(i) = jabsifotsintpar(auxi,9,indexint)
-
-         !NO tabulated coefficient
-c         auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-         !H tabulated coefficient
-c         auxjh(i) = jabsifotsintpar(auxi,12,indexint)
-         !NO2 tabulated coefficient
-c         auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-
-         !Tabulated column
-         auxcoltab(i) = c1_16(auxi,indexint)
-      enddo
-      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,nz2
-c            auxi = nz2-i+1
-c            !N tabulated coefficient
-c            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
-c            !NO tabulated coefficient
-c            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-c            !NO2 tabulated coefficient
-c            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-c         enddo
-c      endif
-
-      call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
-         ind=auxind(i)
-         auxi = klev-i+1
-         !O2 interpolated coefficient
-c         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
-c     $            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)  
-         !N interpolated coefficient
-c          jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
-c     $           wp(i)*auxjn(ind) 
-
-c         !H interpolated coefficient
-c         jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
-c     $            wp(i)*auxjh(ind)         
-      enddo
-      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,klev
-c            ind=auxind(i)
-c            auxi = klev-i+1
-c            !N interpolated coefficient
-c            jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
-c     $           wp(i)*auxjn(ind) 
-c            !NO interpolated coefficient
-c            jfotsout(indexint,10,auxi) = wm(i)*auxjno(ind+1) +
-c     $           wp(i)*auxjno(ind)
-c            !NO2 interpolated coefficient
-c            jfotsout(indexint,13,auxi) = wm(i)*auxjno2(ind+1) +
-c     $           wp(i)*auxjno2(ind)
-c         enddo
-c      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,klev
-         auxcolinp(klev-i+1) = co2colx(i) +  n2colx(i) +
-     $         + cocolx(i) 
-      end do
-
-      do indexint=17,24
-
-c     Interpolations
-
-         do i=1,nz2
-            auxi = nz2-i+1
-            !CO2 tabulated coefficient
-            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
-            !O2 tabulated coefficient
-c            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
-c         if(chemthermod.ge.2) then
-c            do i=1,nz2
-c               auxi = nz2-i+1
-c               !NO tabulated coefficient
-c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-c               !NO2 tabulated coefficient
-c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-c            enddo
-c         endif
-
-         call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-         !Correction to include T variation of CO2 cross section
-         if(indexint.eq.24) then
-            do i=1,klev
-               auxi = klev-i+1
-               if(sigma(indexint,auxi)*
-     $              alfa(indexint,auxi)*coltemp(auxi)
-     $              .lt.60.) then
-                  cortemp(i)=exp(-sigma(indexint,auxi)*
-     $                alfa(indexint,auxi)*coltemp(auxi))
-               else 
-                  cortemp(i)=0.
-               end if
-            enddo
-         else
-            do i=1,klev
-               cortemp(i)=1.
-            enddo
-         end if
-         do i=1,klev           
-            ind=auxind(i)
-            auxi = klev-i+1
-            !O2 interpolated coefficient
-c            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
-c     $           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
-c         if(chemthermod.ge.2) then
-c            do i=1,klev
-c               ind=auxind(i)
-c               auxi = klev-i+1
-c               !NO interpolated coefficient
-c               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
-c     $              wp(i)*auxjno(ind)) * cortemp(i)
-c               !NO2 interpolated coefficient
-c               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
-c     $              wp(i)*auxjno2(ind)) * cortemp(i)
-c            enddo
-c         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,klev
-c         auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
-c     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
-         auxcolinp(klev-i+1) = co2colx(i)  + cocolx(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
-c            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
-            !H2O tabulated coefficient
-c            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
-            !H2O2 tabulated coefficient
-c            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
-c         if(chemthermod.ge.2) then
-c            do i=1,nz2
-c               auxi = nz2-i+1
-c               !NO tabulated coefficient
-c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-c               !NO2 tabulated coefficient
-c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-c            enddo
-c         endif
-         call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-         do i=1,klev
-            ind=auxind(i)
-            auxi = klev-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
-c            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
-c     $           wp(i)*auxjo2(ind)) * cortemp(i)
-            !H2O interpolated coefficient
-c            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
-c     $           wp(i)*auxjh2o(ind)) * cortemp(i)
-            !H2O2 interpolated coefficient
-c            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
-c     $           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
-c         if(chemthermod.ge.2) then
-c            do i=1,klev
-c               ind=auxind(i)
-c               auxi = klev-i+1
-c               !NO interpolated coefficient
-c               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
-c     $              wp(i)*auxjno(ind)) * cortemp(i)
-c               !NO2 interpolated coefficient
-c               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
-c     $              wp(i)*auxjno2(ind)) * cortemp(i)
-c            enddo
-c         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,klev
-         auxcolinp(klev-i+1) = co2colx(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
-c            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
-            !H2O tabulated coefficient
-c            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
-            !H2O2 tabulated coefficient
-c            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)            
-            !Tabulated column
-c            auxcoltab(i) = c30_31(auxi)
-         enddo
-         !Only if chemthermod.ge.2
-c         if(chemthermod.ge.2) then
-c            do i=1,nz2
-c               auxi = nz2-i+1
-c               !NO tabulated coefficient
-c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-c               !NO2 tabulated coefficient
-c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-c            enddo
-c         endif
-
-         call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-         do i=1,klev
-            ind=auxind(i)
-            auxi = klev-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
-c            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
-c     $            wp(i)*auxjo2(ind)) * cortemp(i)
-            !H2O interpolated coefficient
-c            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
-c     $            wp(i)*auxjh2o(ind)) * cortemp(i)
-            !H2O2 interpolated coefficient
-c            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
-c     $            wp(i)*auxjh2o2(ind)) * cortemp(i)            
-         enddo
-c         !Only if chemthermod.ge.2
-c         if(chemthermod.ge.2) then
-c            do i=1,klev 
-c               ind=auxind(i)
-c               auxi = klev-i+1
-c               !NO interpolated coefficient
-c               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
-c     $              wp(i)*auxjno(ind)) * cortemp(i)
-c               !NO2 interpolated coefficient
-c               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
-c     $              wp(i)*auxjno2(ind)) * cortemp(i)
-c            enddo
-c         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,klev
-         auxcolinp(klev-i+1) =co2colx(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
-c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
-         !H2O2 tabulated coefficient
-c         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)         
-         !Tabulated column
-         auxcoltab(i) = c32(auxi)
-      enddo
-      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,nz2
-c            auxi = nz2-i+1
-c            !NO tabulated coefficient
-c            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
-c            !NO2 tabulated coefficient
-c            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
-c         enddo
-c      endif
-      call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
-         ind=auxind(i)
-         auxi = klev-i+1
-         !Correction to include T variation of CO2 cross section
-         if(sigma(indexint,klev-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
-c         jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
-c     $        wp(i)*auxjo2(ind)) * cortemp(i)
-         !H2O2 interpolated coefficient
-c         jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
-c     $        wp(i)*auxjh2o2(ind)) * cortemp(i)         
-      enddo
-      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,klev
-c            auxi = klev-i+1
-c            ind=auxind(i)
-c            !NO interpolated coefficient
-c            jfotsout(indexint,10,auxi) = (wm(i)*auxjno(ind+1) +
-c     $           wp(i)*auxjno(ind)) * cortemp(i)
-           !NO2 interpolated coefficient
-c            jfotsout(indexint,13,auxi) = (wm(i)*auxjno2(ind+1) +
-c     $           wp(i)*auxjno2(ind)) * cortemp(i)
-c         enddo
-c      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
-
-c      indexint=33
-c      do i=1,klev
-c         auxcolinp(klev-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
-c      end do
-
-c     Interpolation
-
-c      do i=1,nz2
-c         auxi = nz2-i+1
-         !O2 tabulated coefficient
-c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
-         !H2O2 tabulated coefficient
-c         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
-         !Tabulated column
-c        auxcoltab(i) = c33(auxi)
-c     enddo
-      !Only if chemthermod.ge.2
-c     if(chemthermod.ge.2) then
-c        do i=1,nz2
-            !NO2 tabulated coefficient
-c           auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
-c        enddo
-c     endif
-c     call interfast
-c    $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-c     do i=1,klev
-c        ind=auxind(i)
-c        auxi = klev-i+1
-         !O2 interpolated coefficient
-c        jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
-c    $        wp(i)*auxjo2(ind)
-c        !H2O2 interpolated coefficient
-c        jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
-c    $        wp(i)*auxjh2o2(ind)         
-c     enddo
-      !Only if chemthermod.ge.2
-c     if(chemthermod.ge.2) then
-c        do i=1,klev
-c           ind=auxind(i)
-c           !NO2 interpolated coefficient
-c           jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
-c    $           wp(i)*auxjno2(ind)
-c        enddo
-c     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
-
-c     indexint=34
-c     do i=1,klev
-c        auxcolinp(klev-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
-c    $        no2colx(i)
-c     end do
-
-c     Interpolation
-
-c     do i=1,nz2
-c        auxi = nz2-i+1
-         !O2 tabulated coefficient
-c        auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
-         !H2O2 tabulated coefficient
-c        auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
-         !O3 tabulated coefficient
-c        auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
-         !Tabulated column
-c        auxcoltab(i) = c34(nz2-i+1)
-c     enddo
-      !Only if chemthermod.ge.2
-c     if(chemthermod.ge.2) then
-c        do i=1,nz2
-            !NO2 tabulated coefficient
-c           auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
-c        enddo
-c     endif
-c     call interfast
-c    $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-c     do i=1,klev
-c        ind=auxind(i)
-c        auxi = klev-i+1
-         !O2 interpolated coefficient
-c        jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
-c    $        wp(i)*auxjo2(ind)
-         !H2O2 interpolated coefficient
-c        jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
-c    $        wp(i)*auxjh2o2(ind)
-         !O3 interpolated coefficient
-c        jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
-c    $        wp(i)*auxjo3(ind)         
-c     enddo
-      !Only if chemthermod.ge.2
-c     if(chemthermod.ge.2) then
-c        do i=1,klev
-c           ind=auxind(i)
-            !NO2 interpolated coefficient
-c           jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
-c    $           wp(i)*auxjno2(ind)
-c        enddo
-c     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
-c      do i=1,klev
-c         auxcolinp(klev-i+1) = o3colx(i) 
-c      end do
-c
-c     Interpolation
-
-c      do i=1,nz2
-c         auxi = nz2-i+1
-         !H2O2 tabulated coefficient
-c         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
-         !O3 tabulated coefficient
-c         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
-         !Tabulated column
-c         auxcoltab(i) = c35(auxi)
-c      enddo
-      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,nz2
-c            !NO2 tabulated coefficient
-c            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
-c         enddo
-c      endif
-c      call interfast
-c     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-c      do i=1,klev
-c         ind=auxind(i)
-c         auxi = klev-i+1
-c         !H2O2 interpolated coefficient
-c         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
-c     $        wp(i)*auxjh2o2(ind)
-c         !O3 interpolated coefficient
-c         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
-c     $        wp(i)*auxjo3(ind)         
-c      enddo
-c      if(chemthermod.ge.2) then
-c         do i=1,klev
-c            ind=auxind(i)
-c            !NO2 interpolated coefficient
-c            jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
-c     $           wp(i)*auxjno2(ind)
-c         enddo
-c      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
-c      do i=1,klev
-c         auxcolinp(klev-i+1) = o3colx(i) 
-c      end do
-
-c     Interpolation
-
-c      do i=1,nz2
-c         auxi = nz2-i+1
-         !O3 tabulated coefficient
-c         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
-         !Tabulated column
-c         auxcoltab(i) = c36(auxi)
-c      enddo
-      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,nz2
-c            !NO2 tabulated coefficient
-c            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
-c         enddo
-c      endif
-c      call interfast
-c     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-c      do i=1,klev
-c         ind=auxind(i)
-c         !O3 interpolated coefficient
-c         jfotsout(indexint,7,klev-i+1) = wm(i)*auxjo3(ind+1) +
-c     $        wp(i)*auxjo3(ind)         
-c      enddo
-c      !Only if chemthermod.ge.2
-c      if(chemthermod.ge.2) then
-c         do i=1,klev
-c            ind=auxind(i)
-c           !NO2 interpolated coefficient
-c            jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
-c     $           wp(i)*auxjno2(ind)
-c         enddo
-c      endif
-
-c     End of interval 36
-
-c     End of interpolation to obtain photoabsorption rates
-
-	
-      return
-      end
-
-
-
-c**********************************************************************
-c**********************************************************************
-
-      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
-     $     co2colx,o3pcolx, n2colx, cocolx)
-
-c     mar 2014        gg            adapted to Venus GCM
-c     nov 2002        fgg           first version
-
-c**********************************************************************
-      use dimphy
-      use conc
-      implicit none
-
-
-c     common variables and constants
-c#include "tracer.h"
-#include "param.h"
-#include "param_v4.h"
-#include "clesphys.h"
-#include "mmol.h"
-
-
-c    local parameters and variables
-
-c     input and output variables
-
-      integer    ig
-      integer    chemthermod
-      integer    nesptherm                   !# of species undergoing chemistry, input
-      real       rm(klev,nesptherm)         !densities (cm-3), input
-      real       tx(klev)                   !temperature profile, input
-      real       iz(klev+1)                 !height profile, input
-      real       zenit                      !SZA, input
-      real       co2colx(klev)              !column density of CO2 (cm^-2), output
-      real       o3pcolx(klev)              !column density of O(3P)(cm^-2), output
-      real       n2colx(klev)               !N2 column density (cm-2), output
-      real       cocolx(klev)               !CO column density (cm-2), output
-
-c      real       o2colx(klev)               !column density of O2(cm^-2), output
-c      real       h2colx(klev)               !H2 column density (cm-2), output
-c      real       h2ocolx(klev)              !H2O column density (cm-2), output
-c      real       h2o2colx(klev)             !column density of H2O2(cm^-2), output
-c      real       o3colx(klev)               !O3 column density (cm-2), output
-c      real       nocolx(klev)               !NO column density (cm-2), output
-c      real       hcolx(klev)                !H column density (cm-2), output
-c      real       no2colx(klev)              !NO2 column density (cm-2), output
-      
-
-c     local variables
-
-      real       xx
-      real       grav(klev)
-      real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
-      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
-
-      real       co2x(klev)
-      real       o3px(klev)
-      real       cox(klev)
-      real       n2x(klev)
-      real       nx(klev)
-
-c      real       o2x(klev)
-c      real       o3x(klev)
-c      real       hx(klev)
-c      real       h2x(klev)
-c      real       h2ox(klev)
-c      real       h2o2x(klev)
-c      real       nox(klev)
-c      real       no2x(klev)
-
-      integer    i,j,k,icol,indexint          !indexes
-
-c     variables for optical path calculation
-
-      integer    nz3
-!      parameter  (nz3=nz*2)
-
-      integer    jj
-      real*8      esp(klev*2)
-      real*8      ilayesp(klev*2)
-      real*8      szalayesp(klev*2)
-      integer     nlayesp
-      real*8      zmini
-      real*8      depth
-      real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
-      real*8      espn2,espn,espno,espco,esph,espno2
-      real*8      rcmnz, rcmmini
-      real*8      szadeg
-
-! Tracer indexes in the thermospheric chemistry:
-!!! ATTENTION. These values have to be identical to those in euvheat.F90
-!!! If the values are changed there, the same has to be done here  !!!
-
-      integer,parameter :: ix_co2=1
-      integer,parameter :: ix_n2=13
-      integer,parameter :: ix_o=3
-      integer,parameter :: ix_co=4
-
-c*************************PROGRAM STARTS*******************************
-
-      nz3 = klev*2
-      do i=1,klev
-        xx = ( radio + iz(i) ) * 1.e5    ! conversion [km] ---> [cm]
-         grav(i) = gg * masa /(xx**2)    ! [cm/s2]
-      end do
-
-      !Scale heights  H = kT /Mg  -->  [cm]
-      xx = kboltzman * tx(klev) * n_avog / grav(klev)   ! g cm mol-1
-
-      Ho3p  = xx / mmolo
-      Hco2  = xx / mmolco2
-      Hco   = xx / mmolco
-      Hn2   = xx / mmoln2
-      Hn    = xx / mmoln
-
-      !Only if O3 chem. required
-c      if(chemthermod.ge.1) 
-!     $     Ho3   = xx / mmol(igcm_o3)
-c     $     Ho3   = xx / mmolo3
-c      !Only if N or ion chem.
-c      if(chemthermod.ge.2) then
-c         Hn2   = xx / mmoln2
-c         Hn    = xx / mmoln
-c         Hno   = xx / mmolno
-c         Hno2  = xx / mmolno2
-c      endif
-      ! first loop in altitude : initialisation
-      do i=klev,1,-1
-         !Column initialisation
-         co2colx(i)  = 0.
-         o3pcolx(i)  = 0.
-         n2colx(i)   = 0.
-         cocolx(i)   = 0.
-
-         !--Densities [cm-3]
-         co2x(i)  = rm(i,ix_co2)
-         o3px(i)  = rm(i,ix_o)
-         cox(i)   = rm(i,ix_co)
-         n2x(i)   = rm(i,ix_n2)
-
-c         write(*,*), '--jthermcalc--', co2x(i)
-
-         !Only if O3 chem. required
-c         if(chemthermod.ge.1) 
-c     $        o3x(i)   = rm(i,i_o3)
-c         !Only if Nitrogen of ion chemistry requested
-c         if(chemthermod.ge.2) then
-c            n2x(i)   = rm(i,i_n2)
-c            nx(i)    = rm(i,i_n)
-c            nox(i)   = rm(i,i_no)
-c            no2x(i)  = rm(i,i_no2)
-c         endif
-      enddo     ! end first loop
-      ! second loop in altitude : column calculations
-      do i=klev,1,-1
-         !Routine to calculate the geometrical length of each layer
-         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
-     $         szalayesp,nlayesp, zmini)
-         if(ilayesp(nlayesp).eq.-1) then
-            co2colx(i)=1.e25
-            o3pcolx(i)=1.e25
-            n2colx(i)=1.e25
-            cocolx(i)=1.e25
-
-c            o2colx(i)=1.e25
-c            o3pcolx(i)=1.e25
-c            h2colx(i)=1.e25
-c            h2ocolx(i)=1.e25
-c            h2o2colx(i)=1.e25
-c            o3colx(i)=1.e25
-c            ncolx(i)=1.e25
-c            nocolx(i)=1.e25
-c            cocolx(i)=1.e25
-c            hcolx(i)=1.e25
-c            no2colx(i)=1.e25
-         else
-            rcmnz = ( radio + iz(klev) ) * 1.e5    ! km --> cm
-            rcmmini = ( radio + zmini ) * 1.e5
-            !Column calculation taking into account the geometrical depth 
-            !calculated before
-            do j=1,nlayesp
-               jj=ilayesp(j)
-               !Top layer
-               if(jj.eq.klev) then
-                  if(zenit.le.60.) then 
-                     o3pcolx(i)=o3pcolx(i)+o3px(klev)*Ho3p*esp(j)
-     $                    *1.e-5
-                     co2colx(i)=co2colx(i)+co2x(klev)*Hco2*esp(j)
-     $                    *1.e-5
-                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
-     $                    *1.e-5
-                     n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
-     $                    *1.e-5
-
-c                     h2o2colx(i)=h2o2colx(i)+
-c     $                    h2o2x(klev)*Hh2o2*esp(j)*1.e-5
-c                     o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)
-c     $                    *1.e-5
-c                     h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)
-c     $                    *1.e-5
-c                     h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)
-c     $                    *1.e-5                     
-c                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
-c     $                    *1.e-5
-c                     hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)
-c     $                    *1.e-5
-                     !Only if O3 chemistry required
-c                     if(chemthermod.ge.1) o3colx(i)=
-c     $                    o3colx(i)+o3x(klev)*Ho3*esp(j)
-c     $                    *1.e-5
-                     !Only if N or ion chemistry requested
-c                     if(chemthermod.ge.2) then
-c                        n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
-c     $                    *1.e-5
-
-c                     endif
-                  else if(zenit.gt.60.) then
-                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
-                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
-                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
-                     espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
-                     espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
-
-c                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
-c                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
-c                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
-c                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
-c                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
-                     !Only if O3 chemistry required
-c                     if(chemthermod.ge.1)                     
-c     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
-c                     !Only if N or ion chemistry requested
-c                     if(chemthermod.ge.2) then
-c                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
-c                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
-c                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
-c                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
-c                     endif
-
-                     co2colx(i) = co2colx(i) + espco2*co2x(klev)
-                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev)
-                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
-                     n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
-
-c                     o2colx(i)  = o2colx(i)  + espo2*o2x(klev)
-c                     h2colx(i)  = h2colx(i)  + esph2*h2x(klev)
-c                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev)
-c                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev)
-c                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
-c                     hcolx(i)   = hcolx(i)   + esph*hx(klev)
-                     !Only if O3 chemistry required
-c                     if(chemthermod.ge.1)                      
-c     $                  o3colx(i) = o3colx(i)  + espo3*o3x(klev)
-c                     !Only if N or ion chemistry requested
-c                     if(chemthermod.ge.2) then
-c                        n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
-c                        ncolx(i)   = ncolx(i)   + espn*nx(klev)
-c                        nocolx(i)  = nocolx(i)  + espno*nox(klev)
-c                        no2colx(i) = no2colx(i) + espno2*no2x(klev)
-c                     endif
-                  endif !Of if zenit.lt.60
-               !Other layers
-               else 
-                  co2colx(i)  = co2colx(i) + 
-     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
-                  o3pcolx(i)  = o3pcolx(i) + 
-     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
-                  cocolx(i)   = cocolx(i) + 
-     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
-                  n2colx(i)   = n2colx(i) + 
-     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
-
-c
-c                  o2colx(i)   = o2colx(i) + 
-c     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
-c                  h2colx(i)   = h2colx(i) + 
-c     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
-c                  h2ocolx(i)  = h2ocolx(i) + 
-c     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
-c                  h2o2colx(i) = h2o2colx(i) + 
-c     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
-c                  hcolx(i)    = hcolx(i) + 
-c     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
-                  !Only if O3 chemistry required
-c                  if(chemthermod.ge.1) 
-c     $                 o3colx(i) = o3colx(i) + 
-c     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
-c                  !Only if N or ion chemistry requested
-c                  if(chemthermod.ge.2) then
-c                     n2colx(i)   = n2colx(i) + 
-c     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
-c                     ncolx(i)    = ncolx(i) + 
-c     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
-c                     nocolx(i)   = nocolx(i) + 
-c     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
-c                     no2colx(i)  = no2colx(i) + 
-c     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
-c                  endif
-
-               endif  !Of if jj.eq.klev
-            end do    !Of do j=1,nlayesp
-         endif        !Of ilayesp(nlayesp).eq.-1
-      enddo           !Of do i=klev,1,-1
-      return
-
-
-      end
-
-
-c**********************************************************************
-c**********************************************************************
-
-      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
-C
-C subroutine to perform linear interpolation in pressure from 1D profile 
-C escin(nl) sampled on pressure grid pin(nl) to profile
-C escout(nlayer) on pressure grid p(nlayer).
-C
-      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
-            wm(n1) = 0.d0
-            wp(n1) = 0.d0
-         else
-            do n = nini,nl-1
-               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
-                  nm(n1)=n
-                  np=n+1
-                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
-                  wp(n1)=1.d0 - wm(n1)
-                  nini = n
-                  exit
-               endif
-            enddo
-         endif
-      enddo
-      return
-      end
-
-
-c**********************************************************************
-c**********************************************************************
-
-      subroutine espesor_optico_A (ig,capa, szadeg,z,
-     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
-
-c     fgg              nov 03      Adaptation to Martian model
-c     malv             jul 03      Corrected z grid. Split in alt & frec codes
-c     fgg              feb 03      first version
-*************************************************************************
-      use dimphy
-      implicit none
-
-
-c     common variables and constants
-#include "param.h"
-#include "param_v4.h"
-
-c     arguments
-
-      real        szadeg                ! I. SZA [rad]
-      real        z                     ! I. altitude of interest [km]
-      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
-                                        !  (=2*klev= max# of layers in ray path)
-      real     iz(klev+1)              ! I. Altitude of each layer
-      real*8        esp(nz3)            ! O. layer widths after geometrically 
-                                        !    amplified; in [cm] except at TOA
-                                        !    where an auxiliary value is used
-      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
-      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
-      integer       nlayesp
-!      real*8        nlayesp             ! O. # layers along ray path at this z
-      real*8        zmini               ! O. Minimum altitud of ray path [km]
-
-
-c     local variables and constants
-
-        integer     j,i,capa
-        integer     jmin                  ! index of min.altitude along ray path
-        real*8      szarad                ! SZA [deg]
-        real*8      zz
-        real*8      diz(klev+1)
-        real*8      rkmnz                 ! distance TOA to center of Planet [km]
-        real*8      rkmmini               ! distance zmini to center of P [km] 
-        real*8      rkmj                  ! intermediate distance to C of P [km]
-
-c external function
-        external  grid_R8          ! Returns index of layer containing the altitude
-                                ! of interest, z; for example, if 
-                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i 
-        integer   grid_R8
-
-*************************************************************************     
-        szarad = dble(szadeg)*3.141592d0/180.d0
-        zz=dble(z)
-        do i=1,klev
-           diz(i)=dble(iz(i))
-        enddo
-        do j=1,nz3 
-          esp(j) = 0.d0
-          szalayesp(j) = 777.d0
-          ilayesp(j) = 0
-        enddo
-        nlayesp = 0
-
-        ! First case: szadeg<60
-        ! The optical thickness will be given by  1/cos(sza)
-        ! We deal with 2 different regions:
-        !   1: First, all layers between z and ztop ("upper part of ray")
-        !   2: Second, the layer at ztop
-        if(szadeg.lt.60.d0) then
-
-           zmini = zz
-           if(abs(zz-diz(klev)).lt.1.d-3) goto 1357
-           ! 1st Zone: Upper part of ray
-           !
-           do j=grid_R8(zz,diz,klev),klev-1
-             nlayesp = nlayesp + 1 
-             ilayesp(nlayesp) = j
-             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
-             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
-             szalayesp(nlayesp) = szadeg
-           end do
-
-           ! 
-           ! 2nd Zone: Top layer
- 1357      continue
-           nlayesp = nlayesp + 1 
-           ilayesp(nlayesp) = klev
-           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
-           szalayesp(nlayesp) = szadeg
-
-
-        ! Second case:  60 < szadeg < 90
-        ! The optical thickness is evaluated.
-        !    (the magnitude of the effect of not using cos(sza) is 3.e-5 
-        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
-        ! We deal with 2 different regions:
-        !   1: First, all layers between z and ztop ("upper part of ray")
-        !   2: Second, the layer at ztop ("uppermost layer")
-        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
-
-           zmini=(radio+zz)*sin(szarad)-radio
-           rkmmini = radio + zmini
-
-           if(abs(zz-diz(klev)).lt.1.d-4) goto 1470
-
-           ! 1st Zone: Upper part of ray
-           !
-           do j=grid_R8(zz,diz,klev),klev-1
-              nlayesp = nlayesp + 1 
-              ilayesp(nlayesp) = j
-              esp(nlayesp) = 
-     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
-     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
-              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
-              rkmj = radio+diz(j)
-              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
-              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
-           end do
- 1470      continue
-           ! 2nd Zone:  Uppermost layer of ray.
-           !
-           nlayesp = nlayesp + 1 
-           ilayesp(nlayesp) = klev
-           rkmnz = radio+diz(klev)
-           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
-           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
-           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
-           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
-
-
-        ! Third case:  szadeg > 90
-        ! The optical thickness is evaluated.
-        ! We deal with 5 different regions:
-        !   1: all layers between z and ztop ("upper part of ray")
-        !   2: the layer at ztop ("uppermost layer")
-        !   3: the lowest layer, at zmini
-        !   4: the layers increasing from zmini to z (here SZA<90)
-        !   5: the layers decreasing from z to zmini (here SZA>90)
-        else if(szadeg.gt.90.d0) then
-
-           zmini=(radio+zz)*sin(szarad)-radio
-           rkmmini = radio + zmini
-
-           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
-             nlayesp = nlayesp + 1 
-             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
-!             esp(nlayesp) = 1.e30
-
-           else
-              jmin=grid_R8(zmini,diz,klev)+1
-              
-
-              if(abs(zz-diz(klev)).lt.1.d-4) goto 9876
-
-              ! 1st Zone: Upper part of ray
-              !
-              do j=grid_R8(zz,diz,klev),klev-1
-                nlayesp = nlayesp + 1 
-                ilayesp(nlayesp) = j
-                esp(nlayesp) = 
-     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
-     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
-                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
-                rkmj = radio+diz(j)
-                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
-                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
-              end do
-
- 9876         continue
-              ! 2nd Zone:  Uppermost layer of ray.
-              !
-              nlayesp = nlayesp + 1 
-              ilayesp(nlayesp) = klev
-              rkmnz = radio+diz(klev)
-              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
-              esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
-              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
-              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
-
-              ! 3er Zone: Lowestmost layer of ray
-              !
-              if ( jmin .ge. 2 ) then      ! If above the planet's surface
-                j=jmin-1
-                nlayesp = nlayesp + 1 
-                ilayesp(nlayesp) = j
-                esp(nlayesp) = 2. * 
-     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
-                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
-                rkmj = radio+diz(j+1)
-                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
-                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
-              endif
-
-              ! 4th zone: Lower part of ray, increasing from zmin to z
-              !    ( layers with SZA < 90 deg )
-              do j=jmin,grid_R8(zz,diz,klev)-1
-                nlayesp = nlayesp + 1 
-                ilayesp(nlayesp) = j
-                esp(nlayesp) = 
-     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
-     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       !  [km]
-                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
-                rkmj = radio+diz(j)
-                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
-                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
-              end do
-
-              ! 5th zone: Lower part of ray, decreasing from z to zmin
-              !    ( layers with SZA > 90 deg )
-              do j=grid_R8(zz,diz,klev)-1, jmin, -1
-                nlayesp = nlayesp + 1 
-                ilayesp(nlayesp) = j
-                esp(nlayesp) = 
-     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
-     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
-                esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
-                rkmj = radio+diz(j)
-                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
-                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
-              end do
-
-           end if
-
-        end if
-
-        return
-
-        end
-
-
-c**********************************************************************
-c***********************************************************************
-
-        function grid_R8 (z, zgrid, nz)
-
-c Returns the index where z is located within vector zgrid
-c The vector zgrid must be monotonously increasing, otherwise program stops.
-c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops. 
-c
-c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le. 
-c MALV    Jul-2003
-c***********************************************************************
-
-        implicit none
-
-c Arguments 
-        integer   nz
-        real*8      z
-        real*8      zgrid(nz)
-        integer   grid_R8
-
-c Local  
-        integer   i, nz1, nznew
-
-c*** CODE START 
-
-        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then 
-           write (*,*) ' GRID/ z outside bounds of zgrid '
-           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
-           stop ' Serious error in GRID.F '
-        endif
-        if ( nz .lt. 2 ) then 
-           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
-           stop ' Serious error in GRID.F '
-        endif
-        if ( zgrid(1) .ge. zgrid(nz) ) then 
-           write (*,*) ' GRID/ zgrid must increase with index'
-           stop ' Serious error in GRID.F '
-        endif
-
-        nz1 = 1
-        nznew = nz/2
-        if ( z .gt. zgrid(nznew) ) then
-           nz1 = nznew
-           nznew = nz
-        endif
-        do i=nz1+1,nznew
-           if ( z. eq. zgrid(i) ) then
-              grid_R8=i
-              return
-              elseif ( z .le. zgrid(i) ) then
-              grid_R8 = i-1
-              return
-           endif
-        enddo
-        grid_R8 = nz
-        return
-
-        end
-
-
-
-!c***************************************************
-!c***************************************************
-
-      subroutine flujo(date)
-
-
-!c     fgg           nov 2002     first version
-!c***************************************************
-      use dimphy
-      use conc
-      implicit none
-
-
-!     common variables and constants
-      include "param.h"
-      include 'param_v4.h'
-      include "clesphys.h"
-
-!     Arguments
-
-       real date
-c       integer, parameter :: dateyr = 2006
-
-!     Local variable and constants
-!    dist_sol : distance venus - soleil
-
-      real, parameter :: dist_sol=0.72333
-      integer i
-      integer inter
-      real    nada
-
-!c*************************************************
-
-      if(date.lt.1985.) date=1985.
-      if(date.gt.2001.) date=2001.
-      
-      do i=1,ninter
-         fluxtop(i)=1.
-         !Variation of solar flux with 11 years solar cycle
-         !For more details, see Gonzalez-Galindo et al. 2005
-         !To be improved in next versions
-        if(i.le.24 .and.solvarmod.eq.0) then
-          fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                  
-     $        *sin(2.*3.1416/11.*(date-1985.-3.1416))          
-     $        +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
-          
-        end if
- !     The solar flux calculated  
- !      is  corrected for
- !     the actual Venus-Sun dist
-        fluxtop(i)=fluxtop(i)*(1./dist_sol)**2
-
-
-       end do
-     
-      return
-      end
Index: trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F	(revision 2464)
@@ -2,5 +2,5 @@
 
       subroutine jthermcalc_e107
-     $     (ig,chemthermod,rm,nesptherm,tx,iz,zenit,zday)
+     $     (ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit)
 
 
@@ -11,44 +11,50 @@
 c MAC July 2003
 c**********************************************************************
-	use dimphy
-	use conc
+
+      use param_v4_h, only: jfotsout,crscabsi2,
+     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
+     .    co2crsc195,co2crsc295,t0,
+     .    jabsifotsintpar,ninter,nz2,
+     .    nabs,e107,date_e107,e107_tab,
+     .    coefit0,coefit1,coefit2,coefit3,coefit4
+
       implicit none
 
-c     common variables and constants
-#include "param.h"
-#include "param_v4.h"
+      include "clesphys.h"
+
 
 c     input and output variables
 
-      integer    ig
-      integer    chemthermod
-      integer    nesptherm                      !Number of species considered
-      real       rm(klev,nesptherm)         !Densities (cm-3)
-      real       tx(klev)                   !temperature
-      real       zenit                          !SZA
-      real       iz(klev)                   !Local altitude
-      real       zday                           !Martian day after Ls=0
-
+      integer,intent(in) :: ig,nlayer
+      integer,intent(in) :: chemthermod
+      integer,intent(in) :: nesptherm          !Number of species considered
+      real,intent(in) :: rm(nlayer,nesptherm)  !Densities (cm-3)
+      real,intent(in) :: tx(nlayer)            !temperature
+      real,intent(in) :: zenit         !SZA
+      real,intent(in) :: iz(nlayer)    !Local altitude
+
+! NB: no output variable! (computed jfotsout() is in module param_v4_h)
 
 c    local parameters and variables
 
-      real       co2colx(klev)              !column density of CO2 (cm^-2)
-      real       o2colx(klev)               !column density of O2(cm^-2)
-      real       o3pcolx(klev)              !column density of O(3P)(cm^-2)
-      real       h2colx(klev)               !H2 column density (cm-2)
-      real       h2ocolx(klev)              !H2O column density (cm-2)
-      real       h2o2colx(klev)             !column density of H2O2(cm^-2)
-      real       o3colx(klev)               !O3 column density (cm-2)
-      real       n2colx(klev)               !N2 column density (cm-2)
-      real       ncolx(klev)                !N column density (cm-2)
-      real       nocolx(klev)               !NO column density (cm-2)
-      real       cocolx(klev)               !CO column density (cm-2)
-      real       hcolx(klev)                !H column density (cm-2)
-      real       no2colx(klev)              !NO2 column density (cm-2)
-      real       t2(klev)
-      real       coltemp(klev)
-      real       sigma(ninter,klev)
-      real       alfa(ninter,klev)
-      real       realday
+      real, parameter :: dist_sol=0.72333
+
+      real       co2colx(nlayer)              !column density of CO2 (cm^-2)
+      real       o2colx(nlayer)               !column density of O2(cm^-2)
+      real       o3pcolx(nlayer)              !column density of O(3P)(cm^-2)
+      real       h2colx(nlayer)               !H2 column density (cm-2)
+      real       h2ocolx(nlayer)              !H2O column density (cm-2)
+      real       h2o2colx(nlayer)             !column density of H2O2(cm^-2)
+      real       o3colx(nlayer)               !O3 column density (cm-2)
+      real       n2colx(nlayer)               !N2 column density (cm-2)
+      real       ncolx(nlayer)                !N column density (cm-2)
+      real       nocolx(nlayer)               !NO column density (cm-2)
+      real       cocolx(nlayer)               !CO column density (cm-2)
+      real       hcolx(nlayer)                !H column density (cm-2)
+      real       no2colx(nlayer)              !NO2 column density (cm-2)
+      real       t2(nlayer)
+      real       coltemp(nlayer)
+      real       sigma(ninter,nlayer)
+      real       alfa(ninter,nlayer)
       
       integer    i,j,k,indexint                 !indexes
@@ -74,21 +80,25 @@
       real*8      auxjh(nz2)
       real*8      auxjno2(nz2)
-      real*8      wp(klev),wm(klev)
-      real*8      auxcolinp(klev)
-      integer     auxind(klev)
+      real*8      wp(nlayer),wm(nlayer)
+      real*8      auxcolinp(nlayer)
+      integer     auxind(nlayer)
       integer     auxi
       integer     ind
-      real*8      cortemp(klev)
+      real*8      cortemp(nlayer)
 
       real*8      limdown                      !limits for interpolation
       real*8      limup                        !        ""
 
-      !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
+      !!!ATTENTION. Here i_co2 has to have the same value than in euvheat.F90
       !!!If the value is changed there, if has to be changed also here !!!!
       integer,parameter :: i_co2=1
 
+      character*20 modname
+      character*80 abort_message
 
 c*************************PROGRAM STARTS*******************************
       
+      modname = 'jthermcalc_e107'
+
       if(zenit.gt.140.) then
          dn='n'
@@ -105,5 +115,5 @@
       !Auxiliar temperature to take into account the temperature dependence
       !of CO2 cross section
-      do i=1,klev
+      do i=1,nlayer
          t2(i)=tx(i)
          if(t2(i).lt.195.0) t2(i)=195.0
@@ -113,11 +123,10 @@
       !Calculation of column amounts 
       call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
-     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
-     $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
+     $     co2colx,o3pcolx,n2colx,cocolx)
 
       !Auxiliar column to include the temperature dependence 
       !of CO2 cross section
-      coltemp(klev)=co2colx(klev)*abs(t2(klev)-t0(klev))
-      do i=klev-1,1,-1
+      coltemp(nlayer)=co2colx(nlayer)*abs(t2(nlayer)-t0(nlayer))
+      do i=nlayer-1,1,-1
         coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR 
      $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5 
@@ -126,5 +135,5 @@
       
       !Calculation of CO2 cross section at temperature t0(i)
-      do i=1,klev
+      do i=1,nlayer
          do indexint=24,32
            sigma(indexint,i)=co2crsc195(indexint-23)
@@ -134,27 +143,15 @@
       end do
 
-      !E10.7 for the day: linear interpolation to tabulated values
-      realday=mod(zday,669.)
-      if(realday.lt.date_e107(1)) then
-         e107=e107_tab(1)
-      else if(realday.ge.date_e107(669)) then
-         e107=e107_tab(669)   
-      else if(realday.ge.date_e107(1).and.
-     $        realday.lt.date_e107(669)) then
-         do i=1,668
-            if(realday.ge.date_e107(i).and.
-     $           realday.lt.date_e107(i+1)) then
-               tinf=i
-               tsup=i+1
-               e107=e107_tab(tinf)+(realday-date_e107(tinf))*
-     $              (e107_tab(tsup)-e107_tab(tinf))
-            endif
-         enddo
-      endif
+      if (solvarmod==0) then
+        e107=fixed_euv_value
+      else
+            abort_message='solvarmod should be 0...'
+            call abort_physic(modname,abort_message,1)
+      endif ! of if (solvarmod==0)
 
       !Photoabsorption coefficients at TOA as a function of E10.7
       do j=1,nabs
          do indexint=1,ninter
-            jfotsout(indexint,j,klev)=coefit0(indexint,j)+
+            jfotsout(indexint,j,nlayer)=coefit0(indexint,j)+
      $           coefit1(indexint,j)*e107+coefit2(indexint,j)*e107**2+
      $           coefit3(indexint,j)*e107**3+coefit4(indexint,j)*e107**4
@@ -178,6 +175,6 @@
 c     Input atmospheric column
       indexint=1
-      do i=1,klev
-         auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint) +
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint) +
      $        o2colx(i)*crscabsi2(2,indexint) + 
      $        o3pcolx(i)*crscabsi2(3,indexint) + 
@@ -213,16 +210,16 @@
 
       call interfast 
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi=klev-i+1
+         auxi=nlayer-i+1
          !CO2 interpolated coefficient
-         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,klev) *
+         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
      $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
      $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
          !O3p interpolated coefficient
-         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,klev) *
+         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
      $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
          !H2 interpolated coefficient
@@ -233,8 +230,8 @@
       !N interpolated coefficient
       if(chemthermod.ge.2) then
-         do i=1,klev
+         do i=1,nlayer
             ind=auxind(i)
-            jfotsout(indexint,9,klev-i+1) =  
-     $           jfotsout(indexint,9,klev) * 
+            jfotsout(indexint,9,nlayer-i+1) =  
+     $           jfotsout(indexint,9,nlayer) * 
      $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
          enddo
@@ -255,6 +252,6 @@
 c     Input atmospheric column
       do indexint=2,15
-         do i=1,klev
-            auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
+         do i=1,nlayer
+            auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
      $           o2colx(i)*crscabsi2(2,indexint)+
      $           o3pcolx(i)*crscabsi2(3,indexint)+
@@ -302,54 +299,54 @@
          endif
 
-          call interfast(wm,wp,auxind,auxcolinp,klev,
+          call interfast(wm,wp,auxind,auxcolinp,nlayer,
      $        auxcoltab,nz2,limdown,limup)
-          do i=1,klev
+          do i=1,nlayer
              ind=auxind(i)
-             auxi = klev-i+1
+             auxi = nlayer-i+1
              !O2 interpolated coefficient
              jfotsout(indexint,2,auxi) = 
-     $            jfotsout(indexint,2,klev) *
+     $            jfotsout(indexint,2,nlayer) *
      $            (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
              !O3p interpolated coefficient
              jfotsout(indexint,3,auxi) = 
-     $            jfotsout(indexint,3,klev) *
+     $            jfotsout(indexint,3,nlayer) *
      $            (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
              !CO2 interpolated coefficient
              jfotsout(indexint,1,auxi) = 
-     $            jfotsout(indexint,1,klev) *
+     $            jfotsout(indexint,1,nlayer) *
      $            (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
              !H2 interpolated coefficient
              jfotsout(indexint,5,auxi) = 
-     $            jfotsout(indexint,5,klev) *
+     $            jfotsout(indexint,5,nlayer) *
      $            (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
              !N2 interpolated coefficient
              jfotsout(indexint,8,auxi) = 
-     $            jfotsout(indexint,8,klev) *
+     $            jfotsout(indexint,8,nlayer) *
      $            (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
              !CO interpolated coefficient
              jfotsout(indexint,11,auxi) = 
-     $            jfotsout(indexint,11,klev) *
+     $            jfotsout(indexint,11,nlayer) *
      $            (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
              !H interpolated coefficient
              jfotsout(indexint,12,auxi) = 
-     $            jfotsout(indexint,12,klev) *
+     $            jfotsout(indexint,12,nlayer) *
      $            (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
           enddo
           !Only if chemthermod.ge.2
           if(chemthermod.ge.2) then
-             do i=1,klev
+             do i=1,nlayer
                 ind=auxind(i)
-                auxi = klev-i+1
+                auxi = nlayer-i+1
                 !N interpolated coefficient
                 jfotsout(indexint,9,auxi) = 
-     $               jfotsout(indexint,9,klev) *
+     $               jfotsout(indexint,9,nlayer) *
      $               (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
                 !NO interpolated coefficient
                 jfotsout(indexint,10,auxi)=
-     $               jfotsout(indexint,10,klev) *
+     $               jfotsout(indexint,10,nlayer) *
      $               (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
                 !NO2 interpolated coefficient
                 jfotsout(indexint,13,auxi)=
-     $               jfotsout(indexint,13,klev) *
+     $               jfotsout(indexint,13,nlayer) *
      $               (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
              enddo
@@ -370,6 +367,6 @@
 c     Input atmospheric column
       indexint=16
-      do i=1,klev
-         auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
      $        o2colx(i)*crscabsi2(2,indexint)+
      $        o3pcolx(i)*crscabsi2(3,indexint)+
@@ -417,45 +414,45 @@
 
       call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = klev-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
      $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
          !CO2 interpolated coefficient
-         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,klev) *
+         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
      $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
          !O3p interpolated coefficient
-         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,klev) *
+         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
      $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
          !N2 interpolated coefficient
-         jfotsout(indexint,8,auxi) = jfotsout(indexint,8,klev) *
+         jfotsout(indexint,8,auxi) = jfotsout(indexint,8,nlayer) *
      $        (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
          !CO interpolated coefficient
          jfotsout(indexint,11,auxi) = 
-     $        jfotsout(indexint,11,klev) *
+     $        jfotsout(indexint,11,nlayer) *
      $        (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
          !H interpolated coefficient
          jfotsout(indexint,12,auxi) = 
-     $        jfotsout(indexint,12,klev) *
+     $        jfotsout(indexint,12,nlayer) *
      $        (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
       enddo
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,klev
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = klev-i+1
+            auxi = nlayer-i+1
             !N interpolated coefficient
             jfotsout(indexint,9,auxi) = 
-     $           jfotsout(indexint,9,klev) *
+     $           jfotsout(indexint,9,nlayer) *
      $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
             !NO interpolated coefficient
             jfotsout(indexint,10,auxi) = 
-     $           jfotsout(indexint,10,klev) * 
+     $           jfotsout(indexint,10,nlayer) * 
      $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
             !NO2 interpolated coefficient
             jfotsout(indexint,13,auxi) = 
-     $           jfotsout(indexint,13,klev) *
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -473,6 +470,6 @@
 c     Input column
 
-      do i=1,klev
-         auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
      $        nocolx(i) + cocolx(i) + no2colx(i)
       end do
@@ -507,9 +504,9 @@
 
          call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
          !Correction to include T variation of CO2 cross section
          if(indexint.eq.24) then
-            do i=1,klev
-               auxi = klev-i+1
+            do i=1,nlayer
+               auxi = nlayer-i+1
                if(sigma(indexint,auxi)*
      $              alfa(indexint,auxi)*coltemp(auxi)
@@ -522,19 +519,19 @@
             enddo
          else
-            do i=1,klev
+            do i=1,nlayer
                cortemp(i)=1.
             enddo
          end if
-         do i=1,klev           
+         do i=1,nlayer           
             ind=auxind(i)
-            auxi = klev-i+1
+            auxi = nlayer-i+1
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = 
-     $           jfotsout(indexint,2,klev) *
+     $           jfotsout(indexint,2,nlayer) *
      $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
      $           cortemp(i)
             !CO2 interpolated coefficient
             jfotsout(indexint,1,auxi) = 
-     $           jfotsout(indexint,1,klev) *
+     $           jfotsout(indexint,1,nlayer) *
      $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) 
      $           * cortemp(i)
@@ -545,10 +542,10 @@
             !N2 interpolated coefficient
             jfotsout(indexint,8,auxi) = 
-     $           jfotsout(indexint,8,klev) *
+     $           jfotsout(indexint,8,nlayer) *
      $           (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind)) * 
      $           cortemp(i)            
             !CO interpolated coefficient
             jfotsout(indexint,11,auxi) = 
-     $           jfotsout(indexint,11,klev) *
+     $           jfotsout(indexint,11,nlayer) *
      $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) * 
      $           cortemp(i)            
@@ -556,15 +553,15 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,klev
+            do i=1,nlayer
                ind=auxind(i)
-               auxi = klev-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=
-     $              jfotsout(indexint,10,klev) *
+     $              jfotsout(indexint,10,nlayer) *
      $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) * 
      $              cortemp(i)
                !NO2 interpolated coefficient
                jfotsout(indexint,13,auxi)=
-     $              jfotsout(indexint,13,klev) *
+     $              jfotsout(indexint,13,nlayer) *
      $              (wm(i)*auxjno2(ind+1)+ wp(i)*auxjno2(ind)) * 
      $              cortemp(i)
@@ -585,6 +582,6 @@
 c     Input atmospheric column
 
-      do i=1,klev
-         auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
      $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
       end do
@@ -620,8 +617,8 @@
          endif
          call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-         do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = klev-i+1
+            auxi = nlayer-i+1
             !Correction to include T variation of CO2 cross section
             if(sigma(indexint,auxi)*alfa(indexint,auxi)*
@@ -634,5 +631,5 @@
             !CO2 interpolated coefficient
             jfotsout(indexint,1,auxi) = 
-     $           jfotsout(indexint,1,klev) *
+     $           jfotsout(indexint,1,nlayer) *
      $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) * 
      $           cortemp(i) *
@@ -641,20 +638,20 @@
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = 
-     $           jfotsout(indexint,2,klev) *
+     $           jfotsout(indexint,2,nlayer) *
      $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
      $           cortemp(i)
             !H2O interpolated coefficient
             jfotsout(indexint,4,auxi) = 
-     $           jfotsout(indexint,4,klev) *
+     $           jfotsout(indexint,4,nlayer) *
      $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) * 
      $           cortemp(i)
             !H2O2 interpolated coefficient
             jfotsout(indexint,6,auxi) = 
-     $           jfotsout(indexint,6,klev) *
+     $           jfotsout(indexint,6,nlayer) *
      $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
      $           cortemp(i)            
             !CO interpolated coefficient
             jfotsout(indexint,11,auxi) = 
-     $           jfotsout(indexint,11,klev) *
+     $           jfotsout(indexint,11,nlayer) *
      $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) * 
      $           cortemp(i)
@@ -662,15 +659,15 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,klev
+            do i=1,nlayer
                ind=auxind(i)
-               auxi = klev-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=
-     $              jfotsout(indexint,10,klev) *
+     $              jfotsout(indexint,10,nlayer) *
      $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) * 
      $              cortemp(i)
                !NO2 interpolated coefficient
                jfotsout(indexint,13,auxi)=
-     $              jfotsout(indexint,13,klev) *
+     $              jfotsout(indexint,13,nlayer) *
      $              (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) * 
      $              cortemp(i)
@@ -693,6 +690,6 @@
 c     Input atmospheric column
 
-      do i=1,klev
-         auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) + 
      $        h2o2colx(i) + nocolx(i) + no2colx(i)
       end do
@@ -727,8 +724,8 @@
 
          call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-         do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = klev-i+1
+            auxi = nlayer-i+1
             !Correction to include T variation of CO2 cross section
             if(sigma(indexint,auxi)*alfa(indexint,auxi)*
@@ -741,5 +738,5 @@
             !CO2 interpolated coefficient
             jfotsout(indexint,1,auxi) = 
-     $           jfotsout(indexint,1,klev) *
+     $           jfotsout(indexint,1,nlayer) *
      $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) * 
      $           cortemp(i) *
@@ -748,15 +745,15 @@
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = 
-     $           jfotsout(indexint,2,klev) *
+     $           jfotsout(indexint,2,nlayer) *
      $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
      $           cortemp(i)
             !H2O interpolated coefficient
             jfotsout(indexint,4,auxi) = 
-     $           jfotsout(indexint,4,klev) *
+     $           jfotsout(indexint,4,nlayer) *
      $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) * 
      $           cortemp(i)
             !H2O2 interpolated coefficient
             jfotsout(indexint,6,auxi) = 
-     $           jfotsout(indexint,6,klev) *
+     $           jfotsout(indexint,6,nlayer) *
      $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
      $           cortemp(i)            
@@ -764,10 +761,10 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,klev 
+            do i=1,nlayer 
                ind=auxind(i)
-               auxi = klev-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=
-     $              jfotsout(indexint,10,klev) *
+     $              jfotsout(indexint,10,nlayer) *
      $              (wm(i)*auxjno(ind+1) +wp(i)*auxjno(ind)) * 
      $              cortemp(i)
@@ -795,6 +792,6 @@
 
       indexint=32
-      do i=1,klev
-         auxcolinp(klev-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) + 
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) + 
      $        nocolx(i) + no2colx(i)
       end do
@@ -824,10 +821,10 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = klev-i+1
+         auxi = nlayer-i+1
          !Correction to include T variation of CO2 cross section
-         if(sigma(indexint,klev-i+1)*alfa(indexint,auxi)*
+         if(sigma(indexint,nlayer-i+1)*alfa(indexint,auxi)*
      $        coltemp(auxi).lt.60.) then
             cortemp(i)=exp(-sigma(indexint,auxi)*
@@ -838,5 +835,5 @@
          !CO2 interpolated coefficient
          jfotsout(indexint,1,auxi) = 
-     $        jfotsout(indexint,1,klev) *
+     $        jfotsout(indexint,1,nlayer) *
      $        (wm(i)*auxjco2(ind+1)+wp(i)*auxjco2(ind)) * 
      $        cortemp(i) *
@@ -845,10 +842,10 @@
          !O2 interpolated coefficient
          jfotsout(indexint,2,auxi) = 
-     $        jfotsout(indexint,2,klev) *
+     $        jfotsout(indexint,2,nlayer) *
      $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) * 
      $        cortemp(i)
          !H2O2 interpolated coefficient
          jfotsout(indexint,6,auxi) = 
-     $        jfotsout(indexint,6,klev) *
+     $        jfotsout(indexint,6,nlayer) *
      $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
      $        cortemp(i)         
@@ -856,15 +853,15 @@
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,klev
-            auxi = klev-i+1
+         do i=1,nlayer
+            auxi = nlayer-i+1
             ind=auxind(i)
             !NO interpolated coefficient
             jfotsout(indexint,10,auxi) = 
-     $           jfotsout(indexint,10,klev) * 
+     $           jfotsout(indexint,10,nlayer) * 
      $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) * 
      $           cortemp(i)
            !NO2 interpolated coefficient
             jfotsout(indexint,13,auxi) = 
-     $           jfotsout(indexint,13,klev) *
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) * 
      $           cortemp(i)
@@ -885,6 +882,6 @@
 
       indexint=33
-      do i=1,klev
-         auxcolinp(klev-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
       end do
 
@@ -908,22 +905,22 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = klev-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
      $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
          !H2O2 interpolated coefficient
-         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,klev) *
+         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
      $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
       enddo
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,klev
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,klev-i+1) = 
-     $           jfotsout(indexint,13,klev) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -943,6 +940,6 @@
 
       indexint=34
-      do i=1,klev
-         auxcolinp(klev-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) + 
      $        no2colx(i)
       end do
@@ -969,25 +966,25 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = klev-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) * 
+         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) * 
      $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
          !H2O2 interpolated coefficient
-         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,klev) *
+         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
      $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
          !O3 interpolated coefficient
-         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,klev) *
+         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
      $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
       enddo
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,klev
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,klev-i+1) = 
-     $           jfotsout(indexint,13,klev) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -1007,6 +1004,6 @@
 
       indexint=35
-      do i=1,klev
-         auxcolinp(klev-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
       end do
 
@@ -1030,21 +1027,21 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = klev-i+1
+         auxi = nlayer-i+1
          !H2O2 interpolated coefficient
-         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,klev) *
+         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
      $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
          !O3 interpolated coefficient
-         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,klev) *
+         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
      $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
       enddo
       if(chemthermod.ge.2) then
-         do i=1,klev
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,klev-i+1) = 
-     $           jfotsout(indexint,13,klev) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -1063,6 +1060,6 @@
 
       indexint=36
-      do i=1,klev
-         auxcolinp(klev-i+1) = o3colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
       end do
 
@@ -1084,19 +1081,19 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
-      do i=1,klev
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
          !O3 interpolated coefficient
-         jfotsout(indexint,7,klev-i+1) = 
-     $        jfotsout(indexint,7,klev) *
+         jfotsout(indexint,7,nlayer-i+1) = 
+     $        jfotsout(indexint,7,nlayer) *
      $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
       enddo
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,klev
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,klev-i+1) = 
-     $           jfotsout(indexint,13,klev) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -1107,5 +1104,292 @@
 c     End of interpolation to obtain photoabsorption rates
 
-
+c     Coefficients are refered to Sun-Mars distance
+c     Correction to the Sun-Venus distance (fixed)
+
+      jfotsout(:,:,:)=jfotsout(:,:,:)*(1.52/dist_sol)**2
+
+      end
+
+c**********************************************************************
+c**********************************************************************
+
+      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
+     $     co2colx,o3pcolx, n2colx, cocolx)
+
+c     mar 2014        gg            adapted to Venus GCM
+c     nov 2002        fgg           first version
+
+c**********************************************************************
+      use dimphy
+      use param_v4_h
+      implicit none
+
+c    common variables and constants
+#include "clesphys.h"
+#include "mmol.h"
+
+c    local parameters and variables
+
+c     input and output variables
+
+      integer    ig
+      integer    chemthermod
+      integer    nesptherm                   !# of species undergoing
+chemistry, input
+      real       rm(klev,nesptherm)         !densities (cm-3), input
+      real       tx(klev)                   !temperature profile, input
+      real       iz(klev+1)                 !height profile, input
+      real       zenit                      !SZA, input
+      real       co2colx(klev)              !column density of CO2 (cm^-2), output
+      real       o3pcolx(klev)              !column density of O(3P)(cm^-2), output
+      real       n2colx(klev)               !N2 column density (cm-2), output
+      real       cocolx(klev)               !CO column density (cm-2), output
+
+c     local variables
+
+      real       xx
+      real       grav(klev)
+      real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
+      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
+
+      real       co2x(klev)
+      real       o3px(klev)
+      real       cox(klev)
+      real       n2x(klev)
+      real       nx(klev)
+
+      integer    i,j,k,icol,indexint          !indexes
+
+      integer    nz3
+
+      integer    jj
+      real*8      esp(klev*2)
+      real*8      ilayesp(klev*2)
+      real*8      szalayesp(klev*2)
+      integer     nlayesp
+      real*8      zmini
+      real*8      depth
+      real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
+      real*8      espn2,espn,espno,espco,esph,espno2
+      real*8      rcmnz, rcmmini
+      real*8      szadeg
+
+! Tracer indexes in the thermospheric chemistry:
+!!! ATTENTION. These values have to be identical to those in euvheat.F90
+!!! If the values are changed there, the same has to be done here  !!!
+
+      integer,parameter :: ix_co2=1
+      integer,parameter :: ix_n2=13
+      integer,parameter :: ix_o=3
+      integer,parameter :: ix_co=4
+
+c*************************PROGRAM STARTS*******************************
+
+      nz3 = klev*2
+      do i=1,klev
+        xx = ( radio + iz(i) ) * 1.e5    ! conversion [km] ---> [cm]
+         grav(i) = gg * masa /(xx**2)    ! [cm/s2]
+      end do
+
+      !Scale heights  H = kT /Mg  -->  [cm]
+      xx = kboltzman * tx(klev) * n_avog / grav(klev)   ! g cm mol-1
+
+      Ho3p  = xx / mmolo
+      Hco2  = xx / mmolco2
+      Hco   = xx / mmolco
+      Hn2   = xx / mmoln2
+      Hn    = xx / mmoln
+
+      !Only if O3 chem. required
+c      if(chemthermod.ge.1) 
+!     $     Ho3   = xx / mmol(igcm_o3)
+c     $     Ho3   = xx / mmolo3
+c      !Only if N or ion chem.
+c      if(chemthermod.ge.2) then
+c         Hn2   = xx / mmoln2
+c         Hn    = xx / mmoln
+c         Hno   = xx / mmolno
+c         Hno2  = xx / mmolno2
+c      endif
+      ! first loop in altitude : initialisation
+      do i=klev,1,-1
+         !Column initialisation
+         co2colx(i)  = 0.
+         o3pcolx(i)  = 0.
+         n2colx(i)   = 0.
+         cocolx(i)   = 0.
+
+         !--Densities [cm-3]
+         co2x(i)  = rm(i,ix_co2)
+         o3px(i)  = rm(i,ix_o)
+         cox(i)   = rm(i,ix_co)
+         n2x(i)   = rm(i,ix_n2)
+
+c         write(*,*), '--jthermcalc--', co2x(i)
+
+         !Only if O3 chem. required
+c         if(chemthermod.ge.1) 
+c     $        o3x(i)   = rm(i,i_o3)
+c         !Only if Nitrogen of ion chemistry requested
+c         if(chemthermod.ge.2) then
+c            n2x(i)   = rm(i,i_n2)
+c            nx(i)    = rm(i,i_n)
+c            nox(i)   = rm(i,i_no)
+c            no2x(i)  = rm(i,i_no2)
+c         endif
+      enddo     ! end first loop
+
+      ! second loop in altitude : column calculations
+      do i=klev,1,-1
+         !Routine to calculate the geometrical length of each layer
+         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
+     $         szalayesp,nlayesp, zmini)
+         if(ilayesp(nlayesp).eq.-1) then
+            co2colx(i)=1.e25
+            o3pcolx(i)=1.e25
+            n2colx(i)=1.e25
+            cocolx(i)=1.e25
+
+c            o2colx(i)=1.e25
+c            o3pcolx(i)=1.e25
+c            h2colx(i)=1.e25
+c            h2ocolx(i)=1.e25
+c            h2o2colx(i)=1.e25
+c            o3colx(i)=1.e25
+c            ncolx(i)=1.e25
+c            nocolx(i)=1.e25
+
+c            cocolx(i)=1.e25
+c            hcolx(i)=1.e25
+c            no2colx(i)=1.e25
+         else
+            rcmnz = ( radio + iz(klev) ) * 1.e5    ! km --> cm
+            rcmmini = ( radio + zmini ) * 1.e5
+            !Column calculation taking into account the geometrical
+            !depth 
+            !calculated before
+            do j=1,nlayesp
+               jj=ilayesp(j)
+               !Top layer
+               if(jj.eq.klev) then
+                  if(zenit.le.60.) then
+                     o3pcolx(i)=o3pcolx(i)+o3px(klev)*Ho3p*esp(j)
+     $                    *1.e-5
+                     co2colx(i)=co2colx(i)+co2x(klev)*Hco2*esp(j)
+     $                    *1.e-5
+                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
+     $                    *1.e-5
+                     n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
+     $                    *1.e-5
+
+c                     h2o2colx(i)=h2o2colx(i)+
+c     $                    h2o2x(klev)*Hh2o2*esp(j)*1.e-5
+c                     o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)
+c     $                    *1.e-5
+c                     h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)
+c     $                    *1.e-5
+c                     h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)
+c     $                    *1.e-5                     
+c                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
+c     $                    *1.e-5
+c                     hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)
+c     $                    *1.e-5
+                     !Only if O3 chemistry required
+c                     if(chemthermod.ge.1) o3colx(i)=
+c     $                    o3colx(i)+o3x(klev)*Ho3*esp(j)
+c     $                    *1.e-5
+                     !Only if N or ion chemistry requested
+c                     if(chemthermod.ge.2) then
+c                        n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
+c     $                    *1.e-5
+
+c                     endif
+                  else if(zenit.gt.60.) then
+                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
+                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
+                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
+                     espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
+                     espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
+
+c                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
+c                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
+c                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
+c                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
+c                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
+                     !Only if O3 chemistry required
+c                     if(chemthermod.ge.1)                     
+c     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
+c                     !Only if N or ion chemistry requested
+c                     if(chemthermod.ge.2) then
+c                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
+c                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
+c                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
+c                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
+c                     endif
+
+                     co2colx(i) = co2colx(i) + espco2*co2x(klev)
+                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev)
+                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
+                     n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
+
+c                     o2colx(i)  = o2colx(i)  + espo2*o2x(klev)
+c                     h2colx(i)  = h2colx(i)  + esph2*h2x(klev)
+c                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev)
+c                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev)
+c                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
+c                     hcolx(i)   = hcolx(i)   + esph*hx(klev)
+                     !Only if O3 chemistry required
+c                     if(chemthermod.ge.1)                      
+c     $                  o3colx(i) = o3colx(i)  + espo3*o3x(klev)
+c                     !Only if N or ion chemistry requested
+c                     if(chemthermod.ge.2) then
+c                        n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
+c                        ncolx(i)   = ncolx(i)   + espn*nx(klev)
+c                        nocolx(i)  = nocolx(i)  + espno*nox(klev)
+c                        no2colx(i) = no2colx(i) + espno2*no2x(klev)
+c                     endif
+                  endif !Of if zenit.lt.60
+               !Other layers
+               else
+                  co2colx(i)  = co2colx(i) +
+     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
+                  o3pcolx(i)  = o3pcolx(i) +
+     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
+                  cocolx(i)   = cocolx(i) +
+     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
+                  n2colx(i)   = n2colx(i) +
+     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
+
+c
+c                  o2colx(i)   = o2colx(i) + 
+c     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
+c                  h2colx(i)   = h2colx(i) + 
+c     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
+c                  h2ocolx(i)  = h2ocolx(i) + 
+c     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
+c                  h2o2colx(i) = h2o2colx(i) + 
+c     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
+c                  hcolx(i)    = hcolx(i) + 
+c     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
+                  !Only if O3 chemistry required
+c                  if(chemthermod.ge.1) 
+c     $                 o3colx(i) = o3colx(i) + 
+c     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
+c                  !Only if N or ion chemistry requested
+c                  if(chemthermod.ge.2) then
+c                     n2colx(i)   = n2colx(i) + 
+c     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
+c                     ncolx(i)    = ncolx(i) + 
+c     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
+c                     nocolx(i)   = nocolx(i) + 
+c     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
+c                     no2colx(i)  = no2colx(i) + 
+c     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
+c                  endif
+
+               endif  !Of if jj.eq.klev
+            end do    !Of do j=1,nlayesp
+         endif        !Of ilayesp(nlayesp).eq.-1
+      enddo           !Of do i=klev,1,-1
       return
 
@@ -1113,6 +1397,340 @@
 
 
-
-
-
-
+c**********************************************************************
+c**********************************************************************
+
+      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
+C
+C subroutine to perform linear interpolation in pressure from 1D profile 
+C escin(nl) sampled on pressure grid pin(nl) to profile
+C escout(nlayer) on pressure grid p(nlayer).
+C
+      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
+            wm(n1) = 0.d0
+            wp(n1) = 0.d0
+         else
+            do n = nini,nl-1
+               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
+                  nm(n1)=n
+                  np=n+1
+                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
+                  wp(n1)=1.d0 - wm(n1)
+                  nini = n
+                  exit
+               endif
+            enddo
+         endif
+      enddo
+      return
+      end
+
+c**********************************************************************
+c**********************************************************************
+
+      subroutine espesor_optico_A (ig,capa, szadeg,z,
+     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
+
+c     fgg              nov 03      Adaptation to Martian model
+c     malv             jul 03      Corrected z grid. Split in alt & frec
+codes
+c     fgg              feb 03      first version
+*************************************************************************
+      use dimphy
+      use param_v4_h
+      implicit none
+
+c     arguments
+
+      real        szadeg                ! I. SZA [rad]
+      real        z                     ! I. altitude of interest [km]
+      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
+                                        !  (=2*klev= max# of layers in
+                                        !  ray path)
+      real     iz(klev+1)              ! I. Altitude of each layer
+      real*8        esp(nz3)            ! O. layer widths after geometrically 
+                                        !    amplified; in [cm] except
+                                        !    at TOA
+                                        !    where an auxiliary value is
+                                        !    used
+      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
+      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "     "
+      integer       nlayesp
+!      real*8        nlayesp             ! O. # layers along ray path at
+!      this z
+      real*8        zmini               ! O. Minimum altitud of ray path [km]
+
+c     local variables and constants
+
+        integer     j,i,capa
+        integer     jmin                  ! index of min.altitude along ray path
+        real*8      szarad                ! SZA [deg]
+        real*8      zz
+        real*8      diz(klev+1)
+        real*8      rkmnz                 ! distance TOA to center of Planet [km]
+        real*8      rkmmini               ! distance zmini to center of P [km] 
+        real*8      rkmj                  ! intermediate distance to C of P [km]
+
+c external function
+        external  grid_R8          ! Returns index of layer containing the altitude
+                                ! of interest, z; for example, if 
+                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) =>
+                                ! grid(z)=i 
+        integer   grid_R8
+
+*************************************************************************     
+        szarad = dble(szadeg)*3.141592d0/180.d0
+        zz=dble(z)
+        do i=1,klev
+           diz(i)=dble(iz(i))
+        enddo
+        do j=1,nz3
+          esp(j) = 0.d0
+          szalayesp(j) = 777.d0
+          ilayesp(j) = 0
+        enddo
+        nlayesp = 0
+
+        ! First case: szadeg<60
+        ! The optical thickness will be given by  1/cos(sza)
+        ! We deal with 2 different regions:
+        !   1: First, all layers between z and ztop ("upper part of
+        !   ray")
+        !   2: Second, the layer at ztop
+        if(szadeg.lt.60.d0) then
+
+           zmini = zz
+           if(abs(zz-diz(klev)).lt.1.d-3) goto 1357
+           ! 1st Zone: Upper part of ray
+           !
+           do j=grid_R8(zz,diz,klev),klev-1
+             nlayesp = nlayesp + 1
+             ilayesp(nlayesp) = j
+             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
+             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
+             szalayesp(nlayesp) = szadeg
+           end do
+
+           ! 
+           ! 2nd Zone: Top layer
+ 1357      continue
+           nlayesp = nlayesp + 1
+           ilayesp(nlayesp) = klev 
+           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
+           szalayesp(nlayesp) = szadeg
+
+        ! Second case:  60 < szadeg < 90
+        ! The optical thickness is evaluated.
+        !    (the magnitude of the effect of not using cos(sza) is 3.e-5 
+        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60,
+        !     approximately)
+        ! We deal with 2 different regions:
+        !   1: First, all layers between z and ztop ("upper part of
+        !   ray")
+        !   2: Second, the layer at ztop ("uppermost layer")
+        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
+
+           zmini=(radio+zz)*sin(szarad)-radio
+           rkmmini = radio + zmini
+
+           if(abs(zz-diz(klev)).lt.1.d-4) goto 1470
+
+           ! 1st Zone: Upper part of ray
+           !
+           do j=grid_R8(zz,diz,klev),klev-1
+              nlayesp = nlayesp + 1 
+              ilayesp(nlayesp) = j
+              esp(nlayesp) = 
+     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
+     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
+              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
+              rkmj = radio+diz(j)
+              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
+              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592
+! [deg]
+           end do
+ 1470      continue
+           ! 2nd Zone:  Uppermost layer of ray.
+           !
+           nlayesp = nlayesp + 1 
+           ilayesp(nlayesp) = klev
+           rkmnz = radio+diz(klev)
+           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
+           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
+           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
+           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
+
+
+        ! Third case:  szadeg > 90
+        ! The optical thickness is evaluated.
+        ! We deal with 5 different regions:
+        !   1: all layers between z and ztop ("upper part of ray")
+        !   2: the layer at ztop ("uppermost layer")
+        !   3: the lowest layer, at zmini
+        !   4: the layers increasing from zmini to z (here SZA<90)
+        !   5: the layers decreasing from z to zmini (here SZA>90)
+        else if(szadeg.gt.90.d0) then
+
+           zmini=(radio+zz)*sin(szarad)-radio
+           rkmmini = radio + zmini
+
+           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
+             nlayesp = nlayesp + 1 
+             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
+!             esp(nlayesp) = 1.e30
+
+           else
+              jmin=grid_R8(zmini,diz,klev)+1
+              
+
+              if(abs(zz-diz(klev)).lt.1.d-4) goto 9876
+
+              ! 1st Zone: Upper part of ray
+              !
+              do j=grid_R8(zz,diz,klev),klev-1
+                nlayesp = nlayesp + 1 
+                ilayesp(nlayesp) = j
+                esp(nlayesp) = 
+     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
+     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
+                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
+                rkmj = radio+diz(j)
+                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
+                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
+! [deg]
+              end do
+
+ 9876         continue
+              ! 2nd Zone:  Uppermost layer of ray.
+              !
+              nlayesp = nlayesp + 1 
+              ilayesp(nlayesp) = klev
+              rkmnz = radio+diz(klev)
+              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 ) !aux.factor[km]
+              esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm]
+              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
+              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
+
+              ! 3er Zone: Lowestmost layer of ray
+              !
+              if ( jmin .ge. 2 ) then      ! If above the planet's surface
+                j=jmin-1
+                nlayesp = nlayesp + 1 
+                ilayesp(nlayesp) = j
+                esp(nlayesp) = 2. * 
+     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
+                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
+                rkmj = radio+diz(j+1)
+                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
+                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
+! [deg]
+              endif
+
+              ! 4th zone: Lower part of ray, increasing from zmin to z
+              !    ( layers with SZA < 90 deg )
+              do j=jmin,grid_R8(zz,diz,klev)-1
+                nlayesp = nlayesp + 1 
+                ilayesp(nlayesp) = j
+                esp(nlayesp) = 
+     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
+     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
+                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
+                rkmj = radio+diz(j)
+                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
+                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
+! [deg]
+              end do
+
+              ! 5th zone: Lower part of ray, decreasing from z to zmin
+              !    ( layers with SZA > 90 deg )
+              do j=grid_R8(zz,diz,klev)-1, jmin, -1
+                nlayesp = nlayesp + 1 
+                ilayesp(nlayesp) = j
+                esp(nlayesp) = 
+     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
+     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )
+! [km]
+                esp(nlayesp) = esp(nlayesp) * 1.d5
+! [cm]
+                rkmj = radio+diz(j)
+                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )
+! [rad]
+                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592
+! [deg]
+              end do
+
+           end if
+
+        end if
+
+        return
+
+        end
+
+
+c**********************************************************************
+c***********************************************************************
+
+        function grid_R8 (z, zgrid, nz)
+
+c Returns the index where z is located within vector zgrid
+c The vector zgrid must be monotonously increasing, otherwise program stops.
+c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops. 
+c
+c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le. 
+c MALV    Jul-2003
+c***********************************************************************
+
+        implicit none
+
+c Arguments 
+        integer   nz
+        real*8      z
+        real*8      zgrid(nz)
+        integer   grid_R8
+
+c Local  
+        integer   i, nz1, nznew
+
+c*** CODE START 
+
+        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then 
+           write (*,*) ' GRID/ z outside bounds of zgrid '
+           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
+           stop ' Serious error in GRID.F '
+        endif
+        if ( nz .lt. 2 ) then 
+           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
+           stop ' Serious error in GRID.F '
+        endif
+        if ( zgrid(1) .ge. zgrid(nz) ) then 
+           write (*,*) ' GRID/ zgrid must increase with index'
+           stop ' Serious error in GRID.F '
+        endif
+
+        nz1 = 1
+        nznew = nz/2
+        if ( z .gt. zgrid(nznew) ) then
+           nz1 = nznew
+           nznew = nz
+        endif
+        do i=nz1+1,nznew
+           if ( z. eq. zgrid(i) ) then
+              grid_R8=i
+              return
+              elseif ( z .le. zgrid(i) ) then
+              grid_R8 = i-1
+              return
+           endif
+        enddo
+        grid_R8 = nz
+        return
+
+        end
+
Index: trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F	(revision 2464)
@@ -44,4 +44,5 @@
       real   lambdamin(nnuve),lambdamax(nnuve) ! in microns
       real   dlambda                  ! in microns 
+      logical upper
 
       nlve = klev
@@ -87,11 +88,21 @@
         read(10,*)
         read(10,*) m,Nb
-cc      GG modif below  
-        if (nlve.le.78.and.m.ne.nlve) then
-         write(*,*) 'This is subroutine load_ksi'
-         print*,'Dimension problem between ksi.txt and nlve'
-         print*,'N levels = ',m,nlve
-         stop
+CC vertical axis check
+        upper=.false.
+cc the number of levels read for the matrix (m) should be the number of levels of the GCM (nlve=klev)
+        if (m.ne.nlve) then
+CC When GCM axis has more levels than the matrix, we must be in the case
+cc with upper atmosphere, ie matrix limited to 78 levels and upper levels set to zero
+cc otherwise there is a problem
+         if ((m.eq.78).and.(nlve.gt.78)) then
+          upper=.true.
+         else
+          write(*,*) 'This is subroutine load_ksi'
+          print*,'Dimension problem between ksi.txt and nlve'
+          print*,'N levels = ',m,nlve
+          stop
+         endif
         endif
+cc wavelength check
         if (Nb.ne.nnuve) then
          write(*,*) 'This is subroutine load_ksi'
@@ -113,4 +124,13 @@
             read(10,format_lect) (ksive(i,j,band,mat),j=0,m+1) ! no unit
          enddo                  ! i
+cc case upper atmosphere: the level for space (m+1) should be raised to nlve+1
+         if (upper) then
+           do i=0,m
+             ksive(i,klev+1,band,mat)=ksive(i,m+1,band,mat)
+             ksive(i,m+1,band,mat)=0.
+           enddo
+             ksive(klev+1,:,band,mat)=ksive(m+1,:,band,mat)
+             ksive(m+1,:,band,mat)=0.
+         endif
         enddo                     ! band
 c       print*,"Matrice ",mat," lue"
Index: trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90	(revision 2464)
@@ -1,4 +1,4 @@
 subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pq,&
-                       ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
+                       ptimestep,pdteuv,pdtconduc,pdqdiff)
 
 USE chemparam_mod
@@ -8,11 +8,5 @@
 implicit none
 
-!#include "dimphys.h"
 #include "comcstfi.h"
-!#include "callkeys.h"
-!#include "comdiurn.h"
-!#include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 #include "diffusion.h"
 
@@ -26,5 +20,4 @@
       real ptimestep
       real pplay(ngrid,nlayer)
-      real zzlay(ngrid,nlayer)
       real pplev(ngrid,nlayer+1)
       real pq(ngrid,nlayer,nq)
@@ -67,5 +60,6 @@
       real*8,dimension(:),allocatable,save ::  wi,Wad,Uthermal,Lambdaexo,Hspecie
       real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
-      character(len=20),dimension(14) :: ListeDiff
+      integer,parameter :: nb_esp_diff=15
+      character(len=20),dimension(nb_esp_diff) :: ListeDiff
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 !     tracer numbering in the molecular diffusion
@@ -130,4 +124,5 @@
 	ListeDiff(13)='o3'
 	ListeDiff(14)='n'
+	ListeDiff(15)='he'
 	i_esc_h=1000
 	i_esc_h2=1000
@@ -140,5 +135,5 @@
 	do nn=1,nq
 
-	do n=1,14
+	do n=1,nb_esp_diff
 	if (ListeDiff(n) .eq. tname(nn)) then
 	indic_diff(nn)=1
@@ -624,4 +619,7 @@
         enddo
         enddo
+
+! DEBUG
+!       print*,"rapport MASSE: ",Mtot2(:)/Mtot1(:)
 
 ! Check mass  conservation of each specie on column
Index: trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F	(revision 2464)
@@ -189,5 +189,10 @@
                   if(pq(ig,l,ico2) .gt. 1.e-6) then
                      oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
-
+                     ! handle the rare cases when pq(ig,l,io)<0
+                     if (pq(ig,l,io).lt.0) then
+                       write(*,*) "nirco2abs: warning ig=",ig," l=",l,
+     &                            " pq(ig,l,io)=",pq(ig,l,io)
+                       oco2gcm=1.e6
+                     endif
                   else
                      oco2gcm=1.e6
Index: trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F	(revision 2464)
@@ -177,7 +177,7 @@
                                 !! k19 & k20
 
-      k20x = 3.d-12
+c      k20x = 3.d-12
 c  TEST GG: double the values of Kvv as recently found by Sharma et al.2014
-c      k20x = 6.d-12   
+      k20x = 6.d-12   
 c  TEST GG: use the minimum value of the experimental bracket's values [1-6]
 c      k20x = 1.d-12 
Index: trunk/LMDZ.VENUS/libf/phyvenus/param.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/param.h	(revision 2463)
+++ 	(revision )
@@ -1,42 +1,0 @@
-
-!***********************************************
-
-!	param.par
-
-!	Parameters for paramhr.f
-!***********************************************
-
-	integer    ninter
-	parameter  (ninter=36)
-
-	integer    nabs
-	parameter  (nabs=13)
-!	parameter  (nabs=5)
-
-	integer    nz2
-	parameter  (nz2=253)
-
-	integer    ninter2
-        parameter  (ninter2=16)
-
-	real       kboltzman                  !cte Boltzman
-	parameter  (kboltzman = 1.381e-16)
-	
-	real       n_avog                    !# de Avogadro
-	parameter  (n_avog = 6.023e23)
-
-	real       gg                        !cte gravitacion
-	parameter  (gg = 6.67259e-8)
-
-	real       masa                      !masa de Venus(g)
-	parameter  (masa = 4.8676e27)
-
-	real       radio
-	parameter  (radio = 6051.8)           !radio de Venus(km)
-
-	integer	   nreact
-	parameter  (nreact=93)
-
-
-
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/param_read.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/param_read.F	(revision 2463)
+++ 	(revision )
@@ -1,443 +1,0 @@
-      subroutine param_read
-
-      use dimphy
-      use conc
-      implicit none
-
-
-c     common variables and constants
-      include 'param.h'
-      include 'param_v4.h'
-c      include 'datafile.h'
- 
- 
-c     local variables
-
-      integer    i,j,k,inter                          !indexes
-      integer ierr
-      real       nada
-      character (len=100) :: datafile="HIGHATM"
-      
-c*************************+PROGRAM STARTS**************************
-
-
-c     data for the UV heating tabulation
-
-      data (crscabsi2(1,j),j=1,16) /5.61031E-19,1.59677E-18,4.7072E-18,
-     $     1.48254e-17,2.07445e-17,2.573e-17,2.901e-17,3.083e-17,
-     $     3.217e-17,3.539e-17,3.658e-17,3.63e-17,3.41239e-17,
-     $     2.71019e-17,4.93677e-17,1.64e-17/
-
-      data (crscabsi2(2,j),j=1,16) /0.27250E-18,0.11650E-17,0.39250E-17,
-     $     0.10630E-16,0.15590E-16,0.17180E-16,0.19270E-16,0.22860E-16,
-     $     0.24270E-16,0.24440E-16,0.25020E-16,0.26600E-16,0.25400E-16,
-     $     0.35800E-16,0.25590E-16,0.16740E-16/
-
-      data (crscabsi2(3,j),j=1,16) /0.2776E-18,0.9792E-18,0.3313E-17,
-     $     0.6621E-17,0.8481E-17,0.9146E-17,0.9414E-17,0.1039E-16,
-     $     0.1012E-16,0.1033E-16,0.1033E-16,0.1033E-16,0.8268E-17,
-     $     0.6563E-17,0.3506E-17,0.3470E-17/
-
-      data (crscabsi2(5,j),j=1,16) /.5E-20,.1077607E-19,.5670491E-19,
-     $     .3322716E-18,.1054509E-17,.1700005E-17,.3171188E-17,
-     $     .4734241E-17,.5108741E-17,.6022236E-17,.6741537E-17,
-     $     .7277079E-17,.9070787E-17,.9708916E-17,.4026281E-17,0.0/
-
-      data (crscabsi2(8,j),j=1,16) /0.0, 7.44175e-19, 2.23167e-18,
-     $    8.46200e-18,1.18275e-17,1.54900e-17,2.32475e-17,2.41373e-17,
-     $     2.55482e-17,2.38431e-17,2.28600e-17,2.35067e-17,2.56000e-17,
-     $     2.64636e-17,2.86260e-17,3.26561e-17/
-
-      data(crscabsi2(9,j),j=1,16) /3.48182e-20,3.37038e-19,1.03077e-18,
-     $     4.01364e-18,6.45e-18,7.8e-18,1.0e-17,1.13500e-17,1.15500e-17,
-     $     1.18000e-17,1.17500e-17,1.16000e-17,1.28667e-17,1.18500e-17,
-     $     1.11000e-17,9.50000e-18/
-
-      data(crscabsi2(10,j),j=1,16) /0.0,9.39833e-19,2.87714e-18,
-     $     9.66900e-18,1.37063e-17,1.61820e-17,2.30450e-17,2.63373e-17,
-     $     2.63773e-17,2.67677e-17,2.64100e-17,2.53000e-17,2.18100e-17,
-     $     2.04941e-17,2.28160e-17,2.93550e-17/
-
-      data(crscabsi2(11,j),j=1,16) /0.0,9.58555e-19,2.52767e-18,
-     $     8.29700e-18,1.21850e-17,1.40500e-17,1.97025e-17,2.12018e-17,
-     $     2.14673e-17,2.20331e-17,2.21500e-17,2.21600e-17,2.33200e-17,
-     $     2.67800e-17,2.56400e-17,3.58561e-17/
-
-      data(crscabsi2(12,j),j=1,16) /0.0,1.0e-20,2.5e-20,1.30400e-19,
-     $     2.93800e-19,4.36000e-19,8.49400e-19,1.29400e-18,1.40500e-18,
-     $     1.67600e-18,1.93400e-18,2.12200e-18,2.75800e-18,3.48400e-18,
-     $     4.17200e-18,5.26000e-18/
-
-      data(crscabsi2(13,j),j=1,16) /0.0,1.60e-18,4.99111e-18,1.48496e-17
-     $     ,2.17395e-17,2.55857e-17,2.87754e-17,3.65571e-17,3.85691e-17,
-     $     4.16286e-17,4.15117e-17,4.05901e-17,3.64000e-17,2.99670e-17,
-     $     2.46796e-17,2.51789e-17/
-
-      data freccen /3.4,7.5,14.5,23.0,30.3,34.1,49.6,50.5,52.5,56.0,
-     $59.0,61.5,68.7,73.1,78.4,83.1,92.4,97.5,99.3,100.1,100.7,102.1,
-     $104.5,116.8,121.3,127.0,130.6,153.7,162.8,171.4
-     $,195.6,206.3,222.0,236.0,289.0,600./
-
-      data co2crsc195/2.05864e-17,5.90557e-20,3.1027e-19,6.70653e-19,
-     $4.55132e-19,8.87122e-20,1.32138e-20,7.22244e-23,2.88002e-26/
-
-      data co2crsc295/2.05897e-17,6.71104e-20,3.45509e-19,7.45711e-19,
-     $4.82752e-19,1.11594e-19,1.98308e-20,1.3853e-22,2.1414e-25/
-
-c     Reads tabulated functions
-
-      !Tabulated column amount
-      open(210, status = 'old',
-c    $file=trim(datafile)//'/EUVDAT/coln.dat',iostat=ierr)
-     $file=trim(datafile)//'/EUVDAT/param_v5/coln.dat',iostat=ierr)
-
-      IF (ierr.NE.0) THEN 
-       write(*,*)'cant find directory EUVDAT containing param_v5 subdir'
-       write(*,*)'(in aeronomars/param_read.F)'
-       write(*,*)'It should be in :', trim(datafile),'/'
-       write(*,*)'1) You can change this directory address in '
-       write(*,*)'   callphys.def with datadir=/path/to/dir'
-       write(*,*)'2) If necessary, EUVDAT (and other datafiles)'
-       write(*,*)'   can be obtained online on:'
-       write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
-       STOP
-      ENDIF
- 
-      !Tabulated photoabsorption coefficients
-      open(220,file=trim(datafile)//'/EUVDAT/param_v5/j2_an.dat')
-      open(230,file=trim(datafile)//'/EUVDAT/param_v5/j3_an.dat')
-      open(240,file=trim(datafile)//'/EUVDAT/param_v5/j1_an.dat')
-      open(250,file=trim(datafile)//'/EUVDAT/param_v5/j2_bn.dat')
-      open(260,file=trim(datafile)//'/EUVDAT/param_v5/j2_cn.dat')
-      open(300,file=trim(datafile)//'/EUVDAT//param_v5/j2_dn.dat')
-      open(270,file=trim(datafile)//'/EUVDAT//param_v5/j1_bn.dat')
-      open(280,file=trim(datafile)//'/EUVDAT//param_v5/j1_cn.dat')
-      open(290,file=trim(datafile)//'/EUVDAT//param_v5/j1_dn.dat')
-      open(150,file=trim(datafile)//'/EUVDAT//param_v5/j4n.dat')
-      open(160,file=trim(datafile)//'/EUVDAT//param_v5/j5n.dat')
-      open(170,file=trim(datafile)//'/EUVDAT//param_v5/j6n.dat')
-      open(180,file=trim(datafile)//'/EUVDAT//param_v5/j7n.dat')
-      open(390,file=trim(datafile)//'/EUVDAT//param_v5/j8_an.dat')
-      open(400,file=trim(datafile)//'/EUVDAT//param_v5/j8_bn.dat')
-      open(410,file=trim(datafile)//'/EUVDAT//param_v5/j9n.dat')
-      open(420,file=trim(datafile)//'/EUVDAT//param_v5/j10_an.dat')
-      open(430,file=trim(datafile)//'/EUVDAT//param_v5/j10_bn.dat')
-      open(440,file=trim(datafile)//'/EUVDAT//param_v5/j10_cn.dat')
-      open(450,file=trim(datafile)//'/EUVDAT//param_v5/j11_an.dat')
-      open(460,file=trim(datafile)//'/EUVDAT//param_v5/j11_bn.dat')
-      open(470,file=trim(datafile)//'/EUVDAT//param_v5/j11_cn.dat')
-      open(480,file=trim(datafile)//'/EUVDAT//param_v5/j12n.dat')
-      open(490,file=trim(datafile)//'/EUVDAT//param_v5/j13_an.dat')
-      open(500,file=trim(datafile)//'/EUVDAT//param_v5/j13_bn.dat')
-      open(510,file=trim(datafile)//'/EUVDAT//param_v5/j13_cn.dat')
-
-      
-      do i=210,300,10
-         read(i,*)
-         read(i,*)
-      end do
-
-      do i=150,180,10
-         read(i,*)
-         read(i,*)
-      end do
-
-      do i=390,510,10
-         read(i,*)
-         read(i,*)
-      enddo
-
-      do i=nz2,1,-1
-         read(210,*) (c1_16(i,j),j=1,16),c17_24(i),c25_29(i),c30_31(i),
-     $        c32(i),c33(i),c34(i),c35(i),c36(i)
-      end do
-
-      do i=nz2,1,-1
-         read(220,*) (jabsifotsintpar(i,2,j),j=1,16)
-      end do
-      
-      do i=nz2,1,-1
-         read(230,*) (jabsifotsintpar(i,3,j),j=1,16)
-      end do
-
-      do i=nz2,1,-1
-         read(240,*) (jabsifotsintpar(i,1,j),j=1,16)
-      end do
-
-      do i=nz2,1,-1
-         read(250,*) (jabsifotsintpar(i,2,j),j=17,24)
-      end do
-
-
-      do i=nz2,1,-1
-         read(260,*) (jabsifotsintpar(i,2,j),j=25,31)
-      end do
-
-      do i=nz2,1,-1
-         read(270,*) (jabsifotsintpar(i,1,j),j=17,24)
-      end do
-
-      do i=nz2,1,-1
-         read(280,*) (jabsifotsintpar(i,1,j),j=25,31)
-      end do
-
-      do i=nz2,1,-1
-         read(290,*) jabsifotsintpar(i,1,32)
-      end do
-
-      do i=nz2,1,-1
-         read(300,*) (jabsifotsintpar(i,2,j),j=32,34)
-      end do
-
-      do i=nz2,1,-1
-         read(160,*) (jabsifotsintpar(i,5,j),j=1,15)
-      end do
-
-      do i=nz2,1,-1
-         read(150,*) (jabsifotsintpar(i,4,j),j=25,31)
-      end do
-
-      do i=nz2,1,-1
-         read(170,*) (jabsifotsintpar(i,6,j),j=25,35)
-      end do
-
-      do i=nz2,1,-1
-         read(180,*) (jabsifotsintpar(i,7,j),j=34,36)
-      end do
-
-      do i=nz2,1,-1
-         read(390,*) (jabsifotsintpar(i,8,j),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(400,*) (jabsifotsintpar(i,8,j),j=17,24)
-      enddo
-
-      do i=nz2,1,-1
-         read(410,*) (jabsifotsintpar(i,9,j),j=1,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(420,*) (jabsifotsintpar(i,10,j),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(430,*) (jabsifotsintpar(i,10,j),j=17,24)
-      enddo
-
-      do i=nz2,1,-1
-         read(440,*) (jabsifotsintpar(i,10,j),j=25,32)
-      enddo
-
-      do i=nz2,1,-1
-         read(450,*) (jabsifotsintpar(i,11,j),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(460,*) (jabsifotsintpar(i,11,j),j=17,24)
-      enddo
-
-      do i=nz2,1,-1
-         read(470,*) (jabsifotsintpar(i,11,j),j=25,29)
-      enddo
-      
-      do i=nz2,1,-1
-         read(480,*) (jabsifotsintpar(i,12,j),j=2,16)
-      enddo
-
-      do i=nz2,1,-1
-         read(490,*) (jabsifotsintpar(i,13,j),j=2,16)
-      enddo
-      
-      do i=nz2,1,-1
-         read(500,*) (jabsifotsintpar(i,13,j),j=17,24)
-      enddo
-      
-      do i=nz2,1,-1
-         read(510,*) (jabsifotsintpar(i,13,j),j=25,36)
-      enddo
-
-      do i=210,300,10
-         close(i)
-      end do
-
-      do i=150,180,10
-         close(i)
-      end do
-
-      do i=390,510,10
-         close(i)
-      enddo
-
-
-c     set t0
-
-      do i=1,nz2
-         t0(i)=195.
-      end do
-
-
-      do i=1,ninter
-         fluxtop(i)=1.
-      end do
-
-      !Parameters for the variation of the solar flux with 11 years cycle
-      open(100,file=trim(datafile)//'/EUVDAT/param_v5/varflujo.dat')
-      read(100,*)
-      do i=1,24
-         read(100,*) inter,ct1(i),p1(i),ct2(i),p2(i),nada
-      end do
-      close(100)
-
-c     dissociation and ionization efficiencies
-
-!      do inter=1,ninter
-!         efdisco2(inter)=0.
-!         efdiso2(inter)=0.
-!         efdish2(inter)=0.
-!         efdish2o(inter)=0.
-!         efdish2o2(inter)=0.
-!         efdiso3(inter)=0.
-!         efdisco(inter)=0.
-!         efdisn2(inter)=0.
-!         efdisno(inter)=0.
-!         efdisno2(inter)=0.
-!         efionco2(inter,1)=0.
-!         efionco2(inter,2)=0.
-!         efionco2(inter,3)=0.
-!         efionco2(inter,4)=0.
-!         efiono3p(inter)=0.
-!         efionn2(inter,1)=0.
-!         efionn2(inter,2)=0.
-!         efionco(inter,1)=0.
-!         efionco(inter,2)=0.
-!         efionco(inter,3)=0.
-!         efionn(inter)=0.
-!         efionh(inter)=0.
-!         efionno(inter)=0.
-!      enddo
-
-
-c     CO2, O2, NO
-
-!      open(120,file=trim(datafile)//'/EUVDAT/param_v5/efdis_inter.dat')
-!      read(120,*)
-!      do i=1,21
-!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
-!      do inter=8,28
-!         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
-!      enddo
-!      do inter=29,ninter
-!         efdisco2(inter)=1.
-!         efdiso2(inter)=1.
-!         efdisno(inter)=1.
-!      enddo
-
-
-c     N2
-
-!      efdisn2(15)=0.1
-!      do inter=16,ninter
-!         efdisn2(inter)=1.
-!      enddo
-
-
-c     CO
-
-!      efdisco(16)=0.5
-!      do inter=17,ninter
-!         efdisco(inter)=1.
-!      enddo
-
-      
-c     O, N, H
-
-!      do inter=1,ninter
-!         efdiso(inter)=0.
-!         efdisn(inter)=0.
-!         efdish(inter)=0.
-!      enddo
-
-
-c     H2O, H2O2, O3, NO2
-
-!      do inter=25,31
-!         efdish2o(inter)=1.
-!      enddo
-!      do inter=25,35
-!         efdish2o2(inter)=1.
-!      enddo
-!      do inter=34,36
-!         efdiso3(inter)=1.
-!      enddo
-!      do inter=27,36
-!         efdisno2(inter)=1.
-!      enddo
-!      do inter=1,15
-!         efdish2(inter)=1.
-!      enddo
-         
-      !4 possible channels for CO2 ionization
-!      do inter=14,16
-!         efionco2(inter,1)=1.-efdisco2(inter)
-!      enddo
-!      efionco2(13,1)=0.805*(1.-efdisco2(13))
-!      efionco2(13,2)=0.195*(1.-efdisco2(13))
-!      do inter=11,12
-!         efionco2(inter,3)=1.-efdisco2(inter)
-!      enddo
-!      efionco2(10,3)=0.9*(1.-efdisco2(10))
-!      efionco2(10,4)=0.1*(1.-efdisco2(10))
-!      do inter=2,9
-!         efionco2(inter,4)=1.-efdisco2(inter)
-!      enddo
-
-      !For O(3p) total ionization under 91.1 nm
-!      do inter=1,16
-!         efiono3p(inter)=1.d0
-!      enddo
-
-      !2 channels for N2 ionization
-!      do inter=9,15
-!         efionn2(inter,1)=1.-efdisn2(inter)
-!      enddo
-!      do inter=2,8
-!         efionn2(inter,2)=1.-efdisn2(inter)
-!      enddo
-      
-      !3 channels for CO ionization
-!      do inter=11,16
-!         efionco(inter,1)=1.-efdisco(inter)
-!      enddo
-!      efionco(10,1)=0.87*(1.-efdisco(10))
-!      efionco(10,2)=0.13*(1.-efdisco(10))
-!      do inter=8,9
-!         efionco(inter,2)=1.-efdisco(inter)
-!      enddo
-!      efionco(7,2)=0.1*(1.-efdisco(7))
-!      efionco(7,3)=0.9*(1.-efdisco(7))
-!      do inter=2,6
-!         efionco(inter,3)=1.-efdisco(inter)
-!      enddo
-
-      !Total ionization under 85 nm for N
-!      do inter=1,16
-!         efionn(inter)=1.
-!      enddo
-
-      !NO
-!      do inter=2,28
-!         efionno(inter)=1.-efdisno(inter)
-!      enddo
-
-      !Total ionization under 90 nm for H
-!      do inter=3,16
-!         efionh(inter)=1.
-!      enddo
-
-
-      return
-
-
-      end
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/param_read_e107.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/param_read_e107.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/param_read_e107.F	(revision 2464)
@@ -1,13 +1,24 @@
       subroutine param_read_e107
 
-	use dimphy
-	use conc
+      use dimphy
+      use param_v4_h, only: jfotsout,crscabsi2,
+     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
+     .    co2crsc195,co2crsc295,t0,
+     .    jabsifotsintpar,ninter,nz2,
+     .    nabs,e107,date_e107,e107_tab,
+     .    coefit0,coefit1,coefit2,coefit3,coefit4,
+     .    efdisco2,efdiso2,efdish2o,
+     .    efdish2o2,efdish2,efdiso3,
+     .    efdiso,efdisn,efdish,
+     .    efdisno,efdisn2,efdisno2,
+     .    efdisco,efionco2,efiono2,efionn2,
+     .    efionco,efiono3p,efionn,
+     .    efionno,efionh,
+     .    fluxtop,ct1,ct2,p1,p2
+
       implicit none
 
 
 c     common variables and constants
-      include 'param.h'
-      include 'param_v4.h'
-c      include 'datafile.h'
       include "clesphys.h"
  
@@ -19,5 +30,5 @@
       real       nada
       character*13  filename
-      character (len=100) :: datafile="/u/forget/WWW/datagcm/datafile"
+      character (len=100) :: datafile="HIGHATM"
       
 c*************************+PROGRAM STARTS**************************
Index: trunk/LMDZ.VENUS/libf/phyvenus/param_v4.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/param_v4.h	(revision 2463)
+++ 	(revision )
@@ -1,191 +1,0 @@
-!**********************************************************************
-
-!	param.cmn
-	
-!	fgg     mar 02          first version
-!**********************************************************************
-
-
-
-      common/uv/ co2crsc195,co2crsc295,crscabsi2,                                          &
-     &     e107,date_e107,e107_tab,freccen               
-
-      real co2crsc195(9)
-      real co2crsc295(9),crscabsi2(13,16)
-      real freccen(36)          !representative wavelenght
-      real e107,date_e107(669),e107_tab(669)
-
-
-
-      common/fotoquimica/  ch2, ch3, ch4, ch5, ch7,ch9,ch10,ch11,                &
-     &                     ch13,ch14,ch15,ch18,ch19,ch20,ch21,ch22,              &
-     &                     ch23,ch24,ch30,ch31,ch32,ch33,ch34,ch35,              &
-     &                     ch36,ch37,ch38,ch39,ch40,ch41,ch42,ch43,              &
-     &                     ch45,ch46,ch47,ch48,ch49,ch50,ch55,ch56,              &
-     &                     ch57,ch58,ch59,ch62,ch63,ch64,ch65,                   &
-     &                     ch66,ch67,ch68,ch69,ch70,ch71,ch72,                   &
-     &                     ch73,ch74,ch75,ch76,ch85,ch86,ch87                   
-!     &                     jdistot,jdistot_b,jion
-
-      real*8 ch2, ch3, ch4, ch5, ch7,ch9,ch10,ch11,ch13,ch14,ch15,ch18
-      real*8 ch19,ch20,ch21,ch22,ch23,ch24,ch30,ch31,ch32,ch33,ch34
-      real*8 ch35,ch36,ch37,ch38,ch39,ch40,ch41,ch42,ch43,ch45
-      real*8 ch46,ch47,ch48,ch49,ch50,ch55,ch56,ch57,ch58,ch59,ch62 
-      real*8 ch63,ch64,ch65,ch66,ch67,ch68,ch69,ch70,ch71
-      real*8 ch72,ch73,ch74,ch75,ch76,ch85,ch86,ch87
-
-						  !reaction rates
-!      real jdistot(nabs,nlayermx)    !photodissociation rates
-!      real jdistot_b(nabs,nlayermx)
-     
-!      real jion(nabs,nlayermx,4)
-
-
-
-!      common/prod/ Pco2,Po2,Po3p,Pco,Ph,Poh,Pho2,Ph2,Ph2o,Po1d,Ph2o2,          &
-!     &             Po3,Pn,Pno,Pno2,Pn2,Pn2d,                                   &
-!     &             Pco2plus,Poplus,Po2plus,Pelect,Pcoplus,Pcplus,              &
-!     &             Pnplus, Pnoplus, Pn2plus, Phplus, Phco2plus,                &
-!     &  	   Pco2tot,Po2tot,Po3ptot,Pcotot,Phtot,Pohtot,                 &
-!     &             Pho2tot,Ph2tot,Ph2otot,Po1dtot,Ph2o2tot,                    &
-!     &             Po3tot,Pntot,Pnotot,Pno2tot,Pn2tot,Pn2dtot,                 &
-!     &             Pcoplustot,Pcplustot,                                       &
-!     &             Pco2plustot,Poplustot,Po2plustot,Pelecttot,                 &
-!     &             Pnplustot,Pnoplustot,Pn2plustot,Phplustot,                  &
-!     &             Phco2plustot  
-
-	
-!      real*8 Pco2(nlayermx,nreact),Po2(nlayermx,nreact)
-!      real*8 Po3p(nlayermx,nreact)
-!      real*8 Pco(nlayermx,nreact),Ph(nlayermx,nreact)
-!      real*8 Poh(nlayermx,nreact)
-!      real*8 Pho2(nlayermx,nreact),Ph2(nlayermx,nreact)
-!      real*8 Ph2o(nlayermx,nreact)
-!      real*8 Po1d(nlayermx,nreact),Ph2o2(nlayermx,nreact)
-!      real*8 Po3(nlayermx,nreact),Pn(nlayermx,nreact)
-!      real*8 Pno(nlayermx,nreact)
-!      real*8 Pno2(nlayermx,nreact),Pn2(nlayermx,nreact)
-!      real*8 Pn2d(nlayermx,nreact)
-!      real*8 Pco2plus(nlayermx,nreact),Poplus(nlayermx,nreact)
-!      real*8 Po2plus(nlayermx,nreact), Pelect(nlayermx,nreact)
-!      real*8 Pcoplus(nlayermx,nreact),Pcplus(nlayermx,nreact)
-!      real*8 Pnplus(nlayermx,nreact),Pnoplus(nlayermx,nreact)
-!      real*8 Pn2plus(nlayermx,nreact),Phplus(nlayermx,nreact)
-!      real*8 Phco2plus(nlayermx,nreact)
-!      real*8 Pco2tot(nlayermx),Po2tot(nlayermx)
-!      real*8 Po3ptot(nlayermx),Pcotot(nlayermx)
-!      real*8 Phtot(nlayermx)
-!      real*8 Pohtot(nlayermx),Pho2tot(nlayermx)
-!      real*8 Ph2tot(nlayermx),Ph2otot(nlayermx)
-!      real*8 Po1dtot(nlayermx),Ph2o2tot(nlayermx)
-!      real*8 Po3tot(nlayermx),Pntot(nlayermx),Pnotot(nlayermx)
-!      real*8 Pno2tot(nlayermx),Pn2tot(nlayermx),Pn2dtot(nlayermx)
-!      real*8 Pco2plustot(nlayermx), Poplustot(nlayermx)
-!      real*8 Po2plustot(nlayermx),Pelecttot(nlayermx)
-!      real*8 Pcoplustot(nlayermx), Pcplustot(nlayermx)
-!      real*8 Pnplustot(nlayermx), Pnoplustot(nlayermx)
-!      real*8 Pn2plustot(nlayermx),Phplustot(nlayermx)
-!      real*8 Phco2plustot(nlayermx)
-
-
-
-!      common/loss/ Lco2,Lo2,Lo3p,Lco,Lh,Loh,Lho2,Lh2,Lh2o,Lo1d,Lh2o2,         &
-!     &             Lo3,Ln,Lno,Lno2,Ln2,Ln2d,Lelect,                           &
-!     &             Lco2plus,Lo2plus,Loplus, Lcoplus, Lcplus,                  &
-!     &             Lnplus,Lnoplus,Ln2plus, Lhplus, Lhco2plus,                 &
-!     &             Lco2tot,Lo2tot,Lo3ptot,Lcotot,Lhtot,Lohtot,                &
-!     &             Lho2tot,Lh2tot,Lh2otot,Lo1dtot,Lh2o2tot,                   &
-!     &             Lo3tot,Lntot,Lnotot,Lno2tot,Ln2tot,Ln2dtot,                &
-!     &             Lcoplustot,Lcplustot,Lco2plustot,Loplustot,                &
-!     &             Lo2plustot,Lelecttot,Lnplustot,Lnoplustot,                 &
-!     &             Ln2plustot,Lhplustot,Lhco2plustot
-
-	
-!      real*8 Lco2(nlayermx,nreact),Lo2(nlayermx,nreact)
-!      real*8 Lo3p(nlayermx,nreact)
-!      real*8 Lco(nlayermx,nreact),Lh(nlayermx,nreact)
-!      real*8 Loh(nlayermx,nreact)
-!      real*8 Lho2(nlayermx,nreact),Lh2(nlayermx,nreact)
-!      real*8 Lh2o(nlayermx,nreact)
-!      real*8 Lo1d(nlayermx,nreact),Lh2o2(nlayermx,nreact)
-!      real*8 Lo3(nlayermx,nreact),Ln(nlayermx,nreact)
-!      real*8 Lno(nlayermx,nreact)
-!      real*8 Lno2(nlayermx,nreact),Ln2(nlayermx,nreact)
-!      real*8 Ln2d(nlayermx,nreact)
-!      real*8 Lco2plus(nlayermx,nreact), Loplus(nlayermx,nreact)
-!      real*8 Lo2plus(nlayermx,nreact), Lelect(nlayermx,nreact)
-!      real*8 Lcoplus(nlayermx,nreact), Lcplus(nlayermx,nreact)
-!      real*8 Lnplus(nlayermx,nreact), Lnoplus(nlayermx,nreact)
-!      real*8 Ln2plus(nlayermx,nreact), Lhplus(nlayermx,nreact)
-!      real*8 Lhco2plus(nlayermx,nreact)
-!      real*8 Lco2tot(nlayermx),Lo2tot(nlayermx)
-!      real*8 Lo3ptot(nlayermx),Lcotot(nlayermx)
-!      real*8 Lhtot(nlayermx),Lohtot(nlayermx)
-!      real*8 Lho2tot(nlayermx),Lh2tot(nlayermx)
-!      real*8 Lh2otot(nlayermx),Lo1dtot(nlayermx)
-!      real*8 Lh2o2tot(nlayermx)
-!      real*8 Lo3tot(nlayermx),Lntot(nlayermx),Lnotot(nlayermx)
-!      real*8 Lno2tot(nlayermx),Ln2tot(nlayermx)
-!      real*8 Ln2dtot(nlayermx)
-!      real*8 Lco2plustot(nlayermx), Loplustot(nlayermx)
-!      real*8 Lo2plustot(nlayermx), Lelecttot(nlayermx)
-!      real*8 Lcoplustot(nlayermx), Lcplustot(nlayermx)
-!      real*8 Lnplustot(nlayermx), Lnoplustot(nlayermx)
-!      real*8 Ln2plustot(nlayermx), Lhplustot(nlayermx)
-!      real*8 Lhco2plustot(nlayermx)
-
-
-
-!      common/vida/ tminco2,tmino2,tmino3p,tminco,tminh,tminoh,                &
-!     &      tminho2,tminh2,tminh2o,tmino1d,tminh2o2,tmino3,tminn,             &
-!     &      tminno,tminno2,tminn2,tminn2d,                                    &
-!     &      tminco2plus,tminoplus,tmino2plus,tmincoplus,                      &
-!     &      tmincplus,tminnplus,tminnoplus,                                   &
-!     &      tminn2plus,tminhplus,tminhco2plus
-
-
-
-!    	real*8 tminco2(nlayermx),tmino2(nlayermx)
-!        real*8 tmino3p(nlayermx),tminco(nlayermx)
-!    	real*8 tminh(nlayermx),tminoh(nlayermx)
-!        real*8 tminho2(nlayermx),tminh2(nlayermx)
-!    	real*8 tminh2o(nlayermx),tmino1d(nlayermx),tminh2o2(nlayermx)
-!        real*8 tmino3(nlayermx),tminn(nlayermx),tminno(nlayermx)
-!        real*8 tminno2(nlayermx),tminn2(nlayermx),tminn2d(nlayermx)
-!        real*8 tminco2plus(nlayermx),tminoplus(nlayermx)
-!        real*8 tmino2plus(nlayermx)
-!        real*8 tmincoplus(nlayermx)
-!        real*8 tmincplus(nlayermx)
-!        real*8 tminnplus(nlayermx),tminnoplus(nlayermx)
-!        real*8 tminn2plus(nlayermx)
-!        real*8 tminhplus(nlayermx)
-!        real*8 tminhco2plus(nlayermx)
-
-
-!      common/flujoctes/ fluxtophr,ct1,p1,ct2,p2   
- 
-!        real fluxtophr(ninter)
-!        real ct1(ninter),p1(ninter),ct2(ninter),p2(ninter)       
-
-!      common/phdisef/ efdisco2, efdiso2, efdish2, efdish2o,                   &
-!     & efdish2o2,efdiso3,efdiso,efdisn,efdish,efdisno,efdisn2,                &
-!     & efdisno2,efdisco,efionco2,efionn2,efionco,                             &
-!     & efiono3p,efionn,efionno,efionh
-  
-
-!        real efdisco2(ninter), efdiso2(ninter), efdish2o(ninter)
-!        real efdish2o2(ninter), efdish2(ninter), efdiso3(ninter)
-!        real efdiso(ninter), efdisn(ninter), efdish(ninter)
-!        real efdisno(ninter), efdisn2(ninter), efdisno2(ninter)
-!        real efdisco(ninter)
-!        real efionco2(ninter,4)
-!        real efionn2(ninter,2)
-!        real efionco(ninter,3)
-!        real efiono3p(ninter),efionn(ninter)
-!        real efionno(ninter),efionh(ninter)
-
-
-
-
-
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/param_v4_h.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/param_v4_h.F90	(revision 2464)
+++ trunk/LMDZ.VENUS/libf/phyvenus/param_v4_h.F90	(revision 2464)
@@ -0,0 +1,431 @@
+MODULE param_v4_h
+
+   IMPLICIT NONE
+
+   integer, parameter :: ninter=36
+   integer, parameter :: nabs=13
+   integer, parameter :: nz2=253
+   integer, parameter :: ninter2=16
+   real*8, parameter :: kboltzman = 1.381e-16
+   real*8, parameter :: n_avog = 6.023e23
+   real*8, parameter :: gg = 6.67259e-8
+   real*8, parameter :: masa = 6.4163e26
+   real*8, parameter :: radio = 3390.
+   integer, parameter :: nreact=93
+   integer, parameter :: tapas=42
+
+   real crscabsi2(nabs,16)    !cross section
+   real c1_16(nz2,16)   !Col. int. 1 (o2+o+h2+n) (cm^-2)
+   real c17_24(nz2)     !Col. int. 17-24 (co2+o2+n2+no+co+no2) (cm^-2)
+   real c25_29(nz2)     !Col. int. 25-29 (co2+o2+h2o+h2o2+no+co+no2) (cm^-2)
+   real c30_31(nz2)     !Col. int. 30-31 (co2+o2+h2o+h2o2+no+no2)
+   real c32(nz2)        !col. int. 32 (co2+h2o2+o2+no+no2) (cm^-2)
+   real c33(nz2)        !col. int. 33 (h2o2+o2+no2) (cm^-2)
+   real c34(nz2)        !col. int. 34 (h2o2+o2+o3+no2) (cm^-2)
+   real c35(nz2)        !col. int. 35 (h2o2+o3+no2) (cm^-2)
+   real c36(nz2)        !col. int. 36 (o3+no2) (cm^-2)
+   real co2crsc195(9)
+   real co2crsc295(9)
+   real t0(nz2)
+   real fluxtop(ninter)
+   real freccen(ninter)          !representative wavelenght
+   real jabsifotsintpar(nz2,nabs,ninter)
+   real e107,date_e107(669),e107_tab(669)
+   real coefit0(ninter,nabs),coefit1(ninter,nabs)
+   real coefit2(ninter,nabs)
+   real coefit3(ninter,nabs),coefit4(ninter,nabs)
+
+   !reaction rates
+   real*8 ch2, ch3, ch4, ch5, ch7,ch9,ch10,ch11,ch13,ch14,ch15,ch18
+   real*8 ch19,ch20,ch21,ch22,ch23,ch24,ch30,ch31,ch32,ch33,ch34
+   real*8 ch35,ch36,ch37,ch38,ch39,ch40,ch41,ch42,ch43,ch45
+   real*8 ch46,ch47,ch48,ch49,ch50,ch55,ch56,ch57,ch58,ch59,ch62
+   real*8 ch63,ch64,ch65,ch66,ch67,ch68,ch69,ch70,ch71
+   real*8 ch72,ch73,ch74,ch75,ch76,ch85,ch86,ch87
+   real*8 rcoef(61,3)
+
+   real fluxtophr(ninter)
+   real ct1(ninter),p1(ninter),ct2(ninter),p2(ninter)
+
+   real efdisco2(ninter), efdiso2(ninter), efdish2o(ninter)
+   real efdish2o2(ninter), efdish2(ninter), efdiso3(ninter)
+   real efdiso(ninter), efdisn(ninter), efdish(ninter)
+   real efdisno(ninter), efdisn2(ninter), efdisno2(ninter)
+   real efdisco(ninter)
+   real efionco2(ninter,4)
+   real efiono2(ninter,2)
+   real efionn2(ninter,2)
+   real efionco(ninter,3)
+   real efiono3p(ninter),efionn(ninter)
+   real efionno(ninter),efionh(ninter)
+
+   !photodissociation rates
+   REAL,SAVE,ALLOCATABLE :: jfotsout(:,:,:)
+   REAL,SAVE,ALLOCATABLE :: jdistot(:,:)
+   REAL,SAVE,ALLOCATABLE :: jdistot_b(:,:)
+   REAL,SAVE,ALLOCATABLE :: jion(:,:,:)
+
+   REAL*8,SAVE,ALLOCATABLE :: Pco2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Po2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Po3p(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pco(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Poh(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pho2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph2o(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Po1d(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph2o2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Po3(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pno(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pno2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn2d(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pco2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Poplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Po2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pelect(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pcoplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pcplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pnplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pnoplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Phplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Phco2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Pco2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Po2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Po3ptot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pcotot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Phtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pohtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pho2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph2otot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Po1dtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Ph2o2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Po3tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pntot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pnotot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pno2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn2dtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pco2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Poplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Po2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pelecttot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pcoplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pcplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pnplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pnoplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Pn2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Phplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Phco2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lco2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo3p(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lco(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Loh(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lho2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh2o(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo1d(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh2o2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo3(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lno(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lno2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln2(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln2d(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lco2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Loplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lelect(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lcoplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lcplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lnplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lnoplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lhplus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lhco2plus(:,:)
+   REAL*8,SAVE,ALLOCATABLE :: Lco2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo3ptot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lcotot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lhtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lohtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lho2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh2otot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo1dtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lh2o2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo3tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lntot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lnotot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lno2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln2tot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln2dtot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lco2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Loplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lo2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lelecttot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lcoplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lcplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lnplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lnoplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Ln2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lhplustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: Lhco2plustot(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminco2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmino2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmino3p(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminco(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminh(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminoh(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminho2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminh2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminh2o(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmino1d(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminh2o2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmino3(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminn(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminno(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminno2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminn2(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminn2d(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminco2plus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminoplus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmino2plus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmincoplus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tmincplus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminnplus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminnoplus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminn2plus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminhplus(:)
+   REAL*8,SAVE,ALLOCATABLE :: tminhco2plus(:)
+
+   CONTAINS
+
+      SUBROUTINE fill_data_thermos
+
+      IMPLICIT NONE
+
+      ! data for the UV heating tabulation
+
+      crscabsi2(1,1:16) = (/ 5.61031E-19,1.59677E-18,4.7072E-18,&
+          1.48254e-17,2.07445e-17,2.573e-17,2.901e-17,3.083e-17,&
+          3.217e-17,3.539e-17,3.658e-17,3.63e-17,3.41239e-17,&
+          2.71019e-17,4.93677e-17,1.64e-17 /) 
+
+      crscabsi2(2,1:16) = (/ 0.27250E-18,0.11650E-17,0.39250E-17,&
+          0.10630E-16,0.15590E-16,0.17180E-16,0.19270E-16,0.22860E-16,&
+          0.24270E-16,0.24440E-16,0.25020E-16,0.26600E-16,0.25400E-16,&
+          0.35800E-16,0.25590E-16,0.16740E-16 /)
+
+      crscabsi2(3,1:16) = (/ 0.2776E-18,0.9792E-18,0.3313E-17,&
+          0.6621E-17,0.8481E-17,0.9146E-17,0.9414E-17,0.1039E-16,&
+          0.1012E-16,0.1033E-16,0.1033E-16,0.1033E-16,0.8268E-17,&
+          0.6563E-17,0.3506E-17,0.3470E-17 /) 
+
+      crscabsi2(5,1:16) = (/ .5E-20,.1077607E-19,.5670491E-19,&
+          .3322716E-18,.1054509E-17,.1700005E-17,.3171188E-17,&
+          .4734241E-17,.5108741E-17,.6022236E-17,.6741537E-17,&
+          .7277079E-17,.9070787E-17,.9708916E-17,.4026281E-17,0.0 /)
+
+      crscabsi2(8,1:16) = (/ 0.0, 7.44175e-19, 2.23167e-18,&
+          8.46200e-18,1.18275e-17,1.54900e-17,2.32475e-17,2.41373e-17,&
+          2.55482e-17,2.38431e-17,2.28600e-17,2.35067e-17,2.56000e-17,&
+          2.64636e-17,2.86260e-17,3.26561e-17 /)
+
+      crscabsi2(9,1:16) = (/ 3.48182e-20,3.37038e-19,1.03077e-18,&
+          4.01364e-18,6.45e-18,7.8e-18,1.0e-17,1.13500e-17,1.15500e-17,&
+          1.18000e-17,1.17500e-17,1.16000e-17,1.28667e-17,1.18500e-17,&
+          1.11000e-17,9.50000e-18 /)
+
+      crscabsi2(10,1:16) = (/ 0.0,9.39833e-19,2.87714e-18,&
+          9.66900e-18,1.37063e-17,1.61820e-17,2.30450e-17,2.63373e-17,&
+          2.63773e-17,2.67677e-17,2.64100e-17,2.53000e-17,2.18100e-17,&
+          2.04941e-17,2.28160e-17,2.93550e-17 /)
+
+      crscabsi2(11,1:16) = (/ 0.0,9.58555e-19,2.52767e-18,&
+          8.29700e-18,1.21850e-17,1.40500e-17,1.97025e-17,2.12018e-17,&
+          2.14673e-17,2.20331e-17,2.21500e-17,2.21600e-17,2.33200e-17,&
+          2.67800e-17,2.56400e-17,3.58561e-17 /)
+
+      crscabsi2(12,1:16) = (/ 0.0,1.0e-20,2.5e-20,1.30400e-19,&
+          2.93800e-19,4.36000e-19,8.49400e-19,1.29400e-18,1.40500e-18,&
+          1.67600e-18,1.93400e-18,2.12200e-18,2.75800e-18,3.48400e-18,&
+          4.17200e-18,5.26000e-18 /)
+
+      crscabsi2(13,1:16) = (/ 0.0,1.60e-18,4.99111e-18,1.48496e-17,&
+          2.17395e-17,2.55857e-17,2.87754e-17,3.65571e-17,3.85691e-17,&
+          4.16286e-17,4.15117e-17,4.05901e-17,3.64000e-17,2.99670e-17,&
+          2.46796e-17,2.51789e-17 /)
+
+      freccen(1:ninter)=(/ 3.4,7.5,14.5,23.0,30.3,34.1,&
+           49.6,50.5,52.5,56.0,&
+           59.0,61.5,68.7,73.1,78.4,83.1,92.4,97.5,99.3,100.1,100.7,102.1,&
+           104.5,116.8,121.3,127.0,130.6,153.7,162.8,171.4,&
+           195.6,206.3,222.0,236.0,289.0,600. /)
+
+      co2crsc195(1:9)=(/ 2.05864e-17,5.90557e-20,3.1027e-19,6.70653e-19,&
+           4.55132e-19,8.87122e-20,1.32138e-20,7.22244e-23,2.88002e-26 /)
+
+      co2crsc295(1:9)=(/2.05897e-17,6.71104e-20,3.45509e-19,7.45711e-19,&
+           4.82752e-19,1.11594e-19,1.98308e-20,1.3853e-22,2.1414e-25 /)
+
+      END SUBROUTINE fill_data_thermos
+
+      SUBROUTINE allocate_param_thermos(nlayer)
+
+        IMPLICIT NONE
+
+        INTEGER :: nlayer
+        allocate(jdistot(nabs,nlayer))    
+        allocate(jdistot_b(nabs,nlayer))
+        allocate(jion(nabs,nlayer,4))
+        allocate(jfotsout(ninter,nabs,nlayer))
+        allocate(Pco2(nlayer,nreact))
+        allocate(Po2(nlayer,nreact))
+        allocate(Po3p(nlayer,nreact))
+        allocate(Pco(nlayer,nreact))
+        allocate(Ph(nlayer,nreact))
+        allocate(Poh(nlayer,nreact))
+        allocate(Pho2(nlayer,nreact))
+        allocate(Ph2(nlayer,nreact))
+        allocate(Ph2o(nlayer,nreact))
+        allocate(Po1d(nlayer,nreact))
+        allocate(Ph2o2(nlayer,nreact))
+        allocate(Po3(nlayer,nreact))
+        allocate(Pn(nlayer,nreact))
+        allocate(Pno(nlayer,nreact))
+        allocate(Pno2(nlayer,nreact))
+        allocate(Pn2(nlayer,nreact))
+        allocate(Pn2d(nlayer,nreact))
+        allocate(Pco2plus(nlayer,nreact))
+        allocate(Poplus(nlayer,nreact))
+        allocate(Po2plus(nlayer,nreact))
+        allocate(Pelect(nlayer,nreact))
+        allocate(Pcoplus(nlayer,nreact))
+        allocate(Pcplus(nlayer,nreact))
+        allocate(Pnplus(nlayer,nreact))
+        allocate(Pnoplus(nlayer,nreact))
+        allocate(Pn2plus(nlayer,nreact))
+        allocate(Phplus(nlayer,nreact))
+        allocate(Phco2plus(nlayer,nreact))
+        allocate(Pco2tot(nlayer))
+        allocate(Po2tot(nlayer))
+        allocate(Po3ptot(nlayer))
+        allocate(Pcotot(nlayer))
+        allocate(Phtot(nlayer))
+        allocate(Pohtot(nlayer))
+        allocate(Pho2tot(nlayer))
+        allocate(Ph2tot(nlayer))
+        allocate(Ph2otot(nlayer))
+        allocate(Po1dtot(nlayer))
+        allocate(Ph2o2tot(nlayer))
+        allocate(Po3tot(nlayer))
+        allocate(Pntot(nlayer))
+        allocate(Pnotot(nlayer))
+        allocate(Pno2tot(nlayer))
+        allocate(Pn2tot(nlayer))
+        allocate(Pn2dtot(nlayer))
+        allocate(Pco2plustot(nlayer))
+        allocate(Poplustot(nlayer))
+        allocate(Po2plustot(nlayer))
+        allocate(Pelecttot(nlayer))
+        allocate(Pcoplustot(nlayer))
+        allocate(Pcplustot(nlayer))
+        allocate(Pnplustot(nlayer))
+        allocate(Pnoplustot(nlayer))
+        allocate(Pn2plustot(nlayer))
+        allocate(Phplustot(nlayer))
+        allocate(Phco2plustot(nlayer))
+        allocate(Lco2(nlayer,nreact))
+        allocate(Lo2(nlayer,nreact))
+        allocate(Lo3p(nlayer,nreact))
+        allocate(Lco(nlayer,nreact))
+        allocate(Lh(nlayer,nreact))
+        allocate(Loh(nlayer,nreact))
+        allocate(Lho2(nlayer,nreact))
+        allocate(Lh2(nlayer,nreact))
+        allocate(Lh2o(nlayer,nreact))
+        allocate(Lo1d(nlayer,nreact))
+        allocate(Lh2o2(nlayer,nreact))
+        allocate(Lo3(nlayer,nreact))
+        allocate(Ln(nlayer,nreact))
+        allocate(Lno(nlayer,nreact))
+        allocate(Lno2(nlayer,nreact))
+        allocate(Ln2(nlayer,nreact))
+        allocate(Ln2d(nlayer,nreact))
+        allocate(Lco2plus(nlayer,nreact))
+        allocate(Loplus(nlayer,nreact))
+        allocate(Lo2plus(nlayer,nreact))
+        allocate(Lelect(nlayer,nreact))
+        allocate(Lcoplus(nlayer,nreact))
+        allocate(Lcplus(nlayer,nreact))
+        allocate(Lnplus(nlayer,nreact))
+        allocate(Lnoplus(nlayer,nreact))
+        allocate(Ln2plus(nlayer,nreact))
+        allocate(Lhplus(nlayer,nreact))
+        allocate(Lhco2plus(nlayer,nreact))
+        allocate(Lco2tot(nlayer))
+        allocate(Lo2tot(nlayer))
+        allocate(Lo3ptot(nlayer))
+        allocate(Lcotot(nlayer))
+        allocate(Lhtot(nlayer))
+        allocate(Lohtot(nlayer))
+        allocate(Lho2tot(nlayer))
+        allocate(Lh2tot(nlayer))
+        allocate(Lh2otot(nlayer))
+        allocate(Lo1dtot(nlayer))
+        allocate(Lh2o2tot(nlayer))
+        allocate(Lo3tot(nlayer))
+        allocate(Lntot(nlayer))
+        allocate(Lnotot(nlayer))
+        allocate(Lno2tot(nlayer))
+        allocate(Ln2tot(nlayer))
+        allocate(Ln2dtot(nlayer))
+        allocate(Lco2plustot(nlayer))
+        allocate(Loplustot(nlayer))
+        allocate(Lo2plustot(nlayer))
+        allocate(Lelecttot(nlayer))
+        allocate(Lcoplustot(nlayer))
+        allocate(Lcplustot(nlayer))
+        allocate(Lnplustot(nlayer))
+        allocate(Lnoplustot(nlayer))
+        allocate(Ln2plustot(nlayer))
+        allocate(Lhplustot(nlayer))
+        allocate(Lhco2plustot(nlayer))
+        allocate(tminco2(nlayer))
+        allocate(tmino2(nlayer))
+        allocate(tmino3p(nlayer))
+        allocate(tminco(nlayer))
+        allocate(tminh(nlayer))
+        allocate(tminoh(nlayer))
+        allocate(tminho2(nlayer))
+        allocate(tminh2(nlayer))
+        allocate(tminh2o(nlayer))
+        allocate(tmino1d(nlayer))
+        allocate(tminh2o2(nlayer))
+        allocate(tmino3(nlayer))
+        allocate(tminn(nlayer))
+        allocate(tminno(nlayer))
+        allocate(tminno2(nlayer))
+        allocate(tminn2(nlayer))
+        allocate(tminn2d(nlayer))
+        allocate(tminco2plus(nlayer))
+        allocate(tminoplus(nlayer))
+        allocate(tmino2plus(nlayer))
+        allocate(tmincoplus(nlayer))
+        allocate(tmincplus(nlayer))
+        allocate(tminnplus(nlayer))
+        allocate(tminnoplus(nlayer))
+        allocate(tminn2plus(nlayer))
+        allocate(tminhplus(nlayer))
+        allocate(tminhco2plus(nlayer))
+
+      END SUBROUTINE allocate_param_thermos
+
+END MODULE param_v4_h
Index: trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 2464)
@@ -1457,5 +1457,5 @@
    cicol  = 0.
 
-   do i = 1,nztable - 1
+   do i = 1,nztable-1
       if (table_colair(i) < col(iz)) then
          cicol = (log(col(iz)) - log(table_colair(i)))           &
Index: trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2464)
@@ -69,4 +69,5 @@
       USE chemparam_mod
       USE conc
+      USE param_v4_h
       USE compo_hedin83_mod2
 !      use ieee_arithmetic
@@ -257,5 +258,4 @@
       EXTERNAL nlthermeq
       EXTERNAL euvheat
-      EXTERNAL param_read
       EXTERNAL param_read_e107
       EXTERNAL conduction
@@ -356,4 +356,5 @@
 c
       REAL :: tr_seri(klon,klev,nqmax)
+      REAL :: tr_hedin(klon,klev,nqmax)
       REAL :: d_tr(klon,klev,nqmax)
 
@@ -411,6 +412,7 @@
 
 c======================================================================
+c======================================================================
 c INITIALISATIONS
-c================
+c======================================================================
 
       modname = 'physiq'
@@ -437,4 +439,5 @@
       END IF 
       
+c======================================================================
 c PREMIER APPEL SEULEMENT
 c========================
@@ -584,6 +587,7 @@
 
          if (callthermos) then
-            if(solvarmod.eq.0) call param_read
-            if(solvarmod.eq.1) call param_read_e107 
+            call fill_data_thermos
+            call allocate_param_thermos(klev)
+            call param_read_e107 
          endif
 
@@ -597,5 +601,4 @@
             rho(ig,j)=pplay(ig,j)*mmean(ig,j)*1e-3/(rnew(ig,j)*t(ig,j))
           enddo
-c        stop
 
         enddo  
@@ -749,9 +752,12 @@
 
       if(callthermos) then
-         call concentrations2(pplay,t,d_t,qx,nqmax)
+         call concentrations2(pplay,t,qx,nqmax)
       endif
 
+c========================
       ENDIF ! debut
-
+c========================
+
+c======================================================================
 ! ------------------------------------------------------
 !   Initializations done at every physical timestep:
@@ -923,5 +929,5 @@
 c====================================================================
 c Orbite et eclairement
-c====================================================================
+c=======================
 
 c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
@@ -931,5 +937,5 @@
       
       zlongi = 0.0
-      dist   = 0.72  ! en UA
+      dist   = 0.72333  ! en UA
 
 c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite 
@@ -952,8 +958,58 @@
       end if
 
+c======================================================================
+c FIN INITIALISATIONS
+c======================================================================
+c======================================================================
+
+c=======================================================
+! CONDUCTION  and  MOLECULAR VISCOSITY
+c=======================================================
+
+        d_t_conduc(:,:)=0.
+        d_u_molvis(:,:)=0.
+        d_v_molvis(:,:)=0.
+
+        IF (callthermos) THEN
+
+           tsurf(:)=t_seri(:,1)
+           call conduction(klon, klev,pdtphys,
+     $            pplay,paprs,t_seri,
+     $            tsurf,zzlev,zzlay,d_t_conduc)
+
+            call molvis(klon, klev,pdtphys,
+     $            pplay,paprs,t_seri,
+     $            u,tsurf,zzlev,zzlay,d_u_molvis)
+
+            call molvis(klon, klev, pdtphys,
+     $            pplay,paprs,t_seri,
+     $            v,tsurf,zzlev,zzlay,d_v_molvis)
+
+            DO k=1,klev
+               DO ig=1,klon
+                   t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
+                   u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
+                   v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
+               ENDDO
+            ENDDO
+        ENDIF
+
+c====================================================================
+c    Compute mean mass, cp and R :
+c------------------------------------
+
+      if(callthermos) then
+         call concentrations2(pplay,t_seri,tr_seri, nqmax)
+      endif
+
+
+c=======================================================
+! CHEMISTRY AND MICROPHYSICS
+c=======================================================
+
       if (iflag_trac.eq.1) then
 !====================================================================
 ! Case 1: pseudo-chemistry with relaxation toward fixed profile
-!====================================================================
+!=========
        if (tr_scheme.eq.1) then
 
@@ -968,8 +1024,7 @@
 ! However, the variable 'source' could be used in physiq
 ! so the call to phytrac_emiss could be to initialise it.
-!====================================================================
-
-         call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
-     I                   debut,lafin,nqmax,
+!=========
+
+         call phytrac_emiss (debut,lafin,nqmax,
      I                   nlon,nlev,dtime,paprs,
      I                   latitude_deg,longitude_deg,
@@ -985,5 +1040,5 @@
 !         nbapp_chem = 24000 => chempas = 4 => zctime = 420 s
 !         nbapp_chem = 12000 => chempas = 8 => zctime = 840 s
-!====================================================================
+!=========
 
          nbapp_chem = 24000
@@ -1131,6 +1186,6 @@
 
             do iq = 1, nqmax - nmicro
-               tr_seri(:,:,iq) = tr_seri(:,:,iq) 
-     $                         + d_tr_chem(:,:,iq)*zctime
+               tr_seri(:,:,iq) = max(tr_seri(:,:,iq) 
+     $                         + d_tr_chem(:,:,iq)*zctime,1.e-30)
             end do
 
@@ -1144,6 +1199,6 @@
             else if (cl_scheme == 2) then
                do iq = nqmax-nmicro+1,nqmax
-                  tr_seri(:,:,iq) = tr_seri(:,:,iq) 
-     $                            + d_tr_sed(:,:,iq)*zctime
+                  tr_seri(:,:,iq) = max(tr_seri(:,:,iq) 
+     $                            + d_tr_sed(:,:,iq)*zctime,1.e-30)
                end do
             end if  ! cl_scheme
@@ -1152,5 +1207,5 @@
          end if     ! mod(itap,chempas)  <------- end of chemistry supercycling
 
-!====================================================================
+!=========
 ! End Case 3: Full chemistry and/or clouds.
 !====================================================================
@@ -1300,7 +1355,4 @@
 !      call writefield_phy('physiq_d_ts',d_ts,1)
 
-c
-c Appeler l'ajustement sec
-c
 c===================================================================
 c Convection seche 
@@ -1371,12 +1423,4 @@
       END IF 
 
-c------------------------------------
-c    Compute mean mass, cp and R :
-c------------------------------------
-
-      if(callthermos) then
-         call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax)
-      endif
-
 c====================================================================
 c RAYONNEMENT
@@ -1385,6 +1429,15 @@
 
       dtimerad = dtime*REAL(radpas)  ! pas de temps du rayonnement (s)
-c====================================================================
-
+
+! update mmean
+         if (ok_chem) then
+            mmean(:,:) = 0.
+            do iq = 1,nqmax - nmicro
+               mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)*m_tr(iq)
+               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
+            enddo
+         endif
+
+cc---------------------------------------------
       if (callnlte .or. callthermos) then
          if (ok_chem) then
@@ -1411,7 +1464,15 @@
         IF(callnlte) THEN                                 
             if(nltemodel.eq.0.or.nltemodel.eq.1) then
-                CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
-     $                    tr_seri, d_t_nlte)
+c nltecool call not correct...
+c                CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
+c     $                    tr_seri, d_t_nlte)
+                abort_message='nltemodel=0 or 1 should not be used...'
+                call abort_physic(modname,abort_message,1)
             else if(nltemodel.eq.2) then                                
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! HEDIN instead of compo for this calculation
+!              call compo_hedin83_mod(pplay,rmu0,   
+!    $                 co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)     
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                CALL nlte_tcool(klon,klev,pplay*9.869e-6,              
      $               t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
@@ -1436,5 +1497,5 @@
      $        CALL nlthermeq(klon, klev, paprs, pplay) 
 
-c
+cc---------------------------------------------
 c       LTE radiative transfert / solar / IR matrix
 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1455,4 +1516,5 @@
 c       heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
 
+cc---------------------------------------------
 c       CO2 near infrared absorption
 c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1460,9 +1522,23 @@
         d_t_nirco2(:,:)=0.
         if (callnirco2) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! HEDIN instead of compo for this calculation
+!          call compo_hedin83_mod(pplay,rmu0,   
+!    $                 co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
+!          tr_hedin(:,:,:)=tr_seri(:,:,:)
+!          tr_hedin(:,:,i_co2)=co2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_co2)
+!          tr_hedin(:,:,i_co) = covmr_gcm(:,:)/mmean(:,:)*m_tr(i_co)
+!          tr_hedin(:,:,i_o)  =  ovmr_gcm(:,:)/mmean(:,:)*m_tr(i_o)
+!          tr_hedin(:,:,i_n2) = n2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_n2)
+!          call nirco2abs (klon, klev, pplay, dist, nqmax, tr_hedin,
+!    .                 rmu0, fract, d_t_nirco2)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
      .                 rmu0, fract, d_t_nirco2)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         endif
 
 
+cc---------------------------------------------
 c          Net atmospheric radiative heating rate (K.s-1)
 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1484,5 +1560,4 @@
 
 cc---------------------------------------------
-
 c          EUV heating rate (K.s-1)
 c          ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1492,10 +1567,13 @@
         IF (callthermos) THEN
 
-c           call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
-c     $          rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm, 
-c     $          covmr_gcm, ovmr_gcm,d_t_euv )
            call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
-     $         rmu0,dtimerad,gmtime,rjourvrai,
+     $         rmu0,dtimerad,gmtime,
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! HEDIN instead of compo for this calculation
+!! cf nlte_tcool and nirco2abs above !!
+!    $         tr_hedin, d_tr, d_t_euv )
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      $         tr_seri, d_tr, d_t_euv )
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                  
            DO k=1,klev
@@ -1508,5 +1586,4 @@
         ENDIF  ! callthermos
 
-c====================================================================
       ENDIF    ! radpas
 c====================================================================
@@ -1523,38 +1600,10 @@
 
       itap   = itap + 1
-
-! CONDUCTION  and  MOLECULAR VISCOSITY
-c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        d_t_conduc(:,:)=0.
-        d_u_molvis(:,:)=0.
-        d_v_molvis(:,:)=0.
-
-        IF (callthermos) THEN
-
-           tsurf(:)=t_seri(:,1)
-           call conduction(klon, klev,pdtphys,
-     $            pplay,paprs,t_seri,
-     $            tsurf,zzlev,zzlay,d_t_conduc)
-
-            call molvis(klon, klev,pdtphys,
-     $            pplay,paprs,t_seri,
-     $            u,tsurf,zzlev,zzlay,d_u_molvis)
-
-            call molvis(klon, klev, pdtphys,
-     $            pplay,paprs,t_seri,
-     $            v,tsurf,zzlev,zzlay,d_v_molvis)
-
-            DO k=1,klev
-               DO ig=1,klon
-                  t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
-                  u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
-                  v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
-               ENDDO
-            ENDDO
-        ENDIF
-
-
+c====================================================================
+
+
+c==============================
 !  --  MOLECULAR DIFFUSION ---
+c==============================
 
           d_q_moldif(:,:,:)=0
@@ -1564,5 +1613,5 @@
              call moldiff_red(klon, klev, nqmax,
      &                   pplay,paprs,t_seri, tr_seri, pdtphys,
-     &                   zzlay,d_t_euv,d_t_conduc,d_q_moldif)
+     &                   d_t_euv,d_t_conduc,d_q_moldif)
 
 
@@ -1572,11 +1621,10 @@
            DO k=1,klev
               DO ig=1,klon
-                tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+ 
-     &                           d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
+                tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+ 
+     &                           d_q_moldif(ig,k,iq)*dtime,1.e-30)
               ENDDO
             ENDDO
            ENDDO
            
-
          ENDIF  ! callthermos & ok_chem
 
@@ -1878,4 +1926,9 @@
       ENDDO
       
+c mise à jour rho,mmean pour sorties
+      if(callthermos) then
+         call concentrations2(pplay,t_seri,tr_seri, nqmax)
+      endif
+
 c------------------------
 c Calcul moment cinetique
@@ -2005,5 +2058,5 @@
          do iq = 1,nqmax - nmicro
             call send_xios_field(tname(iq),
-     $                           qx(:,:,iq)*mmean(:,:)/m_tr(iq))
+     $                           tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
          end do
 
@@ -2012,7 +2065,7 @@
          if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN
             call send_xios_field(tname(i_h2oliq),
-     $             qx(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
+     $             tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
             call send_xios_field(tname(i_h2so4liq),
-     $             qx(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
+     $             tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
             if (ok_sedim) then
                call send_xios_field("Fsedim",fsedim(:,1:klev))
@@ -2022,5 +2075,5 @@
 ! chemical iterations
 
-         call send_xios_field("iter",real(iter))
+         if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
 
       end if
@@ -2031,4 +2084,8 @@
        CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
       ENDIF
+
+!! DEBUG
+!      if (is_master) print*,"itauphy=",itap
+!      if (itap.eq.10) lafin=.true.
 
       if (lafin.and.is_omp_master) then
Index: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 2464)
@@ -15,5 +15,5 @@
 
       use chemparam_mod
-      use conc, only: mmean
+      use conc, only: mmean,rnew
 
       implicit none
@@ -73,4 +73,5 @@
          if (reinit_trac .and. ok_chem) then
 
+!!! in this reinitialisation, trac is VOLUME mixing ratio
             print*, "Tracers are re-initialised"
             trac(:,:,:) = 1.0e-30
@@ -106,6 +107,13 @@
             end if
         
+
+! update mmean
+            mmean(:,:) = 0.
+            do iq = 1,nqmax - nmicro
+               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
+               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
+            enddo
+
 !           convert volume to mass mixing ratio
-
             do iq = 1,nqmax - nmicro
                trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
@@ -269,4 +277,11 @@
 !===================================================================
 
+! update mmean
+            mmean(:,:) = 0.
+            do iq = 1,nqmax - nmicro
+               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
+               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
+            enddo
+
 !     gas phase
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F	(revision 2464)
@@ -4,14 +4,15 @@
 c
 c
-      SUBROUTINE phytrac_emiss (timesimu,
-     I                    debutphy,
+      SUBROUTINE phytrac_emiss (debutphy,
      I                    lafin,
      I                    nqmax,
      I                    nlon,
-     I                    nlev, 
+     I                    nlev,
      I                    pdtphys,
-     I                    paprs,
+     I                    paprs, 
      I                    xlat,xlon,
      O                    tr_seri)
+     
+     
 
 c======================================================================
@@ -28,5 +29,8 @@
 c
 c======================================================================
-      USE infotrac_phy, ONLY: nqtot
+      use infotrac_phy, only: tname, nqtot
+#ifdef CPP_XIOS      
+      use xios_output_mod, only:  send_xios_field
+#endif
       use dimphy
       USE geometry_mod, only: cell_area
@@ -36,4 +40,5 @@
 #include "YOMCST.h"
 #include "clesphys.h"
+     
 c======================================================================
 
@@ -43,5 +48,4 @@
 c   ==========
 
-      real timesimu   ! duree depuis debut simu (s)
       logical debutphy       ! le flag de l'initialisation de la physique
       logical lafin          ! le flag de la fin de la physique
@@ -63,21 +67,29 @@
 cAA ----------------------------
 
-c pour emission volcan
+c pour emission type effusion
       real :: deltatr(klon,klev,nqtot)
 
-      integer,parameter :: nblat=5,nblon=4,nbz=3
-      integer,parameter :: Nemiss=0     ! duree emission (Ed)
-      integer,save :: Nheight(nbz)      ! layer emission
-      real,save :: so2_quantity         ! quantity so2 (kg)
-      real,save :: lat_volcan(nblat),lon_volcan(nblon)
-      real,save :: area_emiss(nblat,nblon)
-      integer,save :: ig_volcan(nblat,nblon)
+
+
+
+      integer,parameter :: nblat=3,nblon=3,nbaire=3,nbflux=2,maxcell=16 
+      integer,parameter :: Nheight=3 ! layer emission (150m)
+      integer,save :: Ncell(nbaire) 
+      real,save :: flux_surface_co2(nbflux) ! flux de CO2 emis (kg/m2/s) 
+      real :: flux(nlon,nqtot) 
+      real,save :: lat_zone(nblat),lon_zone(nblon) 
+      integer,save :: ig_zone(nblat,nblon,nbaire,nbflux,maxcell) 
+      integer,save :: numcell(nblat,nblon,nbaire,nbflux)
+        
 
       INTEGER i, k, it
-      integer ilat,ilon,iz
+      integer ilat,ilon,iaire,iflux,ipos
       real    deltalat,deltalon
-c======================================================================
-
-c EMISSION TRACEURS 
+      
+
+
+c======================================================================
+
+c EMISSIONS TRACEURS 
 
 c---------
@@ -90,6 +102,6 @@
 
         ALLOCATE(M_tr(nqtot))
-        M_tr(:)=64.                 ! SO2
-	
+        M_tr(:)=44.                ! CO2
+
 C=========================================================================
 c Caracteristiques des traceurs emis:
@@ -97,30 +109,33 @@
 
 c nombre total de traceur
-         if (nbz*nblat*nblon .gt. nqtot) then
-            print*, nbz*nblat*nblon, nqtot
+         if (nblat*nblon*nbaire*nbflux+1 .gt. nqtot) then
+            print*, nblat*nblon*nbaire*nbflux+1, nqtot
             write(*,*) "Attention, pas assez de traceurs"
             write(*,*) "le dernier sera bien le dernier"
-         endif
-
-c quantite en kg
-         so2_quantity = 20.*10.**9.
-
-c height (in layer index)
-         Nheight(1) =  6  ! ~ 1 km
-         Nheight(2) = 16  ! ~ 25 km
-         Nheight(3) = 24  ! ~ 50 km
-
-c localisation volcan
-         lat_volcan(1) =  70.
-         lat_volcan(2) =  35.
-         lat_volcan(3) =   0.
-         lat_volcan(4) = -35.
-         lat_volcan(5) = -70.
-         lon_volcan(1) = -125.
-         lon_volcan(2) =  -35.
-         lon_volcan(3) =   55.
-         lon_volcan(4) =  145.
-	 
-	 ig_volcan(ilat,ilon)= 0
+         endif 
+                 
+
+
+
+c flux de CO2 (kg/s/m2)
+         flux_surface_co2(1) = 5.*10.**-9.
+         flux_surface_co2(2) = 5.*10.**-15.
+
+c nombre de cellule pour le cote du carre d'aire 
+         Ncell(1)= 2
+         Ncell(2)= 3
+         Ncell(3)= 4
+      
+
+
+c localisation zone emission
+         lat_zone(1) = 08.
+         lat_zone(2) = -50.
+         lat_zone(3) = 35.
+         lon_zone(1) = -172.
+         lon_zone(2) = -20.
+         lon_zone(3) = 70.
+
+ 
          if ((nbp_lon*nbp_lat)==1) then ! running a 1D simulation
            deltalat=180.
@@ -131,15 +146,28 @@
          endif
 
+         numcell(:,:,:,:)=0
+         ig_zone(:,:,:,:,:)=0
          do i=1,nlon
           do ilat=1,nblat
            do ilon=1,nblon
-            if ((xlat(i).ge.lat_volcan(ilat))
-     &     .and.((xlat(i)-deltalat).lt.lat_volcan(ilat))
-     &     .and.(xlon(i).le.lon_volcan(ilon))
-     &     .and.((xlon(i)+deltalon).gt.lon_volcan(ilon)) ) then
-             ig_volcan(ilat,ilon)= i
-             area_emiss(ilat,ilon) = cell_area(i)
-	     print*,"Lat,lon=",ilat,ilon," OK"
-            end if
+            do iaire=1,nbaire
+
+             if ((xlat(i).ge.lat_zone(ilat))
+     &      .and.((xlat(i)-Ncell(iaire)*deltalat)
+     &      .lt.lat_zone(ilat))
+     &      .and.(xlon(i).le.lon_zone(ilon))
+     &      .and.((xlon(i)+Ncell(iaire)*deltalon)
+     &      .gt.lon_zone(ilon))) then
+
+              do iflux=1,nbflux
+         numcell(ilat,ilon,iaire,iflux)=numcell(ilat,ilon,iaire,iflux)+1
+        ig_zone(ilat,ilon,iaire,iflux,numcell(ilat,ilon,iaire,iflux))= i
+             print*,"Lat,lon,naire,nflux,nlon=",ilat,ilon,iaire,iflux
+     &        ,i," OK"
+              end do
+
+             end if
+
+            end do
            end do
           end do
@@ -148,16 +176,15 @@
 c Reinit des traceurs si necessaire
          if (reinit_trac) then
-	   tr_seri(:,:,:)=0.
-c CAS N2 TRACEUR PASSIF POUR EVALUER PROFIL SOUS 7 KM
-c J'ai mis Nemiss=0 !
+          tr_seri(:,:,:)=0.
+
+
            do i=1,klon
             do k=1,klev
-              tr_seri(i,k,:)=max(min(0.035,
-     &          0.035*(1.-log(paprs(i,k)/6.e6)/log(9.e6/6.e6))),0.)
-            enddo
-           enddo
-c FIN CAS N2 PASSIF
-	 endif
-	 
+              tr_seri(i,k,:)=1.-28/43.44*max(min(0.035,
+     &           0.035*(1.-log(paprs(i,k)/6.e6)/log(9.e6/6.e6))),0.)
+            end do
+           end do
+         endif
+ 
 C=========================================================================
 C=========================================================================
@@ -168,22 +195,21 @@
 
 c======================================================================
-c Emission d'un traceur pendant un certain temps
+c Emission continue d'un traceur
 c necessite raz_date=1 dans run.def
 c et reinit_trac=y
 c======================================================================
        deltatr(:,:,:) = 0.
-
-c source appliquee pendant Nemiss Ed
-       if (timesimu .lt. 86400*Nemiss) then
+       flux(:,:)=0.
 
 c emet les traceurs qui sont presents sur la grille
-        do ilat  = 1,nblat
-        do ilon  = 1,nblon
-         if (ig_volcan(ilat,ilon).ne.0) then
-	 
-          do iz = 1,nbz
-           it=min( (iz-1)*nblat*nblon+(ilat-1)*nblon+ilon , nqtot )	    
-           i=ig_volcan(ilat,ilon)
-	   
+      do ilat  = 1,nblat
+       do ilon  = 1,nblon
+        do iaire = 1,nbaire
+         do iflux = 1,nbflux
+     
+              it=min( (ilat-1)*nblon*nbflux*nbaire+(iaire-1)*nbflux
+     &         +(ilon-1)*nbaire*nbflux+iflux , nqtot )   
+ 
+
 c injection dans une seule cellule:
 c source en kg/kg/s
@@ -196,22 +222,42 @@
 
 c injection dans toute la colonne (a faire):
-            do k=1,Nheight(iz)
-             deltatr(i,k,it) = so2_quantity/(86400.*Nemiss) ! kg/s
-     $  *RG/( area_emiss(ilat,ilon)
-     $       *(paprs(i,1)-paprs(i,Nheight(iz)+1)) )    ! /kg (masse colonne)
-     
-             tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
-	    end do
-	    
+          do ipos=1,maxcell
+ 
+           if (ig_zone(ilat,ilon,iaire,iflux,ipos).ne.0) then
+
+             i=ig_zone(ilat,ilon,iaire,iflux,ipos)         
+             flux(i,it)=flux_surface_co2(iflux)  
+              
+            do k=1,Nheight
+             deltatr(i,k,it) = flux_surface_co2(iflux) ! kg/s/m2
+     $         *RG/(paprs(i,1)-paprs(i,Nheight+1))    ! /kg (masse colonne)
+     
+               tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
+            end do
+
+           end if 
           end do
-	  
-         endif  ! ig_volcan!=0
-        end do
-        end do
-
-       end if  ! duree emission
+ 
+        end do
+        end do
+        end do
+        end do
+
        
 c======================================================================
 c======================================================================
+
+   
+      
+
+     
+#ifdef CPP_XIOS      
+       do it=1,nqtot
+       CALL  send_xios_field("flux_"//tname(it),
+     &                     flux(:,it)) 
+      end do 
+#endif
+
+
 
       RETURN
Index: trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F	(revision 2463)
+++ trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F	(revision 2464)
@@ -285,10 +285,32 @@
 c         PHEAT(j) = PHEAT(j)*2.
 c        endif
-        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
+c BASE:
+c        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
+c          PHEAT(j) = PHEAT(j)*3
+c        endif
+c AFTER Tayloring TdeepD
+         if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
+           PHEAT(j) = PHEAT(j)*3.5
+         endif
+         if ((PPB(j).gt.10.).and.(PPB(j).le.50.)) then
+           PHEAT(j) = PHEAT(j)*1.5
+         endif
+c Options:
+c        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
+c        if ((PPB(j).gt.2.0).and.(PPB(j).le.10.)) then
 c        if ((PPB(j).gt.1.4).and.(PPB(j).le.100.)) then
-         PHEAT(j) = PHEAT(j)*3
-        endif
-c        if ((PPB(j).gt.10.).and.(PPB(j).le.120.)) then
-c         PHEAT(j) = PHEAT(j)*2
+c          PHEAT(j) = PHEAT(j)*3.5
+c          PHEAT(j) = PHEAT(j)*3
+c          PHEAT(j) = PHEAT(j)*2.5
+c        endif
+c        if ((PPB(j).gt.10.).and.(PPB(j).le.35.)) then
+c        if ((PPB(j).gt.10.).and.(PPB(j).le.50.)) then
+c          PHEAT(j) = PHEAT(j)*2
+c          PHEAT(j) = PHEAT(j)*1.5
+c          PHEAT(j) = PHEAT(j)*1.3
+c        endif
+c        if ((PPB(j).gt.35.).and.(PPB(j).le.120.)) then
+c          PHEAT(j) = PHEAT(j)*2
+c          PHEAT(j) = PHEAT(j)*1.5
 c        endif
 c----------------
