Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 1266)
@@ -67,10 +67,6 @@
 !=======================================================================
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "chimiedata.h"
-!#include "tracer.h"
 #include "callkeys.h"
-!#include "conc.h"
 
 !     input:
@@ -668,6 +664,6 @@
         
          if (thermochem) then
-            call chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens,  &
-                             zpress,zlocal,szacol,ptimestep,zday)
+            call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
+                             zdens,zpress,zlocal,szacol,ptimestep,zday)
          end if
 
Index: trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90	(revision 1266)
@@ -1,4 +1,4 @@
-      SUBROUTINE chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens,  &
-           zpress,zlocal,zenit,ptimestep,zday)
+      SUBROUTINE chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp, &
+           zdens,zpress,zlocal,zenit,ptimestep,zday)
 
       use tracer_mod, only: nqmx, igcm_co2, igcm_co, igcm_o, igcm_o1d,  &
@@ -30,20 +30,14 @@
 !    ------------------
 !
-#include "dimensions.h"
-#include "dimphys.h"
 #include "callkeys.h"
-!#include "comdiurn.h"
-#include "param.h"
-#include "param_v4.h"
-!#include"tracer.h"
 !-----------------------------------------------------------------------
 !    Input/Output
 !    ------------
-      integer :: lswitch,ig,chemthermod
-      real :: zycol(nlayermx,nqmx)
-      real :: ztemp(nlayermx)
-      real :: zdens(nlayermx)
-      real :: zpress(nlayermx)			! in mbar
-      real :: zlocal(nlayermx)
+      integer :: lswitch,ig,chemthermod,nlayer
+      real :: zycol(nlayer,nqmx)
+      real :: ztemp(nlayer)
+      real :: zdens(nlayer)
+      real :: zpress(nlayer)			! in mbar
+      real :: zlocal(nlayer)
       real :: zenit
       real :: ptimestep
@@ -426,7 +420,7 @@
 
       !Allocate density vector
-      allocate(rm(nlayermx,nesptherm))
-
-      do l=1,nlayermx
+      allocate(rm(nlayer,nesptherm))
+
+      do l=1,nlayer
         rm(l,i_co2)     = max(zycol(l,g_co2)*zdens(l),1.e-30)
         rm(l,i_co)      = max(zycol(l,g_co)*zdens(l),1.e-30) 
@@ -474,7 +468,7 @@
       call flujo(solarcondate)
       if(solvarmod.eq.0) then
-         call jthermcalc(ig,chemthermod,rm,nesptherm,ztemp,zlocal,zenit)
+         call jthermcalc(ig,nlayer,chemthermod,rm,nesptherm,ztemp,zlocal,zenit)
       else if(solvarmod.eq.1) then
-         call jthermcalc_e107(ig,chemthermod,rm,nesptherm,ztemp,zlocal,zenit,zday)
+         call jthermcalc_e107(ig,nlayer,chemthermod,rm,nesptherm,ztemp,zlocal,zenit,zday)
       endif
          
@@ -482,8 +476,8 @@
       !Chemistry
       call paramfoto_compact  &
-           (ig,chemthermod,lswitch,ztemp,ptimestep,zenit,zlocal,rm,nesptherm)
+           (ig,nlayer,chemthermod,lswitch,ztemp,ptimestep,zenit,zlocal,rm,nesptherm)
 
       !Concentrations back to GCM
-      do l=lswitch,nlayermx
+      do l=lswitch,nlayer
         zycol(l,g_co2)     = max(rm(l,i_co2)     / zdens(l) , 1.e-30)
         zycol(l,g_co)      = max(rm(l,i_co)      / zdens(l) , 1.e-30)
Index: trunk/LMDZ.MARS/libf/aeronomars/chemthermos_readini.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/chemthermos_readini.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/chemthermos_readini.F	(revision 1266)
@@ -1,15 +1,6 @@
       subroutine chemthermos_readini
 
+      use param_v4_h, only: rcoef
       implicit none
-
-
-c     common variables and constants
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-      include 'datafile.h'
-      include "callkeys.h"
-
 
 c     local variables
Index: trunk/LMDZ.MARS/libf/aeronomars/concentrations.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/concentrations.F	(revision 1266)
@@ -27,11 +27,6 @@
 !     declarations
  
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "comdiurn.h"
 #include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 
 !     input/output
Index: trunk/LMDZ.MARS/libf/aeronomars/conduction.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/conduction.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/conduction.F	(revision 1266)
@@ -16,10 +16,4 @@
 c   declarations:
 c-----------------------------------------------------------------------
-
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "surfdat.h"
-!#include "chimiedata.h"
-!#include "conc.h"
 
 c   arguments:
Index: trunk/LMDZ.MARS/libf/aeronomars/deposition.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/deposition.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/deposition.F	(revision 1266)
@@ -12,9 +12,4 @@
       implicit none
 c
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "chimiedata.h"
-!#include "conc.h"
-!#include "surfdat.h"
 c
 c     input
Index: trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90	(revision 1266)
@@ -31,13 +31,5 @@
 !    ------------------
 !
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "comdiurn.h"
-!#include "param.h"
-!#include "param_v4.h"
-!#include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 !-----------------------------------------------------------------------
 !    Input/Output
@@ -413,5 +405,5 @@
          enddo
         !Routine to calculate the UV heating
-         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
+         call hrtherm (ig,nlayer,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
 
 !        euveff=0.16    !UV heating efficiency. Following Fox et al. ASR 1996
Index: trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/hrtherm.F	(revision 1266)
@@ -1,5 +1,6 @@
 c**********************************************************************
 
-      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday,jtot)
+      subroutine hrtherm(ig,nlayer,
+     .      euvmod,rm,nespeuv,tx,iz,zenit,zday,jtot)
 
 
@@ -9,13 +10,9 @@
 c**********************************************************************
 
+      use param_v4_h, only: ninter,nabs,jfotsout,fluxtop,freccen
+
       implicit none
 
 c     common variables and constants
-
-
-      include 'dimensions.h'
-      include 'dimphys.h'
-      include 'param.h'
-      include 'param_v4.h'
       include "callkeys.h"
 
@@ -23,6 +20,6 @@
 c    local parameters and variables
 
-      real       xabsi(nabs,nlayermx) 			!densities
-      real       jergs(ninter,nabs,nlayermx)
+      real       xabsi(nabs,nlayer) 			!densities
+      real       jergs(ninter,nabs,nlayer)
       
       integer    i,j,k,indexint          !indexes
@@ -32,11 +29,11 @@
 c     input and output variables
 
-      integer    ig  ,euvmod 
+      integer    ig  ,euvmod,nlayer 
       integer    nespeuv
-      real       rm(nlayermx,nespeuv)              !density matrix (cm^-3)
-      real       jtot(nlayermx)                    !output: heating rate(erg/s)
-      real       tx(nlayermx)                      !temperature
+      real       rm(nlayer,nespeuv)              !density matrix (cm^-3)
+      real       jtot(nlayer)                    !output: heating rate(erg/s)
+      real       tx(nlayer)                      !temperature
       real       zenit
-      real       iz(nlayermx)
+      real       iz(nlayer)
       real       zday
 
@@ -67,5 +64,5 @@
       end if
       if(dn.eq.'n') then
-        do i=1,nlayermx                                    
+        do i=1,nlayer                                    
 	      jtot(i)=0.
         enddo       
@@ -78,5 +75,5 @@
       jtot(:)=0.
       !All number densities to a single array, xabsi(species,layer)
-      do i=1,nlayermx
+      do i=1,nlayer
          xabsi(1,i)  = rm(i,i_co2)
          xabsi(2,i)  = rm(i,i_o2)
@@ -102,11 +99,12 @@
       !Calculation of photoabsortion coefficient
       if(solvarmod.eq.0) then
-         call jthermcalc(ig,euvmod,rm,nespeuv,tx,iz,zenit)
+         call jthermcalc(ig,nlayer,euvmod,rm,nespeuv,tx,iz,zenit)
       else if (solvarmod.eq.1) then
-         call jthermcalc_e107(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday)
+         call jthermcalc_e107(ig,nlayer,euvmod,
+     .           rm,nespeuv,tx,iz,zenit,zday)
       endif
 
       !Total photoabsorption coefficient
-      do i=1,nlayermx
+      do i=1,nlayer
          jtot(i)=0.
         do j=1,nabs
Index: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 1266)
@@ -34,7 +34,5 @@
 
 #include "dimensions.h"
-!#include "dimphys.h"
 #include "paramet.h"
-!#include "tracer.h"
 #include "comvert.h"
 #include "callkeys.h"
