Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 2334)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 2336)
@@ -1,28 +1,46 @@
 PROGRAM ce0l
-  !
-  ! Purpose: Calls etat0, creates initial states and limit_netcdf
-  !
-  ! interbar=.T. for barycentric interpolation inter_barxy
-  ! extrap  =.T. for data extrapolation, like for the SSTs when file does not
-  !                  contain ocean points only.
-  ! oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
-  ! masque is created in etat0, passed to limit to ensure consistancy.
-
-  USE control_mod, only: DAY_STEP, DAYREF, NSPLIT_PHYS
-  USE etat0dyn, only: etat0dyn_netcdf
-  USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
+!
+!-------------------------------------------------------------------------------
+! Purpose: Initial states and boundary conditions files creation:
+!     * start.nc    for dynamics    (using etat0dyn     routine)
+!     * startphy.nc for physics     (using etat0phys    routine)
+!     * limit.nc    for forced runs (using limit_netcdf routine)
+!-------------------------------------------------------------------------------
+! Notes:
+!     * extrap=.T. (default) for data extrapolation, like for the SSTs when file
+!                   does contain ocean points only.
+!     * "masque" can be:
+!       - read from file "o2a.nc"          (for coupled runs).
+!       - created in etat0phys or etat0dyn (for forced  runs).
+!     It is then passed to limit_netcdf to ensure consistancy.
+!-------------------------------------------------------------------------------
   USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo
-
-  USE etat0phys, only: etat0phys_netcdf
-  USE dimphy, only: KLON
-  USE infotrac, only: TYPE_TRAC, infotrac_init
+  USE control_mod,    ONLY: day_step, dayref, nsplit_phys
+  USE etat0dyn,       ONLY: etat0dyn_netcdf
+  USE etat0phys,      ONLY: etat0phys_netcdf
+  USE limit,          ONLY: limit_netcdf
+  USE netcdf,         ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
+  USE infotrac,       ONLY: type_trac, infotrac_init
+  USE dimphy,         ONLY: klon
   USE test_disvert_m, ONLY: test_disvert
+  USE filtreg_mod,    ONLY: inifilr
+#ifdef inca
+  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
+#endif
+#ifdef CPP_PARA
+  USE mod_const_mpi,  ONLY: init_const_mpi
+  USE parallel_lmdz,  ONLY: init_parallel, mpi_rank, omp_rank, mpi_size
+  USE bands,          ONLY: read_distrib, distrib_phys
+  USE mod_hallo,      ONLY: init_mod_hallo
+  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
+#endif
 
   IMPLICIT NONE
 
-  ! Local variables:
+!-------------------------------------------------------------------------------
+! Local variables:
   include "dimensions.h"
   include "paramet.h"
-  include "comgeom.h"
+  include "comgeom2.h"
   include "comconst.h"
   include "comvert.h"
@@ -30,21 +48,20 @@
   include "temps.h"
   include "logic.h"
-  REAL               :: masque(iip1, jjp1)             !--- CONTINENTAL MASK
-  REAL               :: phis  (iip1, jjp1)             !--- GROUND GEOPOTENTIAL
+  REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
+  REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
   CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
   LOGICAL            :: use_filtre_fft
-  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
-
-  !--- Local variables for ocean mask reading:
+  LOGICAL, PARAMETER :: extrap=.FALSE.
+
+!--- Local variables for ocean mask reading:
   INTEGER            :: nid_o2a, iml_omask, jml_omask, j
   INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
-  REAL, ALLOCATABLE  :: lon_omask(:, :), dlon_omask(:), ocemask(:, :)
-  REAL, ALLOCATABLE  :: lat_omask(:, :), dlat_omask(:), ocetmp (:, :)
+  REAL, ALLOCATABLE  :: lon_omask(:,:), dlon_omask(:), ocemask(:,:)
+  REAL, ALLOCATABLE  :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:)
   REAL               :: date, lev(1)
-
-  !----------------------------------------------------------------------
+!-------------------------------------------------------------------------------
   modname="ce0l"
 
-  !--- Constants
+!--- Constants
   pi     = 4. * ATAN(1.)
   rad    = 6371229.
@@ -59,145 +76,146 @@
 
   CALL conf_gcm( 99, .TRUE. )
-
   dtvr = daysec/REAL(day_step)
-  WRITE(lunout, *)'dtvr', dtvr
-
+  WRITE(lunout,*)'dtvr',dtvr
   CALL iniconst()
   CALL inigeom()
 
+!--- Calendar choice
 #ifdef CPP_IOIPSL
   calnd='gregorian'
   SELECT CASE(calend)
