Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 1304)
+++ trunk/LMDZ.GENERIC/README	(revision 1308)
@@ -1032,2 +1032,7 @@
 And some cosmetic cleanup in rain.F90, vdif_kc.F and turbdiff.F90
 
+== 10/07/2014 == EM
+Some cleanup to simplify dynamics/physics interactions by getting rid 
+of dimphys.h (i.e. the nlayermx parameter) and minimizing use of
+dimension.h in the physics.
+
Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 1308)
@@ -30,6 +30,6 @@
 !     pq                Aerosol mixing ratio
 !     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
-!     QREFvis3d(ngrid,nlayermx,naerkind) \ 3d extinction coefficients
-!     QREFir3d(ngrid,nlayermx,naerkind)  / at reference wavelengths
+!     QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
+!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
 !
 !     Output
@@ -40,9 +40,9 @@
 !=======================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "callkeys.h"
 #include "comcstfi.h"
-#include "comvert.h"
+!#include "comvert.h"
 
       INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
@@ -54,9 +54,9 @@
       REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
       REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
-      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayermx,naerkind) ! extinction coefficient in the visible
-      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayermx,naerkind)
+      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
+      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
       REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
       ! BENJAMIN MODIFS
-      real,intent(in) :: cloudfrac(ngrid,nlayermx) ! cloud fraction
+      real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
       real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
       logical,intent(in) :: clearsky
@@ -135,10 +135,10 @@
            iaer=iaero_co2
 !       1. Initialization
-            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+            aerosol(1:ngrid,1:nlayer,iaer)=0.0
 !       2. Opacity calculation
             if (noaero) then ! aerosol set to zero
-             aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+             aerosol(1:ngrid,1:nlayer,iaer)=0.0
             elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
-               aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
+               aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
                !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
             else
@@ -172,8 +172,8 @@
            iaer=iaero_h2o
 !       1. Initialization
-            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+            aerosol(1:ngrid,1:nlayer,iaer)=0.0
 !       2. Opacity calculation
             if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
-               aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9
+               aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
 
                ! put cloud at cloudlvl
@@ -226,5 +226,5 @@
 
                if(CLFvarying)then
-                  call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer))
+                  call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
                   do ig=1, ngrid
                      do l=1,nlayer-1 ! to stop the rad tran bug
@@ -253,5 +253,5 @@
           iaer=iaero_dust
 !         1. Initialization 
-          aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+          aerosol(1:ngrid,1:nlayer,iaer)=0.0
           
           topdust=30.0 ! km  (used to be 10.0 km) LK
@@ -264,5 +264,5 @@
 !             Typical mixing ratio profile
 
-                 zp=(preff/pplay(ig,l))**(70./topdust)
+                 zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
                  expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
 
@@ -277,6 +277,5 @@
 !          Rescaling each layer to reproduce the choosen (or assimilated)
 !          dust extinction opacity at visible reference wavelength, which
-!          is scaled to the "preff" reference surface pressure available in comvert.h
-!          and stored in startfi.nc
+!          is scaled to the surface pressure pplev(ig,1)
 
             taudusttmp(1:ngrid)=0.
@@ -291,5 +290,5 @@
                   aerosol(ig,l,iaer) = max(1E-20, &
                           dusttau &
-                       *  pplev(ig,1) / preff &
+                       *  pplev(ig,1) / pplev(ig,1) &
                        *  aerosol(ig,l,iaer) &
                        /  taudusttmp(ig))
@@ -307,5 +306,5 @@
 
 !       1. Initialization
-         aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+         aerosol(1:ngrid,1:nlayer,iaer)=0.0
 
 
@@ -317,5 +316,5 @@
 !              Typical mixing ratio profile
 
-               zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
+               zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
                expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
 
@@ -335,5 +334,5 @@
                aerosol(ig,l,iaer) = max(1E-20, &
                           1 &
-                       *  pplev(ig,1) / preff &
+                       *  pplev(ig,1) / pplev(ig,1) &
                        *  aerosol(ig,l,iaer) &
                        /  tauh2so4tmp(ig))
@@ -367,5 +366,5 @@
            iaer=iaero_back2lay
 !       1. Initialization
-            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+            aerosol(1:ngrid,1:nlayer,iaer)=0.0
 !       2. Opacity calculation
           DO ig=1,ngrid
Index: trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90	(revision 1308)
@@ -32,6 +32,6 @@
 !     ==============================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "callkeys.h"
 
@@ -60,5 +60,5 @@
       INTEGER :: grid_i,grid_j
 !     Intermediate variable
-      REAL :: var_tmp,var3d_tmp(ngrid,nlayermx)
+      REAL :: var_tmp,var3d_tmp(ngrid,nlayer)
 !     Bilinear interpolation factors
       REAL :: kx,ky,k1,k2,k3,k4
@@ -161,24 +161,24 @@
       INTEGER :: ngrid,nlayer
 !     Aerosol effective radius used for radiative transfer (meter)
-      REAL :: reffrad(ngrid,nlayermx,naerkind)
+      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind)
 !     Aerosol effective variance used for radiative transfer (n.u.)
-      REAL :: nueffrad(ngrid,nlayermx,naerkind)
+      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind)
 
 !     Outputs
 !     -------
 
-      REAL :: QVISsQREF3d(ngrid,nlayermx,L_NSPECTV,naerkind)
-      REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
-      REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
-
-      REAL :: QIRsQREF3d(ngrid,nlayermx,L_NSPECTI,naerkind)
-      REAL :: omegaIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
-      REAL :: gIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
-
-      REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
-      REAL :: QREFir3d(ngrid,nlayermx,naerkind)
-
-      REAL :: omegaREFvis3d(ngrid,nlayermx,naerkind)
-      REAL :: omegaREFir3d(ngrid,nlayermx,naerkind)
+      REAL,INTENT(OUT) :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
+      REAL,INTENT(OUT) :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
+      REAL,INTENT(OUT) :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
+
+      REAL,INTENT(OUT) :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
+      REAL,INTENT(OUT) :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
+      REAL,INTENT(OUT) :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
+
+      REAL,INTENT(OUT) :: QREFvis3d(ngrid,nlayer,naerkind)
+      REAL,INTENT(OUT) :: QREFir3d(ngrid,nlayer,naerkind)
+
+      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
+      REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
 
       DO iaer = 1, naerkind ! Loop on aerosol kind
@@ -724,5 +724,5 @@
                         k2*omegrefIRgrid(grid_i+1,1,iaer)
           ENDIF                  ! --------------------------------
-        ENDDO !nlayermx
+        ENDDO !nlayer
       ENDDO !ngrid
 
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 1308)
@@ -20,6 +20,6 @@
       implicit none
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90	(revision 1308)
@@ -1,3 +1,3 @@
-subroutine calcenergy_kcm(Tsurf,T,Play,Plev,Qsurf,Q,muvar,Eatmtot)
+subroutine calcenergy_kcm(nlayer,Tsurf,T,Play,Plev,Qsurf,Q,muvar,Eatmtot)
 
 
@@ -11,6 +11,6 @@
 !     ----------------------------------------------------------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 !#include "callkeys.h"
@@ -19,8 +19,9 @@
 
   ! inputs
-  real Tsurf,T(1:nlayermx)
-  real Play(1:nlayermx),Plev(1:nlayermx+1)
-  real Qsurf,Q(1:nlayermx)
-  real muvar(1,nlayermx+1)
+  integer,intent(in) :: nlayer
+  real Tsurf,T(1:nlayer)
+  real Play(1:nlayer),Plev(1:nlayer+1)
+  real Qsurf,Q(1:nlayer)
+  real muvar(1,nlayer+1)
 
   ! internal
@@ -29,5 +30,5 @@
   real s_c,rho_v,L
   double precision p_v, s_v, nul
-  real VMR(1:nlayermx)
+  real VMR(1:nlayer)
 
   ! output
@@ -41,5 +42,5 @@
 
 
-  do il=1,nlayermx
+  do il=1,nlayer
      VMR(il)=Q(il)*(muvar(1,il+1)/mH2O)
   end do
@@ -47,5 +48,5 @@
 
   Eatmtot = 0.0
-  do il=1,nlayermx
+  do il=1,nlayer
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1308)
@@ -35,5 +35,5 @@
 !==================================================================
 
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -42,5 +42,5 @@
 !     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
 !     Layer #1 is the layer near the ground. 
-!     Layer #nlayermx is the layer at the top.
+!     Layer #nlayer is the layer at the top.
 
       INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
@@ -52,14 +52,14 @@
       REAL,INTENT(IN) :: emis(ngrid)     ! LW emissivity
       real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle
-      REAL,INTENT(IN) :: pplev(ngrid,nlayermx+1)  ! inter-layer pressure (Pa)
-      REAL,INTENT(IN) :: pplay(ngrid,nlayermx)    ! mid-layer pressure (Pa)
-      REAL,INTENT(IN) :: pt(ngrid,nlayermx)  ! air temperature (K)
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)  ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer)    ! mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: pt(ngrid,nlayer)  ! air temperature (K)
       REAL,INTENT(IN) :: tsurf(ngrid)        ! surface temperature (K)
       REAL,INTENT(IN) :: fract(ngrid)        ! fraction of day
       REAL,INTENT(IN) :: dist_star           ! distance star-planet (AU)
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
-      real,intent(in) :: muvar(ngrid,nlayermx+1)
-      REAL,INTENT(OUT) :: dtlw(ngrid,nlayermx) ! heating rate (K/s) due to LW
-      REAL,INTENT(OUT) :: dtsw(ngrid,nlayermx) ! heating rate (K/s) due to SW
+      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol tau (kg/kg)
+      real,intent(in) :: muvar(ngrid,nlayer+1)
+      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! heating rate (K/s) due to LW
+      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! heating rate (K/s) due to SW
       REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
       REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
@@ -71,5 +71,5 @@
       REAL,INTENT(OUT) :: tau_col(ngrid) ! diagnostic from aeropacity
 !     for H2O cloud fraction in aeropacity
-      real,intent(in) :: cloudfrac(ngrid,nlayermx)
+      real,intent(in) :: cloudfrac(ngrid,nlayer)
       real,intent(out) :: totcloudfrac(ngrid)
       logical,intent(in) :: clearsky
@@ -79,14 +79,14 @@
 !     Globally varying aerosol optical properties on GCM grid
 !     Not needed everywhere so not in radcommon_h
-      REAL :: QVISsQREF3d(ngrid,nlayermx,L_NSPECTV,naerkind)
-      REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
-      REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
-
-      REAL :: QIRsQREF3d(ngrid,nlayermx,L_NSPECTI,naerkind)
-      REAL :: omegaIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
-      REAL :: gIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
-
-!      REAL :: omegaREFvis3d(ngrid,nlayermx,naerkind)
-!      REAL :: omegaREFir3d(ngrid,nlayermx,naerkind) ! not sure of the point of these...
+      REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
+      REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
+      REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
+
+      REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
+      REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
+      REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
+
+!      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
+!      REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
 
       REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
@@ -141,6 +141,6 @@
       real*8,save ::  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
 
-      !REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
-      !REAL :: QREFir3d(ngrid,nlayermx,naerkind)
+      !REAL :: QREFvis3d(ngrid,nlayer,naerkind)
+      !REAL :: QREFir3d(ngrid,nlayer,naerkind)
       !save QREFvis3d, QREFir3d 
       real, dimension(:,:,:), save, allocatable :: QREFvis3d
@@ -177,5 +177,5 @@
 
 !     included by RW for runaway greenhouse 1D study
-      real vtmp(nlayermx)
+      real vtmp(nlayer)
       REAL*8 muvarrad(L_LEVELS)
 
@@ -195,6 +195,6 @@
          !!! ALLOCATED instances are necessary because of CLFvarying
          !!! strategy to call callcorrk twice in physiq...
-         IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayermx,naerkind))
-         IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayermx,naerkind))
+         IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
+         IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
          ! Effective radius and variance of the aerosols
          IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
@@ -207,5 +207,5 @@
             call abort
          endif
-         call su_aer_radii(ngrid,reffrad,nueffrad)
+         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
 	 
 	 
@@ -266,20 +266,20 @@
 
          if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
-	    call co2_reffrad(ngrid,nq,pq,reffrad(1,1,iaero_co2))
-            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
-            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
+	    call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
+            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
+            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
 	 end if
          if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
-	    call h2o_reffrad(ngrid,pq(1,1,igcm_h2o_ice),pt, &
+	    call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
                              reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
-            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
-            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
+            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
+            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
          endif
          if(iaer.eq.iaero_dust)then
-	    call dust_reffrad(ngrid,reffrad(1,1,iaero_dust))
+	    call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
             print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          if(iaer.eq.iaero_h2so4)then
-	    call h2so4_reffrad(ngrid,reffrad(1,1,iaero_h2so4))
+	    call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
             print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
@@ -318,23 +318,23 @@
             do iaer=1,naerkind
                DO nw=1,L_NSPECTV 
-                  do l=1,nlayermx
-
-                     temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)         &
-                         *QREFvis3d(ig,nlayermx+1-l,iaer)
-
-                     temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)    &
-                         *QREFvis3d(ig,max(nlayermx-l,1),iaer)
+                  do l=1,nlayer
+
+                     temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer)         &
+                         *QREFvis3d(ig,nlayer+1-l,iaer)
+
+                     temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
+                         *QREFvis3d(ig,max(nlayer-l,1),iaer)
 
                      qxvaer(2*l,nw,iaer)  = temp1
                      qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
 
-                     temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
-                     temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
+                     temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
+                     temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
 
                      qsvaer(2*l,nw,iaer)  = temp1
                      qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
 
-                     temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
-                     temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
+                     temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
+                     temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
 
                      gvaer(2*l,nw,iaer)  = temp1
@@ -344,11 +344,11 @@
 
                   qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
-                  qxvaer(2*nlayermx+1,nw,iaer)=0.
+                  qxvaer(2*nlayer+1,nw,iaer)=0.
 
                   qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
-                  qsvaer(2*nlayermx+1,nw,iaer)=0.
+                  qsvaer(2*nlayer+1,nw,iaer)=0.
 
                   gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
-                  gvaer(2*nlayermx+1,nw,iaer)=0.
+                  gvaer(2*nlayer+1,nw,iaer)=0.
 
                end do
@@ -356,23 +356,23 @@
 !     longwave
                DO nw=1,L_NSPECTI 
-                  do l=1,nlayermx
-
-                     temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)         &
-                          *QREFir3d(ig,nlayermx+1-l,iaer)
-
-                     temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)    &
-                          *QREFir3d(ig,max(nlayermx-l,1),iaer)
+                  do l=1,nlayer
+
+                     temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer)         &
+                          *QREFir3d(ig,nlayer+1-l,iaer)
+
+                     temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
+                          *QREFir3d(ig,max(nlayer-l,1),iaer)
 
                      qxiaer(2*l,nw,iaer)  = temp1
                      qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
 
-                     temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
-                     temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
+                     temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
+                     temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
 
                      qsiaer(2*l,nw,iaer)  = temp1
                      qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
 
-                     temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
-                     temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
+                     temp1=gir3d(ig,nlayer+1-l,nw,iaer)
+                     temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
 
                      giaer(2*l,nw,iaer)  = temp1
@@ -382,11 +382,11 @@
 
                   qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
-                  qxiaer(2*nlayermx+1,nw,iaer)=0.
+                  qxiaer(2*nlayer+1,nw,iaer)=0.
 
                   qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
-                  qsiaer(2*nlayermx+1,nw,iaer)=0.
+                  qsiaer(2*nlayer+1,nw,iaer)=0.
 
                   giaer(1,nw,iaer)=giaer(2,nw,iaer)
-                  giaer(2*nlayermx+1,nw,iaer)=0.
+                  giaer(2*nlayer+1,nw,iaer)=0.
 
                end do
@@ -481,5 +481,5 @@
       elseif(varfixed)then
 
-         do l=1,nlayermx        ! here we will assign fixed water vapour profiles globally
+         do l=1,nlayer        ! here we will assign fixed water vapour profiles globally
             RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
             if(RH.lt.0.0) RH=0.0
@@ -507,5 +507,5 @@
 !         call watersat(Ttemp,ptemp,qsat)
 