Index: trunk/LMDZ.MARS/libf/aeronomars/iono.h
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/iono.h	(revision 1265)
+++ 	(revision )
@@ -1,47 +1,0 @@
-      common/eqF/ o1d_eq,ho2_eq,oh_eq,h_eq,n2d_eq,no2_eq,
-     $ o3_eq,no_eq,
-     $ cplus_eq,coplus_eq,oplus_eq,n2plus_eq,hplus_eq,co2plus_eq,
-     $ o2plus_eq,noplus_eq,nplus_eq,hco2plus_eq
-
-	character*1 o1d_eq(nlayermx),ho2_eq(nlayermx),oh_eq(nlayermx)
-      character*1 h_eq(nlayermx),n2d_eq(nlayermx),no2_eq(nlayermx)
-      character*1 o3_eq(nlayermx),no_eq(nlayermx),cplus_eq(nlayermx)
-      character*1 coplus_eq(nlayermx),oplus_eq(nlayermx)
-      character*1 n2plus_eq(nlayermx)
-      character*1 hplus_eq(nlayermx),co2plus_eq(nlayermx)
-      character*1 o2plus_eq(nlayermx)
-      character*1 noplus_eq(nlayermx), nplus_eq(nlayermx) 
-      character*1 hco2plus_eq(nlayermx)      
-
-      common/taus/ tauco2,tauo2,tauo3p,tauco,tauh,tauoh,tauho2,tauh2,
-     $ tauh2o,tauo1d,tauh2o2,tauo3,taun,tauno,taun2,taun2d,tauno2,
-     $ tauco2plus,tauoplus,tauo2plus,taucoplus,taucplus,taunplus,
-     $ taunoplus,taun2plus,tauhplus,tauhco2plus
-
-      real*8    tauco2(nreact,nlayermx)
-      real*8    tauo2(nreact,nlayermx)
-      real*8    tauo3p(nreact,nlayermx)
-      real*8    tauco(nreact,nlayermx)
-      real*8    tauh(nreact,nlayermx)
-      real*8    tauoh(nreact,nlayermx)
-      real*8    tauho2(nreact,nlayermx)
-      real*8    tauh2(nreact,nlayermx)
-      real*8    tauh2o(nreact,nlayermx)
-      real*8    tauo1d(nreact,nlayermx)
-      real*8    tauh2o2(nreact,nlayermx)
-      real*8    tauo3(nreact,nlayermx)
-      real*8    taun(nreact,nlayermx)
-      real*8    tauno(nreact,nlayermx)
-      real*8    taun2(nreact,nlayermx)
-      real*8    taun2d(nreact,nlayermx)
-      real*8    tauno2(nreact,nlayermx)
-      real*8    tauco2plus(nreact,nlayermx)
-      real*8    tauoplus(nreact,nlayermx)
-      real*8    tauo2plus(nreact,nlayermx)
-      real*8    taucoplus(nreact,nlayermx)
-      real*8    taucplus(nreact,nlayermx)
-      real*8    taunplus(nreact,nlayermx) 
-      real*8    taunoplus(nreact,nlayermx)
-      real*8    taun2plus(nreact,nlayermx)
-      real*8    tauhplus(nreact,nlayermx)
-      real*8    tauhco2plus(nreact,nlayermx)
Index: trunk/LMDZ.MARS/libf/aeronomars/iono_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/iono_h.F90	(revision 1266)
+++ trunk/LMDZ.MARS/libf/aeronomars/iono_h.F90	(revision 1266)
@@ -0,0 +1,108 @@
+MODULE iono_h
+
+      IMPLICIT NONE
+
+      character*1,save,allocatable :: o1d_eq(:)
+      character*1,save,allocatable :: ho2_eq(:)
+      character*1,save,allocatable :: oh_eq(:)
+      character*1,save,allocatable :: h_eq(:)
+      character*1,save,allocatable :: n2d_eq(:)
+      character*1,save,allocatable :: no2_eq(:)
+      character*1,save,allocatable :: o3_eq(:)
+      character*1,save,allocatable :: no_eq(:)
+      character*1,save,allocatable :: cplus_eq(:)
+      character*1,save,allocatable :: coplus_eq(:)
+      character*1,save,allocatable :: oplus_eq(:)
+      character*1,save,allocatable :: n2plus_eq(:)
+      character*1,save,allocatable :: hplus_eq(:)
+      character*1,save,allocatable :: co2plus_eq(:)
+      character*1,save,allocatable :: o2plus_eq(:)
+      character*1,save,allocatable :: noplus_eq(:)
+      character*1,save,allocatable :: nplus_eq(:)
+      character*1,save,allocatable :: hco2plus_eq(:)
+      real*8,save,allocatable ::   tauco2(:,:)
+      real*8,save,allocatable ::   tauo2(:,:)
+      real*8,save,allocatable ::   tauo3p(:,:)
+      real*8,save,allocatable ::   tauco(:,:)
+      real*8,save,allocatable ::   tauh(:,:)
+      real*8,save,allocatable ::   tauoh(:,:)
+      real*8,save,allocatable ::   tauho2(:,:)
+      real*8,save,allocatable ::   tauh2(:,:)
+      real*8,save,allocatable ::   tauh2o(:,:)
+      real*8,save,allocatable ::   tauo1d(:,:)
+      real*8,save,allocatable ::   tauh2o2(:,:)
+      real*8,save,allocatable ::   tauo3(:,:)
+      real*8,save,allocatable ::   taun(:,:)
+      real*8,save,allocatable ::   tauno(:,:)
+      real*8,save,allocatable ::   taun2(:,:)
+      real*8,save,allocatable ::   taun2d(:,:)
+      real*8,save,allocatable ::   tauno2(:,:)
+      real*8,save,allocatable ::   tauco2plus(:,:)
+      real*8,save,allocatable ::   tauoplus(:,:)
+      real*8,save,allocatable ::   tauo2plus(:,:)
+      real*8,save,allocatable ::   taucoplus(:,:)
+      real*8,save,allocatable ::   taucplus(:,:)
+      real*8,save,allocatable ::   taunplus(:,:) 
+      real*8,save,allocatable ::   taunoplus(:,:)
+      real*8,save,allocatable ::   taun2plus(:,:)
+      real*8,save,allocatable ::   tauhplus(:,:)
+      real*8,save,allocatable ::   tauhco2plus(:,:)
+
+      CONTAINS
+
+         SUBROUTINE allocate_param_iono(nlayer,nreact)
+
+           IMPLICIT NONE
+
+           INTEGER :: nreact
+           INTEGER :: nlayer
+
+           allocate(o1d_eq(nlayer))
+           allocate(ho2_eq(nlayer))
+           allocate(oh_eq(nlayer))
+           allocate(h_eq(nlayer))
+           allocate(n2d_eq(nlayer))
+           allocate(no2_eq(nlayer))
+           allocate(o3_eq(nlayer))
+           allocate(no_eq(nlayer))
+           allocate(cplus_eq(nlayer))
+           allocate(coplus_eq(nlayer))
+           allocate(oplus_eq(nlayer))
+           allocate(n2plus_eq(nlayer))
+           allocate(hplus_eq(nlayer))
+           allocate(co2plus_eq(nlayer))
+           allocate(o2plus_eq(nlayer))
+           allocate(noplus_eq(nlayer))
+           allocate(nplus_eq(nlayer))
+           allocate(hco2plus_eq(nlayer))
+           allocate(tauco2(nreact,nlayer))
+           allocate(tauo2(nreact,nlayer))
+           allocate(tauo3p(nreact,nlayer))
+           allocate(tauco(nreact,nlayer))
+           allocate(tauh(nreact,nlayer))
+           allocate(tauoh(nreact,nlayer))
+           allocate(tauho2(nreact,nlayer))
+           allocate(tauh2(nreact,nlayer))
+           allocate(tauh2o(nreact,nlayer))
+           allocate(tauo1d(nreact,nlayer))
+           allocate(tauh2o2(nreact,nlayer))
+           allocate(tauo3(nreact,nlayer))
+           allocate(taun(nreact,nlayer))
+           allocate(tauno(nreact,nlayer))
+           allocate(taun2(nreact,nlayer))
+           allocate(taun2d(nreact,nlayer))
+           allocate(tauno2(nreact,nlayer))
+           allocate(tauco2plus(nreact,nlayer))
+           allocate(tauoplus(nreact,nlayer))
+           allocate(tauo2plus(nreact,nlayer))
+           allocate(taucoplus(nreact,nlayer))
+           allocate(taucplus(nreact,nlayer))
+           allocate(taunplus(nreact,nlayer)) 
+           allocate(taunoplus(nreact,nlayer))
+           allocate(taun2plus(nreact,nlayer))
+           allocate(tauhplus(nreact,nlayer))
+           allocate(tauhco2plus(nreact,nlayer))
+         
+         END SUBROUTINE allocate_param_iono
+
+END MODULE iono_h
Index: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F	(revision 1266)
@@ -1,5 +1,6 @@
 c**********************************************************************
 
-      subroutine jthermcalc(ig,chemthermod,rm,nesptherm,tx,iz,zenit)
+      subroutine jthermcalc(ig,nlayer,chemthermod,
+     .      rm,nesptherm,tx,iz,zenit)
 
 
@@ -11,42 +12,41 @@
 c**********************************************************************
 
+      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
+
       implicit none
 
-c     common variables and constants
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-
 c     input and output variables
 
-      integer    ig
+      integer    ig,nlayer
       integer    chemthermod
       integer    nesptherm                      !Number of species considered
-      real       rm(nlayermx,nesptherm)         !Densities (cm-3)
-      real       tx(nlayermx)                   !temperature
+      real       rm(nlayer,nesptherm)         !Densities (cm-3)
+      real       tx(nlayer)                   !temperature
       real       zenit                          !SZA
-      real       iz(nlayermx)                   !Local altitude
+      real       iz(nlayer)                   !Local altitude
 
 
 c    local parameters and variables
 
-      real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
-      real       o2colx(nlayermx)               !column density of O2(cm^-2)
-      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2)
-      real       h2colx(nlayermx)               !H2 column density (cm-2)
-      real       h2ocolx(nlayermx)              !H2O column density (cm-2)
-      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2)
-      real       o3colx(nlayermx)               !O3 column density (cm-2)
-      real       n2colx(nlayermx)               !N2 column density (cm-2)
-      real       ncolx(nlayermx)                !N column density (cm-2)
-      real       nocolx(nlayermx)               !NO column density (cm-2)
-      real       cocolx(nlayermx)               !CO column density (cm-2)
-      real       hcolx(nlayermx)                !H column density (cm-2)
-      real       no2colx(nlayermx)              !NO2 column density (cm-2)
-      real       t2(nlayermx)
-      real       coltemp(nlayermx)
-      real       sigma(ninter,nlayermx)
-      real       alfa(ninter,nlayermx)
+      real       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
@@ -71,10 +71,10 @@
       real*8      auxjh(nz2)
       real*8      auxjno2(nz2)
-      real*8      wp(nlayermx),wm(nlayermx)
-      real*8      auxcolinp(nlayermx)
-      integer     auxind(nlayermx)
+      real*8      wp(nlayer),wm(nlayer)
+      real*8      auxcolinp(nlayer)
+      integer     auxind(nlayer)
       integer     auxi
       integer     ind
-      real*8      cortemp(nlayermx)
+      real*8      cortemp(nlayer)
 
       real*8      limdown                      !limits for interpolation
@@ -102,5 +102,5 @@
       !Auxiliar temperature to take into account the temperature dependence
       !of CO2 cross section
-      do i=1,nlayermx
+      do i=1,nlayer
          t2(i)=tx(i)
          if(t2(i).lt.195.0) t2(i)=195.0
@@ -109,5 +109,5 @@
 
       !Calculation of column amounts 
-      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
+      call column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,
      $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
      $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
@@ -115,6 +115,6 @@
       !Auxiliar column to include the temperature dependence 
       !of CO2 cross section
-      coltemp(nlayermx)=co2colx(nlayermx)*abs(t2(nlayermx)-t0(nlayermx))
-      do i=nlayermx-1,1,-1
+      coltemp(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 
@@ -123,5 +123,5 @@
       
       !Calculation of CO2 cross section at temperature t0(i)
-      do i=1,nlayermx
+      do i=1,nlayer
          do indexint=24,32
            sigma(indexint,i)=co2crsc195(indexint-23)
@@ -148,6 +148,6 @@
 c     Input atmospheric column
       indexint=1
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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) + 
@@ -159,5 +159,4 @@
       limdown=1.e-20
       limup=1.e26
-
 
 c     Interpolations
@@ -185,18 +184,18 @@
 
       call interfast 
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi=nlayermx-i+1
+         auxi=nlayer-i+1
          ! Ehouarn: test
           if ((ind+1.gt.nz2).or.(ind.le.0)) then
             write(*,*) "jthercalc error: ind=",ind,ig,zenit
-            write(*,*) " auxind(1:nlayermx)=",auxind
-            write(*,*) " auxcolinp(:nlayermx)=",auxcolinp
-            write(*,*) " co2colx(:nlayermx)=",co2colx
-            write(*,*) " o2colx(:nlayermx)=",o2colx
-            write(*,*) " o3pcolx(:nlayermx)=",o3pcolx
-            write(*,*) " h2colx(:nlayermx)=",h2colx
-            write(*,*) " ncolx(:nlayermx)=",ncolx
+            write(*,*) " auxind(1:nlayer)=",auxind
+            write(*,*) " auxcolinp(:nlayer)=",auxcolinp
+            write(*,*) " co2colx(:nlayer)=",co2colx
+            write(*,*) " o2colx(:nlayer)=",o2colx
+            write(*,*) " o3pcolx(:nlayer)=",o3pcolx
+            write(*,*) " h2colx(:nlayer)=",h2colx
+            write(*,*) " ncolx(:nlayer)=",ncolx
             write(*,*) " auxcoltab(1:nz2)=",auxcoltab
             write(*,*) " limdown=",limdown
@@ -220,7 +219,7 @@
       !N interpolated coefficient
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
-            jfotsout(indexint,9,nlayermx-i+1) =  wm(i)*auxjn(ind+1) +
+            jfotsout(indexint,9,nlayer-i+1) =  wm(i)*auxjn(ind+1) +
      $         wp(i)*auxjn(ind)
          enddo
@@ -241,6 +240,6 @@
 c     Input atmospheric column
       do indexint=2,15
-         do i=1,nlayermx
-            auxcolinp(nlayermx-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)+
@@ -288,9 +287,9 @@
          endif
 
-          call interfast(wm,wp,auxind,auxcolinp,nlayermx,
+          call interfast(wm,wp,auxind,auxcolinp,nlayer,
      $        auxcoltab,nz2,limdown,limup)
-          do i=1,nlayermx
+          do i=1,nlayer
              ind=auxind(i)
-             auxi = nlayermx-i+1
+             auxi = nlayer-i+1
              !O2 interpolated coefficient
              jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
@@ -317,7 +316,7 @@
           !Only if chemthermod.ge.2
           if(chemthermod.ge.2) then
-             do i=1,nlayermx
+             do i=1,nlayer
                 ind=auxind(i)
-                auxi = nlayermx-i+1
+                auxi = nlayer-i+1
                 !N interpolated coefficient
                 jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
@@ -346,6 +345,6 @@
 c     Input atmospheric column
       indexint=16
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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)+
@@ -393,8 +392,8 @@
 
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
          jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
@@ -418,7 +417,7 @@
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !N interpolated coefficient
             jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
