Index: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/mod_interface_dyn_phys.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynphy_lonlat/mod_interface_dyn_phys.F90	(revision 1543)
+++ trunk/LMDZ.GENERIC/libf/dynphy_lonlat/mod_interface_dyn_phys.F90	(revision 1543)
@@ -0,0 +1,62 @@
+! 
+! $Id: mod_interface_dyn_phys.F90 2351 2015-08-25 15:14:59Z emillour $
+!
+MODULE mod_interface_dyn_phys
+  INTEGER,SAVE,dimension(:),allocatable :: index_i
+  INTEGER,SAVE,dimension(:),allocatable :: index_j
+  
+  
+CONTAINS
+  
+#ifdef CPP_PARA
+! Interface with parallel physics,
+  SUBROUTINE Init_interface_dyn_phys
+    USE mod_phys_lmdz_mpi_data
+    IMPLICIT NONE
+    include 'dimensions.h'    
+    
+    INTEGER :: i,j,k
+    
+    ALLOCATE(index_i(klon_mpi))
+    ALLOCATE(index_j(klon_mpi))
+    
+    k=1
+    IF (is_north_pole) THEN
+      index_i(k)=1
+      index_j(k)=1
+      k=2
+    ELSE
+      DO i=ii_begin,iim
+	index_i(k)=i
+	index_j(k)=jj_begin
+	k=k+1
+       ENDDO
+    ENDIF
+    
+    DO j=jj_begin+1,jj_end-1
+      DO i=1,iim
+	index_i(k)=i
+	index_j(k)=j
+	k=k+1
+      ENDDO
+    ENDDO
+    
+    IF (is_south_pole) THEN
+      index_i(k)=1
+      index_j(k)=jj_end
+    ELSE
+      DO i=1,ii_end
+	index_i(k)=i
+	index_j(k)=jj_end
+	k=k+1
+       ENDDO
+    ENDIF
+  
+  END SUBROUTINE Init_interface_dyn_phys 
+#else
+  SUBROUTINE Init_interface_dyn_phys
+  ! dummy routine for seq case
+  END SUBROUTINE Init_interface_dyn_phys 
+#endif
+! of #ifdef CPP_PARA
+END MODULE mod_interface_dyn_phys
Index: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/iniphysiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/iniphysiq_mod.F90	(revision 1542)
+++ trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/iniphysiq_mod.F90	(revision 1543)
@@ -3,11 +3,14 @@
 CONTAINS
 
-subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
-               rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn,  &
+subroutine iniphysiq(ii,jj,nlayer, &
+                     nbp, communicator, &
+                     punjours, pdayref,ptimestep, &
+                     rlatudyn,rlatvdyn,rlonudyn,rlonvdyn, &
+                     airedyn,cudyn,cvdyn, &
                      prad,pg,pr,pcpp,iflag_phys)
 
-use dimphy, only : klev ! number of atmospheric levels
-use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns
-                                       ! (on full grid)
+use dimphy, only : init_dimphy
+use mod_grid_phy_lmdz, only : klon_glo, & ! number of atmospheric columns (on full grid)
+                              regular_lonlat  ! regular longitude-latitude grid type
 use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid)
                                klon_omp_begin, & ! start index of local omp subgrid
@@ -15,10 +18,14 @@
                                klon_mpi_begin ! start indes of columns (on local mpi grid)
 use control_mod, only: nday 
-use comgeomphy, only : initcomgeomphy, &
-                       cell_area, & ! physics grid area (m2)
-                       dx, & ! cu coeff. (u_covariant = cu * u)
-                       dy, & ! cv coeff. (v_covariant = cv * v)
-                       longitude, & ! longitudes (rad)
-                       latitude ! latitudes (rad)
+use geometry_mod, only: init_geometry, &
+                        cell_area, & ! physics grid area (m2)
+                        longitude, & ! longitudes (rad)
+                        latitude ! latitudes (rad)
+!use comgeomphy, only : initcomgeomphy, &
+!                       cell_area, & ! physics grid area (m2)
+!                       dx, & ! cu coeff. (u_covariant = cu * u)
+!                       dy, & ! cv coeff. (v_covariant = cv * v)
+!                       longitude, & ! longitudes (rad)
+!                       latitude ! latitudes (rad)
 use surf_heat_transp_mod, only: ini_surf_heat_transp
 use infotrac, only : nqtot ! number of advected tracers