-         qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
+         qvar(2*nlayer+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
  
       else
@@ -541,5 +541,5 @@
       
             muvarrad(1) = muvarrad(2)
-            muvarrad(2*nlayermx+1)=muvar(ig,1)
+            muvarrad(2*nlayer+1)=muvar(ig,1)
 
             print*,'Recalculating qvar with VARIABLE epsi for kastprof'
@@ -579,5 +579,5 @@
       
          muvarrad(1) = muvarrad(2)
-         muvarrad(2*nlayermx+1)=muvar(ig,1)         
+         muvarrad(2*nlayer+1)=muvar(ig,1)         
       endif
       
@@ -607,5 +607,5 @@
 
       tlevrad(1) = tlevrad(2)
-      tlevrad(2*nlayermx+1)=tsurf(ig)
+      tlevrad(2*nlayer+1)=tsurf(ig)
       
       tmid(1) = tlevrad(2)
Index: trunk/LMDZ.GENERIC/libf/phystd/callsedim.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callsedim.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/callsedim.F	(revision 1308)
@@ -28,6 +28,6 @@
 c   -------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -65,5 +65,5 @@
       real epaisseur (ngrid,nlay) ! Layer thickness (m)
       real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
-c      real dens(ngrid,nlayermx) ! Mean density of the ice part. accounting for dust core
+c      real dens(ngrid,nlay) ! Mean density of the ice part. accounting for dust core
 
 
@@ -125,5 +125,5 @@
           if (water.and.(iq.eq.igcm_h2o_ice)) then
             ! compute radii for h2o_ice 
-             call h2o_reffrad(ngrid,zqi(1,1,igcm_h2o_ice),zt,
+             call h2o_reffrad(ngrid,nlay,zqi(1,1,igcm_h2o_ice),zt,
      &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
             ! call sedimentation for h2o_ice
Index: trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 1308)
@@ -32,5 +32,5 @@
 !     ptsrf(ngrid)          Surface temperature
 !     
-!     pdt(ngrid,nlayermx)   Time derivative before condensation/sublimation of pt
+!     pdt(ngrid,nlayer)   Time derivative before condensation/sublimation of pt
 !     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf
 !     pqsurf(ngrid,nq)      Sedimentation flux at the surface (kg.m-2.s-1)
@@ -38,7 +38,7 @@
 !     Outputs
 !     -------
-!     pdpsrf(ngrid)         \  Contribution of condensation/sublimation
-!     pdtc(ngrid,nlayermx)  /  to the time derivatives of Ps, pt, and ptsrf
-!     pdtsrfc(ngrid)       /
+!     pdpsrf(ngrid)       \  Contribution of condensation/sublimation
+!     pdtc(ngrid,nlayer)  /  to the time derivatives of Ps, pt, and ptsrf
+!     pdtsrfc(ngrid)     /
 !     
 !     Both
@@ -56,8 +56,8 @@
 !==================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
-#include "comvert.h"
+!#include "comvert.h"
 #include "callkeys.h"
 
@@ -99,30 +99,30 @@
 
       REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles
-      REAL*8 zt(ngrid,nlayermx)
-      REAL zq(ngrid,nlayermx,nq)
+      REAL*8 zt(ngrid,nlayer)
+      REAL zq(ngrid,nlayer,nq)
       REAL zcpi
-      REAL ztcond (ngrid,nlayermx)
-      REAL ztnuc (ngrid,nlayermx)
+      REAL ztcond (ngrid,nlayer)
+      REAL ztnuc (ngrid,nlayer)
       REAL ztcondsol(ngrid) 
       REAL zdiceco2(ngrid)
-      REAL zcondicea(ngrid,nlayermx), zcondices(ngrid)
+      REAL zcondicea(ngrid,nlayer), zcondices(ngrid)
       REAL zfallice(ngrid), Mfallice(ngrid) 
-      REAL zmflux(nlayermx+1)
-      REAL zu(nlayermx),zv(nlayermx)
+      REAL zmflux(nlayer+1)
+      REAL zu(nlayer),zv(nlayer)
       REAL ztsrf(ngrid)
-      REAL ztc(nlayermx), ztm(nlayermx+1) 
-      REAL zum(nlayermx+1) , zvm(nlayermx+1)
+      REAL ztc(nlayer), ztm(nlayer+1) 
+      REAL zum(nlayer+1) , zvm(nlayer+1)
       LOGICAL condsub(ngrid)
       REAL subptimestep
       Integer Ntime
-      real masse (ngrid,nlayermx), w(ngrid,nlayermx,nq)
-      real wq(ngrid,nlayermx+1)
+      real masse (ngrid,nlayer), w(ngrid,nlayer,nq)
+      real wq(ngrid,nlayer+1)
       real vstokes,reff
 
 !     Special diagnostic variables
-      real tconda1(ngrid,nlayermx)
-      real tconda2(ngrid,nlayermx)
-      real zdtsig (ngrid,nlayermx)
-      real zdt (ngrid,nlayermx)
+      real tconda1(ngrid,nlayer)
+      real tconda2(ngrid,nlayer)
+      real zdtsig (ngrid,nlayer)
+      real zdt (ngrid,nlayer)
 
 !-----------------------------------------------------------------------
@@ -205,5 +205,5 @@
 !     zcondices(ngrid)     condensation rate on the ground  (kg/m2/s)
 !     zfallice(ngrid)      flux of ice falling on surface   (kg/m2/s)
-!     pdtc(ngrid,nlayermx) dT/dt due to phase changes       (K/s)
+!     pdtc(ngrid,nlayer) dT/dt due to phase changes       (K/s)
      
 
@@ -301,5 +301,5 @@
             
 !     sedimentation computed from radius computed from q in module radii_mod
-	 call co2_reffrad(ngrid,nq,zq,reffrad)
+	 call co2_reffrad(ngrid,nlayer,nq,zq,reffrad)
 	 
          do  ilay=1,nlayer
@@ -321,5 +321,5 @@
 !     Computing q after sedimentation
 
-         call vlz_fi(ngrid,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
+         call vlz_fi(ngrid,nlayer,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/convadj.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/convadj.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/convadj.F	(revision 1308)
@@ -29,6 +29,6 @@
 !     ------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -57,9 +57,9 @@
       INTEGER jcnt, jadrs(ngrid)
 
-      REAL sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
-      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
-      REAL zh(ngrid,nlayermx)
-      REAL zu2(ngrid,nlayermx),zv2(ngrid,nlayermx)
-      REAL zh2(ngrid,nlayermx), zhc(ngrid,nlayermx)
+      REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
+      REAL zu(ngrid,nlay),zv(ngrid,nlay)
+      REAL zh(ngrid,nlay)
+      REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
+      REAL zh2(ngrid,nlay), zhc(ngrid,nlay)
       REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
 
@@ -67,5 +67,5 @@
       INTEGER iq,ico2
       save ico2
-      REAL zq(ngrid,nlayermx,nq), zq2(ngrid,nlayermx,nq)
+      REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
       REAL zqm(nq),zqco2m
       real m_co2, m_noco2, A , B
Index: trunk/LMDZ.GENERIC/libf/phystd/dimphys.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/dimphys.h	(revision 1304)
+++ 	(revision )
@@ -1,14 +1,0 @@
-!-----------------------------------------------------------------------
-!   INCLUDE 'dimphys.h'
-
-! ngridmx : number of horizontal grid points
-! note: the -1/jjm term will be 0; unless jj=1
-      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)   
-! nlayermx : number of atmospheric layers
-      integer, parameter :: nlayermx = llm 
-! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
-      !integer, parameter :: nsoilmx = 4 ! for a test 
-      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
-      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
-!-----------------------------------------------------------------------
-
Index: trunk/LMDZ.GENERIC/libf/phystd/evap.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/evap.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/evap.F	(revision 1308)
@@ -1,3 +1,3 @@
-      subroutine evap(ngrid,nq,dtime,pt, pq, pdq, pdt,  
+      subroutine evap(ngrid,nlayer,nq,dtime,pt, pq, pdq, pdt,  
      $     dqevap,dtevap, qevap, tevap) 
         
@@ -7,6 +7,6 @@
       implicit none
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 
@@ -24,18 +24,18 @@
 !==================================================================
 
-      INTEGER ngrid,nq
+      INTEGER ngrid,nlayer,nq
 
 ! Arguments:
-      REAL pt(ngrid,nlayermx) 
-      REAL pq(ngrid,nlayermx,nq) 
-      REAL pdt(ngrid,nlayermx) 
-      REAL pdq(ngrid,nlayermx,nq) 
-      REAL dqevap(ngrid,nlayermx)
-      REAL dtevap(ngrid,nlayermx)
-      REAL qevap(ngrid,nlayermx,nq)
+      REAL pt(ngrid,nlayer) 
+      REAL pq(ngrid,nlayer,nq) 
+      REAL pdt(ngrid,nlayer) 
+      REAL pdq(ngrid,nlayer,nq) 
+      REAL dqevap(ngrid,nlayer)
+      REAL dtevap(ngrid,nlayer)
+      REAL qevap(ngrid,nlayer,nq)
       REAL dtime
  
 ! Local:
-      REAL tevap(ngrid,nlayermx) 
+      REAL tevap(ngrid,nlayer) 
       REAL zlvdcp
       REAL zlsdcp
@@ -47,5 +47,5 @@
 !
 
-      DO l=1,nlayermx
+      DO l=1,nlayer
         DO ig=1,ngrid
          qevap(ig,l,igcm_h2o_vap)=pq(ig,l,igcm_h2o_vap)	
@@ -57,5 +57,5 @@
       ENDDO
 
-      DO l = 1, nlayermx  
+      DO l = 1, nlayer
       	DO ig = 1, ngrid
          zlvdcp=RLVTT/RCPD!/(1.0+RVTMP2*qevap(ig,l,igcm_h2o_vap))
@@ -77,4 +77,3 @@
       ENDDO
 
-      RETURN
       END
Index: trunk/LMDZ.GENERIC/libf/phystd/forceWCfn.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/forceWCfn.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/forceWCfn.F	(revision 1308)
@@ -1,3 +1,3 @@
-      subroutine forceWCfn(ngrid,nq,pplev,pt,dq,dqs)
+      subroutine forceWCfn(ngrid,nlayer,nq,pplev,pt,dq,dqs)
 
       USE tracer_h
@@ -18,17 +18,17 @@
 !==================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 
-      INTEGER ngrid,nq
+      INTEGER ngrid,nlayer,nq
 
       real masse, Wtot, Wdiff
 
-      real pplev(ngrid,nlayermx+1)
+      real pplev(ngrid,nlayer+1)
       real pt(ngrid)
 
       real dqs(ngrid,nq) 
-      real dq(ngrid,nlayermx,nq)
+      real dq(ngrid,nlayer,nq)
 
       integer iq, ig, ilay
@@ -37,5 +37,5 @@
         do ig=1,ngrid
            Wtot = 0.0
-           do ilay=1,nlayermx
+           do ilay=1,nlayer
               masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
               Wtot  = Wtot + masse*dq(ig,ilay,iq)
Index: trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 1308)
@@ -39,6 +39,6 @@
 !==================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/ini_archive.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/ini_archive.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/ini_archive.F	(revision 1308)
@@ -40,5 +40,5 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "paramet.h"
 #include "comconst.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 1308)
@@ -10,4 +10,5 @@
       use comsoil_h, only: ini_comsoil_h
       use control_mod, only: ecritphy
+      use planete_mod, only: nres
       use planetwide_mod, only: planetwide_sumval
 
@@ -51,7 +52,7 @@
       USE ioipsl_getincom 
       IMPLICIT NONE
-#include "dimensions.h"
-#include "dimphys.h"
-#include "planete.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
+!#include "planete.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -88,14 +89,4 @@
       avocado = 6.02214179e23   ! added by RW
 
-! --------------------------------------------------------
-!     The usual Tests
-!     --------------
-      IF (nlayer.NE.nlayermx) THEN
-         PRINT*,'STOP in inifis'
-         PRINT*,'Probleme de dimensions :'
-         PRINT*,'nlayer     = ',nlayer
-         PRINT*,'nlayermx   = ',nlayermx
-         STOP
-      ENDIF
 
       ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
Index: trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F	(revision 1308)
@@ -1,4 +1,7 @@
       SUBROUTINE iniorbit
      $     (papoastr,pperiastr,pyear_day,pperi_day,pobliq)
+      
+      USE planete_mod, only: apoastr, periastr, year_day, obliquit,
+     &                       peri_day, e_elips, p_elips, timeperi
       IMPLICIT NONE
 
@@ -19,5 +22,5 @@
 c   ----------
 c   - Doit etre appele avant d'utiliser orbite.
-c   - initialise une partie du common planete.h
+c   - initialise une partie du common planete_mod
 c
 c   Arguments:
@@ -35,5 +38,5 @@
 c   -------------
 
-#include "planete.h"
+!#include "planete.h"
 #include "comcstfi.h"
 
Index: trunk/LMDZ.GENERIC/libf/phystd/iniphysiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniphysiq.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniphysiq.F90	(revision 1308)
@@ -16,6 +16,9 @@
                        rlatd ! latitudes
 use infotrac, only : nqtot ! number of advected tracers
+use planete_mod, only: ini_planete_mod
 
 implicit none
+include "dimensions.h"
+include "comvert.h"
 
 real,intent(in) :: prad ! radius of the planet (m)
@@ -67,4 +70,7 @@
 rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
 
+! copy over preff , ap() and bp() 
+call ini_planete_mod(nlayer,preff,ap,bp)
+
 ! copy some fundamental parameters to physics 
 ! and do some initializations 
Index: trunk/LMDZ.GENERIC/libf/phystd/initracer.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/initracer.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/initracer.F	(revision 1308)
@@ -24,5 +24,5 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 1308)
@@ -34,5 +34,5 @@
 #include "netcdf.inc"
 #include "serre.h"
-#include"dimphys.h"
+!#include"dimphys.h"
 
 c   Arguments:
Index: trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 1308)
@@ -34,5 +34,5 @@
 #include "netcdf.inc"
 #include "serre.h"
-#include"dimphys.h"
+!#include"dimphys.h"
 
 c   Arguments:
Index: trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 1308)
@@ -7,4 +7,5 @@
   use ioipsl_getincom, only: getin
   use comsaison_h, only: mu0, fract, dist_star
+  use planete_mod
 !  use control_mod
 
@@ -30,8 +31,8 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "callkeys.h"
 #include "comcstfi.h"
-#include "planete.h"
+!#include "planete.h"
 !#include "control.h"
 
@@ -42,8 +43,8 @@
   integer nlayer,nlevel,nq
   integer ilay,ilev,iq,iw,iter
-  real play(nlayermx)     ! Pressure at the middle of the layers [Pa]
-  real zlay(nlayermx)     ! Altitude at middle of the layers [km]
-  real plev(nlayermx+1)   ! Intermediate pressure levels [Pa]
-  real temp(nlayermx)     ! temperature at the middle of the layers [K]
+  real play(llm)     ! Pressure at the middle of the layers [Pa]
+  real zlay(llm)     ! Altitude at middle of the layers [km]
+  real plev(llm+1)   ! Intermediate pressure levels [Pa]
+  real temp(llm)     ! temperature at the middle of the layers [K]
   real,allocatable :: q(:,:)   ! tracer mixing ratio [kg/kg]
   real,allocatable :: vmr(:,:) ! tracer mixing ratio [mol/mol]
@@ -53,8 +54,8 @@
   real emis, albedo
 
-  real muvar(nlayermx+1)
-
-  real dtsw(nlayermx) ! heating rate (K/s) due to SW
-  real dtlw(nlayermx) ! heating rate (K/s) due to LW
+  real muvar(llm+1)
+
+  real dtsw(llm) ! heating rate (K/s) due to SW
+  real dtlw(llm) ! heating rate (K/s) due to LW
   real fluxsurf_lw   ! incident LW flux to surf (W/m2)
   real fluxtop_lw    ! outgoing LW flux to space (W/m2)
@@ -64,12 +65,12 @@
 
   ! not used
-  real reffrad(nlayermx,naerkind)
-  real nueffrad(nlayermx,naerkind)
-  real cloudfrac(nlayermx)
+  real reffrad(llm,naerkind)
+  real nueffrad(llm,naerkind)
+  real cloudfrac(llm)
   real totcloudfrac
   real tau_col
 
   real dTstrat
-  real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
+  real aerosol(llm,naerkind) ! aerosol tau (kg/kg)
   real OLR_nu(1,L_NSPECTI)
   real OSR_nu(1,L_NSPECTV)
@@ -93,5 +94,5 @@
 
 
-  nlayer=nlayermx
+  nlayer=llm
   nlevel=nlayer+1
 
@@ -225,6 +226,6 @@
         ! allocate arrays which depend on number of tracers
         allocate(nametrac(nq))
-        allocate(q(nlayermx,nq))
-        allocate(vmr(nlayermx,nq))
+        allocate(q(nlayer,nq))
+        allocate(vmr(nlayer,nq))
         allocate(qsurf(nq))
 
@@ -268,5 +269,5 @@
 
      !    Create vertical profiles