-  CASE('earth_360d')
-     CALL ioconf_calendar('360d')
-     calnd='with 360 days/year'
-  CASE('earth_365d')
-     CALL ioconf_calendar('noleap')
-     calnd='with no leap year'
-  CASE('earth_366d')
-     CALL ioconf_calendar('366d')
-     calnd='with leap years only'
-  CASE('gregorian')
-     CALL ioconf_calendar('gregorian')
-  CASE('standard')
-     CALL ioconf_calendar('gregorian')
-  CASE('julian')
-     CALL ioconf_calendar('julian')
-     calnd='julian'
-  CASE('proleptic_gregorian')
-     CALL ioconf_calendar('gregorian')
-     !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
-  CASE DEFAULT
-     CALL abort_gcm('ce0l', 'Bad choice for calendar', 1)
+    CASE('earth_360d');CALL ioconf_calendar('360d');   calnd='with 360 days/year'
+    CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year'
+    CASE('earth_366d');CALL ioconf_calendar('366d');   calnd='with leap years only'
+    CASE('gregorian'); CALL ioconf_calendar('gregorian')
+    CASE('standard');  CALL ioconf_calendar('gregorian')
+    CASE('julian');    CALL ioconf_calendar('julian'); calnd='julian'
+    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
+  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
+    CASE DEFAULT
+      CALL abort_gcm('ce0l','Bad choice for calendar',1)
   END SELECT
-  WRITE(lunout, *)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
-#endif
-
-  use_filtre_fft=.FALSE.
-  CALL getin('use_filtre_fft', use_filtre_fft)
-  IF(use_filtre_fft) THEN
-     WRITE(lunout, *)"FFT filter not available for sequential dynamics."
-     WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used."
-  ENDIF
-
-  !--- LAND MASK. TWO CASES:
-  !   1) read from ocean model    file "o2a.nc"    (coupled runs)
-  !   2) computed from topography file "Relief.nc" (masque(:, :)=-99999.)
-  ! Coupled simulations (case 1) use the ocean model mask to compute the
-  ! weights to ensure ocean fractions are the same for atmosphere and ocean.
-
-  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
-     WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file found'
-     WRITE(lunout, *)'Forced run.'
-     masque(:, :)=-99999.
-  ELSE
-     iret=NF90_CLOSE(nid_o2a)
-     WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found'
-     WRITE(lunout, *)'Coupled run.'
-     CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
-     IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
-        WRITE(lunout, *)'Mismatching dimensions for ocean mask'
-        WRITE(lunout, *)'iim  = ', iim , ' iml_omask = ', iml_omask
-        WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask
-        CALL abort_gcm(modname, '', 1)
-     END IF
-     ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim ))
-     ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1))
-     CALL flinopen("o2a.nc", .FALSE., iim, jjp1, llm_tmp, lon_omask, &
-          lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
-     CALL flinget(fid, "OceMask",    iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp)
-     CALL flinclo(fid)
-     dlon_omask(1:iim ) = lon_omask(1:iim, 1)
-     dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1)
-     ocemask = ocetmp
-     IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
-        DO j=1, jjp1
-           ocemask(:, j) = ocetmp(:, jjp1-j+1)
-        END DO
-     END IF
-     DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask)
-     IF(prt_level>=1) THEN
-        WRITE(fmt, "(i4, 'i1)')")iim 
-        fmt='('//ADJUSTL(fmt)
-        WRITE(lunout, *)'OCEAN MASK :'
-        WRITE(lunout, fmt) NINT(ocemask)
-     END IF
-     masque(1:iim, :)=1.-ocemask(:, :)
-     masque(iip1 , :)=masque(1, :)
-     DEALLOCATE(ocemask)
-  END IF
-  phis(:, :)=-99999.
-
-  CALL Init_Phys_lmdz(iim, jjp1, llm, 1, (/(jjm-1)*iim+2/))
-  WRITE(lunout, *)'---> klon=', klon
-
-  call infotrac_init
-  CALL iniphysiq(iim, jjm, llm, daysec, dayref, dtphys / nsplit_phys, rlatu, &
-       rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys)
-
-  IF(pressure_exner) CALL test_disvert
-
+  WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
+#endif
+
+#ifdef CPP_PARA
+!--- Physical grid + parallel initializations
+  CALL init_const_mpi()
+  CALL init_parallel()
+  CALL read_distrib()
+  CALL init_mod_hallo()
+  CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
+  CALL init_interface_dyn_phys()
+#else
+  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+#endif
+  WRITE(lunout,*)'---> klon=',klon
+
+!--- Tracers initializations
   IF (type_trac == 'inca') THEN
 #ifdef INCA