@@ -26,8 +33,10 @@
 USE comvert_mod, ONLY: ap,bp,preff
 use inifis_mod, only: inifis
+use physics_distribution_mod, only: init_physics_distribution
 use regular_lonlat_mod, only: init_regular_lonlat, &
                               east, west, north, south, &
                               north_east, north_west, &
                               south_west, south_east
+use mod_interface_dyn_phys, only: init_interface_dyn_phys
 use ioipsl_getin_p_mod, only: getin_p
 
@@ -47,4 +56,6 @@
 integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
 integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
+integer,intent(in) :: nbp ! number of physics columns for this MPI process
+integer,intent(in) :: communicator ! MPI communicator
 real,intent(in) :: rlatudyn(jj+1) ! latitudes of the physics grid
 real,intent(in) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
@@ -59,5 +70,5 @@
 
 integer :: ibegin,iend,offset
-integer :: i,j
+integer :: i,j,k
 character(len=20) :: modname='iniphysiq'
 character(len=80) :: abort_message
@@ -71,31 +82,28 @@
 
 ! global array, on full physics grid:
-real,allocatable :: latfi(:)
-real,allocatable :: lonfi(:)
-real,allocatable :: cufi(:)
-real,allocatable :: cvfi(:)
-real,allocatable :: airefi(:)
+real,allocatable :: latfi_glo(:)
+real,allocatable :: lonfi_glo(:)
+real,allocatable :: cufi_glo(:)
+real,allocatable :: cvfi_glo(:)
+real,allocatable :: airefi_glo(:)
+real,allocatable :: boundslonfi_glo(:,:)
+real,allocatable :: boundslatfi_glo(:,:)
+
+! local arrays, on given MPI/OpenMP domain:
+real,allocatable,save :: latfi(:)
+real,allocatable,save :: lonfi(:)
+real,allocatable,save :: cufi(:)
+real,allocatable,save :: cvfi(:)
+real,allocatable,save :: airefi(:)
+real,allocatable,save :: boundslonfi(:,:)
+real,allocatable,save :: boundslatfi(:,:)
+!$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)
 
 pi=2.*asin(1.0)
 
-IF (nlayer.NE.klev) THEN
-  write(*,*) 'STOP in ',trim(modname)
-  write(*,*) 'Problem with dimensions :'
-  write(*,*) 'nlayer     = ',nlayer
-  write(*,*) 'klev   = ',klev
-  abort_message = ''
-  CALL abort_gcm (modname,abort_message,1)
-ENDIF
-
-!IF (ngrid.NE.klon_glo) THEN
-!  write(*,*) 'STOP in ',trim(modname)
-!  write(*,*) 'Problem with dimensions :'
-!  write(*,*) 'ngrid     = ',ngrid
-!  write(*,*) 'klon   = ',klon_glo
-!  abort_message = ''
-!  CALL abort_gcm (modname,abort_message,1)
-!ENDIF
-
-!call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
+! Initialize Physics distibution and parameters and interface with dynamics
+CALL init_physics_distribution(regular_lonlat,4, &
+                                 nbp,ii,jj+1,nlayer,communicator)
+CALL init_interface_dyn_phys
 
 ! init regular global longitude-latitude grid points and boundaries
@@ -122,32 +130,63 @@
 
 ! Generate global arrays on full physics grid
-allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
-latfi(1)=rlatudyn(1)
-lonfi(1)=0.
-cufi(1) = cudyn(1)
-cvfi(1) = cvdyn(1)
+allocate(latfi_glo(klon_glo),lonfi_glo(klon_glo))
+allocate(cufi_glo(klon_glo),cvfi_glo(klon_glo))
+allocate(airefi_glo(klon_glo))
+allocate(boundslonfi_glo(klon_glo,4))
+allocate(boundslatfi_glo(klon_glo,4))
+
+! North pole
+latfi_glo(1)=rlatudyn(1)
+lonfi_glo(1)=0.
+cufi_glo(1) = cudyn(1)
+cvfi_glo(1) = cvdyn(1)
+boundslonfi_glo(1,north_east)=0
+boundslatfi_glo(1,north_east)=PI/2
+boundslonfi_glo(1,north_west)=2*PI
+boundslatfi_glo(1,north_west)=PI/2
+boundslonfi_glo(1,south_west)=2*PI
+boundslatfi_glo(1,south_west)=rlatvdyn(1)
+boundslonfi_glo(1,south_east)=0
+boundslatfi_glo(1,south_east)=rlatvdyn(1)
 DO j=2,jj
   DO i=1,ii