-     call kcmprof_fn(psurf,qsurf(1),tsurf,    &
+     call kcmprof_fn(nlayer,psurf,qsurf(1),tsurf,    &
           tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
 
@@ -316,5 +317,5 @@
   ! Calculate total atmospheric energy
   Eatmtot=0.0
-  !  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
+  !  call calcenergy_kcm(nlayer,tsurf,temp,play,plev,qsurf,&
   !     q(:,1),muvar,Eatmtot)
 
Index: trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 1308)
@@ -1,3 +1,3 @@
-subroutine kcmprof_fn(psurf_rcm,qsurf_rcm,Tsurf_rcm,Tstra_rcm,P_rcm,Pl_rcm,z_rcm,T_rcm,q_rcm,m_rcm)
+subroutine kcmprof_fn(nlayer,psurf_rcm,qsurf_rcm,Tsurf_rcm,Tstra_rcm,P_rcm,Pl_rcm,z_rcm,T_rcm,q_rcm,m_rcm)
 
 use params_h
@@ -12,6 +12,6 @@
 !     ----------------------------------------------------------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -21,13 +21,14 @@
 
   ! rcm inputs
+  integer nlayer
   real Tsurf_rcm,Tstra_rcm
 
   ! rcm outputs
   real psurf_rcm,qsurf_rcm
-  real P_rcm(1:nlayermx)
-  real Pl_rcm(1:nlayermx+1)
-  real z_rcm(1:nlayermx)
-  real T_rcm(1:nlayermx),q_rcm(1:nlayermx)
-  real m_rcm(1:nlayermx+1)
+  real P_rcm(1:nlayer)
+  real Pl_rcm(1:nlayer+1)
+  real z_rcm(1:nlayer)
+  real T_rcm(1:nlayer),q_rcm(1:nlayer)
+  real m_rcm(1:nlayer+1)
 
   ! rcm for interpolation (should really use log coords?)
@@ -161,13 +162,13 @@
 
   ! define fine pressure grid
-  dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayermx
+  dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayer
 
   P_rcm(1)  = psurf_rcm*exp(dlogp_rcm)
-  do ilay_rcm=1,nlayermx-1
+  do ilay_rcm=1,nlayer-1
      P_rcm(ilay_rcm+1) = P_rcm(ilay_rcm)*exp(dlogp_rcm)
   enddo
 
   Pl_rcm(1) = psurf_rcm
-  do ilev_rcm=2,nlayermx
+  do ilev_rcm=2,nlayer
      ! log-linear interpolation
      Pl_rcm(ilev_rcm) = exp( log( P_rcm(ilev_rcm)*P_rcm(ilev_rcm-1) )/2 )
@@ -318,5 +319,5 @@
   do ilay=2,nlay
 
-     if(ilay_rcm.le.nlayermx)then
+     if(ilay_rcm.le.nlayer)then
      ! interpolate rcm variables
 
@@ -349,10 +350,10 @@
 
   ifinal_rcm=ilay_rcm-1
-  if(ifinal_rcm.lt.nlayermx)then
+  if(ifinal_rcm.lt.nlayer)then
      if(verbose)then
         print*,'Interpolation in kcmprof stopped at layer',ilay_rcm,'!'
      endif
 
-     do ilay_rcm=ifinal_rcm+1,nlayermx
+     do ilay_rcm=ifinal_rcm+1,nlayer
 
         z_rcm(ilay_rcm) = z_rcm(ilay_rcm-1)
@@ -364,5 +365,5 @@
   endif
 
-  do ilay=2,nlayermx
+  do ilay=2,nlayer
      if(T_rcm(ilay).lt.Ttop)then
         T_rcm(ilay)=Ttop
@@ -373,5 +374,5 @@
   if(co2cond)then
      print*,'CO2 condensation haircut - assumes CO2-dominated atmosphere!'
-     do ilay=2,nlayermx
+     do ilay=2,nlayer
         if(P_rcm(ilay).lt.518000.)then
            TCO2cond = (-3167.8)/(log(.01*P_rcm(ilay))-23.23) ! Fanale's formula
Index: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 1308)
@@ -1,4 +1,4 @@
-      subroutine largescale(ngrid,nq,ptimestep, pplev, pplay, pt, pq,   &
-                        pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
+      subroutine largescale(ngrid,nlayer,nq,ptimestep, pplev, pplay,    &
+                    pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
 
 
@@ -23,24 +23,24 @@
 !==================================================================
 
-#include "dimensions.h" 
-#include "dimphys.h"
+!#include "dimensions.h" 
+!#include "dimphys.h"
 #include "comcstfi.h"
 
 #include "callkeys.h"
 
-      INTEGER ngrid,nq
+      INTEGER ngrid,nlayer,nq
 
 !     Arguments
       REAL ptimestep                 ! intervalle du temps (s)
-      REAL pplev(ngrid,nlayermx+1) ! pression a inter-couche
-      REAL pplay(ngrid,nlayermx)   ! pression au milieu de couche
-      REAL pt(ngrid,nlayermx)      ! temperature (K)
-      REAL pq(ngrid,nlayermx,nq) ! tracer mixing ratio (kg/kg)
-      REAL pdt(ngrid,nlayermx)     ! physical temperature tenedency (K/s)
-      REAL pdq(ngrid,nlayermx,nq)! physical tracer tenedency (K/s)
-      REAL pdtlsc(ngrid,nlayermx)  ! incrementation de la temperature (K)
-      REAL pdqvaplsc(ngrid,nlayermx) ! incrementation de la vapeur d'eau
-      REAL pdqliqlsc(ngrid,nlayermx) ! incrementation de l'eau liquide
-      REAL rneb(ngrid,nlayermx)    ! fraction nuageuse
+      REAL pplev(ngrid,nlayer+1) ! pression a inter-couche
+      REAL pplay(ngrid,nlayer)   ! pression au milieu de couche
+      REAL pt(ngrid,nlayer)      ! temperature (K)
+      REAL pq(ngrid,nlayer,nq) ! tracer mixing ratio (kg/kg)
+      REAL pdt(ngrid,nlayer)     ! physical temperature tenedency (K/s)
+      REAL pdq(ngrid,nlayer,nq)! physical tracer tenedency (K/s)
+      REAL pdtlsc(ngrid,nlayer)  ! incrementation de la temperature (K)
+      REAL pdqvaplsc(ngrid,nlayer) ! incrementation de la vapeur d'eau
+      REAL pdqliqlsc(ngrid,nlayer) ! incrementation de l'eau liquide
+      REAL rneb(ngrid,nlayer)    ! fraction nuageuse
 
 
@@ -63,7 +63,7 @@
       
 ! evaporation calculations
-      REAL dqevap(ngrid,nlayermx),dtevap(ngrid,nlayermx)     
-      REAL qevap(ngrid,nlayermx,nq)
-      REAL tevap(ngrid,nlayermx)
+      REAL dqevap(ngrid,nlayer),dtevap(ngrid,nlayer)     
+      REAL qevap(ngrid,nlayer,nq)
+      REAL tevap(ngrid,nlayer)
 
       DOUBLE PRECISION zx_q(ngrid)
@@ -83,13 +83,13 @@
 !     GCM -----> subroutine variables, initialisation of outputs
 
-      pdtlsc(1:ngrid,1:nlayermx)  = 0.0
-      pdqvaplsc(1:ngrid,1:nlayermx)  = 0.0
-      pdqliqlsc(1:ngrid,1:nlayermx) = 0.0
-      rneb(1:ngrid,1:nlayermx) = 0.0
+      pdtlsc(1:ngrid,1:nlayer)  = 0.0
+      pdqvaplsc(1:ngrid,1:nlayer)  = 0.0
+      pdqliqlsc(1:ngrid,1:nlayer) = 0.0
+      rneb(1:ngrid,1:nlayer) = 0.0
       Lcp=RLVTT/RCPD
 
 
       ! Evaporate cloud water/ice
-      call evap(ngrid,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
+      call evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
       ! note: we use qevap but not tevap in largescale/moistadj
             ! otherwise is a big mess
@@ -97,5 +97,5 @@
 
 !  Boucle verticale (du haut vers le bas)
-   DO k = nlayermx, 1, -1
+   DO k = nlayer, 1, -1
 
       zt(1:ngrid)=pt(1:ngrid,k)+(pdt(1:ngrid,k)+dtevap(1:ngrid,k))*ptimestep
@@ -171,7 +171,6 @@
          pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*real(Lcp)
 
-   Enddo ! k= nlayermx, 1, -1
+   Enddo ! k= nlayer, 1, -1
  
 
-      return
       end
Index: trunk/LMDZ.GENERIC/libf/phystd/lect_start_archive.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/lect_start_archive.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/lect_start_archive.F	(revision 1308)
@@ -1,3 +1,4 @@
-      SUBROUTINE lect_start_archive(date,tsurf,tsoil,emis,q2,
+      SUBROUTINE lect_start_archive(ngrid,nlayer,
+     &     date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,h,phisold_newgrid,
      &     q,qsurf,surfith,nid,
@@ -30,5 +31,5 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 !#include "planete.h"
 #include "paramet.h"
@@ -47,4 +48,6 @@
 c=======================================================================
 
+      INTEGER,INTENT(IN) :: ngrid, nlayer
+
 c Old variables dimensions (from file)
 c------------------------------------
@@ -99,13 +102,13 @@
 c variable physique
 c------------------
-      REAL tsurf(ngridmx) ! surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
-      REAL co2ice(ngridmx) ! CO2 ice layer
-      REAL emis(ngridmx)
-      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqtot)
-      REAL tslab(ngridmx,noceanmx)
-      REAL rnat(ngridmx),pctsrf_sic(ngridmx)
-      REAL tsea_ice(ngridmx),sea_ice(ngridmx)
-c     REAL phisfi(ngridmx)
+      REAL tsurf(ngrid) ! surface temperature
+      REAL tsoil(ngrid,nsoilmx) ! soil temperature
+      REAL co2ice(ngrid) ! CO2 ice layer
+      REAL emis(ngrid)
+      REAL q2(ngrid,llm+1),qsurf(ngrid,nqtot)
+      REAL tslab(ngrid,noceanmx)
+      REAL rnat(ngrid),pctsrf_sic(ngrid)
+      REAL tsea_ice(ngrid),sea_ice(ngrid)
+c     REAL phisfi(ngrid)
 
       INTEGER i,j,l
@@ -322,5 +325,5 @@
       allocate(varp1 (imold+1,jmold+1,llm+1))
 
-      write(*,*) 'q2',ngridmx,nlayermx+1
+      write(*,*) 'q2',ngrid,llm+1
       write(*,*) 'q2S',iip1,jjp1,llm+1
       write(*,*) 'q2old',imold+1,jmold+1,lmold+1
@@ -1052,5 +1055,5 @@
       call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsurfs,tsurf)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs,tsurf)
 c     write(44,*) 'tsurf', tsurf
 
@@ -1059,5 +1062,5 @@
 !     &                  imold,jmold,iim,jjm,nsoilmx,
 !     &                   rlonuold,rlatvold,rlonu,rlatv)
-!      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoils,tsoil)
+!      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoils,tsoil)
 c     write(45,*) 'tsoil',tsoil
 
@@ -1065,5 +1068,5 @@
       call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss,emis)
 c     write(46,*) 'emis',emis
 
@@ -1201,5 +1204,5 @@
 
       ! Reshape inertiedatS to scalar grid as inertiedat
-      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,
+      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,
      &                  inertiedatS,inertiedat)
       
@@ -1287,5 +1290,5 @@
 
       ! Reshape tsoilS to scalar grid as tsoil
-       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoilS,tsoil)
+       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS,tsoil)
 
 c-----------------------------------------------------------------------
@@ -1294,21 +1297,21 @@
       call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,noceanmx,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tslabs,tslab)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tslabs,tslab)
 
       call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,rnats,rnat)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,rnats,rnat)
 
       call interp_horiz (pctsrf_sicold,pctsrf_sics,imold,jmold,iim,
      &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,pctsrf_sics,pctsrf_sic)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,pctsrf_sics,pctsrf_sic)
 
       call interp_horiz (tsea_iceold,tsea_ices,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsea_ices,tsea_ice)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsea_ices,tsea_ice)
 
       call interp_horiz (sea_iceold,sea_ices,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,sea_ices,sea_ice)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,sea_ices,sea_ice)
 
 c-----------------------------------------------------------------------
@@ -1334,5 +1337,5 @@
      &                   rlonuold,rlatvold,rlonu,rlatv)
       write (*,*) 'lect_start_archive: q2s ', q2s (1,2,1)  ! INFO
-      call gr_dyn_fi (llm+1,iim+1,jjm+1,ngridmx,q2s,q2)
+      call gr_dyn_fi (llm+1,iim+1,jjm+1,ngrid,q2s,q2)
       write (*,*) 'lect_start_archive: q2 ', q2 (1,2)  ! INFO
 c     write(47,*) 'q2',q2
@@ -1384,5 +1387,5 @@
       enddo
 
-      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
+      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs,qsurf)
 
 c traceurs 3D
@@ -1435,5 +1438,5 @@
       enddo
       
-!      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,co2ices,co2ice)
+!      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,co2ice)
 ! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
 
Index: trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 1308)
@@ -9,4 +9,5 @@
        USE comgeomfi_h
        USE tracer_h
+       USE planete_mod, only: bp
 
        IMPLICIT NONE
@@ -50,9 +51,9 @@
 !    ------------------
 !
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
-#include "comvert.h"
-#include "paramet.h"
+!#include "comvert.h"
+!#include "paramet.h"
 #include "callkeys.h"
 
@@ -84,17 +85,18 @@
 
 !    vertical reorganization of sigma levels
-      REAL zzu(nlayermx),zzv(nlayermx)
-      REAL zzq(nlayermx,nq),zzt(nlayermx)
+      REAL zzu(nlayer),zzv(nlayer)
+      REAL zzq(nlayer,nq),zzt(nlayer)
 !    Dummy variables      
       INTEGER n,l,ig,iq
-      REAL zdtsig(ngrid,nlayermx)
-      REAL zmass(ngrid,nlayermx),zzmass(nlayermx),w(nlayermx+1)
-      REAL zdmass_sum(ngrid,nlayermx+1)
-      REAL zmflux(nlayermx+1)
-      REAL zq1(nlayermx)
+      REAL zdtsig(ngrid,nlayer)
+      REAL zmass(ngrid,nlayer),zzmass(nlayer),w(nlayer+1)
+      REAL zdmass_sum(ngrid,nlayer+1)
+      REAL zmflux(nlayer+1)
+      REAL zq1(nlayer)
       REAL ztsrf(ngrid)
-      REAL ztm(nlayermx+1) 
-      REAL zum(nlayermx+1) , zvm(nlayermx+1)
-      REAL zqm(nlayermx+1,nq),zqm1(nlayermx+1)
+      REAL ztm(nlayer+1) 
+      REAL zum(nlayer+1) , zvm(nlayer+1)
+      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
+      REAL sigma(nlayer+1)
 
 !   local saved variables
@@ -127,6 +129,6 @@
 
       DO ig=1,ngrid
-         zdmass_sum(ig,nlayermx+1)=0.
-         DO l = nlayermx, 1, -1
+         zdmass_sum(ig,nlayer+1)=0.
+         DO l = nlayer, 1, -1
            zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/glat(ig)
 	   zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
@@ -181,26 +183,32 @@
 !  Correction to account for redistribution between sigma or hybrid 
 !  layers when changing surface pressure
-!  zzx quantities have dimension (nlayermx) to speed up calculation
+!  zzx quantities have dimension (nlayer) to speed up calculation
 ! *************************************************************
 
       DO ig=1,ngrid
-         zzt(1:nlayermx)  = pt(ig,1:nlayermx) + pdt(ig,1:nlayermx) * ptimestep
-         zzu(1:nlayermx)  = pu(ig,1:nlayermx) + pdu(ig,1:nlayermx) * ptimestep
-         zzv(1:nlayermx)  = pv(ig,1:nlayermx) + pdv(ig,1:nlayermx) * ptimestep
-         zzq(1:nlayermx,1:nq)=pq(ig,1:nlayermx,1:nq)+pdq(ig,1:nlayermx,1:nq)*ptimestep ! must add the water that has fallen???
+         zzt(1:nlayer)  = pt(ig,1:nlayer) + pdt(ig,1:nlayer) * ptimestep
+         zzu(1:nlayer)  = pu(ig,1:nlayer) + pdu(ig,1:nlayer) * ptimestep
+         zzv(1:nlayer)  = pv(ig,1:nlayer) + pdv(ig,1:nlayer) * ptimestep
+         zzq(1:nlayer,1:nq)=pq(ig,1:nlayer,1:nq)+pdq(ig,1:nlayer,1:nq)*ptimestep ! must add the water that has fallen???
 
 !  Mass fluxes of air through the sigma levels (kg.m-2.s-1)  (>0 when up)
 !  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-
          zmflux(1) = zmassboil(ig)
+         sigma(1)=1
          DO l=1,nlayer
+           ! Ehouarn: shouldn't we rather compute sigma levels than use bp()?
+!           sigma(l+1)=pplev(ig,l+1)/pplev(ig,1)
+!           zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - &
+!                        (sigma(l)-sigma(l+1))*(zdmass_sum(ig,1)+zmflux(1))
+!            if (abs(zmflux(l+1)).lt.1E-13.OR.sigma(l+1).eq.0.) zmflux(l+1)=0.
+           ! Ehouarn: but for now leave things as before
             zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
-            if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
+            if (abs(zmflux(l+1)).lt.1E-13.OR.sigma(l+1).eq.0.) zmflux(l+1)=0.
          END DO
 
 ! Mass of each layer
 ! ------------------ 
-         zzmass(1:nlayermx)=zmass(ig,1:nlayermx)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1))
+         zzmass(1:nlayer)=zmass(ig,1:nlayer)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1))
 
 
@@ -219,10 +227,10 @@
 	 
 !        Van Leer scheme:
-         w(1:nlayermx+1)=-zmflux(1:nlayermx+1)*ptimestep
+         w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
          call vl1d(zzt,2.,zzmass,w,ztm) 
          call vl1d(zzu ,2.,zzmass,w,zum) 
          call vl1d(zzv ,2.,zzmass,w,zvm) 
          do iq=1,nq
-           zq1(1:nlayermx)=zzq(1:nlayermx,iq)
+           zq1(1:nlayer)=zzq(1:nlayer,iq)
            zqm1(1)=zqm(1,iq)
 !		print*,iq
Index: trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 1308)
@@ -1,6 +1,6 @@
-subroutine moistadj(ngrid, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
+subroutine moistadj(ngrid, nlayer, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
 
   use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, RCPV, Psat_water, Lcpdqsat_water
-  USE tracer_h
+  USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
 
   implicit none
@@ -20,28 +20,27 @@
 !=====================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 
-      INTEGER ngrid, nq
-
-      REAL pt(ngrid,nlayermx)            ! temperature (K)
-      REAL pq(ngrid,nlayermx,nq)       ! tracer (kg/kg)
-      REAL pdq(ngrid,nlayermx,nq)
-
-      REAL pdqmana(ngrid,nlayermx,nq)  ! tendency of tracers (kg/kg.s-1)
-      REAL pdtmana(ngrid,nlayermx)       ! temperature increment
+      INTEGER,INTENT(IN) :: ngrid, nlayer, nq
+
+      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K)
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg)
+      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer)   ! mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
+      REAL,INTENT(OUT) :: pdqmana(ngrid,nlayer,nq)  ! tracer tendencies (kg/kg.s-1)
+      REAL,INTENT(OUT) :: pdtmana(ngrid,nlayer) ! temperature increment(K/s)
+      REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction 
 
 !     local variables
-      REAL zt(ngrid,nlayermx)      ! temperature (K)
-      REAL zq(ngrid,nlayermx)      ! humidite specifique (kg/kg)
-      REAL pplev(ngrid,nlayermx+1) ! pression a inter-couche (Pa)
-      REAL pplay(ngrid,nlayermx)   ! pression au milieu de couche (Pa)
-
-      REAL d_t(ngrid,nlayermx)     ! temperature increment
-      REAL d_q(ngrid,nlayermx)     ! incrementation pour vapeur d'eau
-      REAL d_ql(ngrid,nlayermx)    ! incrementation pour l'eau liquide
-      REAL rneb(ngrid,nlayermx) ! cloud fraction 
-      REAL ptimestep
+      REAL zt(ngrid,nlayer)      ! temperature (K)
+      REAL zq(ngrid,nlayer)      ! humidite specifique (kg/kg)
+
+      REAL d_t(ngrid,nlayer)     ! temperature increment
+      REAL d_q(ngrid,nlayer)     ! incrementation pour vapeur d'eau
+      REAL d_ql(ngrid,nlayer)    ! incrementation pour l'eau liquide
 
 !      REAL t_coup
@@ -55,12 +54,12 @@
       INTEGER k1, k1p, k2, k2p
       LOGICAL itest(ngrid)
-      REAL delta_q(ngrid, nlayermx)
-      DOUBLE PRECISION :: cp_new_t(nlayermx), v_cptt(ngrid,nlayermx)
-      REAL cp_delta_t(nlayermx)
-      DOUBLE PRECISION :: v_cptj(nlayermx), v_cptjk1, v_ssig
+      REAL delta_q(ngrid, nlayer)
+      DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer)
+      REAL cp_delta_t(nlayer)
+      DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig
       REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
-      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx),zpsat(ngrid,nlayermx),zdlnpsat(ngrid,nlayermx)
+      REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer)
       REAL zq1(ngrid), zq2(ngrid)
-      DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayermx)
+      DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer)
       DOUBLE PRECISION :: zdp, zdpm
 
@@ -68,5 +67,5 @@
       REAL zflo ! flotabilite
 
-      DOUBLE PRECISION :: local_q(ngrid,nlayermx),local_t(ngrid,nlayermx)
+      DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer)
 
       REAL zdelta, zcor, zcvm5