-     CALL init_const_lmdz(nbtr, anneeref, dayref, iphysiq, day_step, nday)
-     CALL init_inca_para(iim, jjm+1, klon, 1, klon_mpi_para_nb, 0)
-     WRITE(lunout, *)'nbtr =' , nbtr 
-#endif
-  END IF
+    CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,&
+                         nbsrf,is_oce,is_sic,is_ter,is_lic,calend)
+    CALL init_inca_para(iim,jjp1,llm,klon_glo,mpi_size,distrib_phys,&
+                        COMM_LMDZ)
+    WRITE(lunout,*)'nbtr =' , nbtr 
+#endif
+  END IF
+  CALL infotrac_init()
+
+  CALL inifilr()
+  CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
+                 rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+  IF(pressure_exner) CALL test_disvert
+
+#ifdef CPP_PARA
+  IF (mpi_rank==0.AND.omp_rank==0) THEN
+#endif
+  use_filtre_fft=.FALSE.
+  CALL getin('use_filtre_fft',use_filtre_fft)
+  IF(use_filtre_fft) THEN
+     WRITE(lunout,*)"FFT filter not available for sequential dynamics."
+     WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used."
+  ENDIF
+
+!--- LAND MASK. TWO CASES:
+!   1) read from ocean model    file "o2a.nc"    (coupled runs)
+!   2) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
+! Coupled simulations (case 1) use the ocean model mask to compute the
+! weights to ensure ocean fractions are the same for atmosphere and ocean.
+!*******************************************************************************
+  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
+    WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Forced run.'
+    masque(:,:)=-99999.
+  ELSE
+    iret=NF90_CLOSE(nid_o2a)
+    WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Coupled run.'
+    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
+    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
+      WRITE(lunout,*)'Mismatching dimensions for ocean mask'
+      WRITE(lunout,*)'iim  = ',iim ,' iml_omask = ',iml_omask
+      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
+      CALL abort_gcm(modname,'',1)
+    END IF
+    ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim ))
+    ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1))
+    CALL flinopen("o2a.nc", .FALSE.,iim,jjp1,llm_tmp,lon_omask,lat_omask,lev,  &
+     &                                       ttm_tmp,itaul,date,dt,fid)
+    CALL flinget(fid, "OceMask",    iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp)
+    CALL flinclo(fid)
+    dlon_omask(1:iim ) = lon_omask(1:iim,1)
+    dlat_omask(1:jjp1) = lat_omask(1,1:jjp1)
+    ocemask = ocetmp
+    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
+      DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO
+    END IF
+    DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
+    IF(prt_level>=1) THEN
+      WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt)
+      WRITE(lunout,*)'OCEAN MASK :'
+      WRITE(lunout,fmt) NINT(ocemask)
+    END IF
+    masque(1:iim,:)=1.-ocemask(:,:)
+    masque(iip1 ,:)=masque(1,:)
+    DEALLOCATE(ocemask)
+  END IF
+  phis(:,:)=-99999.
+
   IF(ok_etat0) THEN
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, *) '  ***  etat0phy_netcdf ***  '
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) ' interbar = ', interbar
-     CALL etat0phys_netcdf(interbar, masque, phis)
-  END IF
-
-  IF(ok_etat0) THEN
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, *) '  ***  etat0dyn_netcdf ***  '
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) ' interbar = ', interbar
-     CALL etat0dyn_netcdf(interbar, masque, phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0phy_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    CALL etat0phys_netcdf(masque,phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0dyn_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    CALL etat0dyn_netcdf(masque,phis)
   END IF
 
   IF(ok_limit) THEN
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) '  *********************  '
-     WRITE(lunout, *) '  ***  Limit_netcdf ***  '
-     WRITE(lunout, *) '  *********************  '
-     WRITE(lunout, '(//)')
-     CALL limit_netcdf(interbar, extrap, oldice, masque)
-  END IF
-
-  WRITE(lunout, '(//)')
-  WRITE(lunout, *) '  ***************************  '
-  WRITE(lunout, *) '  ***  grilles_gcm_netcdf ***  '
-  WRITE(lunout, *) '  ***************************  '
-  WRITE(lunout, '(//)')
-  CALL grilles_gcm_netcdf_sub(masque, phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,*) '  ***  Limit_netcdf ***  '
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,'(//)')
+    CALL limit_netcdf(masque,phis,extrap)
+  END IF
+
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
+
+#ifdef CPP_PARA
+  END IF
+#endif
 
 END PROGRAM ce0l
+!
+!-------------------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 2334)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 2336)
@@ -12,5 +12,5 @@
 !  routine (to be called after restget):
 !    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
-!                          champ, lon_in2, lat_in2, ibar)
+!                          champ, lon_in2, lat_in2)
 !
 !    *  Variables should have the following names in the NetCDF files:
@@ -36,7 +36,5 @@
   USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
   USE assert_eq_m,    ONLY: assert_eq
-#ifdef CPP_PHYS
   USE indice_sol_mod, ONLY: epsfra
-#endif
   IMPLICIT NONE
 
@@ -54,11 +52,7 @@
   include "serre.h"
   REAL, SAVE :: deg2rad
-#ifndef CPP_PHYS
-  REAL, SAVE :: epsfra= 1.E-5
-#endif
   INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
   REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
   CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
-  CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'
 
 CONTAINS
@@ -66,5 +60,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0dyn_netcdf(ib, masque, phis)
+SUBROUTINE etat0dyn_netcdf(masque, phis)
 !
 !-------------------------------------------------------------------------------
@@ -73,20 +67,17 @@
 ! Notes:  1) This routine is designed to work for Earth
 !         2) If masque(:,:)/=-99999., masque and phis are already known.
-!         Otherwise: read it from ocean model file (coupled run) or compute it.
+!         Otherwise: compute it.
 !-------------------------------------------------------------------------------
   USE control_mod
-#ifdef CPP_PHYS
   USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
   USE regr_pr_o3_m,   ONLY: regr_pr_o3
   USE press_coefoz_m, ONLY: press_coefoz
-#endif
   USE exner_hyb_m,    ONLY: exner_hyb
   USE exner_milieu_m, ONLY: exner_milieu
-  USE infotrac, only: NQTOT, TNAME
+  USE infotrac,       ONLY: nqtot, tname
   USE filtreg_mod
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
-  LOGICAL, INTENT(IN)    :: ib                  !--- Barycentric interpolation
   REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
   REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
@@ -111,12 +102,7 @@
   deg2rad = pi/180.0
 
-! Initializations for tracers and filter
-!*******************************************************************************
-  CALL inifilr()
-
 ! Compute ground geopotential and possibly the mask.
 !*******************************************************************************
   masque_tmp(:,:)=masque(:,:)
-  CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)
   WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
   IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
@@ -132,5 +118,5 @@
 ! Compute psol AND tsol, knowing phis.
 !*******************************************************************************
-  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib, phis, psol)
+  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
 
 ! Mid-levels pressure computation
@@ -147,10 +133,10 @@
 !*******************************************************************************
   uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
-  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
+  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
   CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
- &                           rlonu,rlatu(:jjm),ib)
-  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
+ &                           rlonu,rlatu(:jjm))
+  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
   tpot(:,:,:)=t3d(:,:,:)
-  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
+  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
 
   WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
@@ -165,9 +151,8 @@
 !  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
   qd (:,:,:) = 0.0
-  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
+  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
   ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
   CALL flinclo(fid_dyn)
 
-#ifdef CPP_PHYS
 ! Parameterization of ozone chemistry:
 !*******************************************************************************
@@ -180,20 +165,5 @@
     q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
   END IF
-
-#endif
   q3d(iip1,:,:,:)=q3d(1,:,:,:)
-
-! Intermediate computation
-!*******************************************************************************
-  CALL massdair(p3d,masse)
-  WRITE(lunout,*)' ALPHAX ',alphax
-  DO l=1,llm
-    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
-    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
-    xpn=SUM(xppn)/apoln
-    xps=SUM(xpps)/apols
-    masse(:,1   ,l)=xpn
-    masse(:,jjp1,l)=xps
-  END DO
 
 ! Writing
@@ -215,10 +185,18 @@
   CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
   WRITE(lunout,*)'sortie geopot'
-  CALL caldyn0( itau, uvent, vvent, tpot, psol, pk, phis,               &
-                phi,  w, time+iday-dayref)
+  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
+                phi,  w, pbaru, pbarv, time+iday-dayref)
   WRITE(lunout,*)'sortie caldyn0'     
+#ifdef CPP_PARA
+  CALL dynredem0_loc( "start.nc", dayref, phis)
+#else
   CALL dynredem0( "start.nc", dayref, phis)
+#endif
   WRITE(lunout,*)'sortie dynredem0'
+#ifdef CPP_PARA
+  CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+#else
   CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+#endif
   WRITE(lunout,*)'sortie dynredem1' 
   CALL histclo()
@@ -232,5 +210,5 @@
 !
 SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
-                        champ, lon_in2, lat_in2, ibar)
+                        champ, lon_in2, lat_in2)
 !-------------------------------------------------------------------------------
   IMPLICIT NONE
@@ -252,5 +230,4 @@
   REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
   REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
-  LOGICAL,          INTENT(IN)    :: ibar
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -287,5 +264,5 @@
 !--- INTERPOLATE 3D FIELD IF NEEDED
   IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