-    latfi((j-2)*ii+1+i)= rlatudyn(j)
-    lonfi((j-2)*ii+1+i)= rlonvdyn(i)
-    cufi((j-2)*ii+1+i) = cudyn((j-1)*(ii+1)+i)
-    cvfi((j-2)*ii+1+i) = cvdyn((j-1)*(ii+1)+i)
+    k=(j-2)*ii+1+i
+    latfi_glo((j-2)*ii+1+i)= rlatudyn(j)
+    lonfi_glo((j-2)*ii+1+i)= rlonvdyn(i)
+    cufi_glo((j-2)*ii+1+i) = cudyn((j-1)*(ii+1)+i)
+    cvfi_glo((j-2)*ii+1+i) = cvdyn((j-1)*(ii+1)+i)
+    boundslonfi_glo(k,north_east)=rlonudyn(i)
+    boundslatfi_glo(k,north_east)=rlatvdyn(j-1)
+    boundslonfi_glo(k,north_west)=rlonudyn(i+1)
+    boundslatfi_glo(k,north_west)=rlatvdyn(j-1)
+    boundslonfi_glo(k,south_west)=rlonudyn(i+1)
+    boundslatfi_glo(k,south_west)=rlatvdyn(j)
+    boundslonfi_glo(k,south_east)=rlonudyn(i)
+    boundslatfi_glo(k,south_east)=rlatvdyn(j)
   ENDDO
 ENDDO
-latfi(klon_glo)= rlatudyn(jj+1)
-lonfi(klon_glo)= 0.
-cufi(klon_glo) = cudyn((ii+1)*jj+1)
-cvfi(klon_glo) = cvdyn((ii+1)*jj-ii)
+! South pole
+latfi_glo(klon_glo)= rlatudyn(jj+1)
+lonfi_glo(klon_glo)= 0.
+cufi_glo(klon_glo) = cudyn((ii+1)*jj+1)
+cvfi_glo(klon_glo) = cvdyn((ii+1)*jj-ii)
+boundslonfi_glo(klon_glo,north_east)= 0
+boundslatfi_glo(klon_glo,north_east)= rlatvdyn(jj)
+boundslonfi_glo(klon_glo,north_west)= 2*PI
+boundslatfi_glo(klon_glo,north_west)= rlatvdyn(jj)
+boundslonfi_glo(klon_glo,south_west)= 2*PI
+boundslatfi_glo(klon_glo,south_west)= -PI/2
+boundslonfi_glo(klon_glo,south_east)= 0
+boundslatfi_glo(klon_glo,south_east)= -Pi/2
 
 ! build airefi(), mesh area on physics grid
-allocate(airefi(klon_glo))
-CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,airedyn,airefi)
+CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,airedyn,airefi_glo)
 ! Poles are single points on physics grid
-airefi(1)=sum(airedyn(1:ii,1))
-airefi(klon_glo)=sum(airedyn(1:ii,jj+1))
+airefi_glo(1)=sum(airedyn(1:ii,1))
+airefi_glo(klon_glo)=sum(airedyn(1:ii,jj+1))
 
 ! Sanity check: do total planet area match between physics and dynamics?
 total_area_dyn=sum(airedyn(1:ii,1:jj+1))
-total_area_phy=sum(airefi(1:klon_glo))
+total_area_phy=sum(airefi_glo(1:klon_glo))
 IF (total_area_dyn/=total_area_phy) THEN
   WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
@@ -164,15 +203,24 @@
 !$OMP PARALLEL 
 ! Now generate local lon/lat/cu/cv/area arrays