@@ -444,6 +443,6 @@
 c     Input column
 
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -478,9 +477,9 @@
 
          call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,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,nlayermx
-               auxi = nlayermx-i+1
+            do i=1,nlayer
+               auxi = nlayer-i+1
                if(sigma(indexint,auxi)*
      $              alfa(indexint,auxi)*coltemp(auxi)
@@ -493,11 +492,11 @@
             enddo
          else
-            do i=1,nlayermx
+            do i=1,nlayer
                cortemp(i)=1.
             enddo
          end if
-         do i=1,nlayermx           
+         do i=1,nlayer           
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
@@ -519,7 +518,7 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,nlayermx
+            do i=1,nlayer
                ind=auxind(i)
-               auxi = nlayermx-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
@@ -544,6 +543,6 @@
 c     Input atmospheric column
 
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -579,8 +578,8 @@
          endif
          call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-         do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !Correction to include T variation of CO2 cross section
             if(sigma(indexint,auxi)*alfa(indexint,auxi)*
@@ -611,7 +610,7 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,nlayermx
+            do i=1,nlayer
                ind=auxind(i)
-               auxi = nlayermx-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
@@ -638,6 +637,6 @@
 c     Input atmospheric column
 
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -672,8 +671,8 @@
 
          call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-         do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !Correction to include T variation of CO2 cross section
             if(sigma(indexint,auxi)*alfa(indexint,auxi)*
@@ -701,7 +700,7 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,nlayermx 
+            do i=1,nlayer 
                ind=auxind(i)
-               auxi = nlayermx-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
@@ -728,6 +727,6 @@
 
       indexint=32
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -757,10 +756,10 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !Correction to include T variation of CO2 cross section
-         if(sigma(indexint,nlayermx-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)*
@@ -783,6 +782,6 @@
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
-            auxi = nlayermx-i+1
+         do i=1,nlayer
+            auxi = nlayer-i+1
             ind=auxind(i)
             !NO interpolated coefficient
@@ -808,6 +807,6 @@
 
       indexint=33
-      do i=1,nlayermx
-         auxcolinp(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
       end do
 
@@ -831,8 +830,8 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
          jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
@@ -844,8 +843,8 @@
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
+            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
      $           wp(i)*auxjno2(ind)
          enddo
@@ -865,6 +864,6 @@
 
       indexint=34
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -891,8 +890,8 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
          jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
@@ -907,8 +906,8 @@
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
+            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
      $           wp(i)*auxjno2(ind)
          enddo
@@ -928,6 +927,6 @@
 
       indexint=35
-      do i=1,nlayermx
-         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
       end do
 
@@ -951,8 +950,8 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !H2O2 interpolated coefficient
          jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
@@ -963,8 +962,8 @@
       enddo
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
+            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
      $           wp(i)*auxjno2(ind)
          enddo
@@ -983,6 +982,6 @@
 
       indexint=36
-      do i=1,nlayermx
-         auxcolinp(nlayermx-i+1) = o3colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
       end do
 
@@ -1004,17 +1003,17 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
          !O3 interpolated coefficient
-         jfotsout(indexint,7,nlayermx-i+1) = wm(i)*auxjo3(ind+1) +
+         jfotsout(indexint,7,nlayer-i+1) = wm(i)*auxjo3(ind+1) +
      $        wp(i)*auxjo3(ind)         
       enddo
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
+            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
      $           wp(i)*auxjno2(ind)
          enddo
@@ -1035,5 +1034,5 @@
 c**********************************************************************
 
-      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
+      subroutine column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,
      $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
      $     n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
@@ -1047,13 +1046,10 @@
      &                      igcm_o3, igcm_n2, igcm_n, igcm_no, igcm_no2,
      &                      mmol
+      use param_v4_h, only: radio,gg,masa,kboltzman,n_avog
+
       implicit none
 
 
 c     common variables and constants
-      include "dimensions.h"
-      include "dimphys.h"
-!      include "tracer.h"
-      include 'param.h'
-      include 'param_v4.h'
       include 'callkeys.h'
 
@@ -1066,46 +1062,46 @@
 c     input and output variables
 
-      integer    ig
+      integer    ig,nlayer
       integer    chemthermod
       integer    nesptherm                      !# of species undergoing chemistry, input
-      real       rm(nlayermx,nesptherm)         !densities (cm-3), input
-      real       tx(nlayermx)                   !temperature profile, input
-      real       iz(nlayermx+1)                 !height profile, input
+      real       rm(nlayer,nesptherm)         !densities (cm-3), input
+      real       tx(nlayer)                   !temperature profile, input
+      real       iz(nlayer+1)                 !height profile, input
       real       zenit                          !SZA, input
-      real       co2colx(nlayermx)              !column density of CO2 (cm^-2), output
-      real       o2colx(nlayermx)               !column density of O2(cm^-2), output
-      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2), output
-      real       h2colx(nlayermx)               !H2 column density (cm-2), output
-      real       h2ocolx(nlayermx)              !H2O column density (cm-2), output
-      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2), output
-      real       o3colx(nlayermx)               !O3 column density (cm-2), output
-      real       n2colx(nlayermx)               !N2 column density (cm-2), output
-      real       ncolx(nlayermx)                !N column density (cm-2), output
-      real       nocolx(nlayermx)               !NO column density (cm-2), output
-      real       cocolx(nlayermx)               !CO column density (cm-2), output
-      real       hcolx(nlayermx)                !H column density (cm-2), output
-      real       no2colx(nlayermx)              !NO2 column density (cm-2), output
-      
+      real       co2colx(nlayer)              !column density of CO2 (cm^-2), output
+      real       o2colx(nlayer)               !column density of O2(cm^-2), output
+      real       o3pcolx(nlayer)              !column density of O(3P)(cm^-2), output
+      real       h2colx(nlayer)               !H2 column density (cm-2), output
+      real       h2ocolx(nlayer)              !H2O column density (cm-2), output
+      real       h2o2colx(nlayer)             !column density of H2O2(cm^-2), output
+      real       o3colx(nlayer)               !O3 column density (cm-2), output
+      real       n2colx(nlayer)               !N2 column density (cm-2), output
+      real       ncolx(nlayer)                !N column density (cm-2), output
+      real       nocolx(nlayer)               !NO column density (cm-2), output
+      real       cocolx(nlayer)               !CO column density (cm-2), output
+      real       hcolx(nlayer)                !H column density (cm-2), output
+      real       no2colx(nlayer)              !NO2 column density (cm-2), output
+
 
 c     local variables
 
       real       xx
-      real       grav(nlayermx)
+      real       grav(nlayer)
       real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
       real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
 
-      real       co2x(nlayermx)
-      real       o2x(nlayermx)
-      real       o3px(nlayermx)
-      real       cox(nlayermx)
-      real       hx(nlayermx)
-      real       h2x(nlayermx)
-      real       h2ox(nlayermx)
-      real       h2o2x(nlayermx)
-      real       o3x(nlayermx)
-      real       n2x(nlayermx)
-      real       nx(nlayermx)
-      real       nox(nlayermx)
-      real       no2x(nlayermx)
+      real       co2x(nlayer)
+      real       o2x(nlayer)
+      real       o3px(nlayer)
+      real       cox(nlayer)
+      real       hx(nlayer)
+      real       h2x(nlayer)
+      real       h2ox(nlayer)
+      real       h2o2x(nlayer)
+      real       o3x(nlayer)
+      real       n2x(nlayer)
+      real       nx(nlayer)
+      real       nox(nlayer)
+      real       no2x(nlayer)
 
       integer    i,j,k,icol,indexint          !indexes
@@ -1117,7 +1113,7 @@
 
       integer    jj
-      real*8      esp(nlayermx*2)
-      real*8      ilayesp(nlayermx*2)
-      real*8      szalayesp(nlayermx*2)
+      real*8      esp(nlayer*2)
+      real*8      ilayesp(nlayer*2)
+      real*8      szalayesp(nlayer*2)
       integer     nlayesp
       real*8      zmini
@@ -1149,6 +1145,6 @@
 c*************************PROGRAM STARTS*******************************
 
-      nz3 = nlayermx*2
-      do i=1,nlayermx
+      nz3 = nlayer*2
+      do i=1,nlayer
          xx = ( radio + iz(i) ) * 1.e5
          grav(i) = gg * masa /(xx**2)
@@ -1156,5 +1152,5 @@
 
       !Scale heights
-      xx = kboltzman * tx(nlayermx) * n_avog / grav(nlayermx)
+      xx = kboltzman * tx(nlayer) * n_avog / grav(nlayer)
       Ho3p  = xx / mmol(igcm_o) 
       Hco2  = xx / mmol(igcm_co2)
@@ -1176,5 +1172,5 @@
       endif
       ! first loop in altitude : initialisation
-      do i=nlayermx,1,-1
+      do i=nlayer,1,-1
          !Column initialisation
          co2colx(i)  = 0.
@@ -1212,8 +1208,8 @@
       enddo
       ! second loop in altitude : column calculations
-      do i=nlayermx,1,-1
+      do i=nlayer,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)
+         call espesor_optico_A(ig,i,nlayer,zenit,iz(i),nz3,iz,esp,
+     $         ilayesp,szalayesp,nlayesp, zmini)
          if(ilayesp(nlayesp).eq.-1) then
             co2colx(i)=1.e25
@@ -1231,5 +1227,5 @@
             no2colx(i)=1.e25
          else
-            rcmnz = ( radio + iz(nlayermx) ) * 1.e5
+            rcmnz = ( radio + iz(nlayer) ) * 1.e5
             rcmmini = ( radio + zmini ) * 1.e5
             !Column calculation taking into account the geometrical depth 
@@ -1238,35 +1234,35 @@
                jj=ilayesp(j)
                !Top layer
-               if(jj.eq.nlayermx) then
+               if(jj.eq.nlayer) then
                   if(zenit.le.60.) then 
-                     o3pcolx(i)=o3pcolx(i)+o3px(nlayermx)*Ho3p*esp(j)
+                     o3pcolx(i)=o3pcolx(i)+o3px(nlayer)*Ho3p*esp(j)
      $                    *1.e-5
-                     co2colx(i)=co2colx(i)+co2x(nlayermx)*Hco2*esp(j)
+                     co2colx(i)=co2colx(i)+co2x(nlayer)*Hco2*esp(j)
      $                    *1.e-5
                      h2o2colx(i)=h2o2colx(i)+
-     $                    h2o2x(nlayermx)*Hh2o2*esp(j)*1.e-5
-                     o2colx(i)=o2colx(i)+o2x(nlayermx)*Ho2*esp(j)
+     $                    h2o2x(nlayer)*Hh2o2*esp(j)*1.e-5
+                     o2colx(i)=o2colx(i)+o2x(nlayer)*Ho2*esp(j)
      $                    *1.e-5
-                     h2colx(i)=h2colx(i)+h2x(nlayermx)*Hh2*esp(j)
+                     h2colx(i)=h2colx(i)+h2x(nlayer)*Hh2*esp(j)
      $                    *1.e-5
-                     h2ocolx(i)=h2ocolx(i)+h2ox(nlayermx)*Hh2o*esp(j)
+                     h2ocolx(i)=h2ocolx(i)+h2ox(nlayer)*Hh2o*esp(j)
      $                    *1.e-5                     
-                     cocolx(i)=cocolx(i)+cox(nlayermx)*Hco*esp(j)
+                     cocolx(i)=cocolx(i)+cox(nlayer)*Hco*esp(j)
      $                    *1.e-5
-                     hcolx(i)=hcolx(i)+hx(nlayermx)*Hh*esp(j)
+                     hcolx(i)=hcolx(i)+hx(nlayer)*Hh*esp(j)
      $                    *1.e-5
                      !Only if O3 chemistry required
                      if(chemthermod.ge.1) o3colx(i)=
-     $                    o3colx(i)+o3x(nlayermx)*Ho3*esp(j)
+     $                    o3colx(i)+o3x(nlayer)*Ho3*esp(j)
      $                    *1.e-5
                      !Only if N or ion chemistry requested
                      if(chemthermod.ge.2) then
-                        n2colx(i)=n2colx(i)+n2x(nlayermx)*Hn2*esp(j)
+                        n2colx(i)=n2colx(i)+n2x(nlayer)*Hn2*esp(j)
      $                    *1.e-5
-                        ncolx(i)=ncolx(i)+nx(nlayermx)*Hn*esp(j)
+                        ncolx(i)=ncolx(i)+nx(nlayer)*Hn*esp(j)
      $                       *1.e-5
-                        nocolx(i)=nocolx(i)+nox(nlayermx)*Hno*esp(j)
+                        nocolx(i)=nocolx(i)+nox(nlayer)*Hno*esp(j)
      $                       *1.e-5
-                        no2colx(i)=no2colx(i)+no2x(nlayermx)*Hno2*esp(j)
+                        no2colx(i)=no2colx(i)+no2x(nlayer)*Hno2*esp(j)
      $                       *1.e-5
                      endif
@@ -1290,21 +1286,21 @@
                         espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
                      endif
-                     co2colx(i) = co2colx(i) + espco2*co2x(nlayermx)
-                     o2colx(i)  = o2colx(i)  + espo2*o2x(nlayermx)
-                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayermx)
-                     h2colx(i)  = h2colx(i)  + esph2*h2x(nlayermx)
-                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayermx)
-                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayermx)
-                     cocolx(i)  = cocolx(i)  + espco*cox(nlayermx)
-                     hcolx(i)   = hcolx(i)   + esph*hx(nlayermx)
+                     co2colx(i) = co2colx(i) + espco2*co2x(nlayer)
+                     o2colx(i)  = o2colx(i)  + espo2*o2x(nlayer)
+                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayer)
+                     h2colx(i)  = h2colx(i)  + esph2*h2x(nlayer)
+                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayer)
+                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayer)
+                     cocolx(i)  = cocolx(i)  + espco*cox(nlayer)
+                     hcolx(i)   = hcolx(i)   + esph*hx(nlayer)
                      !Only if O3 chemistry required
                      if(chemthermod.ge.1)                      