-                                                  lat_in2,pls,champ,ibar)
+                                                  lat_in2,pls,champ)
 
 !--- COMPUTE THE REQUIRED FILED
@@ -317,212 +294,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
-!
-!-------------------------------------------------------------------------------
-  USE conf_dat_m, ONLY: conf_dat2d
-  IMPLICIT NONE
-!===============================================================================
-! Purpose:  Compute "phis" just like it would be in start_init_orog.
-!===============================================================================
-! Arguments:
-  REAL, INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
-  REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_orog0"
-  CHARACTER(LEN=256) :: title="RELIEF"
-  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
-  REAL               :: lev(1), date, dt
-  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
-  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
-!-------------------------------------------------------------------------------
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
-  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
-  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
-  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
-
-!--- HIGH RESOLUTION OROGRAPHY
-  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
-
-  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
-  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
-                lev, ttm_tmp, itau, date, dt, fid)
-  ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
-  CALL flinclo(fid)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
-  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
-  DEALLOCATE(lon_ini,lat_ini)
-
-!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
-  WRITE(lunout,*)
-  WRITE(lunout,*)'*** Compute surface geopotential ***'
-
-!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
-  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
-  phis = phis * 9.81
-  phis(iml,:) = phis(1,:)
-  DEALLOCATE(relief_hi,lon_rad,lat_rad)
-
-END SUBROUTINE start_init_orog0
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
-!
-!===============================================================================
-! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
-!          without any call to physics subroutines.
-!===============================================================================
-  IMPLICIT NONE 
-!-------------------------------------------------------------------------------
-! Arguments:
-  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
-  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
-  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
-  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
-  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname="grid_noro0"
-  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
-  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
-  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
-  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
-  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
-  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
-  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
-  LOGICAL :: masque_lu
-  INTEGER :: i, ii, imdp, imar, iext
-  INTEGER :: j, jj, jmdp, jmar, nn
-  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
-  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
-!-------------------------------------------------------------------------------
-  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
-  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
-  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
-  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
-  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
-  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
-  iext=imdp/10
-  xpi = ACOS(-1.)
-  rad = 6371229.
-
-!--- ARE WE USING A READ MASK ?
-  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
-  WRITE(lunout,*)'Masque lu: ',masque_lu
-
-!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
-  ALLOCATE(xusn(imdp+2*iext))
-  xusn(1     +iext:imdp  +iext)=xd(:)
-  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
-  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
-
-  ALLOCATE(yusn(jmdp+2))
-  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
-  yusn(2:jmdp+1)=yd(:)
-  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
-
-  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
-  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
-  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
-  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
-  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
-  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
-  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
-  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
-
-!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
-  ALLOCATE(a(imar+1),b(imar+1))
-  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
-  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
-  a(1)=x(1)-(x(2)-x(1))/2.0
-  a(2:imar+1)= b(1:imar)
-
-  ALLOCATE(c(jmar),d(jmar))
-  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
-  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
-  c(1)=y(1)-(y(2)-y(1))/2.0
-  c(2:jmar)=d(1:jmar-1)
-
-!--- INITIALIZATIONS:
-  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
-  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
-
-!--- SUMMATION OVER GRIDPOINT AREA
-  zleny=xpi/REAL(jmdp)*rad
-  xincr=xpi/REAL(jmdp)/2.
-  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
-  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
-  DO ii = 1, imar+1
-    DO jj = 1, jmar
-      DO j = 2,jmdp+1 
-        zlenx  =zleny  *COS(yusn(j))
-        zbordnor=(xincr+c(jj)-yusn(j))*rad
-        zbordsud=(xincr-d(jj)+yusn(j))*rad
-        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
-        IF(weighy/=0) THEN
-          DO i = 2, imdp+2*iext-1
-            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
-            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
-            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
-            IF(weighx/=0)THEN
-              num_tot(ii,jj)=num_tot(ii,jj)+1.0
-              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
-              weight(ii,jj)=weight(ii,jj)+weighx*weighy
-              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
-            END IF
-          END DO
-        END IF
-      END DO
-    END DO
-  END DO
-
-!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
-  IF(.NOT.masque_lu) THEN
-    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
-  END IF
-  nn=COUNT(weight(:,1:jmar-1)==0.0)
-  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
-  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
-
-!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
-  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
-  WHERE(mask>=0.1) mask_tmp = 1.
-  WHERE(weight(:,:)/=0.0)
-    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
-    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
-  END WHERE
-  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
-
-!--- Values at poles
-  zphi(imar+1,:)=zphi(1,:)
-
-  zweinor=SUM(weight(1:imar,   1),DIM=1)
-  zweisud=SUM(weight(1:imar,jmar),DIM=1)
-  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
-  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
-  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
-
-END SUBROUTINE grid_noro0
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol)
+SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
 !
 !-------------------------------------------------------------------------------