-call initcomgeomphy
-
-!!!!$OMP PARALLEL PRIVATE(ibegin,iend) & 
-!!!	!$OMP SHARED(airefi,cufi,cvfi,lonfi,latfi)
+allocate(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
+allocate(airefi(klon_omp))
+allocate(boundslonfi(klon_omp,4))
+allocate(boundslatfi(klon_omp,4))
+!call initcomgeomphy
 
 offset=klon_mpi_begin-1
-cell_area(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end)
-dx(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end)
-dy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end)
-longitude(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end)
-latitude(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end)
+airefi(1:klon_omp)=airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+cufi(1:klon_omp)=cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+cvfi(1:klon_omp)=cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+lonfi(1:klon_omp)=lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+latfi(1:klon_omp)=latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
+boundslonfi(1:klon_omp,:)=boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
+boundslatfi(1:klon_omp,:)=boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
+
+! copy over local grid longitudes and latitudes
+CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
+                   airefi,cufi,cvfi)
+
+call init_dimphy(klon_omp,nlayer) ! Initialize dimphy module
 
 ! copy over preff , ap() and bp() 
Index: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F	(revision 1542)
+++ trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F	(revision 1543)
@@ -15,4 +15,7 @@
 c=======================================================================
 
+      use mod_phys_lmdz_para, only: is_parallel, is_sequential,
+     &                              is_mpi_root, is_omp_root,
+     &                              is_master
       use infotrac, only: infotrac_init, nqtot, tname
       USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
@@ -20,5 +23,4 @@
       USE surfdat_h, ONLY: phisfi, albedodat,
      &                     zmea, zstd, zsig, zgam, zthe
-      USE comgeomfi_h, ONLY: lati, long, area
       use datafile_mod, only: datadir, surfdir
       use ioipsl_getin_p_mod, only: getin_p
@@ -26,26 +28,23 @@
       use phyredem, only: physdem0, physdem1
       use iostart, only: open_startphy
-      use comgeomphy, only: initcomgeomphy
       use slab_ice_h, only:noceanmx
       use filtreg_mod, only: inifilr
+      USE mod_const_mpi, ONLY: COMM_LMDZ
       USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
       USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
-     .			rad,omeg,g,r,pi
+     .                        rad,omeg,g,r,pi
       USE serre_mod, ONLY: alphax
       USE temps_mod, ONLY: day_ini
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-      use inifis_mod, only: inifis
+      use iniphysiq_mod, only: iniphysiq
       implicit none
 
-#include "dimensions.h"
-!#include "dimphys.h"
+      include "dimensions.h"
       integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
-!#include "planete.h"
-#include "paramet.h"
-#include "comgeom2.h"
-!#include "control.h"
-#include "comdissnew.h"
-#include "netcdf.inc"
-!#include "advtrac.h"
+      include "paramet.h"
+      include "comgeom2.h"
+      include "comdissnew.h"
+      include "netcdf.inc"
+
 c=======================================================================
 c   Declarations
@@ -54,5 +53,5 @@
 c Variables dimension du fichier "start_archive"
 c------------------------------------
-      CHARACTER	relief*3
+      CHARACTER        relief*3
 
 
@@ -94,8 +93,8 @@
 c variable physique
 c------------------
-      REAL tsurf(ngridmx)	! surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
-!      REAL co2ice(ngridmx)	! CO2 ice layer
-      REAL emis(ngridmx)	! surface emissivity
+      REAL tsurf(ngridmx)        ! surface temperature
+      REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
+!      REAL co2ice(ngridmx)        ! CO2 ice layer
+      REAL emis(ngridmx)        ! surface emissivity
       real emisread             ! added by RW
       REAL,ALLOCATABLE :: qsurf(:,:)
@@ -103,5 +102,5 @@
 !      REAL rnaturfi(ngridmx)
       real alb(iip1,jjp1),albfi(ngridmx) ! albedos
-      real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
+      real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
       real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
       REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
@@ -113,5 +112,5 @@
 
       REAL rnat(ngridmx)
-      REAL tslab(ngridmx,nsoilmx) ! slab ocean temperature
+      REAL,ALLOCATABLE :: tslab(:,:) ! slab ocean temperature
       REAL pctsrf_sic(ngridmx) ! sea ice cover
       REAL tsea_ice(ngridmx) ! temperature sea_ice
@@ -145,5 +144,5 @@
       character*80      fichnom
       character*250  filestring