@@ -79,8 +78,5 @@
       INTEGER,SAVE :: i_ice=0  ! water ice
 
-      LOGICAL firstcall
-      SAVE firstcall
-
-      DATA firstcall /.TRUE./
+      LOGICAL,SAVE :: firstcall=.TRUE.
 
       IF (firstcall) THEN
@@ -96,9 +92,9 @@
 
 !     GCM -----> subroutine variables
-      zq(1:ngrid,1:nlayermx)    = pq(1:ngrid,1:nlayermx,i_h2o)+ pdq(1:ngrid,1:nlayermx,i_h2o)*ptimestep
-      zt(1:ngrid,1:nlayermx)    = pt(1:ngrid,1:nlayermx)
-      pdqmana(1:ngrid,1:nlayermx,1:nq)=0.0
-
-      DO k = 1, nlayermx
+      zq(1:ngrid,1:nlayer)    = pq(1:ngrid,1:nlayer,i_h2o)+ pdq(1:ngrid,1:nlayer,i_h2o)*ptimestep
+      zt(1:ngrid,1:nlayer)    = pt(1:ngrid,1:nlayer)
+      pdqmana(1:ngrid,1:nlayer,1:nq)=0.0
+
+      DO k = 1, nlayer
        DO i = 1, ngrid
          if(zq(i,k).lt.0.)then
@@ -108,13 +104,13 @@
       ENDDO
       
-      local_q(1:ngrid,1:nlayermx) = zq(1:ngrid,1:nlayermx)
-      local_t(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx)
-      rneb(1:ngrid,1:nlayermx) = 0.0
-      d_ql(1:ngrid,1:nlayermx) = 0.0
-      d_t(1:ngrid,1:nlayermx)  = 0.0
-      d_q(1:ngrid,1:nlayermx)  = 0.0
+      local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer)
+      local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
+      rneb(1:ngrid,1:nlayer) = 0.0
+      d_ql(1:ngrid,1:nlayer) = 0.0
+      d_t(1:ngrid,1:nlayer)  = 0.0
+      d_q(1:ngrid,1:nlayer)  = 0.0
 
 !     Calculate v_cptt
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
          DO i = 1, ngrid
             v_cptt(i,k) = RCPD * local_t(i,k)
@@ -128,5 +124,5 @@
 
 !     Calculate Gamma * Cp * dz: (gamma is the critical gradient)
-      DO k = 2, nlayermx
+      DO k = 2, nlayer
          DO i = 1, ngrid
             zdp = pplev(i,k)-pplev(i,k+1)
@@ -159,5 +155,5 @@
   810 CONTINUE ! look for k1, the base of the column
       k2 = k2 + 1
-      IF (k2 .GT. nlayermx) GOTO 9999
+      IF (k2 .GT. nlayer) GOTO 9999
       zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
       zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
@@ -169,5 +165,5 @@
 
   820 CONTINUE !! look for k2, the top of the column
-      IF (k2 .EQ. nlayermx) GOTO 821
+      IF (k2 .EQ. nlayer) GOTO 821
       k2p = k2 + 1
       zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
@@ -227,5 +223,5 @@
 !      ENDDO
 
-      DO k = 2, nlayermx
+      DO k = 2, nlayer
          zdpm = pplev(i,k-1) - pplev(i,k)
          zdp = pplev(i,k) - pplev(i,k+1)
@@ -272,5 +268,5 @@
 ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
 
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
       DO i = 1, ngrid
          IF (itest(i)) THEN
@@ -291,5 +287,5 @@
          ENDIF
       ENDDO
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
       DO i = 1, ngrid
          IF (itest(i)) THEN
@@ -300,5 +296,5 @@
       ENDDO
       ENDDO
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
       DO i = 1, ngrid
          IF (itest(i)) THEN
@@ -308,5 +304,5 @@
       ENDDO
 
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
       DO i = 1, ngrid
           local_q(i, k) = MAX(local_q(i, k), seuil_vap)
@@ -314,5 +310,5 @@
       ENDDO
 
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
       DO i = 1, ngrid
          d_t(i,k) = local_t(i,k) - zt(i,k)
@@ -322,5 +318,5 @@
 
 !     now subroutine -----> GCM variables
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
          DO i = 1, ngrid
             
@@ -333,4 +329,3 @@
 
 
-      RETURN
    END
Index: trunk/LMDZ.GENERIC/libf/phystd/newsedim.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/newsedim.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/newsedim.F	(revision 1308)
@@ -16,6 +16,6 @@
 !   ------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 
@@ -47,7 +47,7 @@
 c    Traceurs :
 c    ~~~~~~~~ 
-      real traversee (ngrid,nlayermx)
-      real vstokes(ngrid,nlayermx)
-      real w(ngrid,nlayermx)
+      real traversee (ngrid,nlay)
+      real vstokes(ngrid,nlay)
+      real w(ngrid,nlay)
       real ptop, dztop, Ep, Stra
 
@@ -191,8 +191,7 @@
       end do
 
-      call vlz_fi(ngrid,pqi,2.,masse,w,wq)
+      call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)
 c         write(*,*) ' newsed: wq(6), wq(7), q(6)',
 c    &                wq(1,6),wq(1,7),pqi(1,6)
 
-      RETURN
       END
Index: trunk/LMDZ.GENERIC/libf/phystd/newstart.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/newstart.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/newstart.F	(revision 1308)
@@ -32,6 +32,7 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
-#include "planete.h"
+!#include "dimphys.h"
+      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
+!#include "planete.h"
 #include "paramet.h"
 #include "comconst.h"
@@ -98,5 +99,5 @@
       real emisread             ! added by RW
       REAL,ALLOCATABLE :: qsurf(:,:)
-      REAL q2(ngridmx,nlayermx+1)
+      REAL q2(ngridmx,llm+1)
 !      REAL rnaturfi(ngridmx)
       real alb(iip1,jjp1),albfi(ngridmx) ! albedos
@@ -171,5 +172,5 @@
 
 !     added by BC for cloud fraction setup
-      REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
+      REAL hice(ngridmx),cloudfrac(ngridmx,llm)
       REAL totalfrac(ngridmx)
 
@@ -341,6 +342,6 @@
         write(*,*) 'Reading file STARTFI'
         fichnom = 'startfi.nc'
-        CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqtot,
-     .        day_ini,time,
+        CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
+     .        nqtot,day_ini,time,
      .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
      .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
@@ -534,5 +535,6 @@
 
         write(*,*) 'Reading file START_ARCHIVE'
-        CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
+        CALL lect_start_archive(ngridmx,llm,
+     &   date,tsurf,tsoil,emis,q2,
      &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
      &   surfith,nid,
@@ -956,5 +958,5 @@
                  ucov(1:iip1,1:jjp1,1:llm)=0.
                  vcov(1:iip1,1:jjm,1:llm)=0.
-                 q2(1:ngridmx,1:nlayermx+1)=0.
+                 q2(1:ngridmx,1:llm+1)=0.
                else
                  write(*,*)'problem reading file ',trim(txt),' !'
@@ -1296,5 +1298,5 @@
           ucov(1:iip1,1:jjp1,1:llm)=0
           vcov(1:iip1,1:jjm,1:llm)=0
-          q2(1:ngridmx,1:nlayermx+1)=0
+          q2(1:ngridmx,1:llm+1)=0
 
 c       radequi : Radiative equilibrium profile of temperatures and no winds
@@ -1326,5 +1328,5 @@
           ucov(1:iip1,1:jjp1,1:llm)=0
           vcov(1:iip1,1:jjm,1:llm)=0
-          q2(1:ngridmx,1:nlayermx+1)=0
+          q2(1:ngridmx,1:llm+1)=0
 
 c       coldstart : T set 1K above CO2 frost point and no winds
@@ -1362,5 +1364,5 @@
           ucov(1:iip1,1:jjp1,1:llm)=0
           vcov(1:iip1,1:jjm,1:llm)=0
-          q2(1:ngridmx,1:nlayermx+1)=0
+          q2(1:ngridmx,1:llm+1)=0
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/newtrelax.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/newtrelax.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/newtrelax.F90	(revision 1308)
@@ -1,8 +1,8 @@
-subroutine newtrelax(ngrid,mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall) 
+subroutine newtrelax(ngrid,nlayer,mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall) 
         
   implicit none
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -21,16 +21,17 @@
 !==================================================================
  
-  integer ngrid
  
   ! Input
-  real mu0(ngrid)                    ! cosine of sun incident angle
-  real sinlat(ngrid)                 ! sine of latitude
-  real temp(ngrid,nlayermx)          ! temperature at each layer (K)
-  real pplay(ngrid,nlayermx)         ! pressure at each layer (Pa)
-  real pplev(ngrid,nlayermx+1)       ! pressure at each level (Pa)
-  real popsk(ngrid,nlayermx)         ! pot. T to T converter
+  integer,intent(in) :: ngrid, nlayer
+  logical,intent(in) :: firstcall
+  real,intent(in) :: mu0(ngrid)            ! cosine of sun incident angle
+  real,intent(in) :: sinlat(ngrid)         ! sine of latitude
+  real,intent(in) :: temp(ngrid,nlayer)    ! temperature at each layer (K)
+  real,intent(in) :: pplay(ngrid,nlayer)   ! pressure at each layer (Pa)
+  real,intent(in) :: pplev(ngrid,nlayer+1) ! pressure at each level (Pa)
+  real,intent(in) :: popsk(ngrid,nlayer)   ! pot. T to T converter
 
   ! Output
-  real dtrad(ngrid,nlayermx) 
+  real,intent(out) :: dtrad(ngrid,nlayer) 
 
   ! Internal
@@ -44,5 +45,4 @@
   real sig, f_sig, sig_trop
   integer l,ig
-  logical firstcall
 
 
@@ -53,5 +53,5 @@
   if(firstcall) then
 
-     ALLOCATE(Trelax(ngrid,nlayermx))
+     ALLOCATE(Trelax(ngrid,nlayer))
 
      print*,'-----------------------------------------------------'
@@ -66,5 +66,5 @@
            T_surf = 126. + 239.*mu0(ig)
            T_trop = 140. + 52.*mu0(ig)
-           do l=1,nlayermx
+           do l=1,nlayer
 
               if(mu0(ig).le.0.0)then ! night side
@@ -86,5 +86,5 @@
         sig_trop=(T_trop/T_surf)**(1./rcp)
 
-        do l=1,nlayermx
+        do l=1,nlayer
            do ig=1,ngrid
 
@@ -109,10 +109,8 @@
      endif
 
-     firstcall=.false.
- 
   endif
 
   ! Calculate radiative forcing
-  do l=1,nlayermx
+  do l=1,nlayer
      do ig=1,ngrid
         dtrad(ig,l) = -(temp(ig,l) - Trelax(ig,l)) / tau_relax
@@ -129,4 +127,3 @@
   !call writediagfi(ngrid,'ThetaZ','stellar zenith angle','deg',2,mu0)
 
-  return
 end subroutine newtrelax
Index: trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90	(revision 1308)
@@ -15,6 +15,6 @@
 
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -67,6 +67,6 @@
 
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -210,6 +210,6 @@
       use slab_ice_h
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -451,6 +451,6 @@
     use slab_ice_h 
 
-#include "dimensions.h"
-#include "comcstfi.h"
+!#include "dimensions.h"
+!#include "comcstfi.h"
 !    INCLUDE "iniprint.h"
 #include "callkeys.h"
@@ -534,5 +534,5 @@
 !  Transport diffusif
 !         IF (ok_soil_hdiff) THEN
-             CALL divgrad_phy(noceanmx,tmp_tslab_loc,dtdiff_loc)
+             CALL divgrad_phy(ngrid,noceanmx,tmp_tslab_loc,dtdiff_loc)
              dtdiff_loc=dtdiff_loc*soil_hdiff*50./SUM(slabh)!*100.
  !           dtdiff_loc(:,1)=dtdiff_loc(:,1)*soil_hdiff*50./SUM(slabh)*0.8
@@ -544,5 +544,5 @@
 ! Calcul  de transport par Ekman
 
-          CALL slab_ekman2(taux_slab,tauy_slab,tslab,dtekman_loc)
+          CALL slab_ekman2(ngrid,taux_slab,tauy_slab,tslab,dtekman_loc)
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/orbite.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/orbite.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/orbite.F	(revision 1308)
@@ -1,3 +1,5 @@
       subroutine orbite(pls,pdist_star,pdecli)
+      
+      use planete_mod, only: p_elips, e_elips, timeperi, obliquit
       implicit none
 
@@ -23,5 +25,5 @@
 c   -------------
 
-#include "planete.h"
+!#include "planete.h"
 #include "comcstfi.h"
 
Index: trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F90	(revision 1308)
@@ -1,3 +1,3 @@
-subroutine phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq, &
+subroutine phyetat0 (ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
                      day_ini,time,tsurf,tsoil, &
                      emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
@@ -20,7 +20,7 @@
 !======================================================================
 #include "netcdf.inc"
-#include "dimensions.h"
-#include "dimphys.h"
-#include "planete.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
+!#include "planete.h"
 #include "comcstfi.h"
 
@@ -33,4 +33,5 @@
 !  inputs:
   integer,intent(in) :: ngrid
+  integer,intent(in) :: nlayer
   character*(*),intent(in) :: fichnom ! "startfi.nc" file
   integer,intent(in) :: tab0
@@ -45,8 +46,8 @@
   real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
   real,intent(out) :: emis(ngrid) ! surface emissivity
-  real,intent(out) :: q2(ngrid, llm+1) ! 
+  real,intent(out) :: q2(ngrid,nlayer+1) ! 
   real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
 ! real co2ice(ngrid) ! co2 ice cover
-  real,intent(out) :: cloudfrac(ngrid,nlayermx)
+  real,intent(out) :: cloudfrac(ngrid,nlayer)
   real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
   real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,noceanmx)  
Index: trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90	(revision 1308)
@@ -17,7 +17,9 @@
                       put_var, put_field, length
   use mod_grid_phy_lmdz, only : klon_glo
+  use planete_mod, only: year_day, periastr, apoastr, peri_day, &
+                         obliquit, z0, lmixmin, emin_turb
 
   implicit none
-#include "planete.h"
+!#include "planete.h"
 #include "comcstfi.h"
   character(len=*), intent(in) :: filename
@@ -147,5 +149,5 @@
 !#include "temps.h"
 #include "comcstfi.h"
-#include "planete.h"
+!#include "planete.h"
 #include "callkeys.h"
 !======================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1308)
@@ -12,5 +12,5 @@
       use gases_h, only: gnom, gfrac
       use radcommon_h, only: sigma, eclipse, glat, grav
-      use radii_mod, only: h2o_reffrad, co2_reffrad, h2o_cloudrad
+      use radii_mod, only: h2o_reffrad, co2_reffrad
       use aerosol_mod, only: iaero_co2, iaero_h2o
       use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe, &
@@ -32,5 +32,6 @@
       use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
       use mod_phys_lmdz_para, only : is_master
-
+      use planete_mod, only: apoastr, periastr, year_day, peri_day, &
+                            obliquit, nres, z0
 
       implicit none
@@ -102,8 +103,8 @@
 !   ------
 !
-!    pdu(ngrid,nlayermx)        \
-!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
-!    pdt(ngrid,nlayermx)         /  variables due to physical processes.
-!    pdq(ngrid,nlayermx)        /
+!    pdu(ngrid,nlayer)        \
+!    pdv(ngrid,nlayer)         \  Temporal derivative of the corresponding
+!    pdt(ngrid,nlayer)         /  variables due to physical processes.
+!    pdq(ngrid,nlayer)        /
 !    pdpsrf(ngrid)             /
 !    tracerdyn                 call tracer in dynamical part of GCM ?
@@ -131,9 +132,9 @@
 !    ------------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "callkeys.h"
 #include "comcstfi.h"
-#include "planete.h"
+!#include "planete.h"
 !#include "control.h"
 #include "netcdf.inc"
@@ -179,5 +180,5 @@
 !     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
 !      aerosol optical properties:
-!      real aerosol(ngrid,nlayermx,naerkind) 
+!      real aerosol(ngrid,nlayer,naerkind) 
 !     this is now internal to callcorrk and hence no longer needed here
 
@@ -207,6 +208,6 @@
 !     aerosol (dust or ice) extinction optical depth  at reference wavelength 
 !     for the "naerkind" optically active aerosols:
-      real aerosol(ngrid,nlayermx,naerkind) 
-      real zh(ngrid,nlayermx)      ! potential temperature (K)
+      real aerosol(ngrid,nlayer,naerkind) 
+      real zh(ngrid,nlayer)      ! potential temperature (K)
 
       character*80 fichier 
@@ -228,56 +229,56 @@
       real zls                       ! solar longitude (rad)
       real zday                      ! date (time since Ls=0, in martian days)
-      real zzlay(ngrid,nlayermx)   ! altitude at the middle of the layers
-      real zzlev(ngrid,nlayermx+1) ! altitude at layer boundaries
+      real zzlay(ngrid,nlayer)   ! altitude at the middle of the layers
+      real zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
       real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
 
 !     Tendencies due to various processes:
       real dqsurf(ngrid,nq)
-      real cldtlw(ngrid,nlayermx)                           ! (K/s) LW heating rate for clear areas
-      real cldtsw(ngrid,nlayermx)                           ! (K/s) SW heating rate for clear areas
+      real cldtlw(ngrid,nlayer)                           ! (K/s) LW heating rate for clear areas
+      real cldtsw(ngrid,nlayer)                           ! (K/s) SW heating rate for clear areas
       real zdtsurf(ngrid)                                   ! (K/s)
-      real dtlscale(ngrid,nlayermx)
-      real zdvdif(ngrid,nlayermx),zdudif(ngrid,nlayermx)  ! (m.s-2)
-      real zdhdif(ngrid,nlayermx), zdtsdif(ngrid)         ! (K/s)
-      real zdtdif(ngrid,nlayermx)			      ! (K/s)
-      real zdvadj(ngrid,nlayermx),zduadj(ngrid,nlayermx)  ! (m.s-2)
-      real zdhadj(ngrid,nlayermx)                           ! (K/s)
-      real zdtgw(ngrid,nlayermx)                            ! (K/s)
-      real zdtmr(ngrid,nlayermx)
-      real zdugw(ngrid,nlayermx),zdvgw(ngrid,nlayermx)    ! (m.s-2)
-      real zdtc(ngrid,nlayermx),zdtsurfc(ngrid)
-      real zdvc(ngrid,nlayermx),zduc(ngrid,nlayermx)
-      real zdumr(ngrid,nlayermx),zdvmr(ngrid,nlayermx)
+      real dtlscale(ngrid,nlayer)
+      real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
+      real zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (K/s)
+      real zdtdif(ngrid,nlayer)			      ! (K/s)
+      real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
+      real zdhadj(ngrid,nlayer)                           ! (K/s)
+      real zdtgw(ngrid,nlayer)                            ! (K/s)
+      real zdtmr(ngrid,nlayer)
+      real zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
+      real zdtc(ngrid,nlayer),zdtsurfc(ngrid)
+      real zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
+      real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)
       real zdtsurfmr(ngrid)
       
