Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 1527)
+++ /trunk/LMDZ.MARS/README	(revision 1528)
@@ -2237,2 +2237,8 @@
   of relevent variables for outputs in time_phylmdz_mod.
 
+== 02/04/2016 == EM
+- Got rid of references to "dimensions.h" from physics packages:
+  use nbp_lon (=iim), nbp_lat (==jjp1) and nbp_lev 
+  from module mod_grid_phy_lmdz (in phy_common) instead.
+- Added "ioipsl_getin_p_mod.F90" (getin_p routine) in phy_common to
+  correctly read in parameters from *.def files in a parallel environment.
Index: /trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90	(revision 1528)
@@ -4,4 +4,5 @@
       use tracer_mod
       USE comvert_mod, ONLY: aps,bps
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       implicit none
 
@@ -26,7 +27,7 @@
 !  ----------
 !
-!  pq(iip1,jjp1,llm,nq)  Advected fields, ie chemical species here
+!  pq(nbp_lon+1,nbp_lat,nbp_lev,nq)  Advected fields, ie chemical species here
 !  qsurf(ngrid,nq)     Amount of tracer on the surface (kg/m2)
-!  ps(iip1,jjp1)           Surface pressure (Pa)
+!  ps(nbp_lon+1,nbp_lat)           Surface pressure (Pa)
 !  flagh2o                 flag for initialisation of h2o (1: yes / 0: no)
 !  flagthermo              flag for initialisation of thermosphere only (1: yes / 0: no)
@@ -34,8 +35,6 @@
 !=======================================================================
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "callkeys.h"
-#include "datafile.h"
+      include "callkeys.h"
+      include "datafile.h"
 
 ! inputs :
@@ -43,5 +42,5 @@
       integer,intent(in) :: ngrid         ! number of atmospheric columns in the physics
       integer,intent(in) :: nq                    ! number of tracers
-      real,intent(in) :: ps(iip1,jjp1)            ! surface pressure in the gcm (Pa)   
+      real,intent(in) :: ps(nbp_lon+1,nbp_lat)            ! surface pressure in the gcm (Pa)   
       integer,intent(in) :: flagh2o               ! flag for h2o initialisation
       integer,intent(in) :: flagthermo            ! flag for thermosphere initialisation only
@@ -49,5 +48,5 @@
 ! outputs :
 
-      real,intent(out) :: pq(iip1,jjp1,llm,nq)  ! advected fields, ie chemical species
+      real,intent(out) :: pq(nbp_lon+1,nbp_lat,nbp_lev,nq)  ! advected fields, ie chemical species
       real,intent(out) :: qsurf(ngrid,nq)     ! surface values (kg/m2) of tracers
 
@@ -56,5 +55,5 @@
       integer :: iq, i, j, l, n, nbqchem
       integer :: count, ierr, dummy
-      real    :: mmean(iip1,jjp1,llm)             ! mean molecular mass (g)
+      real    :: mmean(nbp_lon+1,nbp_lat,nbp_lev)             ! mean molecular mass (g)
       real    :: pgcm                             ! pressure at each layer in the gcm (Pa)
 
@@ -547,7 +546,7 @@
 ! 3. initialization of tracers
 
-      do i = 1,iip1
-         do j = 1,jjp1
-            do l = 1,llm
+      do i = 1,nbp_lon+1
+         do j = 1,nbp_lat
+            do l = 1,nbp_lev
 
                pgcm = aps(l) + bps(l)*ps(i,j)  ! gcm pressure
@@ -596,7 +595,7 @@
       if (igcm_ch4 /= 0) then
          vmr = 10.e-9       
-         do i = 1,iip1
-            do j = 1,jjp1
-               do l = 1,llm
+         do i = 1,nbp_lon+1
+            do j = 1,nbp_lat
+               do l = 1,nbp_lev
                   pq(i,j,l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(i,j,l)
                end do
@@ -623,7 +622,7 @@
          end if
 
-         do i = 1,iip1
-            do j = 1,jjp1
-               do l = 1,llm
+         do i = 1,nbp_lon+1
+            do j = 1,nbp_lat
+               do l = 1,nbp_lev
                   ! all ions to 0     
                   pq(i,j,l,igcm_co2plus)  = 0.
Index: /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1528)
@@ -7,13 +7,5 @@
 implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
-!#include "comcstfi.h"
-!#include "callkeys.h"
-!#include "comdiurn.h"
-!#include "chimiedata.h"
-!#include "tracer.h"
-!#include "conc.h"
-#include "diffusion.h"
+include "diffusion.h"
 
 ! July 2014 JYC ADD BALISTIC Transport coupling to compute wup for H and H2
@@ -937,5 +929,4 @@
         use tracer_mod, only: nqmx
         IMPLICIT NONE
-!#include "dimensions.h"
 
         INTEGER,INTENT(IN) :: nl,nq
@@ -973,5 +964,4 @@
 	use tracer_mod, only: nqmx
         IMPLICIT NONE
-!#include "dimensions.h"
 
 	INTEGER :: nl,nq,l
@@ -1058,5 +1048,4 @@
         use tracer_mod, only: nqmx
 	IMPLICIT NONE
-!#include "dimensions.h"
 	
 	INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
@@ -1400,5 +1389,4 @@
 	use tracer_mod, only: nqmx
         IMPLICIT NONE
-!#include "dimensions.h"
 	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
 	INTEGER,DIMENSION(1) :: indP
@@ -1502,5 +1490,4 @@
         use tracer_mod, only: nqmx
         IMPLICIT NONE
-!#include "dimensions.h"
         INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
         INTEGER,DIMENSION(1) :: indP
Index: /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 1528)
@@ -19,7 +19,6 @@
 !==========================================================================
 
-#include "dimensions.h"
-#include "callkeys.h"
-#include "chimiedata.h"
+      include "callkeys.h"
+      include "chimiedata.h"
 
 ! input