-      integer Lmodif,iq,thermo
+      integer Lmodif,iq
       character modif*20
       real z_reel(iip1,jjp1)
@@ -155,9 +154,7 @@
       real val, val2, val3, val4 ! to store temporary variables
       real :: iceith=2000 ! thermal inertia of subterranean ice
-      integer iref,jref
 
       INTEGER :: itau
       
-      INTEGER :: nq,numvanle
       character(len=20) :: txt ! to store some text
       character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
@@ -166,7 +163,4 @@
       real :: profile(llm+1) ! to store an atmospheric profile + surface value
 
-!     added by RW for test
-      real pmean, phi0
-
 !     added by BC for equilibrium temperature startup
       real teque
@@ -193,9 +187,15 @@
                   ! make deliberate choice of these values elsewhere.
 
+      planet_type="generic"
+
 ! initialize "serial/parallel" related stuff
-      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-      call initcomgeomphy
+! (required because we call tabfi() below, before calling iniphysiq)
+      is_sequential=.true.
+      is_parallel=.false.
+      is_mpi_root=.true.
+      is_omp_root=.true.
+      is_master=.true.
+      
 ! Load tracer number and names:
-!      call iniadvtrac(nqtot,numvanle)
       call infotrac_init
 ! allocate arrays
@@ -203,6 +203,12 @@
       allocate(qsurf(ngridmx,nqtot))
       
-      planet_type="generic"
-
+! get value of nsoilmx and allocate corresponding arrays
+      nsoilmx=18 ! default value
+      call getin_p("nsoilmx",nsoilmx)
+      
+      allocate(tsoil(ngridmx,nsoilmx))
+      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
+      allocate(tslab(ngridmx,nsoilmx))
+      
 c=======================================================================
 c   Choice of the start file(s) to use
@@ -284,15 +290,4 @@
       IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
       IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
-      ! ALLOCATE ARRAYS in comsoil_h (if not already done)
-      IF (.not.ALLOCATED(layer))
-     . ALLOCATE(layer(nsoilmx))
-      IF (.not.ALLOCATED(mlayer))
-     . ALLOCATE(mlayer(0:nsoilmx-1))
-      IF (.not.ALLOCATED(inertiedat))
-     . ALLOCATE(inertiedat(ngridmx,nsoilmx))
-      ! ALLOCATE ARRAYS IN comgeomfi_h (done in inifis usually)
-      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx))
-      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx))
-      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx))
 
 c-----------------------------------------------------------------------
@@ -337,5 +332,5 @@
           write(*,*) i,tab_cntrl(i)
         enddo
-	
+        
         ! Lmodif set to 0 to disable modifications possibility in phyeta0                           
         write(*,*) 'Reading file START'
@@ -355,7 +350,7 @@
         do i=1,ngridmx
           albfi(i) = albedodat(i)
-	  do j=1,nsoilmx
+          do j=1,nsoilmx
            ithfi(i,j) = inertiedat(i,j)
-	  enddo
+          enddo
         ! build a surfithfi(:) using 1st layer of ithfi(:), which might
         ! be needed later on if reinitializing soil thermal inertia
@@ -371,5 +366,5 @@
       endif
 c-----------------------------------------------------------------------
-c		Initialisation des constantes dynamique
+c                Initialisation des constantes dynamique
 c-----------------------------------------------------------------------
 
@@ -451,30 +446,11 @@
       idum=0
 