-     $                  o3colx(i) = o3colx(i)  + espo3*o3x(nlayermx)
+     $                  o3colx(i) = o3colx(i)  + espo3*o3x(nlayer)
                      !Only if N or ion chemistry requested
                      if(chemthermod.ge.2) then
-                        n2colx(i)  = n2colx(i)  + espn2*n2x(nlayermx)
-                        ncolx(i)   = ncolx(i)   + espn*nx(nlayermx)
-                        nocolx(i)  = nocolx(i)  + espno*nox(nlayermx)
-                        no2colx(i) = no2colx(i) + espno2*no2x(nlayermx)
+                        n2colx(i)  = n2colx(i)  + espn2*n2x(nlayer)
+                        ncolx(i)   = ncolx(i)   + espn*nx(nlayer)
+                        nocolx(i)  = nocolx(i)  + espno*nox(nlayer)
+                        no2colx(i) = no2colx(i) + espno2*no2x(nlayer)
                      endif
                   endif !Of if zenit.lt.60
@@ -1342,9 +1338,8 @@
      $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
                   endif
-               endif  !Of if jj.eq.nlayermx
+               endif  !Of if jj.eq.nlayer
             end do    !Of do j=1,nlayesp
          endif        !Of ilayesp(nlayesp).eq.-1
-
-      enddo           !Of do i=nlayermx,1,-1
+      enddo           !Of do i=nlayer,1,-1
 
 
@@ -1394,5 +1389,5 @@
 c**********************************************************************
 
-      subroutine espesor_optico_A (ig,capa, szadeg,z,
+      subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z,
      @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
 
@@ -1402,13 +1397,6 @@
 *************************************************************************
 
+      use param_v4_h, only: radio
       implicit none
-
-
-c     common variables and constants
-
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
 
 c     arguments
@@ -1416,7 +1404,7 @@
       real        szadeg                ! I. SZA [rad]
       real        z                     ! I. altitude of interest [km]
-      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
-                                        !  (=2*nlayermx= max# of layers in ray path)
-      real     iz(nlayermx+1)              ! I. Altitude of each layer
+      integer     nz3,ig,nlayer         ! I. dimension of esp, ylayesp, etc...
+                                        !  (=2*nlayer= max# of layers in ray path)
+      real     iz(nlayer+1)              ! I. Altitude of each layer
       real*8        esp(nz3)            ! O. layer widths after geometrically 
                                         !    amplified; in [cm] except at TOA
@@ -1435,5 +1423,5 @@
         real*8      szarad                ! SZA [deg]
         real*8      zz
-        real*8      diz(nlayermx+1)
+        real*8      diz(nlayer+1)
         real*8      rkmnz                 ! distance TOA to center of Planet [km]
         real*8      rkmmini               ! distance zmini to center of P [km] 
@@ -1448,5 +1436,5 @@
         szarad = dble(szadeg)*3.141592d0/180.d0
         zz=dble(z)
-        do i=1,nlayermx
+        do i=1,nlayer
            diz(i)=dble(iz(i))
         enddo
@@ -1466,8 +1454,8 @@
 
            zmini = zz
-           if(abs(zz-diz(nlayermx)).lt.1.d-3) goto 1357
+           if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357
            ! 1st Zone: Upper part of ray
            !
-           do j=grid_R8(zz,diz,nlayermx),nlayermx-1
+           do j=grid_R8(zz,diz,nlayer),nlayer-1
              nlayesp = nlayesp + 1 
              ilayesp(nlayesp) = j
@@ -1481,5 +1469,5 @@
  1357      continue
            nlayesp = nlayesp + 1 
-           ilayesp(nlayesp) = nlayermx
+           ilayesp(nlayesp) = nlayer
            esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
            szalayesp(nlayesp) = szadeg
@@ -1498,9 +1486,9 @@
            rkmmini = radio + zmini
 
-           if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 1470
+           if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470
 
            ! 1st Zone: Upper part of ray
            !
-           do j=grid_R8(zz,diz,nlayermx),nlayermx-1
+           do j=grid_R8(zz,diz,nlayer),nlayer-1
               nlayesp = nlayesp + 1 
               ilayesp(nlayesp) = j
@@ -1517,6 +1505,6 @@
            !
            nlayesp = nlayesp + 1 
-           ilayesp(nlayesp) = nlayermx
-           rkmnz = radio+diz(nlayermx)
+           ilayesp(nlayesp) = nlayer
+           rkmnz = radio+diz(nlayer)
            esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
            esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
@@ -1549,18 +1537,18 @@
 
            else
-              jmin=grid_R8(zmini,diz,nlayermx)+1
+              jmin=grid_R8(zmini,diz,nlayer)+1
               !Correction for possible rounding errors when SZA very close 
               !to 90 degrees
-              if(jmin.gt.grid_R8(zz,diz,nlayermx)) then
+              if(jmin.gt.grid_R8(zz,diz,nlayer)) then
                  write(*,*)'jthermcalc warning: possible rounding error'
                  write(*,*)'point,sza,layer:',ig,szadeg,capa
-                 jmin=grid_R8(zz,diz,nlayermx)
+                 jmin=grid_R8(zz,diz,nlayer)
               endif
 
-              if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 9876
+              if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876
 
               ! 1st Zone: Upper part of ray
               !
-              do j=grid_R8(zz,diz,nlayermx),nlayermx-1
+              do j=grid_R8(zz,diz,nlayer),nlayer-1
                 nlayesp = nlayesp + 1 
                 ilayesp(nlayesp) = j
@@ -1578,6 +1566,6 @@
               !
               nlayesp = nlayesp + 1 
-              ilayesp(nlayesp) = nlayermx
-              rkmnz = radio+diz(nlayermx)
+              ilayesp(nlayesp) = nlayer
+              rkmnz = radio+diz(nlayer)
               esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
               esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
@@ -1601,5 +1589,5 @@
               ! 4th zone: Lower part of ray, increasing from zmin to z
               !    ( layers with SZA < 90 deg )
-              do j=jmin,grid_R8(zz,diz,nlayermx)-1
+              do j=jmin,grid_R8(zz,diz,nlayer)-1
                 nlayesp = nlayesp + 1 
                 ilayesp(nlayesp) = j
@@ -1615,5 +1603,5 @@
               ! 5th zone: Lower part of ray, decreasing from z to zmin
               !    ( layers with SZA > 90 deg )
-              do j=grid_R8(zz,diz,nlayermx)-1, jmin, -1
+              do j=grid_R8(zz,diz,nlayer)-1, jmin, -1
                 nlayesp = nlayesp + 1 
                 ilayesp(nlayesp) = j
@@ -1721,13 +1709,10 @@
 
       use comsaison_h, only: dist_sol
+      use param_v4_h, only: ninter, 
+     .                      fluxtop, ct1, ct2, p1, p2
       implicit none
 
 
 !     common variables and constants
-      include "dimensions.h"
-      include "dimphys.h"
-!      include "comsaison.h"
-      include 'param.h'
-      include 'param_v4.h'
       include "callkeys.h"
 
Index: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F	(revision 1266)
@@ -2,5 +2,5 @@
 
       subroutine jthermcalc_e107
-     $     (ig,chemthermod,rm,nesptherm,tx,iz,zenit,zday)
+     $     (ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,zday)
 
 
@@ -12,21 +12,22 @@
 c**********************************************************************
 
+      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 "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-
 c     input and output variables
 
-      integer    ig
+      integer    ig,nlayer
       integer    chemthermod
       integer    nesptherm                      !Number of species considered
-      real       rm(nlayermx,nesptherm)         !Densities (cm-3)
-      real       tx(nlayermx)                   !temperature
+      real       rm(nlayer,nesptherm)         !Densities (cm-3)
+      real       tx(nlayer)                   !temperature
       real       zenit                          !SZA
-      real       iz(nlayermx)                   !Local altitude
+      real       iz(nlayer)                   !Local altitude
       real       zday                           !Martian day after Ls=0
 
@@ -34,21 +35,21 @@
 c    local parameters and variables
 
-      real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
-      real       o2colx(nlayermx)               !column density of O2(cm^-2)
-      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2)
-      real       h2colx(nlayermx)               !H2 column density (cm-2)
-      real       h2ocolx(nlayermx)              !H2O column density (cm-2)
-      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2)
-      real       o3colx(nlayermx)               !O3 column density (cm-2)
-      real       n2colx(nlayermx)               !N2 column density (cm-2)
-      real       ncolx(nlayermx)                !N column density (cm-2)
-      real       nocolx(nlayermx)               !NO column density (cm-2)
-      real       cocolx(nlayermx)               !CO column density (cm-2)
-      real       hcolx(nlayermx)                !H column density (cm-2)
-      real       no2colx(nlayermx)              !NO2 column density (cm-2)
-      real       t2(nlayermx)
-      real       coltemp(nlayermx)
-      real       sigma(ninter,nlayermx)
-      real       alfa(ninter,nlayermx)
+      real       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)
       real       realday
       
@@ -75,10 +76,10 @@
       real*8      auxjh(nz2)
       real*8      auxjno2(nz2)
-      real*8      wp(nlayermx),wm(nlayermx)
-      real*8      auxcolinp(nlayermx)
-      integer     auxind(nlayermx)
+      real*8      wp(nlayer),wm(nlayer)
+      real*8      auxcolinp(nlayer)
+      integer     auxind(nlayer)
       integer     auxi
       integer     ind
-      real*8      cortemp(nlayermx)
+      real*8      cortemp(nlayer)
 
       real*8      limdown                      !limits for interpolation
@@ -106,5 +107,5 @@
       !Auxiliar temperature to take into account the temperature dependence
       !of CO2 cross section
-      do i=1,nlayermx
+      do i=1,nlayer
          t2(i)=tx(i)
          if(t2(i).lt.195.0) t2(i)=195.0
@@ -113,5 +114,5 @@
 
       !Calculation of column amounts 
-      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
+      call column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,
      $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
      $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
@@ -119,6 +120,6 @@
       !Auxiliar column to include the temperature dependence 
       !of CO2 cross section