Index: /trunk/LMDZ.MARS/libf/phy_common/ioipsl_getin_p_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phy_common/ioipsl_getin_p_mod.F90	(revision 1528)
+++ /trunk/LMDZ.MARS/libf/phy_common/ioipsl_getin_p_mod.F90	(revision 1528)
@@ -0,0 +1,169 @@
+!
+! $Id$
+!
+MODULE ioipsl_getin_p_mod
+! To use getin in a parallel context
+!---------------------------------------------------------------------
+#ifdef CPP_IOIPSL
+USE ioipsl, ONLY: getin
+#else
+USE ioipsl_getincom, ONLY: getin
+#endif
+USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+use mod_phys_lmdz_para, only : bcast
+!-
+IMPLICIT NONE
+!-
+PRIVATE
+PUBLIC :: getin_p
+!-
+INTERFACE getin_p
+
+  MODULE PROCEDURE getinrs_p, getinr1d_p, getinr2d_p, &
+ &                 getinis_p, getini1d_p, getini2d_p, &
+ &                 getincs_p, 		              &
+ &                 getinls_p, getinl1d_p, getinl2d_p
+END INTERFACE
+!-
+CONTAINS
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!   Definition des getin -> bcast      !!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!! -- Les chaines de caracteres -- !!
+  
+  SUBROUTINE getincs_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    CHARACTER(LEN=*),INTENT(INOUT) :: VarOut    
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getincs_p
+
+!! -- Les entiers -- !!
+  
+  SUBROUTINE getinis_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    INTEGER,INTENT(INOUT) :: VarOut    
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinis_p
+
+  SUBROUTINE getini1d_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    INTEGER,INTENT(INOUT) :: VarOut(:)
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getini1d_p
+
+  SUBROUTINE getini2d_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    INTEGER,INTENT(INOUT) :: VarOut(:,:)
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getini2d_p
+
+!! -- Les flottants -- !!
+  
+  SUBROUTINE getinrs_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    REAL,INTENT(INOUT) :: VarOut
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinrs_p
+
+  SUBROUTINE getinr1d_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    REAL,INTENT(INOUT) :: VarOut(:)
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinr1d_p
+
+  SUBROUTINE getinr2d_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    REAL,INTENT(INOUT) :: VarOut(:,:)
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinr2d_p
+
+!! -- Les Booleens -- !!
+  
+  SUBROUTINE getinls_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    LOGICAL,INTENT(INOUT) :: VarOut
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinls_p
+
+  SUBROUTINE getinl1d_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    LOGICAL,INTENT(INOUT) :: VarOut(:)
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinl1d_p
+
+  SUBROUTINE getinl2d_p(VarIn,VarOut)
+    IMPLICIT NONE    
+    CHARACTER(LEN=*),INTENT(IN) :: VarIn
+    LOGICAL,INTENT(INOUT) :: VarOut(:,:)
+
+!$OMP BARRIER
+    IF (is_mpi_root .AND. is_omp_root) THEN
+    	CALL getin(VarIn,VarOut)
+    ENDIF
+    CALL bcast(VarOut)
+  END SUBROUTINE getinl2d_p
+!-
+!-----------------------------
+!-----------------------------
+!-----------------------------
+
+END MODULE ioipsl_getin_p_mod
+
Index: /trunk/LMDZ.MARS/libf/phymars/albedocaps.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/albedocaps.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/albedocaps.F90	(revision 1528)
@@ -12,5 +12,5 @@
 implicit none
 
-#include"callkeys.h"
+include"callkeys.h"
 
 ! arguments:
@@ -93,6 +93,5 @@
                   
 implicit none
-#include"dimensions.h"
-#include"datafile.h"
+include"datafile.h"
 
 ! arguments:
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/abort_gcm.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/abort_gcm.F	(revision 1528)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/abort_gcm.F	(revision 1528)
@@ -0,0 +1,1 @@
+link ../../dyn3d/abort_gcm.F
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/disvert.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/disvert.F	(revision 1528)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/disvert.F	(revision 1528)
@@ -0,0 +1,1 @@
+link ../../dyn3d/disvert.F
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1528)
@@ -827,4 +827,2 @@
 c***********************************************************************
 
-#include "../dyn3d/disvert.F"
-#include "../dyn3d/abort_gcm.F"
Index: /trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90	(revision 1528)
@@ -1,4 +1,6 @@
 module eofdump_mod
 ! this module controls the production of data for EOFs
+! it won't work if run in parallel (but it's OK, we don't use it anymore...)
+! Mainly kept for reference.
 implicit none
 ! Dump profiles for EOFs every ieofs physics timesteps,
@@ -15,4 +17,5 @@
       subroutine eofdump(ngrid,nlayer,u,v,t,rho,ps)
 
+      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
       implicit none
 !
@@ -21,5 +24,4 @@
 !     Corrected small bug in sampling rate/count, EM 11/2007
 !
-#include "dimensions.h"
 !
 
@@ -52,7 +54,7 @@
       if (mod(count+1,ieofs).eq.0) then
 !        write(*,*)'eofdump: dump --> ps(1)=',ps(1)
-        do i=1,iim,eofskip
-          do j=1+eofskip/2,jjm+1,eofskip
-            ig = 1+ (j-2)*iim +i
+        do i=1,nbp_lon,eofskip
+          do j=1+eofskip/2,nbp_lat,eofskip
+            ig = 1+ (j-2)*nbp_lon +i
 #ifdef NC_DOUBLE
             write(uedata) (real(u(ig,l)),l=1,nlayer)
@@ -82,9 +84,9 @@
       use time_phylmdz_mod, only: daysec, dtphys
       USE comvert_mod, ONLY: aps,bps
+      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
       implicit none
 !
 !     Initialise dumping of profiles for EOF calculations
 !
-#include "dimensions.h"
 
       integer,intent(in) :: ngrid ! total number of physics grid points
@@ -96,5 +98,5 @@
 
       if (firstcall) then
-         npgrid=ngrid+2*(iim-1)
+         npgrid=ngrid+2*(nbp_lon-1)
          firstcall=.false.
       endif
@@ -111,11 +113,11 @@
       open(uehead,file='profiles.hdr',form='formatted')
       write(uehead,*) 0.E+0,0,0,ieofs,1,0
-      write(uehead,*) iim,npgrid/iim,npgrid,nlayer
+      write(uehead,*) nbp_lon,npgrid/nbp_lon,npgrid,nlayer
 
-      do i=1,iim,eofskip
-        do j=1+eofskip/2,jjm+1,eofskip    
-          ig = 1+ (j-2)*iim +i
+      do i=1,nbp_lon,eofskip
+        do j=1+eofskip/2,nbp_lat,eofskip    
+          ig = 1+ (j-2)*nbp_lon +i
           if(j.eq.1) stop 'Problem in ineofdump.F'