@@ -534,5 +304,4 @@
   REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
   REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
-  LOGICAL, INTENT(IN)  :: ibar
   REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
   REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
@@ -606,7 +375,7 @@
     ALLOCATE(field(iml,jml))
     CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
-    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)
-    CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,       &
-                              lon_in, lat_in, lon_in2, lat_in2, field)
+    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
+    CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana,              &
+                                        lon_in, lat_in, lon_in2, lat_in2, field)
   ELSE IF(SIZE(field)/=SIZE(z)) THEN
     msg='The '//TRIM(msg)//' field we have does not have the right size'
@@ -625,5 +394,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar)
+SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
 !
 !-------------------------------------------------------------------------------
@@ -639,5 +408,4 @@
   REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
   REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
-  LOGICAL, INTENT(IN)  :: ibar
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -667,5 +435,5 @@
   ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
   CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
-                       lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
+                       lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
   DEALLOCATE(lon_ini, lat_ini)
 
@@ -673,5 +441,5 @@
   ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
   DO il = 1,llm_dyn
-    CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),     &
+    CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),          &
                           lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
   END DO
@@ -708,14 +476,13 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
 !
 !-------------------------------------------------------------------------------
   USE inter_barxy_m, ONLY: inter_barxy
-  USE grid_atob_m,   ONLY: grille_m
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*), INTENT(IN)  :: nam
-  LOGICAL,          INTENT(IN)  :: ibar, ibeg
+  LOGICAL,          INTENT(IN)  :: ibeg
   REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
   REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
@@ -735,14 +502,10 @@
   j2=SIZE(lat2)
   ALLOCATE(vtmp(i1-1,j1))
-  IF(ibar) THEN
-    IF(ibeg.AND.prt_level>1) THEN
-      WRITE(lunout,*)"---------------------------------------------------------"
-      WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
-      WRITE(lunout,*)"---------------------------------------------------------"
-    END IF
-    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
-  ELSE
-    CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
-  END IF
+  IF(ibeg.AND.prt_level>1) THEN
+    WRITE(lunout,*)"---------------------------------------------------------"
+    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+    WRITE(lunout,*)"---------------------------------------------------------"
+  END IF
+  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
 
Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 2334)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 2336)
@@ -36,5 +36,6 @@
   USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
   USE assert_eq_m,        ONLY: assert_eq
-  USE dimphy
+  USE dimphy,             ONLY: klon
+  USE conf_dat_m,         ONLY: conf_dat2d
   USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
     rlon, solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
@@ -51,18 +52,15 @@
   include "paramet.h"
   include "comgeom2.h"
-  include "comvert.h"
   include "comconst.h"
   include "dimsoil.h"
   include "temps.h"
-  include "comdissnew.h"
-  include "serre.h"
   include "clesphys.h"
   REAL, SAVE :: deg2rad
   REAL, SAVE, ALLOCATABLE :: tsol(:)
-!  REAL, SAVE, ALLOCATABLE :: rugo(:,:)  ! ??? COMPUTED BUT NOT USED ???
   INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
   REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
-  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc"
-  CHARACTER(LEN=256), PARAMETER :: title="RELIEF"
+  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
+  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
+  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
 
 
@@ -72,10 +70,12 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0phys_netcdf(ib, masque, phis)
+SUBROUTINE etat0phys_netcdf(masque, phis)
 !
 !-------------------------------------------------------------------------------
 ! Purpose: Creates initial states
 !-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
+! Notes:  1) This routine is designed to work for Earth
+!         2) If masque(:,:)/=-99999., masque and phis are already known.
+!         Otherwise: compute it.
 !-------------------------------------------------------------------------------
   USE control_mod
@@ -84,13 +84,9 @@
   USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
   USE indice_sol_mod
-  USE conf_phys_m,    ONLY: conf_phys
-  USE exner_hyb_m,    ONLY: exner_hyb
-  USE exner_milieu_m, ONLY: exner_milieu
-  USE test_disvert_m, ONLY: test_disvert
-  USE grid_atob_m,    ONLY: grille_m
+  USE conf_phys_m, ONLY: conf_phys
+  USE init_ssrf_m, ONLY: start_init_subsurf
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
-  LOGICAL, INTENT(IN)    :: ib          !--- Barycentric interpolation
   REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
   REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
@@ -99,18 +95,10 @@
   CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
   INTEGER            :: i, j, l, ji, iml, jml