-      coltemp(nlayermx)=co2colx(nlayermx)*abs(t2(nlayermx)-t0(nlayermx))
-      do i=nlayermx-1,1,-1
+      coltemp(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 
@@ -127,5 +128,5 @@
       
       !Calculation of CO2 cross section at temperature t0(i)
-      do i=1,nlayermx
+      do i=1,nlayer
          do indexint=24,32
            sigma(indexint,i)=co2crsc195(indexint-23)
@@ -157,5 +158,5 @@
       do j=1,nabs
          do indexint=1,ninter
-            jfotsout(indexint,j,nlayermx)=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
@@ -179,6 +180,6 @@
 c     Input atmospheric column
       indexint=1
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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) + 
@@ -214,16 +215,16 @@
 
       call interfast 
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi=nlayermx-i+1
+         auxi=nlayer-i+1
          !CO2 interpolated coefficient
-         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx) *
+         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
      $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
          !H2 interpolated coefficient
@@ -234,8 +235,8 @@
       !N interpolated coefficient
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
-            jfotsout(indexint,9,nlayermx-i+1) =  
-     $           jfotsout(indexint,9,nlayermx) * 
+            jfotsout(indexint,9,nlayer-i+1) =  
+     $           jfotsout(indexint,9,nlayer) * 
      $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
          enddo
@@ -256,6 +257,6 @@
 c     Input atmospheric column
       do indexint=2,15
-         do i=1,nlayermx
-            auxcolinp(nlayermx-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)+
@@ -303,54 +304,54 @@
          endif
 
-          call interfast(wm,wp,auxind,auxcolinp,nlayermx,
+          call interfast(wm,wp,auxind,auxcolinp,nlayer,
      $        auxcoltab,nz2,limdown,limup)
-          do i=1,nlayermx
+          do i=1,nlayer
              ind=auxind(i)
-             auxi = nlayermx-i+1
+             auxi = nlayer-i+1
              !O2 interpolated coefficient
              jfotsout(indexint,2,auxi) = 
-     $            jfotsout(indexint,2,nlayermx) *
+     $            jfotsout(indexint,2,nlayer) *
      $            (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
              !O3p interpolated coefficient
              jfotsout(indexint,3,auxi) = 
-     $            jfotsout(indexint,3,nlayermx) *
+     $            jfotsout(indexint,3,nlayer) *
      $            (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
              !CO2 interpolated coefficient
              jfotsout(indexint,1,auxi) = 
-     $            jfotsout(indexint,1,nlayermx) *
+     $            jfotsout(indexint,1,nlayer) *
      $            (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
              !H2 interpolated coefficient
              jfotsout(indexint,5,auxi) = 
-     $            jfotsout(indexint,5,nlayermx) *
+     $            jfotsout(indexint,5,nlayer) *
      $            (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
              !N2 interpolated coefficient
              jfotsout(indexint,8,auxi) = 
-     $            jfotsout(indexint,8,nlayermx) *
+     $            jfotsout(indexint,8,nlayer) *
      $            (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
              !CO interpolated coefficient
              jfotsout(indexint,11,auxi) = 
-     $            jfotsout(indexint,11,nlayermx) *
+     $            jfotsout(indexint,11,nlayer) *
      $            (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
              !H interpolated coefficient
              jfotsout(indexint,12,auxi) = 
-     $            jfotsout(indexint,12,nlayermx) *
+     $            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,nlayermx
+             do i=1,nlayer
                 ind=auxind(i)
-                auxi = nlayermx-i+1
+                auxi = nlayer-i+1
                 !N interpolated coefficient
                 jfotsout(indexint,9,auxi) = 
-     $               jfotsout(indexint,9,nlayermx) *
+     $               jfotsout(indexint,9,nlayer) *
      $               (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
                 !NO interpolated coefficient
                 jfotsout(indexint,10,auxi)=
-     $               jfotsout(indexint,10,nlayermx) *
+     $               jfotsout(indexint,10,nlayer) *
      $               (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
                 !NO2 interpolated coefficient
                 jfotsout(indexint,13,auxi)=
-     $               jfotsout(indexint,13,nlayermx) *
+     $               jfotsout(indexint,13,nlayer) *
      $               (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
              enddo
@@ -371,6 +372,6 @@
 c     Input atmospheric column
       indexint=16
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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)+
@@ -418,45 +419,45 @@
 
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx) *
+     $        jfotsout(indexint,11,nlayer) *
      $        (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
          !H interpolated coefficient
          jfotsout(indexint,12,auxi) = 
-     $        jfotsout(indexint,12,nlayermx) *
+     $        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,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !N interpolated coefficient
             jfotsout(indexint,9,auxi) = 
-     $           jfotsout(indexint,9,nlayermx) *
+     $           jfotsout(indexint,9,nlayer) *
      $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
             !NO interpolated coefficient
             jfotsout(indexint,10,auxi) = 
-     $           jfotsout(indexint,10,nlayermx) * 
+     $           jfotsout(indexint,10,nlayer) * 
      $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
             !NO2 interpolated coefficient
             jfotsout(indexint,13,auxi) = 
-     $           jfotsout(indexint,13,nlayermx) *
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -474,6 +475,6 @@
 c     Input column
 
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -508,9 +509,9 @@
 
          call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,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,nlayermx
-               auxi = nlayermx-i+1
+            do i=1,nlayer
+               auxi = nlayer-i+1
                if(sigma(indexint,auxi)*
      $              alfa(indexint,auxi)*coltemp(auxi)
@@ -523,19 +524,19 @@
             enddo
          else
-            do i=1,nlayermx
+            do i=1,nlayer
                cortemp(i)=1.
             enddo
          end if
-         do i=1,nlayermx           
+         do i=1,nlayer           
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = 
-     $           jfotsout(indexint,2,nlayermx) *
+     $           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,nlayermx) *
+     $           jfotsout(indexint,1,nlayer) *
      $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) 
      $           * cortemp(i)
@@ -546,10 +547,10 @@
             !N2 interpolated coefficient
             jfotsout(indexint,8,auxi) = 
-     $           jfotsout(indexint,8,nlayermx) *
+     $           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,nlayermx) *
+     $           jfotsout(indexint,11,nlayer) *
      $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) * 
      $           cortemp(i)            
@@ -557,15 +558,15 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,nlayermx
+            do i=1,nlayer
                ind=auxind(i)
-               auxi = nlayermx-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=
-     $              jfotsout(indexint,10,nlayermx) *
+     $              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,nlayermx) *
+     $              jfotsout(indexint,13,nlayer) *
      $              (wm(i)*auxjno2(ind+1)+ wp(i)*auxjno2(ind)) * 
      $              cortemp(i)
@@ -586,6 +587,6 @@
 c     Input atmospheric column
 
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -621,8 +622,8 @@
          endif
          call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-         do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !Correction to include T variation of CO2 cross section
             if(sigma(indexint,auxi)*alfa(indexint,auxi)*
@@ -635,5 +636,5 @@
             !CO2 interpolated coefficient
             jfotsout(indexint,1,auxi) = 
-     $           jfotsout(indexint,1,nlayermx) *
+     $           jfotsout(indexint,1,nlayer) *
      $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) * 
      $           cortemp(i) *
@@ -642,20 +643,20 @@
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = 
-     $           jfotsout(indexint,2,nlayermx) *
+     $           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,nlayermx) *
+     $           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,nlayermx) *
+     $           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,nlayermx) *
+     $           jfotsout(indexint,11,nlayer) *
      $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) * 
      $           cortemp(i)
@@ -663,15 +664,15 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,nlayermx
+            do i=1,nlayer
                ind=auxind(i)
-               auxi = nlayermx-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=
-     $              jfotsout(indexint,10,nlayermx) *
+     $              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,nlayermx) *
+     $              jfotsout(indexint,13,nlayer) *
      $              (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) * 
      $              cortemp(i)
@@ -694,6 +695,6 @@
 c     Input atmospheric column
 
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -728,8 +729,8 @@
 
          call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-         do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+         do i=1,nlayer
             ind=auxind(i)
-            auxi = nlayermx-i+1
+            auxi = nlayer-i+1
             !Correction to include T variation of CO2 cross section
             if(sigma(indexint,auxi)*alfa(indexint,auxi)*
@@ -742,5 +743,5 @@
             !CO2 interpolated coefficient
             jfotsout(indexint,1,auxi) = 
-     $           jfotsout(indexint,1,nlayermx) *
+     $           jfotsout(indexint,1,nlayer) *
      $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) * 
      $           cortemp(i) *
@@ -749,15 +750,15 @@
             !O2 interpolated coefficient
             jfotsout(indexint,2,auxi) = 
-     $           jfotsout(indexint,2,nlayermx) *
+     $           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,nlayermx) *
+     $           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,nlayermx) *
+     $           jfotsout(indexint,6,nlayer) *
      $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
      $           cortemp(i)            
@@ -765,10 +766,10 @@
          !Only if chemthermod.ge.2
          if(chemthermod.ge.2) then
-            do i=1,nlayermx 
+            do i=1,nlayer 
                ind=auxind(i)
-               auxi = nlayermx-i+1
+               auxi = nlayer-i+1
                !NO interpolated coefficient
                jfotsout(indexint,10,auxi)=
-     $              jfotsout(indexint,10,nlayermx) *
+     $              jfotsout(indexint,10,nlayer) *
      $              (wm(i)*auxjno(ind+1) +wp(i)*auxjno(ind)) * 
      $              cortemp(i)
@@ -796,6 +797,6 @@
 
       indexint=32
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -825,10 +826,10 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !Correction to include T variation of CO2 cross section
-         if(sigma(indexint,nlayermx-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)*
@@ -839,5 +840,5 @@
          !CO2 interpolated coefficient
          jfotsout(indexint,1,auxi) = 
-     $        jfotsout(indexint,1,nlayermx) *
+     $        jfotsout(indexint,1,nlayer) *
      $        (wm(i)*auxjco2(ind+1)+wp(i)*auxjco2(ind)) * 
      $        cortemp(i) *
@@ -846,10 +847,10 @@
          !O2 interpolated coefficient
          jfotsout(indexint,2,auxi) = 
-     $        jfotsout(indexint,2,nlayermx) *
+     $        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,nlayermx) *
+     $        jfotsout(indexint,6,nlayer) *
      $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) * 
      $        cortemp(i)         
@@ -857,15 +858,15 @@
       !Only if chemthermod.ge.2
       if(chemthermod.ge.2) then
-         do i=1,nlayermx
-            auxi = nlayermx-i+1
+         do i=1,nlayer
+            auxi = nlayer-i+1
             ind=auxind(i)
             !NO interpolated coefficient
             jfotsout(indexint,10,auxi) = 
-     $           jfotsout(indexint,10,nlayermx) * 
+     $           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,nlayermx) *
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) * 
      $           cortemp(i)
@@ -886,6 +887,6 @@
 
       indexint=33
-      do i=1,nlayermx
-         auxcolinp(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
       end do
 
@@ -909,22 +910,22 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = 
-     $           jfotsout(indexint,13,nlayermx) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -944,6 +945,6 @@
 
       indexint=34
-      do i=1,nlayermx
-         auxcolinp(nlayermx-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
@@ -970,25 +971,25 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !O2 interpolated coefficient
-         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayermx) * 
+         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,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = 
-     $           jfotsout(indexint,13,nlayermx) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -1008,6 +1009,6 @@
 
       indexint=35
-      do i=1,nlayermx
-         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
       end do
 
@@ -1031,21 +1032,21 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
-         auxi = nlayermx-i+1
+         auxi = nlayer-i+1
          !H2O2 interpolated coefficient
-         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayermx) *
+         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,nlayermx) *
+         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,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = 
-     $           jfotsout(indexint,13,nlayermx) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
@@ -1064,6 +1065,6 @@
 
       indexint=36
-      do i=1,nlayermx
-         auxcolinp(nlayermx-i+1) = o3colx(i) + no2colx(i)
+      do i=1,nlayer
+         auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
       end do
 
@@ -1085,19 +1086,19 @@
       endif
       call interfast