-          if(j.eq.jjm+1) stop 'Problem in ineofdump.F'
+          if(j.eq.nbp_lat) stop 'Problem in ineofdump.F'
 #ifdef NC_DOUBLE
           write(uehead,*) real(long(ig)*180./pi),real(lati(ig)*180./pi)
Index: /trunk/LMDZ.MARS/libf/phymars/getslopes.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/getslopes.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/getslopes.F90	(revision 1528)
@@ -3,8 +3,9 @@
 use comgeomfi_h, only: long, lati
 use slope_mod, only: theta_sl, psi_sl
-USE comcstfi_h
+use comcstfi_h, only: g, rad, pi
+use mod_phys_lmdz_para, only: is_parallel
+use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
 implicit none
 
-#include "dimensions.h"
 
 ! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
@@ -16,20 +17,26 @@
 integer,intent(in) :: ngrid ! nnumber of atmospheric columns
 real,intent(in) :: geopot(ngrid)     ! geopotential on phy grid
-real topogrid(iim,jjm+1) ! topography on lat/lon grid with poles and only one -180/180 point
-real latigrid(iim,jjm+1),longgrid(iim,jjm+1) ! meshgrid of latitude and longitude values (radians)
+real topogrid(nbp_lon,nbp_lat) ! topography on lat/lon grid with poles and only one -180/180 point
+real latigrid(nbp_lon,nbp_lat),longgrid(nbp_lon,nbp_lat) ! meshgrid of latitude and longitude values (radians)
 real theta_val ! slope inclination
 real psi_val   ! slope orientation
-real gradx(iim,jjm+1) ! x: latitude-wise topography gradient,  increasing northward
-real grady(iim,jjm+1) ! y: longitude-wise topography gradient, increasing westward
+real gradx(nbp_lon,nbp_lat) ! x: latitude-wise topography gradient,  increasing northward
+real grady(nbp_lon,nbp_lat) ! y: longitude-wise topography gradient, increasing westward
 integer i,j,ig0
 integer id2,idm1 ! a trick to compile testphys1d with debug option
 
+if (is_parallel) then
+  ! This routine only works in serial mode so stop now.
+  write(*,*) "getslopes Error: this routine is not designed to run in parallel"
+  stop
+endif
+
 id2  = 2
-idm1 = iim-1
+idm1 = nbp_lon-1
 
 ! rearrange topography on a 2d array
-do j=2,jjm
-   ig0= 1+(j-2)*iim
-   do i=1,iim
+do j=2,nbp_lat-1
+   ig0= 1+(j-2)*nbp_lon
+   do i=1,nbp_lon
       topogrid(i,j)=geopot(ig0+i)/g
       latigrid(i,j)=lati(ig0+i)
@@ -41,7 +48,7 @@
 latigrid(:,1) = lati(1)
 longgrid(:,1) = long(1)
-topogrid(:,jjm+1) = geopot(ngrid)/g
-latigrid(:,jjm+1) = lati(ngrid)
-longgrid(:,jjm+1) = long(ngrid)
+topogrid(:,nbp_lat) = geopot(ngrid)/g
+latigrid(:,nbp_lat) = lati(ngrid)
+longgrid(:,nbp_lat) = long(ngrid)
 
 
@@ -49,14 +56,14 @@
 ! compute topography gradient
 ! topogrid and rad are both in meters
-do j=2,jjm
-   do i=1,iim
+do j=2,nbp_lat-1
+   do i=1,nbp_lon
      gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1))
      gradx(i,j) = gradx(i,j) / rad
    enddo
-   grady(1,j) = (topogrid(id2,j) - topogrid(iim,j)) / (2*pi+longgrid(id2,j)-longgrid(iim,j)) 
+   grady(1,j) = (topogrid(id2,j) - topogrid(nbp_lon,j)) / (2*pi+longgrid(id2,j)-longgrid(nbp_lon,j)) 
    grady(1,j) = grady(1,j) / rad
-   grady(iim,j) = (topogrid(1,j) - topogrid(idm1,j)) / (2*pi+longgrid(1,j)-longgrid(idm1,j)) 
-   grady(iim,j) = grady(iim,j) / rad
-   do i=2,iim-1
+   grady(nbp_lon,j) = (topogrid(1,j) - topogrid(idm1,j)) / (2*pi+longgrid(1,j)-longgrid(idm1,j)) 
+   grady(nbp_lon,j) = grady(nbp_lon,j) / rad
+   do i=2,nbp_lon-1
      grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j)) 
      grady(i,j) = grady(i,j) / rad
@@ -66,6 +73,6 @@
 gradx(:,1) = 0.
 grady(:,1) = 0.
-gradx(:,jjm+1) = 0.
-grady(:,jjm+1) = 0.
+gradx(:,nbp_lat) = 0.
+grady(:,nbp_lat) = 0.
 
 
@@ -74,8 +81,8 @@
 theta_sl(:) = 0.
 psi_sl(:)   = 0.
-do j=2,jjm
-   do i=1,iim
+do j=2,nbp_lat-1
+   do i=1,nbp_lon
    
-     ig0= 1+(j-2)*iim
+     ig0= 1+(j-2)*nbp_lon
    
      theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 ))
Index: /trunk/LMDZ.MARS/libf/phymars/inistats.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 1528)
@@ -5,9 +5,8 @@
       USE comcstfi_h, ONLY: pi
       USE time_phylmdz_mod, ONLY: daysec,dtphys
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       implicit none
 
-      include "dimensions.h"
-      include "paramet.h"
-      include "comgeom.h"
       include "statto.h"
       include "netcdf.inc"
@@ -16,5 +15,6 @@
       integer :: nid
       integer :: l,nsteppd
-      real, dimension(llm) ::  sig_s
+      real, dimension(nbp_lev) ::  sig_s
+      real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
       integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
       real, dimension(istime) :: lt
@@ -40,8 +40,12 @@
       write (*,*) 
 
-      do l= 1, llm
+      do l= 1, nbp_lev
          sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
          pseudoalt(l)=-10.*log(presnivs(l)/preff)   
       enddo
+      
+      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+      !add extra redundant point (180 degrees, since lon_reg starts at -180
+      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
 
       if (is_master) then
@@ -54,8 +58,8 @@
       endif
 
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_lat)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_lon)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
-      ierr = NF_DEF_DIM (nid, "llmp1", llm+1, idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
+      ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
       ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
 
@@ -69,14 +73,14 @@
      &            "degrees_north",1,idim_lat,nvarid,ierr)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
 #else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
 #endif
       call def_var_stats(nid,"longitude","East longitude",
      &            "degrees_east",1,idim_lon,nvarid,ierr)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
 #else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 