-  REAL               :: phystep
+  LOGICAL            :: read_mask
+  REAL               :: phystep, dummy
   REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp
   REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
   REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
   REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
-
-!--- Local variables for sea-ice reading:
-  LOGICAL           :: read_mask
-  INTEGER           :: iml_lic, jml_lic, isst(klon-2)
-  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1)
-  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)
-  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)
-  REAL              :: date, lev(1), dummy
-  REAL              :: flic_tmp(SIZE(masque,1),SIZE(masque,2))
 
 !--- Arguments for conf_phys
@@ -132,5 +120,5 @@
   jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
 
-! Grid construction and miscellanous initializations.
+! Physics configuration
 !*******************************************************************************
   CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
@@ -145,18 +133,8 @@
                    read_climoz,                                         &
                    alp_offset)
-
   CALL phys_state_var_init(read_climoz)
 
-  co2_ppm0 = co2_ppm  !--- Initial atmospheric CO2 conc. from .def file
-
-  rlat(1) = pi/2.
-  DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
-  rlat(klon) = - pi/2.
-  rlat(:)=rlat(:)*(180.0/pi)
-
-  rlon(1) = 0.0
-  DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
-  rlon(klon) = 0.0
-  rlon(:)=rlon(:)*(180.0/pi)
+!--- Initial atmospheric CO2 conc. from .def file
+  co2_ppm0 = co2_ppm
 
 ! Compute ground geopotential, sub-cells quantities and possibly the mask.
@@ -174,9 +152,9 @@
     WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
   END IF
- CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
+  CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
 
 ! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
 !*******************************************************************************
-  CALL start_init_phys(rlonv, rlatu, rlonu, rlatv, ib, phis)
+  CALL start_init_phys(rlonu, rlatv, phis)
 
 ! Some initializations.
@@ -188,56 +166,7 @@
     CALL regr_lat_time_climoz(read_climoz)
 
-! Sub-surfaces initialization
-!*******************************************************************************
-!--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
-  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
-  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
-  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
-  ALLOCATE( fraclic(iml_lic,jml_lic))
-  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
- &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
-  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
-  CALL flinclo(fid)
-  WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
-  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees
-  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
-  dlon_lic(:)=lon_lic(:,1)
-  dlat_lic(:)=lat_lic(1,:)
-  CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) )
-  flic_tmp(iml,:)=flic_tmp(1,:)
-
-!--- To the physical grid
-  pctsrf(:,:) = 0.
-  CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic))
-
-!--- Adequation with soil/sea mask
-  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
-  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
-  pctsrf(:,is_ter)=zmasq(:)
-  DO ji=1,klon
-    IF(zmasq(ji)>EPSFRA) THEN 
-      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
-        pctsrf(ji,is_lic)=zmasq(ji)
-        pctsrf(ji,is_ter)=0.
-      ELSE
-        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
-        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
-          pctsrf(ji,is_ter)=0.
-          pctsrf(ji,is_lic)=zmasq(ji)
-        END IF 
-      END IF 
-    END IF 
-  END DO 
-
-!--- Sub-surface ocean and sea ice (sea ice set to zero for start).
-  pctsrf(:,is_oce)=(1.-zmasq(:))
-  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
-  IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:)
-  isst=0
-  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
-
-!--- It is checked that the sub-surfaces sum is equal to 1.
-  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
-  IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
+! Sub-surfaces initialization.
+!*******************************************************************************
+  CALL start_init_subsurf(read_mask)
 
 ! Write physical initial state
@@ -321,6 +250,4 @@
 !   described in LOTT & MILLER (1997) and LOTT(1999).
 !===============================================================================
-  USE conf_dat_m,  ONLY: conf_dat2d
-!  USE grid_atob_m, ONLY: rugsoro
   USE grid_noro_m, ONLY: grid_noro
   IMPLICIT NONE
@@ -331,6 +258,5 @@
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_orog"
-  CHARACTER(LEN=256) :: title="RELIEF"
+  CHARACTER(LEN=256) :: modname
   INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
   REAL               :: lev(1), date, dt
@@ -340,4 +266,5 @@
   REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
 !-------------------------------------------------------------------------------
+  modname="start_init_orog"
   iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
   jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
@@ -350,5 +277,5 @@
                 lev, ttm_tmp, itau, date, dt, fid)
   ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
+  CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
   CALL flinclo(fid)
 
@@ -360,5 +287,5 @@
 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
   ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
+  CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
   DEALLOCATE(lon_ini,lat_ini)
 
@@ -378,12 +305,4 @@
   phis = phis * 9.81
   phis(iml,:) = phis(1,:)