-     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
-      do i=1,nlayermx
+     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
+      do i=1,nlayer
          ind=auxind(i)
          !O3 interpolated coefficient
-         jfotsout(indexint,7,nlayermx-i+1) = 
-     $        jfotsout(indexint,7,nlayermx) *
+         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,nlayermx
+         do i=1,nlayer
             ind=auxind(i)
             !NO2 interpolated coefficient
-            jfotsout(indexint,13,nlayermx-i+1) = 
-     $           jfotsout(indexint,13,nlayermx) *
+            jfotsout(indexint,13,nlayer-i+1) = 
+     $           jfotsout(indexint,13,nlayer) *
      $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
          enddo
Index: trunk/LMDZ.MARS/libf/aeronomars/moldiff.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiff.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiff.F	(revision 1266)
@@ -10,13 +10,4 @@
       USE comcstfi_h
       implicit none
-
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "callkeys.h"
-!#include "comdiurn.h"
-!#include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
-
 
 c
Index: trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1266)
@@ -7,13 +7,5 @@
 implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "callkeys.h"
-!#include "comdiurn.h"
-!#include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 #include "diffusion.h"
-
 
 !
Index: trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff.F	(revision 1266)
@@ -16,11 +16,6 @@
 c
 c=======================================================================
-#include "dimensions.h"
-#include "dimphys.h"
 #include "callkeys.h"
-!#include "comdiurn.h"
 #include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 
 c-----------------------------------------------------------------------
Index: trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F	(revision 1266)
@@ -12,11 +12,6 @@
 c
 c=======================================================================
-#include "dimensions.h"
-#include "dimphys.h"
 #include "callkeys.h"
-!#include "comdiurn.h"
 #include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 #include "diffusion.h"
 
Index: trunk/LMDZ.MARS/libf/aeronomars/molvis.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/molvis.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/molvis.F	(revision 1266)
@@ -19,10 +19,4 @@
 c   declarations:
 c-----------------------------------------------------------------------
-
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "surfdat.h"
-!#include "chimiedata.h"
-!#include "conc.h"
 
 c   arguments:
Index: trunk/LMDZ.MARS/libf/aeronomars/param.h
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param.h	(revision 1265)
+++ 	(revision )
@@ -1,40 +1,0 @@
-!***********************************************
-
-!	param.par
-
-!	Parameters for paramhr.f
-!***********************************************
-
-	integer    ninter
-	parameter  (ninter=36)
-
-	integer    nabs
-	parameter  (nabs=13)
-
-	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 Marte(g)
-	parameter  (masa = 6.4163e26)
-
-	real       radio
-	parameter  (radio = 3390.)           !radio de Marte(km)
-
-	integer	   nreact
-	parameter  (nreact=93)
-
-
-
-
Index: trunk/LMDZ.MARS/libf/aeronomars/param_read.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_read.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/param_read.F	(revision 1266)
@@ -1,12 +1,21 @@
       subroutine param_read
 
+      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,
+     .    efdisco2,efdiso2,efdish2o,
+     .    efdish2o2,efdish2,efdiso3,
+     .    efdiso,efdisn,efdish,
+     .    efdisno,efdisn2,efdisno2,
+     .    efdisco,efionco2,efionn2,
+     .    efionco,efiono3p,efionn,
+     .    efionno,efionh,
+     .    fluxtop,ct1,ct2,p1,p2
+
       implicit none
 
 
 c     common variables and constants
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
       include 'datafile.h'
  
@@ -20,68 +29,4 @@
       
 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
Index: trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F	(revision 1266)
@@ -1,12 +1,23 @@
       subroutine param_read_e107
 
+      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,efionn2,
+     .    efionco,efiono3p,efionn,
+     .    efionno,efionh,
+     .    fluxtop,ct1,ct2,p1,p2
+
       implicit none
 
 
 c     common variables and constants
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
       include 'datafile.h'
       include "callkeys.h"
Index: trunk/LMDZ.MARS/libf/aeronomars/param_v4.h
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_v4.h	(revision 1265)
+++ 	(revision )
@@ -1,210 +1,0 @@
-!**********************************************************************
-
-!	param.cmn
-	
-!	fgg     mar 02          first version
-!**********************************************************************
-
-
-      common/uv/ crscabsi2,c1_16,c17_24,c25_29,c30_31,c32,                         &
-     &   c33,c34,c35,c36,jfotsout,co2crsc195,co2crsc295,                           &
-     &   t0,fluxtop,freccen,jabsifotsintpar,e107,date_e107,e107_tab,               &
-     &   coefit0,coefit1,coefit2,coefit3,coefit4
-     
-     
-
-      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 jfotsout(ninter,nabs,nlayermx)
-      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)
-
-
-      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,rcoef,             &
-     &                     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
-      real*8 rcoef(61,3)
-						  !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.MARS/libf/aeronomars/param_v4_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/param_v4_h.F90	(revision 1266)
+++ trunk/LMDZ.MARS/libf/aeronomars/param_v4_h.F90	(revision 1266)
@@ -0,0 +1,430 @@
+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 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.MARS/libf/aeronomars/paramfoto_compact.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F	(revision 1266)
@@ -2,24 +2,21 @@
 
       subroutine paramfoto_compact
-     $(ig,chemthermod,lswitch,tx,timestep,zenit,zx,rm,nesptherm)
+     $(ig,nlayer,chemthermod,lswitch,tx,timestep,zenit,zx,rm,nesptherm)
   
 c     may 2008      FGG+MALV,GG             
 c**********************************************************************
 
+      use iono_h
+      use param_v4_h
       implicit none
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-      include 'iono.h'
 
 c     arguments
 
-      integer lswitch,ig,nesptherm,chemthermod
-      real zdens(nlayermx)
-      real tx(nlayermx)
+      integer lswitch,ig,nesptherm,chemthermod,nlayer
+      real zdens(nlayer)
+      real tx(nlayer)
       real zenit
-      real zx(nlayermx)
-      real rm(nlayermx,nesptherm)
+      real zx(nlayer)
+      real rm(nlayer,nesptherm)
       real timestep
 
@@ -67,7 +64,7 @@
       real*8  delectx,dnplusx
 
-      real*8  jdistot8(nabs,nlayermx)
-      real*8  jdistot8_b(nabs,nlayermx)
-      real*8  jion8(nabs,nlayermx,4)
+      real*8  jdistot8(nabs,nlayer)
+      real*8  jdistot8_b(nabs,nlayer)
+      real*8  jion8(nabs,nlayer,4)
       real*8  tx8
 	  
@@ -76,11 +73,11 @@
       real*8  alfa_laststep, IonMostAbundant
 
-      real*8  tmin(nlayermx)
+      real*8  tmin(nlayer)
       real*8  fmargin1,critere
 
-      integer compmin(nlayermx)
+      integer compmin(nlayer)
       integer i,j,k
       integer numpasos
-      integer n_comp_en_EQ(nlayermx), paso
+      integer n_comp_en_EQ(nlayer), paso
 
 ! Tracer indexes in the thermospheric chemistry:
@@ -124,5 +121,5 @@
       timefrac_sec=dble(timestep)
 
-      do i=nlayermx,lswitch,-1
+      do i=nlayer,lswitch,-1
 c     Concentrations to real*8
 
@@ -179,5 +176,5 @@
          !from photoabsorption rates and ionization-to-dissociation 
          !branching ratios
-         call phdisrate(ig,chemthermod,zenit,i)   
+         call phdisrate(ig,nlayer,chemthermod,zenit,i)   
          ! Conversion to double precision
          do j=1,nabs
@@ -194,5 +191,5 @@
 		 
          !Lifetimes and temporal integration 