Index: /trunk/LMDZ.MARS/libf/phymars/iniwrite.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/iniwrite.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/iniwrite.F	(revision 1528)
@@ -1,3 +1,3 @@
-      SUBROUTINE iniwrite(nid,idayref,phis)
+      SUBROUTINE iniwrite(nid,idayref,phis,area)
 
       use comsoil_h, only: mlayer, nsoilmx
@@ -8,4 +8,6 @@
       USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       IMPLICIT NONE
 
@@ -26,8 +28,5 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 c   Arguments:
@@ -36,5 +35,6 @@
       integer,intent(in) :: nid        ! NetCDF file ID
       INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
-      real,intent(in) :: phis(ip1jmp1) ! surface geopotential
+      real,intent(in) :: phis(nbp_lon+1,nbp_lat) ! surface geopotential
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -44,7 +44,8 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
-
-      integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
-      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
+
+      integer :: nvarid,idim_index,idim_rlonv
+      integer :: idim_rlatu,idim_llmp1,idim_llm
       integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
       integer, dimension(2) :: id  
@@ -54,7 +55,7 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = real(iim)
-      tab_cntrl(2)  = real(jjm)
-      tab_cntrl(3)  = real(llm)
+      tab_cntrl(1)  = real(nbp_lon)
+      tab_cntrl(2)  = real(nbp_lat-1)
+      tab_cntrl(3)  = real(nbp_lev)
       tab_cntrl(4)  = real(idayref)
       tab_cntrl(5)  = rad
@@ -101,10 +102,10 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
-      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
-      ierr = NF_DEF_DIM (nid, "interlayer", (llm+1), idim_llmp1)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+!      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
+!      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
+      ierr = NF_DEF_DIM (nid, "interlayer", (nbp_lev+1), idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
       ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
 c
@@ -131,18 +132,19 @@
 c --------------------------
 c  longitudes and latitudes
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
-     .                       "Longitudes at u nodes")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
-#endif
+!
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+!     .                       "Longitudes at u nodes")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
+!#endif
 c
 c --------------------------
@@ -158,10 +160,14 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+c
+c --------------------------
+      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+      !add extra redundant point (180 degrees, since lon_reg starts at -180
+      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -175,7 +181,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -201,18 +207,18 @@
 c
 c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
-     .                       "Latitudes at v nodes")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
-#endif
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
+!     .                       "Latitudes at v nodes")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
+!#endif
 c
 c --------------------------
@@ -276,40 +282,40 @@
 c  Mesh area and conversion coefficients cov. <-> contra. <--> natural
 
-      id(1)=idim_rlonu
-      id(2)=idim_rlatu
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
-#endif
-c
-      id(1)=idim_rlonv
-      id(2)=idim_rlatv
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
-#endif
+!      id(1)=idim_rlonu
+!      id(2)=idim_rlatu
+c
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+!     .             "Conversion coefficients cov <--> natural")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+!#endif
+c
+!      id(1)=idim_rlonv
+!      id(2)=idim_rlatv
+c
+c --------------------------
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+!     .             "Conversion coefficients cov <--> natural")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+!#endif
 c
       id(1)=idim_rlonv
@@ -327,7 +333,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 c
@@ -352,7 +358,8 @@
 c
 
-      write(*,*)'iniwrite: iim,jjm,llm,idayref',iim,jjm,llm,idayref
+      write(*,*)'iniwrite: nbp_lon,nbp_lat,nbp_lev,idayref',
+     & nbp_lon,nbp_lat,nbp_lev,idayref
       write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',
-     s rad,omeg,g,mugaz,rcp
+     & rad,omeg,g,mugaz,rcp
       write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
 
Index: /trunk/LMDZ.MARS/libf/phymars/iniwritesoil.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/iniwritesoil.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/iniwritesoil.F90	(revision 1528)
@@ -1,3 +1,3 @@
-subroutine iniwritesoil(nid,ngrid)
+subroutine iniwritesoil(nid,ngrid,inertia,area)
 
 ! initialization routine for 'writediagoil'. Here we create/define
@@ -5,17 +5,18 @@
 ! (time-independent) parameters.
 
-use comsoil_h, only: mlayer, inertiedat, nsoilmx
-USE comcstfi_h
+use comsoil_h, only: mlayer, nsoilmx
+USE comcstfi_h, only: pi
+USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+use mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
 
 implicit none
 
-#include"dimensions.h"
-#include"paramet.h"
-#include"comgeom.h"
-#include"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
 integer,intent(in) :: ngrid
 integer,intent(in) :: nid ! NetCDF output file ID
+real,intent(in) :: inertia(nbp_lon+1,nbp_lat,nsoilmx)
+real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 ! Local variables:
@@ -30,6 +31,7 @@
 integer,dimension(3) :: dimids ! to store IDs of dimensions of a variable
 character(len=60) :: text ! to store some text
-real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
+real,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
 integer :: i,j,l,ig0
+real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
 
 ! 1. Define the dimensions
@@ -38,11 +40,9 @@
 
 ! Define the dimensions
-ierr=NF_DEF_DIM(nid,"longitude",iip1,idim_rlonv)
-! iip1 known from paramet.h
+ierr=NF_DEF_DIM(nid,"longitude",nbp_lon+1,idim_rlonv)
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define longitude dimension"
 endif
-ierr=NF_DEF_DIM(nid,"latitude",jjp1,idim_rlatu)
-! jjp1 known from paramet.h
+ierr=NF_DEF_DIM(nid,"latitude",nbp_lat,idim_rlatu)
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define latitude dimension"
@@ -81,11 +81,15 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
+lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+!add extra redundant point (180 degrees, since lon_reg starts at -180
+lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+
 ! Write longitude to file
 ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,rlonv*(180./pi))
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,rlonv*(180./pi))
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,lon_reg_ext*(180./pi))
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,lon_reg_ext*(180./pi))
 #endif
 ! Note: rlonv is known from comgeom.h and pi from comcstfi.h