-      real zdmassmr(ngrid,nlayermx),zdpsrfmr(ngrid)
+      real zdmassmr(ngrid,nlayer),zdpsrfmr(ngrid)
       real zdmassmr_col(ngrid)
 
-      real zdqdif(ngrid,nlayermx,nq), zdqsdif(ngrid,nq)
-      real zdqevap(ngrid,nlayermx)
-      real zdqsed(ngrid,nlayermx,nq), zdqssed(ngrid,nq)
-      real zdqdev(ngrid,nlayermx,nq), zdqsdev(ngrid,nq)
-      real zdqadj(ngrid,nlayermx,nq)
-      real zdqc(ngrid,nlayermx,nq)
-      real zdqmr(ngrid,nlayermx,nq),zdqsurfmr(ngrid,nq)
-      real zdqlscale(ngrid,nlayermx,nq)
+      real zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
+      real zdqevap(ngrid,nlayer)
+      real zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
+      real zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
+      real zdqadj(ngrid,nlayer,nq)
+      real zdqc(ngrid,nlayer,nq)
+      real zdqmr(ngrid,nlayer,nq),zdqsurfmr(ngrid,nq)
+      real zdqlscale(ngrid,nlayer,nq)
       real zdqslscale(ngrid,nq)
-      real zdqchim(ngrid,nlayermx,nq)
+      real zdqchim(ngrid,nlayer,nq)
       real zdqschim(ngrid,nq)
 
-      real zdteuv(ngrid,nlayermx)    ! (K/s)
-      real zdtconduc(ngrid,nlayermx) ! (K/s)
-      real zdumolvis(ngrid,nlayermx)
-      real zdvmolvis(ngrid,nlayermx)
-      real zdqmoldiff(ngrid,nlayermx,nq)
+      real zdteuv(ngrid,nlayer)    ! (K/s)
+      real zdtconduc(ngrid,nlayer) ! (K/s)
+      real zdumolvis(ngrid,nlayer)
+      real zdvmolvis(ngrid,nlayer)
+      real zdqmoldiff(ngrid,nlayer,nq)
 
 !     Local variables for local calculations:
       real zflubid(ngrid)
-      real zplanck(ngrid),zpopsk(ngrid,nlayermx)
-      real zdum1(ngrid,nlayermx)
-      real zdum2(ngrid,nlayermx)
+      real zplanck(ngrid),zpopsk(ngrid,nlayer)
+      real zdum1(ngrid,nlayer)
+      real zdum2(ngrid,nlayer)
       real ztim1,ztim2,ztim3, z1,z2
       real ztime_fin
-      real zdh(ngrid,nlayermx)
+      real zdh(ngrid,nlayer)
       real gmplanet
       real taux(ngrid),tauy(ngrid)
@@ -288,21 +289,21 @@
 ! local variables only used for diagnostics (output in file "diagfi" or "stats")
 ! ------------------------------------------------------------------------------
-      real ps(ngrid), zt(ngrid,nlayermx)
-      real zu(ngrid,nlayermx),zv(ngrid,nlayermx)
-      real zq(ngrid,nlayermx,nq)
+      real ps(ngrid), zt(ngrid,nlayer)
+      real zu(ngrid,nlayer),zv(ngrid,nlayer)
+      real zq(ngrid,nlayer,nq)
       character*2 str2
       character*5 str5
-      real zdtadj(ngrid,nlayermx)
-      real zdtdyn(ngrid,nlayermx)
+      real zdtadj(ngrid,nlayer)
+      real zdtdyn(ngrid,nlayer)
       real,allocatable,dimension(:,:),save :: ztprevious
-      real reff(ngrid,nlayermx) ! effective dust radius (used if doubleq=T)
+      real reff(ngrid,nlayer) ! effective dust radius (used if doubleq=T)
       real qtot1,qtot2            ! total aerosol mass
       integer igmin, lmin
       logical tdiag
 
-      real zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
       real zstress(ngrid), cd
-      real hco2(nq), tmean, zlocal(nlayermx)
-      real vmr(ngrid,nlayermx) ! volume mixing ratio
+      real hco2(nq), tmean, zlocal(nlayer)
+      real vmr(ngrid,nlayer) ! volume mixing ratio
 
       real time_phys
@@ -318,6 +319,6 @@
 
 !     included by RW for H2O precipitation
-      real zdtrain(ngrid,nlayermx)
-      real zdqrain(ngrid,nlayermx,nq)
+      real zdtrain(ngrid,nlayer)
+      real zdqrain(ngrid,nlayer,nq)
       real zdqsrain(ngrid)
       real zdqssnow(ngrid)
@@ -325,12 +326,12 @@
 
 !     included by RW for H2O Manabe scheme
-      real dtmoist(ngrid,nlayermx)
-      real dqmoist(ngrid,nlayermx,nq)
-
-      real qvap(ngrid,nlayermx)
-      real dqvaplscale(ngrid,nlayermx)
-      real dqcldlscale(ngrid,nlayermx)
-      real rneb_man(ngrid,nlayermx)
-      real rneb_lsc(ngrid,nlayermx)
+      real dtmoist(ngrid,nlayer)
+      real dqmoist(ngrid,nlayer,nq)
+
+      real qvap(ngrid,nlayer)
+      real dqvaplscale(ngrid,nlayer)
+      real dqcldlscale(ngrid,nlayer)
+      real rneb_man(ngrid,nlayer)
+      real rneb_lsc(ngrid,nlayer)
 
 !     included by RW to account for surface cooling by evaporation
@@ -339,10 +340,10 @@
 
 !     to test energy conservation (RW+JL)
-      real mass(ngrid,nlayermx),massarea(ngrid,nlayermx)
+      real mass(ngrid,nlayer),massarea(ngrid,nlayer)
       real dEtot, dEtots, AtmToSurf_TurbFlux
       real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
-      real dEzRadsw(ngrid,nlayermx),dEzRadlw(ngrid,nlayermx),dEzdiff(ngrid,nlayermx)
+      real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
       real dEdiffs(ngrid),dEdiff(ngrid)
-      real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayermx), lscaledEz(ngrid,nlayermx)
+      real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayer), lscaledEz(ngrid,nlayer)
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
       real dtmoist_max,dtmoist_min
@@ -351,8 +352,8 @@
 
 !     included by BC for evaporation     
-      real qevap(ngrid,nlayermx,nq)
-      real tevap(ngrid,nlayermx)
-      real dqevap1(ngrid,nlayermx)
-      real dtevap1(ngrid,nlayermx)
+      real qevap(ngrid,nlayer,nq)
+      real tevap(ngrid,nlayer)
+      real dqevap1(ngrid,nlayer)
+      real dtevap1(ngrid,nlayer)
 
 !     included by BC for hydrology
@@ -367,5 +368,5 @@
 
 !     included by RW for RH diagnostic
-      real qsat(ngrid,nlayermx), RH(ngrid,nlayermx), H2Omaxcol(ngrid),psat_tmp
+      real qsat(ngrid,nlayer), RH(ngrid,nlayer), H2Omaxcol(ngrid),psat_tmp
 
 !     included by RW for hydrology
@@ -378,6 +379,6 @@
 !     included by BC for double radiative transfer call
       logical clearsky
-      real zdtsw1(ngrid,nlayermx)
-      real zdtlw1(ngrid,nlayermx)
+      real zdtsw1(ngrid,nlayer)
+      real zdtlw1(ngrid,nlayer)
       real fluxsurf_lw1(ngrid)     
       real fluxsurf_sw1(ngrid)  
@@ -402,20 +403,20 @@
 
 !     included by RW for temp convadj conservation test
-      real playtest(nlayermx)
-      real plevtest(nlayermx)
-      real ttest(nlayermx)
-      real qtest(nlayermx)
+      real playtest(nlayer)
+      real plevtest(nlayer)
+      real ttest(nlayer)
+      real qtest(nlayer)
       integer igtest
 
 !     included by RW for runway greenhouse 1D study
-      real muvar(ngrid,nlayermx+1)
+      real muvar(ngrid,nlayer+1)
 
 !     included by RW for variable H2O particle sizes
       real,allocatable,dimension(:,:,:),save :: reffrad ! aerosol effective radius (m)
       real,allocatable,dimension(:,:,:),save :: nueffrad ! aerosol effective radius variance
-!      real :: nueffrad_dummy(ngrid,nlayermx,naerkind) !! AS. This is temporary. Check below why.
-      real :: reffh2oliq(ngrid,nlayermx) ! liquid water particles effective radius (m)
-      real :: reffh2oice(ngrid,nlayermx) ! water ice particles effective radius (m)
-   !   real reffH2O(ngrid,nlayermx)
+!      real :: nueffrad_dummy(ngrid,nlayer,naerkind) !! AS. This is temporary. Check below why.
+      real :: reffh2oliq(ngrid,nlayer) ! liquid water particles effective radius (m)
+      real :: reffh2oice(ngrid,nlayer) ! water ice particles effective radius (m)
+   !   real reffH2O(ngrid,nlayer)
       real reffcol(ngrid,naerkind)
 
@@ -469,5 +470,5 @@
         ALLOCATE(rnat(ngrid))         
         ALLOCATE(emis(ngrid))   
-        ALLOCATE(dtrad(ngrid,nlayermx))
+        ALLOCATE(dtrad(ngrid,nlayer))
         ALLOCATE(fluxrad_sky(ngrid))
         ALLOCATE(fluxrad(ngrid))    
@@ -475,12 +476,12 @@
         ALLOCATE(fluxgrd(ngrid))  
         ALLOCATE(qsurf(ngrid,nq))  
-        ALLOCATE(q2(ngrid,nlayermx+1)) 
-        ALLOCATE(ztprevious(ngrid,nlayermx))
-        ALLOCATE(cloudfrac(ngrid,nlayermx))
+        ALLOCATE(q2(ngrid,nlayer+1)) 
+        ALLOCATE(ztprevious(ngrid,nlayer))
+        ALLOCATE(cloudfrac(ngrid,nlayer))
         ALLOCATE(totcloudfrac(ngrid))
         ALLOCATE(hice(ngrid))
         ALLOCATE(qsurf_hist(ngrid,nq))
-        ALLOCATE(reffrad(ngrid,nlayermx,naerkind))
-        ALLOCATE(nueffrad(ngrid,nlayermx,naerkind))
+        ALLOCATE(reffrad(ngrid,nlayer,naerkind))
+        ALLOCATE(nueffrad(ngrid,nlayer,naerkind))
         ALLOCATE(ice_initial(ngrid))
         ALLOCATE(ice_min(ngrid)) 
@@ -494,6 +495,6 @@
         ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
         ALLOCATE(sensibFlux(ngrid))
-        ALLOCATE(zdtlw(ngrid,nlayermx))
-        ALLOCATE(zdtsw(ngrid,nlayermx))
+        ALLOCATE(zdtlw(ngrid,nlayer))
+        ALLOCATE(zdtsw(ngrid,nlayer))
         ALLOCATE(tau_col(ngrid))
         ALLOCATE(pctsrf_sic(ngrid))
@@ -541,5 +542,5 @@
 !        read startfi (initial parameters)
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call phyetat0(ngrid,"startfi.nc",0,0,nsoilmx,nq,   &
+         call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,   &
                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
                cloudfrac,totcloudfrac,hice,  &
@@ -587,5 +588,5 @@
 
 
-           CALL init_masquv(zmasq)
+           CALL init_masquv(ngrid,zmasq)
 
 
@@ -721,10 +722,10 @@
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       
-      pdu(1:ngrid,1:nlayermx) = 0.0
-      pdv(1:ngrid,1:nlayermx) = 0.0
+      pdu(1:ngrid,1:nlayer) = 0.0
+      pdv(1:ngrid,1:nlayer) = 0.0
       if ( .not.nearco2cond ) then
-         pdt(1:ngrid,1:nlayermx) = 0.0
+         pdt(1:ngrid,1:nlayer) = 0.0
       endif
-      pdq(1:ngrid,1:nlayermx,1:nq) = 0.0
+      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
       pdpsrf(1:ngrid)       = 0.0
       zflubid(1:ngrid)      = 0.0
@@ -773,6 +774,6 @@
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      zzlay(1:ngrid,1:nlayermx)=pphi(1:ngrid,1:nlayermx)
-      do l=1,nlayermx         
+      zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
+      do l=1,nlayer         
       zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
       enddo
@@ -892,9 +893,9 @@
               endif
 	      if(water) then
-                  muvar(1:ngrid,1:nlayermx)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayermx,igcm_h2o_vap)) 
-		  muvar(1:ngrid,nlayermx+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayermx,igcm_h2o_vap)) 
+                  muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) 
+		  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 
 		  ! take into account water effect on mean molecular weight
 	      else
-                  muvar(1:ngrid,1:nlayermx+1)=mugaz
+                  muvar(1:ngrid,1:nlayer+1)=mugaz
 	      endif 
       
@@ -950,6 +951,6 @@
                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
                     
-                    zdtlw(ig,1:nlayermx) = ntf*zdtlw1(ig,1:nlayermx) + tf*zdtlw(ig,1:nlayermx)
-                    zdtsw(ig,1:nlayermx) = ntf*zdtsw1(ig,1:nlayermx) + tf*zdtsw(ig,1:nlayermx)
+                    zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
+                    zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
 
                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
@@ -976,5 +977,5 @@
 
 !             Net atmospheric radiative heating rate (K.s-1)
-              dtrad(1:ngrid,1:nlayermx)=zdtsw(1:ngrid,1:nlayermx)+zdtlw(1:ngrid,1:nlayermx)
+              dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
 
            elseif(newtonian)then
@@ -982,5 +983,5 @@
 !          b) Call Newtonian cooling scheme
 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-              call newtrelax(ngrid,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
+              call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
 
               zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
@@ -1001,5 +1002,5 @@
 
 
-              dtrad(1:ngrid,1:nlayermx)=0.0
+              dtrad(1:ngrid,1:nlayer)=0.0
               ! hence no atmospheric radiative heating
 
@@ -1015,5 +1016,5 @@
         zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
-        pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+dtrad(1:ngrid,1:nlayermx)
+        pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
 
 !-------------------------
@@ -1048,6 +1049,6 @@
          zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
 
-         zdum1(1:ngrid,1:nlayermx)=0.0
-         zdum2(1:ngrid,1:nlayermx)=0.0
+         zdum1(1:ngrid,1:nlayer)=0.0
+         zdum2(1:ngrid,1:nlayer)=0.0
 
 
@@ -1066,5 +1067,5 @@
 	 else
 	 
-           zdh(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)/zpopsk(1:ngrid,1:nlayermx)
+           zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
  
            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
@@ -1077,12 +1078,12 @@
               taux,tauy,lastcall)
 
-           zdtdif(1:ngrid,1:nlayermx)=zdhdif(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
-	   zdqevap(1:ngrid,1:nlayermx)=0.
+           zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
+	   zdqevap(1:ngrid,1:nlayer)=0.
 
          end if !turbdiff
 
-         pdv(1:ngrid,1:nlayermx)=pdv(1:ngrid,1:nlayermx)+zdvdif(1:ngrid,1:nlayermx)
-         pdu(1:ngrid,1:nlayermx)=pdu(1:ngrid,1:nlayermx)+zdudif(1:ngrid,1:nlayermx)
-         pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtdif(1:ngrid,1:nlayermx)
+         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
+         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
+         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
          zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
 
@@ -1094,5 +1095,5 @@
 
          if (tracer) then 
-           pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq)
+           pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
            dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
          end if ! of if (tracer)
@@ -1171,8 +1172,8 @@
       if(calladj) then
 
-         zdh(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)/zpopsk(1:ngrid,1:nlayermx)
-         zduadj(1:ngrid,1:nlayermx)=0.0
-         zdvadj(1:ngrid,1:nlayermx)=0.0
-         zdhadj(1:ngrid,1:nlayermx)=0.0
+         zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
+         zduadj(1:ngrid,1:nlayer)=0.0
+         zdvadj(1:ngrid,1:nlayer)=0.0
+         zdhadj(1:ngrid,1:nlayer)=0.0
 
 
@@ -1184,11 +1185,11 @@
               zdqadj)
 
-         pdu(1:ngrid,1:nlayermx) = pdu(1:ngrid,1:nlayermx) + zduadj(1:ngrid,1:nlayermx)
-         pdv(1:ngrid,1:nlayermx) = pdv(1:ngrid,1:nlayermx) + zdvadj(1:ngrid,1:nlayermx)
-         pdt(1:ngrid,1:nlayermx)    = pdt(1:ngrid,1:nlayermx) + zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx)
-         zdtadj(1:ngrid,1:nlayermx) = zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
+         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
+         pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
+         pdt(1:ngrid,1:nlayer)    = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
+         zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
 
          if(tracer) then 
-            pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqadj(1:ngrid,1:nlayermx,1:nq) 
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) 
          end if
 
@@ -1240,10 +1241,10 @@
               zdqc)
 
-         pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtc(1:ngrid,1:nlayermx)
-         pdv(1:ngrid,1:nlayermx)=pdv(1:ngrid,1:nlayermx)+zdvc(1:ngrid,1:nlayermx)
-         pdu(1:ngrid,1:nlayermx)=pdu(1:ngrid,1:nlayermx)+zduc(1:ngrid,1:nlayermx)
+         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
+         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
+         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
          zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
 
-         pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqc(1:ngrid,1:nlayermx,1:nq) 
+         pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq) 
          ! Note: we do not add surface co2ice tendency
          ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
@@ -1283,14 +1284,14 @@
 !             ----------------
 