-c Initialisation coordonnees /aires
-c -------------------------------
-! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
-!       rlatu() and rlonv() are given in radians
-      latfi(1)=rlatu(1)
-      lonfi(1)=0.
-      DO j=2,jjm
-         DO i=1,iim
-            latfi((j-2)*iim+1+i)=rlatu(j)
-            lonfi((j-2)*iim+1+i)=rlonv(i)
-         ENDDO
-      ENDDO
-      latfi(ngridmx)=rlatu(jjp1)
-      lonfi(ngridmx)=0.
-       
-      ! build airefi(), mesh area on physics grid
-      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
-      ! Poles are single points on physics grid
-      airefi(1)=sum(aire(1:iim,1))
-      airefi(ngridmx)=sum(aire(1:iim,jjm+1))
-
-! also initialize various physics flags/settings which might be needed
-!    (for instance initracer needs to know about some flags, and/or
-!      'datafile' path may be changed by user)
-      call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
-     &                latfi,lonfi,airefi,rad,g,r,cpp)
+! Initialize the physics
+         CALL iniphysiq(iim,jjm,llm,
+     &                  (jjm-1)*iim+2,comm_lmdz,
+     &                  daysec,day_ini,dtphys,
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
 
 c=======================================================================
@@ -560,7 +536,7 @@
      &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
         write(*,*) "OK, read start_archive file"
-	! copy soil thermal inertia
-	ithfi(:,:)=inertiedat(:,:)
-	
+        ! copy soil thermal inertia
+        ithfi(:,:)=inertiedat(:,:)
+        
         ierr= NF_CLOSE(nid)
 
@@ -727,8 +703,8 @@
           DO j=1,jjp1
              DO i=1,iip1
-                alb(i,j) = alb_bb	! albedo
-		do isoil=1,nsoilmx
-                  ith(i,j,isoil) = ith_bb	! thermal inertia
-		enddo
+                alb(i,j) = alb_bb        ! albedo
+                do isoil=1,nsoilmx
+                  ith(i,j,isoil) = ith_bb        ! thermal inertia
+                enddo
              END DO
           END DO
@@ -867,5 +843,5 @@
          if (yes.eq.'y') then
            write(*,*) 'Value?'
-		read(*,*,iostat=ierr) psea
+                read(*,*,iostat=ierr) psea
              DO i=1,iip1
                DO j=1,jjp1
@@ -874,5 +850,5 @@
                   ENDDO
                ENDDO
-		write(*,*) 'psea=',psea
+                write(*,*) 'psea=',psea
           else
             write(*,*) 'no'
@@ -1467,15 +1443,15 @@
 !       ----------------------------------------------------------------------
 
-	else if (trim(modif) .eq. 'therm_ini_s') then
+        else if (trim(modif) .eq. 'therm_ini_s') then
 !          write(*,*)"surfithfi(1):",surfithfi(1)
-	  do isoil=1,nsoilmx
-	    inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
-	  enddo
+          do isoil=1,nsoilmx
+            inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
+          enddo
           write(*,*)'OK: Soil thermal inertia has been reset to referenc
      &e surface values'
-!	  write(*,*)"inertiedat(1,1):",inertiedat(1,1)
-	  ithfi(:,:)=inertiedat(:,:)
-	 ! recast ithfi() onto ith()
-	 call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+!          write(*,*)"inertiedat(1,1):",inertiedat(1,1)
+          ithfi(:,:)=inertiedat(:,:)
+         ! recast ithfi() onto ith()
+         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
 ! Check:
 !         do i=1,iip1
@@ -1485,5 +1461,5 @@
 !             enddo
 !           enddo
-!	 enddo
+!         enddo
 
 
@@ -1660,5 +1636,5 @@
 
 c=======================================================================
-c	Formats 
+c        Formats 
 c=======================================================================
 
Index: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F	(revision 1542)
+++ trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F	(revision 1543)
@@ -21,25 +21,28 @@
       use infotrac, only: infotrac_init, nqtot, tname
       USE comsoil_h
-      USE comgeomfi_h, ONLY: lati, long, area
+      
+!      USE comgeomfi_h, ONLY: lati, long, area
 !      use control_mod
-      use comgeomphy, only: initcomgeomphy
+!      use comgeomphy, only: initcomgeomphy
       use slab_ice_h, only: noceanmx
 ! to use  'getin'
       USE ioipsl_getincom
       USE planete_mod, only: year_day
+      USE mod_const_mpi, ONLY: COMM_LMDZ
       USE control_mod, only: planet_type
       USE callkeys_mod, ONLY: ok_slab_ocean
       use filtreg_mod, only: inifilr
       USE comvert_mod, ONLY: ap,bp
-      USE comconst_mod, ONLY: cpp,g
+      USE comconst_mod, ONLY: daysec,dtphys,rad,g,r,cpp
       USE logic_mod, ONLY: grireg
       USE temps_mod, ONLY: day_ini
+      USE iniphysiq_mod, ONLY: iniphysiq
       implicit none
 
-#include "dimensions.h"
+      include "dimensions.h"
       integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
-#include "paramet.h"
-#include "comdissip.h"
-#include "comgeom.h"
+      include "paramet.h"
+      include "comdissip.h"
+      include "comgeom.h"
 !#include "control.h"
 
@@ -47,5 +50,5 @@
 !#include "planete.h"
 !#include"advtrac.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 c-----------------------------------------------------------------------
 c   Declarations
@@ -68,7 +71,7 @@
 c Variable Physiques (grille physique)
 c ------------------------------------
-      REAL tsurf(ngridmx)	! Surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
-      REAL co2ice(ngridmx)	! CO2 ice layer
+      REAL tsurf(ngridmx)        ! Surface temperature
+      REAL,ALLOCATABLE :: tsoil(:,:) ! Soil temperature
+      REAL co2ice(ngridmx)        ! CO2 ice layer
       REAL q2(ngridmx,llm+1)
       REAL,ALLOCATABLE :: qsurf(:,:)
@@ -93,6 +96,6 @@
       REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
       REAL tsurfS(ip1jmp1)
-      REAL tsoilS(ip1jmp1,nsoilmx)
-      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
+      REAL,ALLOCATABLE :: tsoilS(:,:)
+      REAL,ALLOCATABLE :: ithS(:,:) ! Soil Thermal Inertia
       REAL co2iceS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
@@ -144,14 +147,4 @@
       grireg   = .TRUE.
 
-! initialize "serial/parallel" related stuff
-      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-      call initcomgeomphy
-
-      ! ALLOCATE ARRAYS IN comgeomfi_h (usually done in inifis)
-      ! this must be here for start2archive to work
-      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx))
-      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx))
-      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx))
-
       planet_type="generic"
 