@@ -117,9 +121,8 @@
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,rlatu*(180./pi))
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,rlatu*(180./pi))
-#endif
-! Note: rlatu is known from comgeom.h and pi from comcstfi.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,lat_reg*(180./pi))
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,lat_reg*(180./pi))
+#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write longitude variable"
@@ -209,9 +212,8 @@
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,aire)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,aire)
-#endif
-! Note: aire is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,area)
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,area)
+#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write area variable"
@@ -240,30 +242,11 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
-! Recast data along 'dynamics' grid
-! Note: inertiedat is known from comsoil_h
-
-do l=1,nsoilmx
-  ! handle the poles
-  do i=1,iip1
-    data3(i,1,l)=inertiedat(1,l)
-    data3(i,jjp1,l)=inertiedat(ngrid,l)
-  enddo
-  ! rest of the grid
-  do j=2,jjm
-    ig0=1+(j-2)*iim
-    do i=1,iim
-      data3(i,j,l)=inertiedat(ig0+i,l)
-    enddo
-    data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude
-  enddo
-enddo ! of do l=1,nsoilmx
-
 ! Write data2 to file
 ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,data3)
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,inertia)
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,inertia)
 #endif
 if (ierr.ne.NF_NOERR) then
Index: /trunk/LMDZ.MARS/libf/phymars/mkstat.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/mkstat.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/mkstat.F90	(revision 1528)
@@ -11,23 +11,21 @@
 
 use mod_phys_lmdz_para, only : is_master
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
 
 implicit none
 
-#include "dimensions.h"
-#include "statto.h"
-#include "netcdf.inc"
+include "statto.h"
+include "netcdf.inc"
 
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
 integer :: ierr,nid,nbvar,i,ndims,lt,nvarid
 integer, dimension(4) :: id,varid,start,size
 integer, dimension(5) :: dimids
 character (len=50) :: name,nameout,units,title
-real, dimension(iip1,jjp1,llm) :: sum3d,square3d,mean3d,sd3d
-real, dimension(iip1,jjp1) :: sum2d,square2d,mean2d,sd2d
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: sum3d,square3d,mean3d,sd3d
+real, dimension(nbp_lon+1,nbp_lat) :: sum2d,square2d,mean2d,sd2d
 real, dimension(istime) :: time
-real, dimension(jjp1) :: lat
-real, dimension(iip1) :: lon
-real, dimension(llm) :: alt
+real, dimension(nbp_lat) :: lat
+real, dimension(nbp_lon+1) :: lon
+real, dimension(nbp_lev) :: alt
 logical :: lcopy=.true.
 !integer :: latid,lonid,altid,timeid
@@ -107,5 +105,5 @@
 !      dimout(4)=timeid
 
-      size=(/iip1,jjp1,llm,1/)
+      size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
       do lt=1,istime
          start=(/1,1,1,lt/)
@@ -137,5 +135,5 @@
 !      dimout(3)=timeid
 
-      size=(/iip1,jjp1,1,0/)
+      size=(/nbp_lon+1,nbp_lat,1,0/)
       do lt=1,istime
          start=(/1,1,lt,0/)
Index: /trunk/LMDZ.MARS/libf/phymars/newcondens.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/newcondens.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/newcondens.F	(revision 1528)
@@ -61,6 +61,5 @@
 c    ------------------
 c
-#include "dimensions.h"
-#include "callkeys.h"
+      include "callkeys.h"
 
 c-----------------------------------------------------------------------
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1528)
@@ -148,5 +148,4 @@
 c    ------------------
 
-#include "dimensions.h"
 #include "callkeys.h"
 #include "comg1d.h"
Index: /trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 1528)
@@ -13,4 +13,5 @@
 #endif
       USE comcstfi_h
+      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
       IMPLICIT NONE
 c=======================================================================
@@ -22,5 +23,4 @@
 c   Declarations:
 c   -------------
-#include "dimensions.h"
 #include "callkeys.h"
 #include "datafile.h"
@@ -55,5 +55,5 @@
       INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
 #endif
-      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
+      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis
 
       REAL step,count,ratiolat
@@ -228,6 +228,8 @@
              zelon = i - 180.
              zelat = 90. - j 
-            if ((abs(lati_glo(ig)*180./pi-zelat).le.90./real(jjm)) .and.
-     &        (abs(long_glo(ig)*180./pi-zelon).le.180./real(iim))) then
+            if ((abs(lati_glo(ig)*180./pi-zelat).le.
+     &           90./real(nbp_lat-1)) .and.
+     &          (abs(long_glo(ig)*180./pi-zelon).le.
+     &           180./real(nbp_lon))) then
               ! count all points in that GCM grid point
               nb_ice(ig,1) = nb_ice(ig,1) + 1
@@ -240,8 +242,8 @@
 
         ! projection of nb_ice on GCM lat and lon axes
-          latice(1+(ig-2)/iim,:) =
-     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
-          lonice(1+mod(ig-2,iim),:) = 
-     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
+          latice(1+(ig-2)/nbp_lon,:) =
+     &     latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
+          lonice(1+mod(ig-2,nbp_lon),:) = 
+     &     lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
 
          enddo ! of do ig=2,klon_glo-1
@@ -253,6 +255,6 @@
          latice(1,:)   = nb_ice(1,:)
          lonice(1,:)   = nb_ice(1,:)
-         latice(jjm,:) = nb_ice(ngrid,:)
-         lonice(iim,:) = nb_ice(ngrid,:)
+         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
+         lonice(nbp_lon,:) = nb_ice(ngrid,:)
       
      
@@ -271,5 +273,5 @@
     
          ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
-         do i=1,jjm/2
+         do i=1,(nbp_lat-1)/2
           step  = 1. ! threshold to add ice cap
           count = 0. ! number of ice GCM caps at this latitude
@@ -280,20 +282,20 @@
           ! put ice caps while there is not enough ice,
           ! as long as the threshold is above 20%
-          do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
+          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
            count = 0.
            ! loop over GCM longitudes
-           do j=1,iim
+           do j=1,nbp_lon
             ! if the detected ice ratio in the GCM grid point 
             ! is more than 'step', then add ice