-	       dqmoist(1:ngrid,1:nlayermx,1:nq)=0.
-	       dtmoist(1:ngrid,1:nlayermx)=0.
-
-	       call moistadj(ngrid,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
-
-               pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap)   &
-	                  +dqmoist(1:ngrid,1:nlayermx,igcm_h2o_vap)
-               pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) =pdq(1:ngrid,1:nlayermx,igcm_h2o_ice)     &
-                          +dqmoist(1:ngrid,1:nlayermx,igcm_h2o_ice)
-               pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+dtmoist(1:ngrid,1:nlayermx)
+	       dqmoist(1:ngrid,1:nlayer,1:nq)=0.
+	       dtmoist(1:ngrid,1:nlayer)=0.
+
+	       call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
+
+               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)   &
+	                  +dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
+                          +dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
+               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
 
                !-------------------------
@@ -1321,9 +1322,9 @@
 !        Large scale condensation/evaporation
 !        --------------------------------
-               call largescale(ngrid,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
-
-               pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+dtlscale(1:ngrid,1:nlayermx)
-               pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayermx)
-               pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) = pdq(1:ngrid,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayermx)
+               call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
+
+               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
 
                !-------------------------
@@ -1363,5 +1364,5 @@
             if(waterrain)then
 
-               zdqrain(1:ngrid,1:nlayermx,1:nq) = 0.0
+               zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0
                zdqsrain(1:ngrid)    = 0.0
                zdqssnow(1:ngrid)    = 0.0
@@ -1370,9 +1371,9 @@
                    zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
 
-               pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) &
-	             +zdqrain(1:ngrid,1:nlayermx,igcm_h2o_vap)
-               pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) = pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) &
-	             +zdqrain(1:ngrid,1:nlayermx,igcm_h2o_ice)
-               pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+zdtrain(1:ngrid,1:nlayermx)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
+	             +zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
+	             +zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
+               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
                dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) ! a bug was here
                dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
@@ -1419,5 +1420,5 @@
 !        ------------- 
         if (sedimentation) then 
-           zdqsed(1:ngrid,1:nlayermx,1:nq) = 0.0
+           zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
            zdqssed(1:ngrid,1:nq)  = 0.0
 
@@ -1455,5 +1456,5 @@
               ! and as in rain.F90, whether it falls as rain or snow depends
               ! only on the surface temperature
-           pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqsed(1:ngrid,1:nlayermx,1:nq)
+           pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
 
@@ -1483,12 +1484,12 @@
          if(mass_redistrib) then
 
-            zdmassmr(1:ngrid,1:nlayermx) = mass(1:ngrid,1:nlayermx) * &
-	       (   zdqevap(1:ngrid,1:nlayermx)                          &
-	         + zdqrain(1:ngrid,1:nlayermx,igcm_h2o_vap)             &
-		 + dqmoist(1:ngrid,1:nlayermx,igcm_h2o_vap)             &
-		 + dqvaplscale(1:ngrid,1:nlayermx) )
+            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * &
+	       (   zdqevap(1:ngrid,1:nlayer)                          &
+	         + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
+		 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
+		 + dqvaplscale(1:ngrid,1:nlayer) )
 
             do ig = 1, ngrid
-	       zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayermx))
+	       zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
             enddo
             
@@ -1503,13 +1504,13 @@
 	 
 
-            pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqmr(1:ngrid,1:nlayermx,1:nq)
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
             dqsurf(1:ngrid,1:nq)         = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
-            pdt(1:ngrid,1:nlayermx)        = pdt(1:ngrid,1:nlayermx) + zdtmr(1:ngrid,1:nlayermx)
-            pdu(1:ngrid,1:nlayermx)        = pdu(1:ngrid,1:nlayermx) + zdumr(1:ngrid,1:nlayermx)
-            pdv(1:ngrid,1:nlayermx)        = pdv(1:ngrid,1:nlayermx) + zdvmr(1:ngrid,1:nlayermx)
+            pdt(1:ngrid,1:nlayer)        = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
+            pdu(1:ngrid,1:nlayer)        = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
+            pdv(1:ngrid,1:nlayer)        = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
 	    pdpsrf(1:ngrid)                = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
             zdtsurf(1:ngrid)               = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
 	    
-!	    print*,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayermx,igcm_h2o_vap) 
+!	    print*,'after mass redistrib, q=',pq(211,1:nlayer,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayer,igcm_h2o_vap) 
 	 endif
 
@@ -1672,14 +1673,14 @@
  
       ! temperature, zonal and meridional wind
-      zt(1:ngrid,1:nlayermx) = pt(1:ngrid,1:nlayermx) + pdt(1:ngrid,1:nlayermx)*ptimestep
-      zu(1:ngrid,1:nlayermx) = pu(1:ngrid,1:nlayermx) + pdu(1:ngrid,1:nlayermx)*ptimestep
-      zv(1:ngrid,1:nlayermx) = pv(1:ngrid,1:nlayermx) + pdv(1:ngrid,1:nlayermx)*ptimestep
+      zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
+      zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
+      zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
 
       ! diagnostic
-      zdtdyn(1:ngrid,1:nlayermx)     = pt(1:ngrid,1:nlayermx)-ztprevious(1:ngrid,1:nlayermx)
-      ztprevious(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx)
+      zdtdyn(1:ngrid,1:nlayer)     = pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)
+      ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
 
       if(firstcall)then
-         zdtdyn(1:ngrid,1:nlayermx)=0.0
+         zdtdyn(1:ngrid,1:nlayer)=0.0
       endif
 
@@ -1690,5 +1691,5 @@
 
       ! tracers
-      zq(1:ngrid,1:nlayermx,1:nq) = pq(1:ngrid,1:nlayermx,1:nq) + pdq(1:ngrid,1:nlayermx,1:nq)*ptimestep
+      zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
 
       ! surface pressure
@@ -1798,5 +1799,5 @@
          do iq=1,nq
            do ig=1,ngrid
-              qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
+              qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
            enddo
          enddo
@@ -1805,14 +1806,14 @@
          reffcol(1:ngrid,1:naerkind)=0.0
          if(co2cond.and.(iaero_co2.ne.0))then
-            call co2_reffrad(ngrid,nq,zq,reffrad(1,1,iaero_co2))
+            call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2))
             do ig=1,ngrid
-               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
+               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)*reffrad(ig,1:nlayer,iaero_co2)*mass(ig,1:nlayer))
             enddo
          endif
          if(water.and.(iaero_h2o.ne.0))then
-            call h2o_reffrad(ngrid,zq(1,1,igcm_h2o_ice),zt, &
+            call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, &
                              reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
             do ig=1,ngrid
-               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
+               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)*reffrad(ig,1:nlayer,iaero_h2o)*mass(ig,1:nlayer))
             enddo
          endif
@@ -1989,5 +1990,5 @@
               end do
             if (water) then
-               vmr=zq(1:ngrid,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
+               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
                call wstats(ngrid,"vmr_h2ovapor",                           &
                           "H2O vapour volume mixing ratio","mol/mol",       & 
Index: trunk/LMDZ.GENERIC/libf/phystd/planete.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/planete.h	(revision 1304)
+++ 	(revision )
@@ -1,15 +1,0 @@
-!-----------------------------------------------------------------------
-! INCLUDE planet.h
-
-      COMMON/planet/apoastr,periastr,year_day,peri_day,                 &
-     &       obliquit,nres,                                             &
-     &       z0,lmixmin,emin_turb,coefvis,coefir,                       &
-     &       timeperi,e_elips,p_elips
-
-      real apoastr,periastr,year_day,peri_day,                          &
-     &     obliquit,nres,                                               &
-     &     z0,lmixmin,emin_turb,coefvis,coefir,                         &
-     &       timeperi,e_elips,p_elips
-      
-
-!-----------------------------------------------------------------------
Index: trunk/LMDZ.GENERIC/libf/phystd/planete_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/planete_mod.F90	(revision 1308)
+++ trunk/LMDZ.GENERIC/libf/phystd/planete_mod.F90	(revision 1308)
@@ -0,0 +1,42 @@
+MODULE planete_mod
+  IMPLICIT NONE
+  
+  REAL :: apoastr ! maximum star-planet distance (AU)
+  REAL :: periastr ! minimum star-planet distance (AU)
+  REAL :: year_day ! length of year (sols)
+  REAL :: peri_day ! date of periastron (sols since N. spring)
+  REAL :: obliquit ! Obliquity of the planet (deg)
+  REAL :: nres ! tidal resonance ratio
+  REAL :: z0 ! surface roughness (m)
+  REAL :: lmixmin ! mixing length
+  REAL :: emin_turb ! minimal energy
+  REAL :: coefvis
+  REAL :: coefir
+  REAL :: timeperi
+  REAL :: e_elips
+  REAL :: p_elips
+  
+  REAL :: preff ! reference surface pressure (Pa)
+  REAL,ALLOCATABLE :: ap(:) ! hybrid coordinate at layer interface
+  REAL,ALLOCATABLE :: bp(:) ! hybrid coordinate at layer interface
+  
+  CONTAINS
+  
+  subroutine ini_planete_mod(nlayer,preff_dyn,ap_dyn,bp_dyn)
+  
+  implicit none
+  integer,intent(in) :: nlayer ! number of atmospheric layers
+  real,intent(in) :: preff_dyn ! reference surface pressure (Pa)
+  real,intent(in) :: ap_dyn(nlayer+1) ! hybrid coordinate at interfaces
+  real,intent(in) :: bp_dyn(nlayer+1) ! hybrid coordinate at interfaces
+  
+  allocate(ap(nlayer+1))
+  allocate(bp(nlayer+1))
+  
+  preff=preff_dyn
+  ap(:)=ap_dyn(:)
+  bp(:)=bp_dyn(:)
+  
+  end subroutine ini_planete_mod
+  
+END MODULE planete_mod
Index: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 1308)
@@ -20,5 +20,5 @@
 
 !==================================================================
-   subroutine su_aer_radii(ngrid,reffrad,nueffrad)
+   subroutine su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
 !==================================================================
 !     Purpose
@@ -39,11 +39,12 @@
 
       include "callkeys.h"
-      include "dimensions.h"
-      include "dimphys.h"
-
-      integer,intent(in) :: ngrid
-
-      real, intent(out) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
-      real, intent(out) :: nueffrad(ngrid,nlayermx,naerkind)     !variance     
+!      include "dimensions.h"
+!      include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+      integer,intent(in) :: nlayer
+
+      real, intent(out) :: reffrad(ngrid,nlayer,naerkind)      !aerosols radii (K)
+      real, intent(out) :: nueffrad(ngrid,nlayer,naerkind)     !variance     
 
       logical, save :: firstcall=.true.
@@ -58,26 +59,26 @@
 
             if(iaer.eq.iaero_co2)then ! CO2 ice
-               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-4
-               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
+               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4
+               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
 
             if(iaer.eq.iaero_h2o)then ! H2O ice
-               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-5
-               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
+               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
+               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
 
             if(iaer.eq.iaero_dust)then ! dust
-               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-5
-               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
+               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
+               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
  
             if(iaer.eq.iaero_h2so4)then ! H2O ice
-               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-6
-               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
+               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-6
+               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
 	    
 	    if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
-               reffrad(1:ngrid,1:nlayermx,iaer) = 2.e-6
-               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
+               reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
+               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
 
@@ -126,5 +127,5 @@
 
 !==================================================================
-   subroutine h2o_reffrad(ngrid,pq,pt,reffrad,nueffrad)
+   subroutine h2o_reffrad(ngrid,nlayer,pq,pt,reffrad,nueffrad)
 !==================================================================
 !     Purpose
@@ -141,14 +142,15 @@
 
       include "callkeys.h"
-      include "dimensions.h"
-      include "dimphys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
       include "comcstfi.h"
 
       integer,intent(in) :: ngrid
-
-      real, intent(in) :: pq(ngrid,nlayermx) !water ice mixing ratios (kg/kg)
-      real, intent(in) :: pt(ngrid,nlayermx) !temperature (K)
-      real, intent(out) :: reffrad(ngrid,nlayermx)      !aerosol radii
-      real, intent(out) :: nueffrad(ngrid,nlayermx) ! dispersion      
+      integer,intent(in) :: nlayer
+
+      real, intent(in) :: pq(ngrid,nlayer) !water ice mixing ratios (kg/kg)
+      real, intent(in) :: pt(ngrid,nlayer) !temperature (K)
+      real, intent(out) :: reffrad(ngrid,nlayer)      !aerosol radii
+      real, intent(out) :: nueffrad(ngrid,nlayer) ! dispersion      
 
       integer :: ig,l
@@ -158,5 +160,5 @@
 
       if (radfixed) then
-         do l=1,nlayermx
+         do l=1,nlayer
             do ig=1,ngrid
                zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
@@ -167,5 +169,5 @@
          enddo
       else
-         do l=1,nlayermx
+         do l=1,nlayer
             do ig=1,ngrid
                zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
@@ -186,5 +188,5 @@
 
 !==================================================================
-   subroutine h2o_cloudrad(ngrid,pql,reffliq,reffice)
+   subroutine h2o_cloudrad(ngrid,nlayer,pql,reffliq,reffice)
 !==================================================================
 !     Purpose
@@ -201,12 +203,13 @@
 
       include "callkeys.h"
-      include "dimensions.h"
-      include "dimphys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
       include "comcstfi.h"
 
       integer,intent(in) :: ngrid
-
-      real, intent(in) :: pql(ngrid,nlayermx) !condensed water mixing ratios (kg/kg)
-      real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (m)
+      integer,intent(in) :: nlayer
+
+      real, intent(in) :: pql(ngrid,nlayer) !condensed water mixing ratios (kg/kg)
+      real, intent(out) :: reffliq(ngrid,nlayer),reffice(ngrid,nlayer)     !liquid and ice water particle radii (m)
 
       real,external :: CBRT            
@@ -214,8 +217,8 @@
 
       if (radfixed) then
-	 reffliq(1:ngrid,1:nlayermx)= rad_h2o
-	 reffice(1:ngrid,1:nlayermx)= rad_h2o_ice
+	 reffliq(1:ngrid,1:nlayer)= rad_h2o
+	 reffice(1:ngrid,1:nlayer)= rad_h2o_ice
       else
-	 do k=1,nlayermx
+	 do k=1,nlayer
 	   do i=1,ngrid
 	     reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
@@ -234,5 +237,5 @@
 
 !==================================================================
-   subroutine co2_reffrad(ngrid,nq,pq,reffrad)
+   subroutine co2_reffrad(ngrid,nlayer,nq,pq,reffrad)
 !==================================================================
 !     Purpose
@@ -249,12 +252,12 @@
 
       include "callkeys.h"
-      include "dimensions.h"
-      include "dimphys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
       include "comcstfi.h"
 
-      integer,intent(in) :: ngrid,nq
-
-      real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
-      real, intent(out) :: reffrad(ngrid,nlayermx)      !co2 ice particles radii (K)
+      integer,intent(in) :: ngrid,nlayer,nq
+
+      real, intent(in) :: pq(ngrid,nlayer,nq) !tracer mixing ratios (kg/kg)
+      real, intent(out) :: reffrad(ngrid,nlayer)      !co2 ice particles radii (m)
 
       integer :: ig,l
@@ -265,7 +268,7 @@
 
       if (radfixed) then
-         reffrad(1:ngrid,1:nlayermx) = 5.e-5 ! CO2 ice
+         reffrad(1:ngrid,1:nlayer) = 5.e-5 ! CO2 ice
       else
-         do l=1,nlayermx
+         do l=1,nlayer
             do ig=1,ngrid
                zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
@@ -281,5 +284,5 @@
 
 !==================================================================
-   subroutine dust_reffrad(ngrid,reffrad)
+   subroutine dust_reffrad(ngrid,nlayer,reffrad)
 !==================================================================
 !     Purpose
@@ -294,13 +297,14 @@
       Implicit none
 
-      include "callkeys.h"
-      include "dimensions.h"
-      include "dimphys.h"
-
-      integer,intent(in) :: ngrid
-
-      real, intent(out) :: reffrad(ngrid,nlayermx)      !dust particles radii (K)
+!      include "callkeys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+      integer,intent(in) :: nlayer
+
+      real, intent(out) :: reffrad(ngrid,nlayer)      !dust particles radii (m)
             
-      reffrad(1:ngrid,1:nlayermx) = 2.e-6 ! dust
+      reffrad(1:ngrid,1:nlayer) = 2.e-6 ! dust
 
    end subroutine dust_reffrad
@@ -309,5 +313,5 @@
 
 !==================================================================
-   subroutine h2so4_reffrad(ngrid,reffrad)
+   subroutine h2so4_reffrad(ngrid,nlayer,reffrad)
 !==================================================================
 !     Purpose
@@ -322,13 +326,14 @@
       Implicit none
 
-      include "callkeys.h"
-      include "dimensions.h"
-      include "dimphys.h"
-
-      integer,intent(in) :: ngrid
-
-      real, intent(out) :: reffrad(ngrid,nlayermx)      !h2so4 particle radii (K)
+!      include "callkeys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+      integer,intent(in) :: nlayer
+
+      real, intent(out) :: reffrad(ngrid,nlayer)      !h2so4 particle radii (m)
                 
-      reffrad(1:ngrid,1:nlayermx) = 1.e-6 ! h2so4
+      reffrad(1:ngrid,1:nlayer) = 1.e-6 ! h2so4
 
    end subroutine h2so4_reffrad
@@ -352,10 +357,10 @@
 
      include "callkeys.h"
-     include "dimensions.h"
-     include "dimphys.h"
-
-      integer,intent(in) :: ngrid
-
-      real, intent(out) :: reffrad(ngrid,nlayermx)      ! particle radii 
+!     include "dimensions.h"
+!     include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+
+      real, intent(out) :: reffrad(ngrid,nlayer)      ! particle radii (m)
       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
       INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
Index: trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 1308)
@@ -1,3 +1,3 @@
-subroutine rain(ngrid,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
+subroutine rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
 
 
@@ -22,30 +22,31 @@
 !==================================================================
 
-  include "dimensions.h"
-  include "dimphys.h"
+!  include "dimensions.h"
+!  include "dimphys.h"
   include "comcstfi.h"
   include "callkeys.h"
 
 !     Arguments
-      integer,intent(in) :: ngrid ! number of atmospherci columns
+      integer,intent(in) :: ngrid ! number of atmospheric columns
+      integer,intent(in) :: nlayer ! number of atmospheric layers
       integer,intent(in) :: nq ! number of tracers
       real,intent(in) :: ptimestep    ! time interval
-      real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa)
-      real,intent(in) :: pplay(ngrid,nlayermx)   ! mid-layer pressure (Pa)
-      real,intent(in) :: t(ngrid,nlayermx) ! input temperature (K)
-      real,intent(in) :: pdt(ngrid,nlayermx) ! input tendency on temperature (K/s)      
-      real,intent(in) :: pq(ngrid,nlayermx,nq)  ! tracers (kg/kg)
-      real,intent(in) :: pdq(ngrid,nlayermx,nq) ! input tendency on tracers
-      real,intent(out) :: d_t(ngrid,nlayermx) ! temperature tendency (K/s)
-      real,intent(out) :: dqrain(ngrid,nlayermx,nq) ! tendency of H2O precipitation (kg/kg.s-1)
+      real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      real,intent(in) :: pplay(ngrid,nlayer)   ! mid-layer pressure (Pa)
+      real,intent(in) :: t(ngrid,nlayer) ! input temperature (K)
+      real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)      
+      real,intent(in) :: pq(ngrid,nlayer,nq)  ! tracers (kg/kg)
+      real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers
+      real,intent(out) :: d_t(ngrid,nlayer) ! temperature tendency (K/s)
+      real,intent(out) :: dqrain(ngrid,nlayer,nq) ! tendency of H2O precipitation (kg/kg.s-1)
       real,intent(out) :: dqsrain(ngrid)  ! rain flux at the surface (kg.m-2.s-1)
       real,intent(out) :: dqssnow(ngrid)  ! snow flux at the surface (kg.m-2.s-1)