-
-!--- COMPUTE SURFACE ROUGHNESS
-!  WRITE(lunout,*)
-!  WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'
-!  ALLOCATE(tmp_var(iml-1,jml))
-!  CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)
-!  ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
-!  DEALLOCATE(tmp_var)
   DEALLOCATE(relief_hi,lon_rad,lat_rad)
 
@@ -405,5 +324,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis)
+SUBROUTINE start_init_phys(lon_in,lat_in,phis)
 !
 !===============================================================================
@@ -413,11 +332,9 @@
 !-------------------------------------------------------------------------------
 ! Arguments:
-  REAL,    INTENT(IN) :: lon_in (:),  lat_in (:)     ! dim (iml) (jml)
-  REAL,    INTENT(IN) :: lon_in2(:),  lat_in2(:)     ! dim (iml) (jml2)
-  LOGICAL, INTENT(IN) :: ibar
+  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
   REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc"
+  CHARACTER(LEN=256) :: modname
   REAL               :: date, dt
   INTEGER            :: iml, jml, jml2, itau(1)
@@ -426,15 +343,15 @@
   REAL, ALLOCATABLE  :: ts(:,:), qs(:,:)
 !-------------------------------------------------------------------------------
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),              TRIM(modname)//" jml")
-  jml2=SIZE(lat_in2)
+  modname="start_init_phys"
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
+  jml=SIZE(phis,2); jml2=SIZE(lat_in)
 
   WRITE(lunout,*)'Opening the surface analysis'
-  CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
-  WRITE(lunout,*) 'Values read: ',   iml_phys, jml_phys, llm_phys, ttm_phys
+  CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
+  WRITE(lunout,*) 'Values read: ',  iml_phys, jml_phys, llm_phys, ttm_phys
 
   ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
   ALLOCATE(levphys_ini(llm_phys))
-  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
+  CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
                 lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
 
@@ -445,6 +362,6 @@
 
   ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
-  CALL get_var_phys('ST'  ,ts)                   !--- SURFACE TEMPERATURE
-  CALL get_var_phys('CDSW',qs)                   !--- SOIL MOISTURE
+  CALL get_var_phys(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
+  CALL get_var_phys(qsolvar,qs)                   !--- SOIL MOISTURE
   CALL flinclo(fid_phys)
   DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
@@ -463,5 +380,4 @@
 !
 !-------------------------------------------------------------------------------
-  USE conf_dat_m, ONLY: conf_dat2d
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -474,13 +390,13 @@
 !-------------------------------------------------------------------------------
   SELECT CASE(title)
-    CASE('SP');        tllm=0
-    CASE('ST','CDSW'); tllm=llm_phys
+    CASE(psrfvar);         tllm=0
+    CASE(tsrfvar,qsolvar); tllm=llm_phys
   END SELECT
   IF(ALLOCATED(field)) RETURN
   ALLOCATE(field(iml,jml)); field(:,:)=0.
   CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
-  CALL conf_dat2d(title, lon_ini, lat_ini,  lon_rad, lat_rad, var_ana, ibar)
-  CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,         &
-                            lon_in, lat_in, lon_in2, lat_in2, field)
+  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
+  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
+                                      lon_in,  lat_in,  field)
 
 END SUBROUTINE get_var_phys
@@ -495,41 +411,35 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
 !
 !-------------------------------------------------------------------------------
   USE inter_barxy_m, ONLY: inter_barxy
-  USE grid_atob_m,   ONLY: grille_m
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*), INTENT(IN)  :: nam
-  LOGICAL,          INTENT(IN)  :: ibar, ibeg
+  LOGICAL,          INTENT(IN)  :: ibeg
   REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
   REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
-  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
   REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
   REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=256) :: modname="interp_startvar"
+  CHARACTER(LEN=256) :: modname
   INTEGER            :: ii, jj, i1, j1, j2
   REAL, ALLOCATABLE  :: vtmp(:,:)
 !-------------------------------------------------------------------------------
-  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
-  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
-  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
-  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
-  j2=SIZE(lat2)
+  modname="interp_startvar"
+  ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
+  jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
+  i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
+  j1=SIZE(varo,2); j2=SIZE(lat2)
   ALLOCATE(vtmp(i1-1,j1))
-  IF(ibar) THEN
-    IF(ibeg.AND.prt_level>1) THEN
-      WRITE(lunout,*)"--------------------------------------------------------"
-      WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
-      WRITE(lunout,*)"--------------------------------------------------------"
-    END IF
-    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
-  ELSE
-    CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
+  IF(ibeg.AND.prt_level>1) THEN
+    WRITE(lunout,*)"--------------------------------------------------------"
+    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+    WRITE(lunout,*)"--------------------------------------------------------"
   END IF
+  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
 