-            if (real(nb_ice((i-1)*iim+1+j,2)) 
-     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
-                  watercaptag_glo((i-1)*iim+1+j) = .true.
+            if (real(nb_ice((i-1)*nbp_lon+1+j,2)) 
+     &        / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
+                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
                   count = count + 1
             endif
-           enddo ! of do j=1,iim
-           !print*, 'step',step,count,ratiolat*iim
+           enddo ! of do j=1,nbp_lon
+           !print*, 'step',step,count,ratiolat*nbp_lon
            step = step - 0.01
           enddo ! of do while
-          !print*, 'step',step,count,ratiolat*iim
+          !print*, 'step',step,count,ratiolat*nbp_lon
 
          enddo ! of do i=1,jjm/2
@@ -304,5 +306,5 @@
          print*,'Surfini: predefined ice caps'
       
-         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
+         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
            
           print*,'water ice caps distribution for 32x24 resolution'
@@ -317,5 +319,5 @@
 !---------------------   OUTLIERS  ----------------------------
 
-         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
+         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
 
           print*,'water ice caps distribution for 64x48 resolution'
@@ -355,5 +357,5 @@
         
           print*,'No predefined ice location for this resolution :',
-     &           iim,jjm
+     &           nbp_lon,nbp_lat-1
           print*,'Please change icelocationmode in surfini.F'
           print*,'Or add some new definitions ...'
Index: /trunk/LMDZ.MARS/libf/phymars/writediagfi.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 1528)
@@ -40,15 +40,14 @@
 !=================================================================
       use surfdat_h, only: phisfi
+      use comgeomphy, only: airephy
       use time_phylmdz_mod, only: ecritphy, day_step, iphysiq, day_ini
       USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
      &                               is_master, gather
-      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat, nbp_lev
       implicit none
 
 ! Commons
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 ! Arguments on input:
@@ -56,21 +55,21 @@
       character (len=*),intent(in) :: nom,titre,unite
       integer,intent(in) :: dim
-      real,intent(in) :: px(ngrid,llm)
+      real,intent(in) :: px(ngrid,nbp_lev)
 
 ! Local variables:
 
-      real*4 dx3(iip1,jjp1,llm) ! to store a 3D data set
-      real*4 dx2(iip1,jjp1)     ! to store a 2D (surface) data set
-      real*4 dx1(llm)           ! to store a 1D (column) data set
+      real*4 dx3(nbp_lon+1,nbp_lat,nbp_lev) ! to store a 3D data set
+      real*4 dx2(nbp_lon+1,nbp_lat)     ! to store a 2D (surface) data set
+      real*4 dx1(nbp_lev)           ! to store a 1D (column) data set
       real*4 dx0
 
       real*4,save :: date
 
-      REAL phis(ip1jmp1)
+      REAL phis((nbp_lon+1),nbp_lat)
+      REAL area((nbp_lon+1),nbp_lat)
 
       integer irythme
       integer ierr,ierr2
-      integer iq
-      integer i,j,l,zmax , ig0
+      integer i,j,l, ig0
 
       integer,save :: zitau=0
@@ -98,14 +97,16 @@
 #ifdef CPP_PARA
 ! Added to work in parallel mode
-      real dx3_glop(klon_glo,llm)
-      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
+      real dx3_glop(klon_glo,nbp_lev)
+      real dx3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set
       real dx2_glop(klon_glo)
-      real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+      real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
       real px2(ngrid)
-!      real dx1_glo(llm)          ! to store a 1D (column) data set
+!      real dx1_glo(nbp_lev)          ! to store a 1D (column) data set
 !      real dx0_glo
       real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
 #else
       real phisfi_glo(ngrid) ! surface geopotential on global physics grid
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -181,6 +182,9 @@
           ! Gather phisfi() geopotential on physics grid
           call Gather(phisfi,phisfi_glo)
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
 #else
          phisfi_glo(:)=phisfi(:)
+         areafi_glo(:)=airephy(:)
 #endif
 
@@ -209,7 +213,25 @@
          ierr = NF_ENDDEF(nid)
 
+         ! Build phis() and area()
+         do i=1,nbp_lon+1 ! poles
+           phis(i,1)=phisfi_glo(1)
+           phis(i,nbp_lat)=phisfi_glo(klon_glo)
+           ! for area, divide at the poles by nbp_lon
+           area(i,1)=areafi_glo(1)/nbp_lon
+           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+         enddo
+         do j=2,nbp_lat-1
+           ig0= 1+(j-2)*nbp_lon
+           do i=1,nbp_lon
+              phis(i,j)=phisfi_glo(ig0+i)
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           phis(nbp_lon+1,j)=phis(1,j)
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+         
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
-         call iniwrite(nid,day_ini,phis)
+         call iniwrite(nid,day_ini,phis,area)
 
          endif ! of if (is_master)
@@ -290,6 +312,6 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -298,15 +320,15 @@
 !         Passage variable physique -->  variable dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-           DO l=1,llm
-             DO i=1,iip1
+           DO l=1,nbp_lev
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -338,7 +360,7 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
-           edges(3)=llm
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
+           edges(3)=nbp_lev
            edges(4)=1
 !#ifdef NC_DOUBLE
@@ -374,6 +396,6 @@
             call Grid1Dto2D_glo(dx2_glop,dx2_glo)
             ! copy dx2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=dx2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -384,14 +406,14 @@
 !         recast (copy) variable from physics grid to dynamics grid
 
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx2(i,1)=px(1,1)
-                dx2(i,jjp1)=px(ngrid,1)
+                dx2(i,nbp_lat)=px(ngrid,1)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx2(i,j)=px(ig0+i,1)
                 ENDDO
-                dx2(iip1,j)=dx2(1,j)
+                dx2(nbp_lon+1,j)=dx2(1,j)
              ENDDO
 #endif
@@ -420,6 +442,6 @@
            corner(2)=1
            corner(3)=ntime
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=1
 
@@ -451,5 +473,5 @@
 !         Passage variable physique -->  physique dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-          do l=1,llm
+          do l=1,nbp_lev
             dx1(l)=px(1,l)
           enddo
@@ -473,5 +495,5 @@
            corner(2)=ntime
            
-           edges(1)=llm
+           edges(1)=nbp_lev
            edges(2)=1
 !#ifdef NC_DOUBLE
Index: /trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90	(revision 1528)
@@ -12,16 +12,14 @@
 ! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
 
-use comsoil_h, only: nsoilmx
+use comsoil_h, only: nsoilmx, inertiedat
+use comgeomphy, only: airephy
 use time_phylmdz_mod, only: ecritphy, day_step, iphysiq
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
+                              nbp_lon, nbp_lat
 
 implicit none
 
-#include"dimensions.h"
-#include"paramet.h"
-!#include"control.h"
-!#include"comsoil.h"
-#include"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
@@ -33,10 +31,8 @@
 integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
 real,dimension(ngrid,nsoilmx),intent(in) :: px ! variable
-! Note: nsoilmx is a parameter set in 'comsoil_h'
 
 ! Local variables:
-real*4,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
-! Note iip1,jjp1 known from paramet.h; nsoilmx known from comsoil_h
-real*4,dimension(iip1,jjp1) :: data2 ! to store 2D data
+real*4,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
+real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
 real*4 :: data0 ! to store 0D data
 integer :: i,j,l ! for loops
@@ -44,4 +40,11 @@
 
 real*4,save :: date ! time counter (in elapsed days)
+
+real :: inertia((nbp_lon+1),nbp_lat,nsoilmx)
+real :: area((nbp_lon+1),nbp_lat)
+
+real :: inertiafi_glo(klon_glo,nsoilmx)
+real :: areafi_glo(klon_glo)
+
 integer,save :: isample ! sample rate at which data is to be written to output
 integer,save :: ntime=0 ! counter to internally store time steps
@@ -61,7 +64,7 @@
 ! Added to work in parallel mode
 real dx3_glop(klon_glo,nsoilmx)
-real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
+real dx3_glo(nbp_lon,nbp_lat,nsoilmx) ! to store a global 3D data set
 real dx2_glop(klon_glo)
-real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
 real px2(ngrid)
 #endif
@@ -92,8 +95,38 @@
     stop
    endif
+
+#ifdef CPP_PARA
+   ! Gather inertiedat() soil thermal inertia on physics grid
+   call Gather(inertiedat,inertiafi_glo)
+   ! Gather airephy() mesh area on physics grid
+   call Gather(airephy,areafi_glo)
+#else
+         inertiafi_glo(:,:)=inertiedat(:,:)
+         areafi_glo(:)=airephy(:)
+#endif
+
+   ! build inertia() and area()
+   do i=1,nbp_lon+1 ! poles
+     inertia(i,1,1:nsoilmx)=inertiafi_glo(1,1:nsoilmx)
+     inertia(i,nbp_lat,1:nsoilmx)=inertiafi_glo(klon_glo,1:nsoilmx)
+     ! for area, divide at the poles by nbp_lon
+     area(i,1)=areafi_glo(1)/nbp_lon
+     area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+   enddo
+   do j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     do i=1,nbp_lon
+        inertia(i,j,1:nsoilmx)=inertiafi_glo(ig0+i,1:nsoilmx)
+        area(i,j)=areafi_glo(ig0+i)
+     enddo
+     ! handle redundant point in longitude
+     inertia(nbp_lon+1,j,1:nsoilmx)=inertia(1,j,1:nsoilmx)
+     area(nbp_lon+1,j)=area(1,j)
+   enddo
+   
+   ! write "header" of file (longitudes, latitudes, geopotential, ...)
+   call iniwritesoil(nid,ngrid,inertia,area)
+
   endif ! of if (is_master)
-
-  ! Define dimensions and axis attributes
-  call iniwritesoil(nid,ngrid)
   
   ! set zitau to -1 to be compatible with zitau incrementation step below
@@ -149,6 +182,6 @@
     call Grid1Dto2D_glo(dx3_glop,dx3_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-    data3(iip1,:,:)=data3(1,:,:)
+    data3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+    data3(nbp_lon+1,:,:)=data3(1,:,:)
   endif
 !$OMP END MASTER
@@ -157,15 +190,15 @@
   do l=1,nsoilmx
     ! handle the poles
-    do i=1,iip1
+    do i=1,nbp_lon+1
       data3(i,1,l)=px(1,l)
-      data3(i,jjp1,l)=px(ngrid,l)
+      data3(i,nbp_lat,l)=px(ngrid,l)
     enddo
     ! rest of the grid
-    do j=2,jjm
-      ig0=1+(j-2)*iim
-      do i=1,iim
+    do j=2,nbp_lat-1
+      ig0=1+(j-2)*nbp_lon
+      do i=1,nbp_lon
         data3(i,j,l)=px(ig0+i,l)
       enddo
-      data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude
+      data3(nbp_lon+1,j,l)=data3(1,j,l) ! extra (modulo) longitude
     enddo
   enddo
@@ -196,6 +229,6 @@
   corners(4)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=nsoilmx
   edges(4)=1
@@ -224,6 +257,6 @@
     call Grid1Dto2D_glo(dx2_glop,dx2_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    data2(1:iim,:)=dx2_glo(1:iim,:)
-    data2(iip1,:)=data2(1,:)
+    data2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+    data2(nbp_lon+1,:)=data2(1,:)
   endif
 !$OMP END MASTER
@@ -231,15 +264,15 @@
 #else
   ! handle the poles
-  do i=1,iip1
+  do i=1,nbp_lon+1
     data2(i,1)=px(1,1)
-    data2(i,jjp1)=px(ngrid,1)
+    data2(i,nbp_lat)=px(ngrid,1)
   enddo
   ! rest of the grid
-  do j=2,jjm
-    ig0=1+(j-2)*iim
-    do i=1,iim
+  do j=2,nbp_lat-1
+    ig0=1+(j-2)*nbp_lon
+    do i=1,nbp_lon
       data2(i,j)=px(ig0+i,1)
     enddo
-    data2(iip1,j)=data2(1,j) ! extra (modulo) longitude
+    data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
   enddo
 #endif
@@ -267,6 +300,6 @@
   corners(3)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=1
   
Index: /trunk/LMDZ.MARS/libf/phymars/wstats.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 1527)
+++ /trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 1528)
@@ -2,20 +2,17 @@
 
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
-use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
+                              nbp_lon, nbp_lat, nbp_lev
 implicit none
 
-#include "dimensions.h"
-#include "statto.h"
-#include "netcdf.inc"
+include "statto.h"
+include "netcdf.inc"
 
 integer,intent(in) :: ngrid
 character (len=*),intent(in) :: nom,titre,unite
 integer,intent(in) :: dim
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
-real,intent(in) :: px(ngrid,llm)
-real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
-real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
+real,intent(in) :: px(ngrid,nbp_lev)
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: mean3d,sd3d,dx3
+real, dimension(nbp_lon+1,nbp_lat) :: mean2d,sd2d,dx2
 character (len=50) :: namebis
 character (len=50), save :: firstvar
@@ -31,16 +28,16 @@
 ! Added to work in parallel mode
 #ifdef CPP_PARA
-real px3_glop(klon_glo,llm) ! to store a 3D data set on global physics grid
-real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
 real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
-real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
 real px2(ngrid)
-real px3(ngrid,llm)
+real px3(ngrid,nbp_lev)
 #else
 ! When not running in parallel mode:
-real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
-real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
 real px2_glop(ngrid) ! to store a 2D data set on global physics grid
-real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
 #endif
 
@@ -64,5 +61,5 @@
 #ifdef CPP_PARA
  if (dim.eq.3) then
-  px3(1:ngrid,1:llm)=px(1:ngrid,1:llm)
+  px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
   ! Gather fieds on a "global" (without redundant longitude) array
   call Gather(px3,px3_glop)
@@ -71,6 +68,6 @@
     call Grid1Dto2D_glo(px3_glop,px3_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    dx3(1:iim,:,:)=px3_glo(1:iim,:,:)
-    dx3(iip1,:,:)=dx3(1,:,:)
+    dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
+    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
   endif
 !$OMP END MASTER
@@ -84,6 +81,6 @@
             call Grid1Dto2D_glo(px2_glop,px2_glo)
             ! copy px2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=px2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -92,14 +89,14 @@
 #else
   if (dim.eq.3) then
-    px3_glop(:,1:llm)=px(:,1:llm)
+    px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
 !  Passage variable physique -->  variable dynamique
-    DO l=1,llm
-      DO i=1,iim
+    DO l=1,nbp_lev
+      DO i=1,nbp_lon
          px3_glo(i,1,l)=px(1,l)
-         px3_glo(i,jjp1,l)=px(ngrid,l)
+         px3_glo(i,nbp_lat,l)=px(ngrid,l)
       ENDDO
-      DO j=2,jjm
-         ig0= 1+(j-2)*iim
-         DO i=1,iim
+      DO j=2,nbp_lat-1
+         ig0= 1+(j-2)*nbp_lon
+         DO i=1,nbp_lon
             px3_glo(i,j,l)=px(ig0+i,l)
          ENDDO
@@ -109,11 +106,11 @@
     px2_glop(:)=px(:,1)
 !    Passage variable physique -->  physique dynamique
-   DO i=1,iim
+   DO i=1,nbp_lon
      px2_glo(i,1)=px(1,1)
-     px2_glo(i,jjp1)=px(ngrid,1)
+     px2_glo(i,nbp_lat)=px(ngrid,1)
    ENDDO
-   DO j=2,jjm
-     ig0= 1+(j-2)*iim
-     DO i=1,iim
+   DO j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     DO i=1,nbp_lon
         px2_glo(i,j)=px(ig0+i,1)
      ENDDO
@@ -195,10 +192,10 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/iip1,jjp1,llm,1/)
+      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
       mean3d(:,:,:)=0
       sd3d(:,:,:)=0
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lev,1,0/)
       mean2d(:,:)=0
       sd2d(:,:)=0