-      real,intent(in) :: rneb(ngrid,nlayermx) ! cloud fraction
-
-      REAL zt(ngrid,nlayermx)         ! working temperature (K)
-      REAL ql(ngrid,nlayermx)         ! liquid water (Kg/Kg)
-      REAL q(ngrid,nlayermx)          ! specific humidity (Kg/Kg)
-      REAL d_q(ngrid,nlayermx)        ! water vapor increment
-      REAL d_ql(ngrid,nlayermx)       ! liquid water / ice increment
+      real,intent(in) :: rneb(ngrid,nlayer) ! cloud fraction
+
+      REAL zt(ngrid,nlayer)         ! working temperature (K)
+      REAL ql(ngrid,nlayer)         ! liquid water (Kg/Kg)
+      REAL q(ngrid,nlayer)          ! specific humidity (Kg/Kg)
+      REAL d_q(ngrid,nlayer)        ! water vapor increment
+      REAL d_ql(ngrid,nlayer)       ! liquid water / ice increment
 
 !     Subroutine options
@@ -73,5 +74,5 @@
 !     Local variables
       INTEGER i, k, n
-      REAL zqs(ngrid,nlayermx),Tsat(ngrid,nlayermx), zdelta, zcor
+      REAL zqs(ngrid,nlayer),Tsat(ngrid,nlayer), zdelta, zcor
       REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
 
@@ -80,10 +81,10 @@
       REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid)
 
-      real reffh2oliq(ngrid,nlayermx),reffh2oice(ngrid,nlayermx)
+      real reffh2oliq(ngrid,nlayer),reffh2oice(ngrid,nlayer)
   
       real ttemp, ptemp, psat_tmp
-      real tnext(ngrid,nlayermx)
-
-      real l2c(ngrid,nlayermx)
+      real tnext(ngrid,nlayer)
+
+      real l2c(ngrid,nlayer)
       real dWtot
 
@@ -151,5 +152,5 @@
 
 !     GCM -----> subroutine variables
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
       DO i = 1, ngrid
 
@@ -172,12 +173,12 @@
 
 !     Initialise the outputs
-      d_t(1:ngrid,1:nlayermx) = 0.0
-      d_q(1:ngrid,1:nlayermx) = 0.0
-      d_ql(1:ngrid,1:nlayermx) = 0.0
+      d_t(1:ngrid,1:nlayer) = 0.0
+      d_q(1:ngrid,1:nlayer) = 0.0
+      d_ql(1:ngrid,1:nlayer) = 0.0
       zrfl(1:ngrid) = 0.0
       zrfln(1:ngrid) = 0.0
 
       ! calculate saturation mixing ratio
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
          DO i = 1, ngrid
             ttemp = zt(i,k) 
@@ -190,5 +191,5 @@
 
       ! get column / layer conversion factor
-      DO k = 1, nlayermx
+      DO k = 1, nlayer
          DO i = 1, ngrid
             l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
@@ -199,5 +200,5 @@
       ! We carry the rain with us and calculate that added by warm/cold precipitation
       ! processes and that subtracted by evaporation at each level.
-      DO k = nlayermx, 1, -1
+      DO k = nlayer, 1, -1
 
          IF (evap_prec) THEN ! note no rneb dependence!
@@ -273,5 +274,5 @@
 
  !recalculate liquid water particle radii
-	   call h2o_cloudrad(ngrid,ql,reffh2oliq,reffh2oice)
+	   call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
 
            SELECT CASE(precip_scheme)
@@ -357,5 +358,5 @@
          endif ! if precip_scheme=1
 
-      ENDDO ! of DO k = nlayermx, 1, -1
+      ENDDO ! of DO k = nlayer, 1, -1
 
 !     Rain or snow on the ground
@@ -376,11 +377,11 @@
 !     now subroutine -----> GCM variables
       if (evap_prec) then
-        dqrain(1:ngrid,1:nlayermx,i_vap)=d_q(1:ngrid,1:nlayermx)/ptimestep
-        d_t(1:ngrid,1:nlayermx)=d_t(1:ngrid,1:nlayermx)/ptimestep
+        dqrain(1:ngrid,1:nlayer,i_vap)=d_q(1:ngrid,1:nlayer)/ptimestep
+        d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
       else
-        dqrain(1:ngrid,1:nlayermx,i_vap)=0.0
-        d_t(1:ngrid,1:nlayermx)=0.0
+        dqrain(1:ngrid,1:nlayer,i_vap)=0.0
+        d_t(1:ngrid,1:nlayer)=0.0
       endif
-      dqrain(1:ngrid,1:nlayermx,i_ice) = d_ql(1:ngrid,1:nlayermx)/ptimestep
+      dqrain(1:ngrid,1:nlayer,i_ice) = d_ql(1:ngrid,1:nlayer)/ptimestep
 
     end subroutine rain
Index: trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 1308)
@@ -16,4 +16,5 @@
       use comgeomphy, only: initcomgeomphy
       use slab_ice_h, only: noceanmx
+      use planete_mod
 
       implicit none
@@ -45,8 +46,8 @@
 #include "dimensions.h"
 #include "paramet.h"
-#include "dimphys.h"
+!include "dimphys.h"
 #include "callkeys.h"
 #include "comcstfi.h"
-#include "planete.h"
+!#include "planete.h"
 !#include "control.h"
 #include "comvert.h"
@@ -67,10 +68,10 @@
       REAL day           ! date durant le run
       REAL time             ! time (0<time<1 ; time=0.5 a midi)
-      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
-      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
+      REAL play(llm)   ! Pressure at the middle of the layers (Pa)
+      REAL plev(llm+1) ! intermediate pressure levels (pa)
       REAL psurf,tsurf(1)      
-      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
+      REAL u(llm),v(llm)  ! zonal, meridional wind
       REAL gru,grv   ! prescribed "geostrophic" background wind
-      REAL temp(nlayermx)   ! temperature at the middle of the layers
+      REAL temp(llm)   ! temperature at the middle of the layers
       REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
       REAL,ALLOCATABLE :: qsurf(:)      ! tracer surface budget (e.g. kg.m-2)
@@ -81,10 +82,10 @@
       integer :: i_h2o_vap=0  ! tracer index of h2o vapor
       REAL emis(1)             ! surface layer
-      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
-      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
+      REAL q2(llm+1)   ! Turbulent Kinetic Energy
+      REAL zlay(llm)   ! altitude estimee dans les couches (km)
 
 c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
-      REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
-      REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
+      REAL du(llm),dv(llm),dtemp(llm)
+      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
       REAL dpsurf    
       REAL,ALLOCATABLE :: dq(:,:)
@@ -94,8 +95,8 @@
 !      INTEGER thermo
       REAL zls
-      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
-      REAL pks, ptif, w(nlayermx)
+      REAL phi(llm),h(llm),s(llm)
+      REAL pks, ptif, w(llm)
       INTEGER ierr, aslun
-      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
+      REAL tmp1(0:llm),tmp2(0:llm)
       Logical  tracerdyn
       integer :: nq !=1 ! number of tracers
@@ -111,8 +112,8 @@
 !     added by RW for autozlevs computation
       real nu, xx, pMIN, zlev, Htop
-      real logplevs(nlayermx)
+      real logplevs(llm)
 
 !     added by BC
-      REAL cloudfrac(1,nlayermx)
+      REAL cloudfrac(1,llm)
       REAL hice(1),totcloudfrac(1)
       REAL qzero1D   !initial water amount on the ground
@@ -253,8 +254,8 @@
           if (nq>0) then
             allocate(tname(nq))
-            allocate(q(nlayermx,nq))
+            allocate(q(llm,nq))
             allocate(qsurf(nq))
-            allocate(dq(nlayermx,nq))
-            allocate(dqdyn(nlayermx,nq))
+            allocate(dq(llm,nq))
+            allocate(dqdyn(llm,nq))
           else
             write(*,*) "rcm1d: Error. You set tracer=.true."
@@ -298,6 +299,6 @@
         allocate(tname(1)) 
         allocate(qsurf(1))
-        allocate(q(nlayermx,1))
-        allocate(dq(nlayermx,1))
+        allocate(q(llm,1))
+        allocate(dq(llm,1))
       
        ! Check that tracer boolean is compliant with number of tracers
@@ -484,5 +485,5 @@
 c  --------------
 c
-      nlayer=nlayermx
+      nlayer=llm
       nlevel=nlayer+1
       nsoil=nsoilmx
@@ -890,5 +891,5 @@
          OPEN(12,file='profile.out',form='formatted')
          write(12,*) tsurf
-         DO ilayer=1,nlayermx
+         DO ilayer=1,llm
             write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
          ENDDO
Index: trunk/LMDZ.GENERIC/libf/phystd/soil.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1308)
@@ -17,6 +17,6 @@
 !-----------------------------------------------------------------------
 
-      include "dimensions.h"
-      include "dimphys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/start2archive.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/start2archive.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/start2archive.F	(revision 1308)
@@ -27,8 +27,10 @@
 ! to use  'getin'
       USE ioipsl_getincom
+      USE planete_mod
 
       implicit none
 
 #include "dimensions.h"
+      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
 #include "paramet.h"
 #include "comconst.h"
@@ -41,6 +43,6 @@
 #include "ener.h"
 
-#include "dimphys.h"
-#include "planete.h"
+!#include "dimphys.h"
+!#include "planete.h"
 !#include"advtrac.h"
 #include "netcdf.inc"
@@ -69,5 +71,5 @@
       REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
       REAL co2ice(ngridmx)	! CO2 ice layer
-      REAL q2(ngridmx,nlayermx+1)
+      REAL q2(ngridmx,llm+1)
       REAL,ALLOCATABLE :: qsurf(:,:)
       REAL emis(ngridmx)
@@ -80,5 +82,5 @@
 !     added by FF for cloud fraction setup
       REAL hice(ngridmx)
-      REAL cloudfrac(ngridmx,nlayermx),totalcloudfrac(ngridmx)
+      REAL cloudfrac(ngridmx,llm),totalcloudfrac(ngridmx)
 
 !     added by BC for slab ocean
@@ -199,6 +201,6 @@
 
 
-      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqtot,day_ini_fi,
-     .      timefi,
+      CALL phyetat0 (ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot,
+     .      day_ini_fi,timefi,
      .      tsurf,tsoil,emis,q2,qsurf,
 !       change FF 05/2011
Index: trunk/LMDZ.GENERIC/libf/phystd/stellarlong.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/stellarlong.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/stellarlong.F	(revision 1308)
@@ -1,3 +1,5 @@
       SUBROUTINE stellarlong(pday,pstellong)
+      
+      USE planete_mod, ONLY: year_day, peri_day, e_elips, timeperi
       IMPLICIT NONE
 
@@ -40,5 +42,5 @@
 c   -------------
 
-#include "planete.h"
+!#include "planete.h"
 #include "comcstfi.h"
 
Index: trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90	(revision 1308)
@@ -7,5 +7,5 @@
 CONTAINS
 
-      SUBROUTINE divgrad_phy(nlevs,temp,delta)
+      SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
 
 
@@ -14,5 +14,5 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "paramet.h"
 #include "comcstfi.h"
@@ -20,12 +20,13 @@
 #include "comhdiff.h"
       
-      INTEGER nlevs, ll
-      REAL temp(ngridmx,nlevs)
+      INTEGER,INTENT(IN) :: ngrid, nlevs
+      REAL,INTENT(IN) :: temp(ngrid,nlevs)
+      REAL,INTENT(OUT) :: delta(ngrid,nlevs)
       REAL delta_2d(ip1jmp1,nlevs)
-      REAL delta(ngridmx,nlevs)
+      INTEGER :: ll
       REAL ghx(ip1jmp1,nlevs), ghy(ip1jm,nlevs)
 
 
-      CALL gr_fi_dyn(nlevs,ngridmx,iip1,jjp1,temp,delta_2d)
+      CALL gr_fi_dyn(nlevs,ngrid,iip1,jjp1,temp,delta_2d)
       CALL grad(nlevs,delta_2d,ghx,ghy)
       DO ll=1,nlevs
@@ -36,5 +37,5 @@
       END DO
       CALL diverg(nlevs,ghx,ghy,delta_2d)
-      CALL gr_dyn_fi(nlevs,iip1,jjp1,ngridmx,delta_2d,delta)
+      CALL gr_dyn_fi(nlevs,iip1,jjp1,ngrid,delta_2d,delta)
 
 
@@ -43,10 +44,10 @@
 
 
-      SUBROUTINE init_masquv(zmasq)
+      SUBROUTINE init_masquv(ngrid,zmasq)
      
       IMPLICIT NONE
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "paramet.h"
 #include "comcstfi.h"   
@@ -55,5 +56,6 @@
 
 
-      REAL zmasq(ngridmx)
+      INTEGER,INTENT(IN) :: ngrid
+      REAL zmasq(ngrid)
       REAL zmasq_2d(ip1jmp1)
       REAL ff(ip1jm)
@@ -66,5 +68,5 @@
       zmasqv=1.
 
-      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zmasq,zmasq_2d)
+      CALL gr_fi_dyn(1,ngrid,iip1,jjp1,zmasq,zmasq_2d)
 
       DO i=1,ip1jmp1-1
@@ -104,5 +106,5 @@
 
 
-      SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy)
+      SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
   
           use slab_ice_h 
@@ -111,5 +113,5 @@
       
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "paramet.h"
 #include "comcstfi.h"
@@ -118,4 +120,5 @@
 #include "comhdiff.h"
 
+      INTEGER,INTENT(IN) :: ngrid
       INTEGER ij
       REAL txv(ip1jm),fluxm(ip1jm),tyv(ip1jm)
@@ -123,6 +126,6 @@
       REAL tyu(ip1jmp1),txu(ip1jmp1),fluxz(ip1jmp1),fluxv(ip1jmp1)
       REAL dt(ip1jmp1,noceanmx),ts(ip1jmp1,noceanmx)
-      REAL tx_phy(ngridmx),ty_phy(ngridmx)
-      REAL dt_phy(ngridmx,noceanmx),ts_phy(ngridmx,noceanmx)
+      REAL tx_phy(ngrid),ty_phy(ngrid)
+      REAL dt_phy(ngrid,noceanmx),ts_phy(ngrid,noceanmx)
 
 
@@ -130,6 +133,6 @@
 
 ! Passage taux,y sur grilles 2d
-      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,tx_phy,txu)
-      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,ty_phy,tyu)
+      CALL gr_fi_dyn(1,ngrid,iip1,jjp1,tx_phy,txu)
+      CALL gr_fi_dyn(1,ngrid,iip1,jjp1,ty_phy,tyu)
 ! Passage grille u,v
 ! Multiplication par f ou eps.
@@ -144,5 +147,5 @@
             
 ! Calcul de T, Tdeep
-      CALL gr_fi_dyn(2,ngridmx,iip1,jjp1,ts_phy,ts)
+      CALL gr_fi_dyn(2,ngrid,iip1,jjp1,ts_phy,ts)
        
 ! Flux de masse
@@ -222,5 +225,5 @@
 
 ! Retour grille physique
-      CALL gr_dyn_fi(2,iip1,jjp1,ngridmx,dt,dt_phy)
+      CALL gr_dyn_fi(2,iip1,jjp1,ngrid,dt,dt_phy)
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/surface_nature.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/surface_nature.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/surface_nature.F	(revision 1308)
@@ -33,6 +33,6 @@
 !==================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/surfini.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/surfini.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/surfini.F	(revision 1308)
@@ -15,6 +15,6 @@
 c   Declarations:
 c   -------------
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "callkeys.h"
 c
Index: trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 1308)
@@ -34,5 +34,5 @@
 c                              comparer avec le day_ini dynamique)
 c
-c      - lmax:    tab_cntrl(tab0+2) (pour test avec nlayermx)
+c      - lmax:    tab_cntrl(tab0+2) (pour test avec nlayer)
 c
 c      - p_rad
@@ -51,9 +51,11 @@
       use iostart, only: get_var
       use mod_phys_lmdz_para, only: is_parallel
+      use planete_mod, only: year_day, periastr, apoastr, peri_day,
+     &                       obliquit, z0, lmixmin, emin_turb
 
       implicit none
  
 #include "comcstfi.h"
-#include "planete.h"
+!#include "planete.h"
 #include "netcdf.inc"
 #include "callkeys.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90	(revision 1308)
@@ -1,3 +1,3 @@
-      subroutine totalcloudfrac(ngrid,nq,rneb,totalrneb,pplev,pq,tau)
+      subroutine totalcloudfrac(ngrid,nlayer,nq,rneb,totalrneb,pplev,pq,tau)
 
       use watercommon_h
@@ -19,18 +19,19 @@
 !==================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
 
       integer,intent(in) :: ngrid        ! number of atmospheric columns
+      integer,intent(in) :: nlayer       ! number of atmospheric layers
       integer,intent(in) :: nq           ! number of tracers
-      real,intent(in) :: rneb(ngrid,nlayermx)    ! cloud fraction     
+      real,intent(in) :: rneb(ngrid,nlayer)    ! cloud fraction     
       real,intent(out) :: totalrneb(ngrid)       ! total cloud fraction 