-         call lifetimes(ig,i,chemthermod,zenit,zx,
+         call lifetimes(ig,i,nlayer,chemthermod,zenit,zx,
      $        jdistot8,jdistot8_b,jion8,
      $        tmin(i),compmin(i), 
@@ -220,7 +217,6 @@
         endif   !Of chemthermod.eq.3
 
-         call timemarching ( ig,i,chemthermod,n_comp_en_EQ,compmin,
-     $       tmin,timefrac_sec, deltat,fmargin1)
-
+         call timemarching ( ig,i,nlayer,chemthermod,n_comp_en_EQ,
+     .       compmin,tmin,timefrac_sec, deltat,fmargin1)
 
          !Number of timesteps
@@ -300,5 +296,5 @@
             end if
             !Calculation of productions and losses
-            call prodsandlosses ( ig, i , chemthermod,zenit, zx,
+            call prodsandlosses (ig,i,nlayer,chemthermod,zenit, zx,
      &                 jdistot8, jdistot8_b, jion8,
      &                              co2xinput, o2xinput, o3pxinput,
@@ -430,5 +426,5 @@
             !Third, those species (among the 16 that can be in PE) that are in PE
             call EF_oscilacion
-     &           ( ig,i, paso,chemthermod,zenit, zx,
+     &           ( ig,i,nlayer, paso,chemthermod,zenit, zx,
      &           jdistot8, jdistot8_b,jion8,
      &           deltat,
@@ -1181,20 +1177,26 @@
 c**********************************************************************
 
-      subroutine phdisrate(ig,chemthermod,zenit,i)
+      subroutine phdisrate(ig,nlayer,chemthermod,zenit,i)
 
 c     apr 2002       fgg           first version
 c**********************************************************************
 
+      use param_v4_h, only: ninter,nabs,
+     .    jfotsout,fluxtop,
+     .    jion,jdistot,jdistot_b,
+     .    efdisco2,efdiso2,efdish2o,
+     .    efdish2o2,efdish2,efdiso3,
+     .    efdiso,efdisn,efdish,
+     .    efdisno,efdisn2,efdisno2,
+     .    efdisco,efionco2,efionn2,
+     .    efionco,efiono3p,efionn,
+     .    efionno,efionh
 
       implicit none
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
 
 c     arguments
 
       integer         i                !altitude
-      integer         ig,chemthermod
+      integer         ig,chemthermod,nlayer
       real            zenit
 
@@ -1203,5 +1205,5 @@
       integer         inter,iz,j
       real            lambda
-      real            jdis(nabs,ninter,nlayermx)
+      real            jdis(nabs,ninter,nlayer)
       character*1     dn
 
@@ -1388,10 +1390,15 @@
 c***************************************************************************
 
-	
+        use param_v4_h, only: rcoef,
+     .  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
+
+
+
 	implicit none
-	include "dimensions.h"
-	include "dimphys.h"
-	include 'param.h'
-	include 'param_v4.h'
 
 c Arguments 	
@@ -2010,5 +2017,5 @@
 
       subroutine lifetimes
-     &   ( ig,i,chemthermod,zenit,zx,jdistot8, jdistot8_b, jion8,
+     &   ( ig,i,nlayer,chemthermod,zenit,zx,jdistot8, jdistot8_b, jion8,
      $     xtmin, xcompmin, xn_comp_en_EQ,
      $     co2xini,o2xini,o3pxini,coxini,hxini,ohxini,ho2xini,h2xini,
@@ -2028,22 +2035,20 @@
 c**********************************************************************
 
+      use iono_h
+      use param_v4_h
+
       implicit none
 
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-      include 'iono.h'
       include 'callkeys.h'
 
 c     arguments
 c
-      integer   i,ig                                  ! I. Layer  
+      integer   i,ig,nlayer                               ! I. Layer  
       integer   chemthermod
       real      zenit
-      real      zx(nlayermx)
-      real*8    jdistot8(nabs,nlayermx)                  ! I.
-      real*8    jdistot8_b(nabs,nlayermx)                ! I.
-      real*8    jion8(nabs,nlayermx,4)                ! I.
+      real      zx(nlayer)
+      real*8    jdistot8(nabs,nlayer)                  ! I.
+      real*8    jdistot8_b(nabs,nlayer)                ! I.
+      real*8    jion8(nabs,nlayer,4)                ! I.
 
       real*8    xtmin                         ! O. 
@@ -2714,5 +2719,5 @@
 c**********************************************************************
 
-      subroutine timemarching( ig,i,chemthermod,n_comp_en_EQ,
+      subroutine timemarching( ig,i,nlayer,chemthermod,n_comp_en_EQ,
      $     compmin,tmin,timefrac_sec, deltat,fmargin1 ) 
 
@@ -2725,19 +2730,20 @@
 c**********************************************************************
 
+      use iono_h
+      use param_v4_h, only: 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
+
       implicit none
-
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-      include 'iono.h'
 
 c     arguments
 c
-      integer   i,ig                              ! I. Layer  
+      integer   i,ig,nlayer                              ! I. Layer  
       integer   chemthermod
-      integer   n_comp_en_EQ(nlayermx)         ! Number of species in PE
-      integer   compmin(nlayermx)              ! Species with minimum lifetime
-      real*8    tmin(nlayermx)                 ! Minimum lifetime
+      integer   n_comp_en_EQ(nlayer)         ! Number of species in PE
+      integer   compmin(nlayer)              ! Species with minimum lifetime
+      real*8    tmin(nlayer)                 ! Minimum lifetime
       real*8    timefrac_sec                   ! I. 
       real*8    deltat                         ! O. TimeMarching step
@@ -2893,5 +2899,5 @@
 c**********************************************************************
 
-      subroutine prodsandlosses ( ig,i,chemthermod,zenit,zx,
+      subroutine prodsandlosses ( ig,i,nlayer,chemthermod,zenit,zx,
      &                 jdistot8, jdistot8_b, jion8,
      &                              co2xinput, o2xinput, o3pxinput,
@@ -2916,20 +2922,17 @@
 c**********************************************************************
 
+      use param_v4_h
       implicit none
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
 
 c     arguments
 c
-      integer   ig
+      integer   ig,nlayer
       integer   i                                 ! I. Layer  
       integer   chemthermod
-      real      zx(nlayermx)
+      real      zx(nlayer)
       real      zenit
-      real*8    jdistot8(nabs,nlayermx)
-      real*8    jdistot8_b(nabs,nlayermx)
-      real*8    jion8(nabs,nlayermx,4)
+      real*8    jdistot8(nabs,nlayer)
+      real*8    jdistot8_b(nabs,nlayer)
+      real*8    jion8(nabs,nlayer,4)
       real*8    co2xinput,o2xinput,o3pxinput,coxinput
       real*8    ho2xinput,h2xinput,hxinput,ohxinput
@@ -2939,5 +2942,4 @@
       real*8    cplusxinput,noplusxinput,n2plusxinput,hplusxinput
       real*8    electxinput,nplusxinput,hco2plusxinput
-
 
 c     local variables
@@ -4128,5 +4130,5 @@
 
       subroutine EF_oscilacion
-     &               (ig,i,paso,chemthermod,zenit, zx,
+     &               (ig,i,nlayer,paso,chemthermod,zenit, zx,
      &                 jdistot8, jdistot8_b,jion8, 
      &                 tminaux, 
@@ -4170,25 +4172,28 @@
 c**********************************************************************
 
+      use iono_h
+      use param_v4_h, only: nabs,
+     .  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
+
+
       implicit none
 
-      include "dimensions.h"
-      include "dimphys.h"
-      include 'param.h'
-      include 'param_v4.h'
-      include 'iono.h'
-
-
 c     arguments
       
-      integer   ig
+      integer   ig,nlayer
       integer   i         ! I. Layer                              
       integer   paso      ! I. paso temporal del timemarching, 1,numpasos
       integer   chemthermod
       real*8    tminaux   ! I.
-      real      zx(nlayermx)
+      real      zx(nlayer)
       real      zenit
-      real*8    jdistot8(nabs,nlayermx)                  ! I.
-      real*8    jdistot8_b(nabs,nlayermx)                ! I.
-      real*8    jion8(nabs,nlayermx,4)
+      real*8    jdistot8(nabs,nlayer)                  ! I.
+      real*8    jdistot8_b(nabs,nlayer)                ! I.
+      real*8    jion8(nabs,nlayer,4)
 
       real*8    co2xoutput,o2xoutput,o3pxoutput,coxoutput,h2xoutput
@@ -4277,4 +4282,14 @@
       external ionsec_o2plus
       real*8   ionsec_o2plus
+
+      external avg
+      real*8   avg
+
+      external dif
+      real*8   dif
+
+      real*8 log1
+      real*8 log2
+      real*8 log3
 
 ccccccccccccccc CODE STARTS 
@@ -4712,16 +4727,14 @@
             
             if (o1d_eq(i).eq.'Y') then
-               avg_pares = ( dlog10(o1dxpares(1)) +
-     &              dlog10(o1dxpares(2)) +
-     &              dlog10(o1dxpares(3)) )*0.333
-               dif_pares = ( abs(dlog10(o1dxpares(1))-avg_pares) + 
-     &              abs(dlog10(o1dxpares(2))-avg_pares) +   
-     &              abs(dlog10(o1dxpares(3))-avg_pares) ) * 0.333   
-               avg_impar = ( dlog10(o1dximpar(1)) + 
-     &              dlog10(o1dximpar(2)) +
-     &              dlog10(o1dximpar(3)) )*0.333
-               dif_impar = ( abs(dlog10(o1dximpar(1))-avg_impar) + 
-     &              abs(dlog10(o1dximpar(2))-avg_impar) +   
-     &              abs(dlog10(o1dximpar(3))-avg_impar) ) * 0.333
+               log1 = dlog10(o1dxpares(1))
+               log2 = dlog10(o1dxpares(2))
+               log3 = dlog10(o1dxpares(3))
+               avg_pares = avg(log1,log2,log3)
+               dif_pares = dif(log1,log2,log3,avg_pares)
+               log1 = dlog10(o1dximpar(1))
+               log2 = dlog10(o1dximpar(2))
+               log3 = dlog10(o1dximpar(3))
+               avg_impar = avg(log1,log2,log3)
+               dif_impar = dif(log1,log2,log3,avg_impar)
                dispersion = dif_pares + dif_impar
                dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4733,16 +4746,14 @@
             
             if (oh_eq(i).eq.'Y') then
-               avg_pares = ( dlog10(ohxpares(1)) +
-     &              dlog10(ohxpares(2)) +
-     &              dlog10(ohxpares(3)) )*0.333
-               dif_pares = ( abs(dlog10(ohxpares(1))-avg_pares) + 
-     &              abs(dlog10(ohxpares(2))-avg_pares) +   
-     &              abs(dlog10(ohxpares(3))-avg_pares) ) * 0.333   
-               avg_impar = ( dlog10(ohximpar(1)) + 
-     &              dlog10(ohximpar(2)) +
-     &              dlog10(ohximpar(3)) )*0.333
-               dif_impar = ( abs(dlog10(ohximpar(1))-avg_impar) + 
-     &              abs(dlog10(ohximpar(2))-avg_impar) +   
-     &              abs(dlog10(ohximpar(3))-avg_impar) ) * 0.333
+               log1 = dlog10(ohxpares(1))
+               log2 = dlog10(ohxpares(2))
+               log3 = dlog10(ohxpares(3))
+               avg_pares = avg(log1,log2,log3)
+               dif_pares = dif(log1,log2,log3,avg_pares)
+               log1 = dlog10(ohximpar(1))
+               log2 = dlog10(ohximpar(2))
+               log3 = dlog10(ohximpar(3))
+               avg_impar = avg(log1,log2,log3)
+               dif_impar = dif(log1,log2,log3,avg_impar)
                dispersion = dif_pares + dif_impar
                dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4755,16 +4766,14 @@
             
             if (ho2_eq(i).eq.'Y') then
-               avg_pares = ( dlog10(ho2xpares(1)) +
-     &              dlog10(ho2xpares(2)) +
-     &              dlog10(ho2xpares(3)) )*0.333
-               dif_pares = ( abs(dlog10(ho2xpares(1))-avg_pares) + 
-     &              abs(dlog10(ho2xpares(2))-avg_pares) +   
-     &              abs(dlog10(ho2xpares(3))-avg_pares) ) * 0.333   
-               avg_impar = ( dlog10(ho2ximpar(1)) + 
-     &              dlog10(ho2ximpar(2)) +
-     &              dlog10(ho2ximpar(3)) )*0.333
-               dif_impar = ( abs(dlog10(ho2ximpar(1))-avg_impar) + 
-     &              abs(dlog10(ho2ximpar(2))-avg_impar) +   
-     &              abs(dlog10(ho2ximpar(3))-avg_impar) ) * 0.333
+               log1 = dlog10(ho2xpares(1))
+               log2 = dlog10(ho2xpares(2))
+               log3 = dlog10(ho2xpares(3))
+               avg_pares = avg(log1,log2,log3)
+               dif_pares = dif(log1,log2,log3,avg_pares)
+               log1 = dlog10(ho2ximpar(1))
+               log2 = dlog10(ho2ximpar(2))
+               log3 = dlog10(ho2ximpar(3))
+               avg_impar = avg(log1,log2,log3)
+               dif_impar = dif(log1,log2,log3,avg_impar)
                dispersion = dif_pares + dif_impar
                dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4776,16 +4785,14 @@
             
             if (h_eq(i).eq.'Y') then
-               avg_pares = ( dlog10(hxpares(1)) +
-     &              dlog10(hxpares(2)) +
-     &              dlog10(hxpares(3)) )*0.333
-               dif_pares = ( abs(dlog10(hxpares(1))-avg_pares) + 
-     &              abs(dlog10(hxpares(2))-avg_pares) +   
-     &              abs(dlog10(hxpares(3))-avg_pares) ) * 0.333   
-               avg_impar = ( dlog10(hximpar(1)) + 
-     &              dlog10(hximpar(2)) + 
-     &              dlog10(hximpar(3)) )*0.333
-               dif_impar = ( abs(dlog10(hximpar(1))-avg_impar) + 
-     &              abs(dlog10(hximpar(2))-avg_impar) +   
-     &              abs(dlog10(hximpar(3))-avg_impar) ) * 0.333
+               log1 = dlog10(hxpares(1))
+               log2 = dlog10(hxpares(2))
+               log3 = dlog10(hxpares(3))
+               avg_pares = avg(log1,log2,log3)
+               dif_pares = dif(log1,log2,log3,avg_pares)
+               log1 = dlog10(hximpar(1))
+               log2 = dlog10(hximpar(2))
+               log3 = dlog10(hximpar(3))
+               avg_impar = avg(log1,log2,log3)
+               dif_impar = dif(log1,log2,log3,avg_impar)
                dispersion = dif_pares + dif_impar
                dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4799,16 +4806,14 @@
             if(chemthermod.ge.2) then
                if (n2d_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(n2dxpares(1)) +
-     &                 dlog10(n2dxpares(2)) +
-     &                 dlog10( n2dxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(n2dxpares(1))-avg_pares) + 
-     &                 abs(dlog10(n2dxpares(2))-avg_pares) +   
-     &                 abs(dlog10(n2dxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(n2dximpar(1)) + 
-     &                 dlog10(n2dximpar(2)) +
-     &                 dlog10( n2dximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(n2dximpar(1))-avg_impar) + 
-     &                 abs(dlog10(n2dximpar(2))-avg_impar) +   
-     &                 abs(dlog10(n2dximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(n2dxpares(1))
+                  log2 = dlog10(n2dxpares(2))
+                  log3 = dlog10(n2dxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(n2dximpar(1))
+                  log2 = dlog10(n2dximpar(2))
+                  log3 = dlog10(n2dximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4820,16 +4825,14 @@
             
                if (no2_eq(i).eq.'Y') then
-                  avg_pares = (dlog10(no2xpares(1)) + 
-     &                 dlog10(no2xpares(2)) +
-     &                 dlog10( no2xpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(no2xpares(1))-avg_pares) + 
-     &                 abs(dlog10(no2xpares(2))-avg_pares) +   
-     &                 abs(dlog10(no2xpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(no2ximpar(1)) + 
-     &                 dlog10(no2ximpar(2)) +
-     &                 dlog10( no2ximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(no2ximpar(1))-avg_impar) + 
-     &                 abs(dlog10(no2ximpar(2))-avg_impar) +   
-     &                 abs(dlog10(no2ximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(no2xpares(1))
+                  log2 = dlog10(no2xpares(2))
+                  log3 = dlog10(no2xpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(no2ximpar(1))
+                  log2 = dlog10(no2ximpar(2))
+                  log3 = dlog10(no2ximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4846,16 +4849,14 @@
             if(chemthermod.eq.3) then
                if (cplus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(cplusxpares(1)) +
-     &                 dlog10(cplusxpares(2)) +
-     &                 dlog10(cplusxpares(3)) ) * 0.333
-                  dif_pares = ( abs(dlog10(cplusxpares(1))-avg_pares) + 
-     &                 abs(dlog10(cplusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(cplusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(cplusximpar(1)) + 
-     &                 dlog10(cplusximpar(2)) +
-     &                 dlog10(cplusximpar(3)) ) * 0.333
-                  dif_impar = ( abs(dlog10(cplusximpar(1))-avg_impar) + 
-     &                 abs(dlog10(cplusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(cplusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(cplusxpares(1))
+                  log2 = dlog10(cplusxpares(2))
+                  log3 = dlog10(cplusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(cplusximpar(1))
+                  log2 = dlog10(cplusximpar(2))
+                  log3 = dlog10(cplusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4867,16 +4868,14 @@
                
                if (coplus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(coplusxpares(1)) +
-     &                 dlog10(coplusxpares(2)) +
-     &                 dlog10(coplusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(coplusxpares(1))-avg_pares)+ 
-     &                 abs(dlog10(coplusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(coplusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(coplusximpar(1)) + 
-     &                 dlog10(coplusximpar(2)) +
-     &                 dlog10(coplusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(coplusximpar(1))-avg_impar)+ 
-     &                 abs(dlog10(coplusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(coplusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(coplusxpares(1))
+                  log2 = dlog10(coplusxpares(2))
+                  log3 = dlog10(coplusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(coplusximpar(1))
+                  log2 = dlog10(coplusximpar(2))
+                  log3 = dlog10(coplusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4888,16 +4887,14 @@
                
                if (oplus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(oplusxpares(1)) +
-     &                 dlog10(oplusxpares(2)) +
-     &                 dlog10( oplusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(oplusxpares(1))-avg_pares) + 
-     &                 abs(dlog10(oplusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(oplusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(oplusximpar(1)) + 
-     &                 dlog10(oplusximpar(2)) +
-     &                 dlog10(oplusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(oplusximpar(1))-avg_impar) + 
-     &                 abs(dlog10(oplusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(oplusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(oplusxpares(1))
+                  log2 = dlog10(oplusxpares(2))
+                  log3 = dlog10(oplusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(oplusximpar(1))
+                  log2 = dlog10(oplusximpar(2))
+                  log3 = dlog10(oplusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4909,16 +4906,14 @@
                
                if (n2plus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(n2plusxpares(1)) +
-     &                 dlog10(n2plusxpares(2)) +
-     &                 dlog10(n2plusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(n2plusxpares(1))-avg_pares)+ 
-     &                 abs(dlog10(n2plusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(n2plusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(n2plusximpar(1)) + 
-     &                 dlog10(n2plusximpar(2))+
-     &                 dlog10(n2plusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(n2plusximpar(1))-avg_impar)+ 
-     &                 abs(dlog10(n2plusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(n2plusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(n2plusxpares(1))
+                  log2 = dlog10(n2plusxpares(2))
+                  log3 = dlog10(n2plusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(n2plusximpar(1))
+                  log2 = dlog10(n2plusximpar(2))
+                  log3 = dlog10(n2plusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4930,16 +4925,14 @@
                
                if (hplus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(hplusxpares(1)) +
-     &                 dlog10(hplusxpares(2)) +
-     &                 dlog10( hplusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(hplusxpares(1))-avg_pares) + 
-     &                 abs(dlog10(hplusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(hplusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(hplusximpar(1)) + 
-     &                 dlog10(hplusximpar(2)) +
-     &                 dlog10(hplusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(hplusximpar(1))-avg_impar) + 
-     &                 abs(dlog10(hplusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(hplusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(hplusxpares(1))
+                  log2 = dlog10(hplusxpares(2))
+                  log3 = dlog10(hplusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(hplusximpar(1))
+                  log2 = dlog10(hplusximpar(2))
+                  log3 = dlog10(hplusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4951,16 +4944,14 @@
                
                if (co2plus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(co2plusxpares(1)) +
-     &                 dlog10(co2plusxpares(2)) +
-     &                 dlog10(co2plusxpares(3)) )*0.333
-                  dif_pares = (abs(dlog10(co2plusxpares(1))-avg_pares)+ 
-     &                 abs(dlog10(co2plusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(co2plusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(co2plusximpar(1)) + 
-     &                 dlog10(co2plusximpar(2)) +
-     &                 dlog10(co2plusximpar(3)) )*0.333
-                  dif_impar = (abs(dlog10(co2plusximpar(1))-avg_impar)+ 
-     &                 abs(dlog10(co2plusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(co2plusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(co2plusxpares(1))
+                  log2 = dlog10(co2plusxpares(2))
+                  log3 = dlog10(co2plusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(co2plusximpar(1))
+                  log2 = dlog10(co2plusximpar(2))
+                  log3 = dlog10(co2plusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4972,16 +4963,14 @@
                
                if (o2plus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(o2plusxpares(1)) +
-     &                 dlog10(o2plusxpares(2)) +
-     &                 dlog10(o2plusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(o2plusxpares(1))-avg_pares)+ 
-     &                 abs(dlog10(o2plusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(o2plusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(o2plusximpar(1)) + 
-     &                 dlog10(o2plusximpar(2)) +
-     &                 dlog10(o2plusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(o2plusximpar(1))-avg_impar)+ 
-     &                 abs(dlog10(o2plusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(o2plusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(o2plusxpares(1))
+                  log2 = dlog10(o2plusxpares(2))
+                  log3 = dlog10(o2plusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(o2plusximpar(1))
+                  log2 = dlog10(o2plusximpar(2))
+                  log3 = dlog10(o2plusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -4993,16 +4982,14 @@
                
                if (noplus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(noplusxpares(1)) +
-     &                 dlog10(noplusxpares(2)) +
-     &                 dlog10(noplusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(noplusxpares(1))-avg_pares)+ 
-     &                 abs(dlog10(noplusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(noplusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(noplusximpar(1)) + 
-     &                 dlog10(noplusximpar(2)) +
-     &                 dlog10(noplusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(noplusximpar(1))-avg_impar)+ 
-     &                 abs(dlog10(noplusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(noplusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(noplusxpares(1))
+                  log2 = dlog10(noplusxpares(2))
+                  log3 = dlog10(noplusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(noplusximpar(1))
+                  log2 = dlog10(noplusximpar(2))
+                  log3 = dlog10(noplusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -5014,16 +5001,14 @@
                
                if (nplus_eq(i).eq.'Y') then
-                  avg_pares = ( dlog10(nplusxpares(1)) +
-     &                 dlog10(nplusxpares(2)) +
-     &                 dlog10( nplusxpares(3)) )*0.333
-                  dif_pares = ( abs(dlog10(nplusxpares(1))-avg_pares) + 
-     &                 abs(dlog10(nplusxpares(2))-avg_pares) +   
-     &                 abs(dlog10(nplusxpares(3))-avg_pares) ) * 0.333   
-                  avg_impar = ( dlog10(nplusximpar(1)) + 
-     &                 dlog10(nplusximpar(2)) +
-     &                 dlog10(nplusximpar(3)) )*0.333
-                  dif_impar = ( abs(dlog10(nplusximpar(1))-avg_impar) + 
-     &                 abs(dlog10(nplusximpar(2))-avg_impar) +   
-     &                 abs(dlog10(nplusximpar(3))-avg_impar) ) * 0.333
+                  log1 = dlog10(nplusxpares(1))
+                  log2 = dlog10(nplusxpares(2))
+                  log3 = dlog10(nplusxpares(3))
+                  avg_pares = avg(log1,log2,log3)
+                  dif_pares = dif(log1,log2,log3,avg_pares)
+                  log1 = dlog10(nplusximpar(1))
+                  log2 = dlog10(nplusximpar(2))
+                  log3 = dlog10(nplusximpar(3))
+                  avg_impar = avg(log1,log2,log3)
+                  dif_impar = dif(log1,log2,log3,avg_impar)
                   dispersion = dif_pares + dif_impar
                   dif_pares_impar = abs(avg_pares-avg_impar)
@@ -5035,16 +5020,14 @@
 
                if (hco2plus_eq(i).eq.'Y') then
-                 avg_pares = ( dlog10(hco2plusxpares(1)) +
-     &                dlog10(hco2plusxpares(2)) +
-     &                dlog10(hco2plusxpares(3)) )*0.333
-                 dif_pares = (abs(dlog10(hco2plusxpares(1))-avg_pares)+ 
-     &                abs(dlog10(hco2plusxpares(2))-avg_pares) +   
-     &                abs(dlog10(hco2plusxpares(3))-avg_pares) ) * 0.333   
-                 avg_impar = ( dlog10(hco2plusximpar(1)) + 
-     &                dlog10(hco2plusximpar(2)) +
-     &                dlog10(hco2plusximpar(3)) )*0.333
-                 dif_impar = (abs(dlog10(hco2plusximpar(1))-avg_impar)+ 
-     &                abs(dlog10(hco2plusximpar(2))-avg_impar) +   
-     &                abs(dlog10(hco2plusximpar(3))-avg_impar) ) * 0.333
+                 log1 = dlog10(hco2plusxpares(1))
+                 log2 = dlog10(hco2plusxpares(2))
+                 log3 = dlog10(hco2plusxpares(3))
+                 avg_pares = avg(log1,log2,log3)
+                 dif_pares = dif(log1,log2,log3,avg_pares)
+                 log1 = dlog10(hco2plusximpar(1))
+                 log2 = dlog10(hco2plusximpar(2))
+                 log3 = dlog10(hco2plusximpar(3))
+                 avg_impar = avg(log1,log2,log3)
+                 dif_impar = dif(log1,log2,log3,avg_impar)
                  dispersion = dif_pares + dif_impar
                  dif_pares_impar = abs(avg_pares-avg_impar)
@@ -5120,4 +5103,24 @@
       end
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        function avg(log1,log2,log3)
+        implicit none
+        real*8 avg
+        real*8 log1,log2,log3
+        avg = (log1+log2+log3)*0.333
+        return
+        end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        function dif(log1,log2,log3,avg)
+        implicit none
+        real*8 dif
+        real*8 avg
+        real*8 log1,log2,log3
+        dif = (abs(log1-avg) +
+     &                abs(log2-avg) +
+     &                abs(log3-avg) ) * 0.333
+        return
+        end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 c**********************************************************************
Index: trunk/LMDZ.MARS/libf/aeronomars/perosat.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/perosat.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/perosat.F	(revision 1266)
@@ -25,9 +25,4 @@
 c   -------------
 
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
 c
 c   arguments:
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 1266)
@@ -17,9 +17,6 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "chimiedata.h" 
 #include "callkeys.h"
-!#include "tracer.h"
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -193,6 +190,4 @@
       implicit none 
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "chimiedata.h" 
 #include "callkeys.h" 
@@ -790,6 +785,4 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "chimiedata.h" 
 
@@ -1163,8 +1156,5 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "tracer.h"
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -1268,8 +1258,5 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "tracer.h"
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -1365,7 +1352,4 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
-
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c     inputs:                                                    c
Index: trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 1266)
@@ -20,10 +20,6 @@
 
 #include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "tracer.h"
-!#include "dimradmars.h"
 #include "chimiedata.h"
-!#include "conc.h"
 
 ! input
Index: trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F	(revision 1265)
+++ trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F	(revision 1266)
@@ -9,13 +9,5 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "comdiurn.h"
-!#include "param.h"
-!#include "param_v4.h"
-!#include "chimiedata.h"
-!#include "conc.h"
-
 
       integer,intent(in) :: ngrid ! number of atmospheric columns