@@ -208,5 +205,5 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/iip1,jjp1,llm,1/)
+      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
@@ -223,5 +220,5 @@
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lat,1,0/)
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
@@ -242,9 +239,9 @@
 
 if (dim.eq.3) then
-  dx3(1:iim,:,:)=px3_glo(:,:,:)
-  dx3(iip1,:,:)=dx3(1,:,:)
+  dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
+  dx3(nbp_lon+1,:,:)=dx3(1,:,:)
 else ! dim.eq.2
-  dx2(1:iim,:)=px2_glo(:,:)
-  dx2(iip1,:)=dx2(1,:)
+  dx2(1:nbp_lon,:)=px2_glo(:,:)
+  dx2(nbp_lon+1,:)=dx2(1,:)
 endif
 
@@ -287,40 +284,37 @@
 !======================================================
 subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
 
 implicit none
 
-include "dimensions.h"
 include "netcdf.inc"
 
 integer, intent(in) :: nid,varid,dim,indx,ngrid
-real, dimension(ngrid,llm), intent(in) :: px
+real, dimension(ngrid,nbp_lev), intent(in) :: px
 integer, intent(out) :: ierr
-
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
 
 integer :: l,i,j,ig0
 integer, dimension(4) :: start,sizes
-real, dimension(iip1,jjp1,llm) :: dx3
-real, dimension(iip1,jjp1) :: dx2
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
+real, dimension(nbp_lon+1,nbp_lat) :: dx2
 
 if (dim.eq.3) then
 
    start=(/1,1,1,indx/)