-      real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa)
-      real,intent(in) :: pq(ngrid,nlayermx,nq)   ! tracers (.../kg_of_air)
-      real,intent(in) :: tau(ngrid,nlayermx)
+      real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      real,intent(in) :: pq(ngrid,nlayer,nq)   ! tracers (.../kg_of_air)
+      real,intent(in) :: tau(ngrid,nlayer)
 
-      real, dimension(nlayermx+1) :: masse
+      real, dimension(nlayer+1) :: masse
       integer, parameter          :: recovery=7
       integer ltau_max
@@ -47,5 +48,5 @@
       real clear,tau_min
       real, parameter ::  tau_c=0.1 !threshold of optical depth for the calculation of total cloud fraction 
-      real rneb2(nlayermx)
+      real rneb2(nlayer)
 
 
@@ -55,5 +56,5 @@
          if (recovery.eq.1) then
             clear = (1.-rneb(ig,1))
-            do l=2,nlayermx       
+            do l=2,nlayer      
                clear = clear*(1.-rneb(ig,l))
             enddo
@@ -62,5 +63,5 @@
          elseif (recovery.eq.2) then
             totalrneb(ig) = rneb(ig,1)
-            do l=2,14 !nlayermx       
+            do l=2,14 !nlayer    
                totalrneb(ig) = max(rneb(ig,l),totalrneb(ig))
             enddo
@@ -68,5 +69,5 @@
          elseif (recovery.eq.3) then
             totalrneb(ig) = rneb(ig,1)
-            do l=2,nlayermx       
+            do l=2,nlayer    
                totalrneb(ig) = min(rneb(ig,l),totalrneb(ig))
             enddo
@@ -77,5 +78,5 @@
          elseif (recovery.eq.5) then
             totalrneb(ig) = rneb(ig,1)            
-            do l=1,nlayermx
+            do l=1,nlayer
                masse(l)=pq(ig,l,igcm_h2o_ice)*(pplev(ig,l)-pplev(ig,l+1))
             enddo
@@ -85,10 +86,10 @@
          elseif (recovery.eq.6) then
             totalrneb(ig) = 0.            
-            do l=1,nlayermx
+            do l=1,nlayer
                masse(l)=pq(ig,l,igcm_h2o_ice)*(pplev(ig,l)-pplev(ig,l+1))
                masse(l)=max(masse(l),0.)
             enddo
             massetot=sum(masse,dim=1)
-            do l=1,nlayermx
+            do l=1,nlayer
                totalrneb(ig) = totalrneb(ig)+rneb(ig,l)*masse(l)/massetot
             enddo
@@ -96,10 +97,10 @@
          elseif (recovery.eq.7) then
 
-            rneb2(:)=rneb(ig,1:nlayermx)
-	    tau_min=MIN(tau_c,MAXVAL(tau(ig,1:nlayermx))/2.)
-            do l=1,nlayermx
+            rneb2(:)=rneb(ig,1:nlayer)
+	    tau_min=MIN(tau_c,MAXVAL(tau(ig,1:nlayer))/2.)
+            do l=1,nlayer
                if(tau(ig,l)<tau_min) rneb2(l)=0.      	
             enddo
-            totalrneb(ig)=maxval(rneb2(1:nlayermx))
+            totalrneb(ig)=maxval(rneb2(1:nlayer))
 
          endif                  ! (recovery=)   
Index: trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1308)
@@ -41,6 +41,6 @@
 !     ------------
 
-      include "dimensions.h"
-      include "dimphys.h"
+!      include "dimensions.h"
+!      include "dimphys.h"
       include "comcstfi.h"
       include "callkeys.h"
@@ -48,5 +48,6 @@
 !     arguments
 !     ---------
-      INTEGER,INTENT(IN) :: ngrid,nlay
+      INTEGER,INTENT(IN) :: ngrid
+      INTEGER,INTENT(IN) :: nlay
       REAL,INTENT(IN) :: ptimestep
       REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
@@ -86,18 +87,18 @@
 
       REAL z4st,zdplanck(ngrid)
-      REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
+      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
       REAL zcdv(ngrid),zcdh(ngrid)
       REAL zcdv_true(ngrid),zcdh_true(ngrid)
-      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
-      REAL zh(ngrid,nlayermx),zt(ngrid,nlayermx)
+      REAL zu(ngrid,nlay),zv(ngrid,nlay)
+      REAL zh(ngrid,nlay),zt(ngrid,nlay)
       REAL ztsrf(ngrid)
       REAL z1(ngrid),z2(ngrid)
-      REAL zmass(ngrid,nlayermx)
-      REAL zfluxv(ngrid,nlayermx),zfluxt(ngrid,nlayermx),zfluxq(ngrid,nlayermx)
-      REAL zb0(ngrid,nlayermx)
-      REAL zExner(ngrid,nlayermx),zovExner(ngrid,nlayermx)
-      REAL zcv(ngrid,nlayermx),zdv(ngrid,nlayermx)  !inversion coefficient for winds
-      REAL zct(ngrid,nlayermx),zdt(ngrid,nlayermx)  !inversion coefficient for temperature
-      REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx)  !inversion coefficient for tracers
+      REAL zmass(ngrid,nlay)
+      REAL zfluxv(ngrid,nlay),zfluxt(ngrid,nlay),zfluxq(ngrid,nlay)
+      REAL zb0(ngrid,nlay)
+      REAL zExner(ngrid,nlay),zovExner(ngrid,nlay)
+      REAL zcv(ngrid,nlay),zdv(ngrid,nlay)  !inversion coefficient for winds
+      REAL zct(ngrid,nlay),zdt(ngrid,nlay)  !inversion coefficient for temperature
+      REAL zcq(ngrid,nlay),zdq(ngrid,nlay)  !inversion coefficient for tracers
       REAL zcst1
       REAL zu2!, a
@@ -110,7 +111,7 @@
 !     -------
       INTEGER iq
-      REAL zq(ngrid,nlayermx,nq)
-      REAL zqnoevap(ngrid,nlayermx) !special for water case to compute where evaporated water goes.
-      REAL pdqevap(ngrid,nlayermx) !special for water case to compute where evaporated water goes.
+      REAL zq(ngrid,nlay,nq)
+      REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
+      REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
       REAL zdmassevap(ngrid)
       REAL rho(ngrid)         ! near-surface air density
@@ -241,5 +242,5 @@
 !     ------------------------------------------------------ 
 
-      call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature 
+      call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature 
       
 !     Adding eddy mixing to mimic 3D general circulation in 1D
@@ -262,5 +263,5 @@
  
 !JL12 calculate the flux coefficients (tables multiplied elements by elements)
-      zfluxv(1:ngrid,1:nlayermx)=zkv(1:ngrid,1:nlayermx)*zb0(1:ngrid,1:nlayermx)
+      zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
       
 !-----------------------------------------------------------------------
@@ -364,6 +365,6 @@
 !     JL12 calculate the flux coefficients (tables multiplied elements by elements)
 !     ---------------------------------------------------------------
-      zfluxq(1:ngrid,1:nlayermx)=zkh(1:ngrid,1:nlayermx)*zb0(1:ngrid,1:nlayermx) !JL12 we save zfluxq which doesn't need the Exner factor
-      zfluxt(1:ngrid,1:nlayermx)=zfluxq(1:ngrid,1:nlayermx)*zExner(1:ngrid,1:nlayermx)
+      zfluxq(1:ngrid,1:nlay)=zkh(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay) !JL12 we save zfluxq which doesn't need the Exner factor
+      zfluxt(1:ngrid,1:nlay)=zfluxq(1:ngrid,1:nlay)*zExner(1:ngrid,1:nlay)
 
       DO ig=1,ngrid
Index: trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F	(revision 1308)
@@ -1,7 +1,7 @@
-      SUBROUTINE vdif_kc(ngrid,dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
+      SUBROUTINE vdif_kc(ngrid,nlay,dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
       IMPLICIT NONE
 c.......................................................................
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 c.......................................................................
 c
@@ -28,14 +28,15 @@
 c.......................................................................
       INTEGER,INTENT(IN) :: ngrid
+      INTEGER,INTENT(IN) :: nlay
       REAL,INTENT(IN) :: dt,g
-      REAL,INTENT(IN) :: zlev(ngrid,nlayermx+1)
-      REAL,INTENT(IN) :: zlay(ngrid,nlayermx)
-      REAL,INTENT(IN) :: u(ngrid,nlayermx)
-      REAL,INTENT(IN) :: v(ngrid,nlayermx)
-      REAL,INTENT(IN) :: teta(ngrid,nlayermx)
+      REAL,INTENT(IN) :: zlev(ngrid,nlay+1)
+      REAL,INTENT(IN) :: zlay(ngrid,nlay)
+      REAL,INTENT(IN) :: u(ngrid,nlay)
+      REAL,INTENT(IN) :: v(ngrid,nlay)
+      REAL,INTENT(IN) :: teta(ngrid,nlay)
       REAL,INTENT(IN) :: cd(ngrid)
-      REAL,INTENT(INOUT) :: q2(ngrid,nlayermx+1)
-      REAL,INTENT(OUT) :: km(ngrid,nlayermx+1)
-      REAL,INTENT(OUT) :: kn(ngrid,nlayermx+1)
+      REAL,INTENT(INOUT) :: q2(ngrid,nlay+1)
+      REAL,INTENT(OUT) :: km(ngrid,nlay+1)
+      REAL,INTENT(OUT) :: kn(ngrid,nlay+1)
 c.......................................................................
 c
@@ -50,9 +51,8 @@
 c
 c.......................................................................
-      INTEGER,PARAMETER :: nlay=nlayermx
-      INTEGER,PARAMETER :: nlev=nlayermx+1
-      REAL unsdz(ngrid,nlayermx)
-      REAL unsdzdec(ngrid,nlayermx+1)
-      REAL q(ngrid,nlayermx+1)
+      INTEGER :: nlev
+      REAL unsdz(ngrid,nlay)
+      REAL unsdzdec(ngrid,nlay+1)
+      REAL q(ngrid,nlay+1)
 c.......................................................................
 c
@@ -64,5 +64,5 @@
 c
 c.......................................................................
-      REAL kmpre(ngrid,nlayermx+1)
+      REAL kmpre(ngrid,nlay+1)
       REAL qcstat
       REAL q2cstat
@@ -72,5 +72,5 @@
 c
 c.......................................................................
-      REAL long(ngrid,nlayermx+1)
+      REAL long(ngrid,nlay+1)
 c.......................................................................
 c
@@ -95,8 +95,8 @@
       REAL mcstat
       REAL m2cstat
-      REAL m(ngrid,nlayermx+1)
-      REAL mpre(ngrid,nlayermx+1)
-      REAL m2(ngrid,nlayermx+1)
-      REAL n2(ngrid,nlayermx+1)
+      REAL m(ngrid,nlay+1)
+      REAL mpre(ngrid,nlay+1)
+      REAL m2(ngrid,nlay+1)
+      REAL n2(ngrid,nlay+1)
 c.......................................................................
 c
@@ -120,9 +120,9 @@
       LOGICAL gnsup
       REAL gm
-c      REAL ri(ngrid,nlayermx+1)
-      REAL sn(ngrid,nlayermx+1)
-      REAL snq2(ngrid,nlayermx+1)
-      REAL sm(ngrid,nlayermx+1)
-      REAL smq2(ngrid,nlayermx+1)
+c      REAL ri(ngrid,nlaye+1)
+      REAL sn(ngrid,nlay+1)
+      REAL snq2(ngrid,nlay+1)
+      REAL sm(ngrid,nlay+1)
+      REAL smq2(ngrid,nlay+1)
 c.......................................................................
 c
@@ -179,4 +179,5 @@
 
 ! initialization of local variables:
+      nlev=nlay+1
       long(:,:)=0.
       n2(:,:)=0.
Index: trunk/LMDZ.GENERIC/libf/phystd/vdifc.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/vdifc.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/vdifc.F	(revision 1308)
@@ -38,6 +38,6 @@
 !     ------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -77,17 +77,17 @@
 
       REAL z4st,zdplanck(ngrid)
-      REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
+      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
       REAL zcdv(ngrid),zcdh(ngrid)
       REAL zcdv_true(ngrid),zcdh_true(ngrid)
-      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
-      REAL zh(ngrid,nlayermx)
+      REAL zu(ngrid,nlay),zv(ngrid,nlay)
+      REAL zh(ngrid,nlay)
       REAL ztsrf2(ngrid)
       REAL z1(ngrid),z2(ngrid)
-      REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
-      REAL zb0(ngrid,nlayermx)
-      REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
+      REAL za(ngrid,nlay),zb(ngrid,nlay)
+      REAL zb0(ngrid,nlay)
+      REAL zc(ngrid,nlay),zd(ngrid,nlay)
       REAL zcst1
       REAL zu2!, a
-      REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx)
+      REAL zcq(ngrid,nlay),zdq(ngrid,nlay)
       REAL evap(ngrid)
       REAL zcq0(ngrid),zdq0(ngrid)
@@ -101,5 +101,5 @@
 !     variables added for CO2 condensation
 !     ------------------------------------
-      REAL hh                   !, zhcond(ngrid,nlayermx)
+      REAL hh                   !, zhcond(ngrid,nlay)
 !     REAL latcond,tcond1mb
 !     REAL acond,bcond
@@ -110,5 +110,5 @@
 !     -------
       INTEGER iq
-      REAL zq(ngrid,nlayermx,nq)
+      REAL zq(ngrid,nlay,nq)
       REAL zq1temp(ngrid)
       REAL rho(ngrid)         ! near-surface air density
@@ -123,6 +123,6 @@
       real z1_T(ngrid), z2_T(ngrid)
       real zb_T(ngrid)
-      real zc_T(ngrid,nlayermx)
-      real zd_T(ngrid,nlayermx)
+      real zc_T(ngrid,nlay)
+      real zd_T(ngrid,nlay)
       real lat1(ngrid), lat2(ngrid)
       real surfh2otot
@@ -249,5 +249,5 @@
 !     ------------------------------------------------------ 
 
-      call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay
+      call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay
      &     ,pu,pv,ph,zcdv_true
      &     ,pq2,zkv,zkh)
Index: trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F	(revision 1308)
@@ -1,3 +1,3 @@
-      SUBROUTINE vlz_fi(ngrid,q,pente_max,masse,w,wq)
+      SUBROUTINE vlz_fi(ngrid,nlayer,q,pente_max,masse,w,wq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -16,6 +16,6 @@
       IMPLICIT NONE
 c
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 
 c
@@ -23,9 +23,9 @@
 c   Arguments:
 c   ----------
-      integer ngrid
-      real masse(ngrid,llm),pente_max
-      REAL q(ngrid,llm)
-      REAL w(ngrid,llm)
-      REAL wq(ngrid,llm+1)
+      integer,intent(in) :: ngrid, nlayer
+      real,intent(in) :: masse(ngrid,nlayer),pente_max
+      REAL,INTENT(INOUT) :: q(ngrid,nlayer)
+      REAL,INTENT(INOUT) :: w(ngrid,nlayer)
+      REAL,INTENT(OUT) :: wq(ngrid,nlayer+1)
 c
 c      Local 
@@ -35,5 +35,6 @@
 c
 
-      real dzq(ngrid,llm),dzqw(ngrid,llm),adzqw(ngrid,llm),dzqmax
+      real dzq(ngrid,nlayer),dzqw(ngrid,nlayer),adzqw(ngrid,nlayer)
+      real dzqmax
       real newmasse
       real sigw, Mtot, MQtot
@@ -47,5 +48,5 @@
 c    sens de W
 
-      do l=2,llm
+      do l=2,nlayer
          do ij=1,ngrid
             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
@@ -54,5 +55,5 @@
       enddo
 
-      do l=2,llm-1
+      do l=2,nlayer-1
          do ij=1,ngrid
 #ifdef CRAY
@@ -73,5 +74,5 @@
       do ij=1,ngrid
          dzq(ij,1)=0.
-         dzq(ij,llm)=0.
+         dzq(ij,nlayer)=0.
       enddo
 c ---------------------------------------------------------------
@@ -83,5 +84,5 @@
 c      No flux at the model top:
        do ij=1,ngrid
-          wq(ij,llm+1)=0.
+          wq(ij,nlayer+1)=0.
        enddo
 
@@ -89,5 +90,5 @@
 c      ===============================
 
-       do l = 1,llm          ! loop different than when w<0
+       do l = 1,nlayer          ! loop different than when w<0
         do ij=1,ngrid
 
@@ -107,13 +108,13 @@
             Mtot = masse(ij,m)
             MQtot = masse(ij,m)*q(ij,m)
-            if(m.ge.llm)goto 88
+            if(m.ge.nlayer)goto 88
             do while(w(ij,l).gt.(Mtot+masse(ij,m+1)))
                 m=m+1
                 Mtot = Mtot + masse(ij,m)
                 MQtot = MQtot + masse(ij,m)*q(ij,m)
-                if(m.ge.llm)goto 88
+                if(m.ge.nlayer)goto 88
             end do
  88         continue
-            if (m.lt.llm) then
+            if (m.lt.nlayer) then
                 sigw=(w(ij,l)-Mtot)/masse(ij,m+1)
                 wq(ij,l)=(MQtot + (w(ij,l)-Mtot)*
@@ -137,5 +138,5 @@
        end do
 
-       do l = 1,llm-1  ! loop different than when w>0
+       do l = 1,nlayer-1  ! loop different than when w>0
         do ij=1,ngrid
          if(w(ij,l+1).le.0)then
@@ -176,5 +177,5 @@
  99    continue
 
-      do l=1,llm
+      do l=1,nlayer
          do ij=1,ngrid
 
@@ -191,5 +192,3 @@
 
 
-
-      return
       end
Index: trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 1308)
@@ -54,5 +54,5 @@
 ! Commons
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "paramet.h"
 !#include "control.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 1308)
@@ -54,5 +54,5 @@
 ! Commons
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "paramet.h"
 !#include "control.h"
Index: trunk/LMDZ.GENERIC/libf/phystd/wstats.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/wstats.F90	(revision 1304)
+++ trunk/LMDZ.GENERIC/libf/phystd/wstats.F90	(revision 1308)
@@ -7,5 +7,5 @@
 
 #include "dimensions.h"
-#include "dimphys.h"
+!#include "dimphys.h"
 #include "comconst.h"
 #include "statto.h"
@@ -293,5 +293,5 @@
 
 include "dimensions.h"
-include "dimphys.h"
+!include "dimphys.h"
 include "netcdf.inc"
 