@@ -160,5 +153,4 @@
 c=======================================================================
 ! Load tracer number and names:
-!      call iniadvtrac(nqtot,numvanle)
       call infotrac_init
 
@@ -168,5 +160,5 @@
       allocate(qsurfS(ip1jmp1,nqtot))
 ! other array allocations:
-      call ini_comsoil_h(ngridmx)
+!      call ini_comsoil_h(ngridmx) ! done via iniphysiq
 
       fichnom = 'start.nc'
@@ -198,5 +190,45 @@
 
       ierr = NF_CLOSE(nid1)
+
+! Get value of the "subsurface_layers" dimension from physics start file
+      fichnom = 'startfi.nc'
+      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
+       IF (ierr.NE.NF_NOERR) THEN
+         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
+        CALL ABORT
+       ENDIF
+      ierr = NF_INQ_DIMID(nid1,"subsurface_layers",varid)
+      IF (ierr .NE. NF_NOERR) THEN
+       PRINT*, "start2archive: No subsurface_layers dimension!!"
+       CALL abort
+      ENDIF
+      ierr = NF_INQ_DIMLEN(nid1,varid,nsoilmx)
+      IF (ierr .NE. NF_NOERR) THEN
+       PRINT*, "start2archive: Failed reading subsurface_layers value!!"
+       CALL abort
+      ENDIF
+      ierr = NF_CLOSE(nid1)
       
+      ! allocate arrays of nsoilmx size
+      allocate(tsoil(ngridmx,nsoilmx))
+      allocate(tsoilS(ip1jmp1,nsoilmx))
+      allocate(ithS(ip1jmp1,nsoilmx))
+
+c-----------------------------------------------------------------------
+c   Initialisations 
+c-----------------------------------------------------------------------
+
+      CALL defrun_new(99, .FALSE. )
+      call iniconst
+      call inigeom
+      call inifilr
+
+! Initialize the physics
+         CALL iniphysiq(iim,jjm,llm,
+     &                  (jjm-1)*iim+2,comm_lmdz,
+     &                  daysec,day_ini,dtphys,
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
 
       fichnom = 'startfi.nc'
@@ -256,12 +288,4 @@
 c *****************************************************************
 
-c-----------------------------------------------------------------------
-c   Initialisations 
-c-----------------------------------------------------------------------
-
-      CALL defrun_new(99, .FALSE. )
-      call iniconst
-      call inigeom
-      call inifilr
       CALL pression(ip1jmp1, ap, bp, ps, p3d)
       call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