-   sizes=(/iip1,jjp1,llm,1/)
+   sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
 
 !  Passage variable physique -->  variable dynamique
 
-   DO l=1,llm
-      DO i=1,iip1
+   DO l=1,nbp_lev
+      DO i=1,nbp_lon+1
          dx3(i,1,l)=px(1,l)
-         dx3(i,jjp1,l)=px(ngrid,l)
+         dx3(i,nbp_lat,l)=px(ngrid,l)
       ENDDO
-      DO j=2,jjm
-         ig0= 1+(j-2)*iim
-         DO i=1,iim
+      DO j=2,nbp_lat-1
+         ig0= 1+(j-2)*nbp_lon
+         DO i=1,nbp_lon
             dx3(i,j,l)=px(ig0+i,l)
          ENDDO
-         dx3(iip1,j,l)=dx3(1,j,l)
+         dx3(nbp_lon+1,j,l)=dx3(1,j,l)
       ENDDO
    ENDDO
@@ -335,18 +329,18 @@
 
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lat,1,0/)
 
 !    Passage variable physique -->  physique dynamique
 
-  DO i=1,iip1
+  DO i=1,nbp_lon+1
      dx2(i,1)=px(1,1)
-     dx2(i,jjp1)=px(ngrid,1)
+     dx2(i,nbp_lat)=px(ngrid,1)
   ENDDO
-  DO j=2,jjm
-     ig0= 1+(j-2)*iim
-     DO i=1,iim
+  DO j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     DO i=1,nbp_lon
         dx2(i,j)=px(ig0+i,1)
      ENDDO
-     dx2(iip1,j)=dx2(1,j)
+     dx2(nbp_lon+1,j)=dx2(1,j)
   ENDDO
 
@@ -376,5 +370,5 @@
 implicit none
 
-#include "netcdf.inc"
+include "netcdf.inc"
 
 integer,intent(in) :: nid ! NetCDF file ID
