Index: /branches/UKMO/restart_datestamp/NEMOGCM/ARCH/arch-XC40_METO.fcm
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/ARCH/arch-XC40_METO.fcm	(revision 5539)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/ARCH/arch-XC40_METO.fcm	(revision 5539)
@@ -0,0 +1,52 @@
+# compiler options for Archer CRAY XC-40 (using crayftn compiler)
+#
+# NCDF_HOME   root directory containing lib and include subdirectories for netcdf4
+# HDF5_HOME   root directory containing lib and include subdirectories for HDF5
+# XIOS_HOME   root directory containing lib for XIOS
+# OASIS_HOME  root directory containing lib for OASIS
+#
+# NCDF_INC    netcdf4 include file
+# NCDF_LIB    netcdf4 library
+# XIOS_INC    xios include file    (taken into accound only if key_iomput is activated)
+# XIOS_LIB    xios library         (taken into accound only if key_iomput is activated)
+# OASIS_INC   oasis include file   (taken into accound only if key_oasis3 is activated)
+# OASIS_LIB   oasis library        (taken into accound only if key_oasis3 is activated)
+#
+# FC          Fortran compiler command
+# FCFLAGS     Fortran compiler flags
+# FFLAGS      Fortran 77 compiler flags
+# LD          linker
+# LDFLAGS     linker flags, e.g. -L<lib dir> if you have libraries
+# FPPFLAGS    pre-processing flags
+# AR          assembler
+# ARFLAGS     assembler flags
+# MK          make
+# USER_INC    complete list of include files
+# USER_LIB    complete list of libraries to pass to the linker
+#
+# Note that:
+#  - unix variables "$..." are accpeted and will be evaluated before calling fcm.
+#  - fcm variables are starting with a % (and not a $)
+#
+%NCDF_HOME           /opt/cray/netcdf-hdf5parallel/4.3.2/CRAY/83
+%HDF5_HOME           /opt/cray/hdf5-parallel/1.8.13/CRAY/83
+%XIOS_HOME           /projects/ocean/nemo/xios/xios_1.0_r618_20150619/XIOS
+
+%NCDF_INC            -I%NCDF_HOME/include -I%HDF5_HOME/include
+%NCDF_LIB            -L%HDF5_HOME/lib -L%NCDF_HOME/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 -lz
+%XIOS_INC            -I%XIOS_HOME/inc 
+%XIOS_LIB            -L%XIOS_HOME/lib -lxios
+
+%CPP                 cpp
+%FC                  ftn
+%FCFLAGS             -emf -s real64 -s integer32  -O2 -hflex_mp=intolerant -Rb
+%FFLAGS              -emf -s real64 -s integer32  -O0 -hflex_mp=strict -Rb
+%LD                  ftn
+%FPPFLAGS            -P -E -traditional-cpp
+%LDFLAGS             -hbyteswapio 
+%AR                  ar 
+%ARFLAGS             -r
+%MK                  gmake
+
+%USER_INC            %XIOS_INC
+%USER_LIB            %XIOS_LIB
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 5539)
@@ -99,4 +99,6 @@
    ln_ssr      = .false.   !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
    nn_fwb      = 0         !  FreshWater Budget: =0 unchecked
+   ln_apr_dyn  = .false.    !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
+
 /
 !-----------------------------------------------------------------------
@@ -158,10 +160,10 @@
 !          !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
 !          !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_apr      = 'patm'    ,         -1        ,'somslpre',    .true.      , .true.  , 'yearly'  ,  ''      ,   '' , ''
+   sn_apr     = 'amm12_mslp'     ,          1        ,  'p_msl'      , .false.      , .false. , 'daily'   ,  ''      ,  ''  , ''
 
-   cn_dir      = './'       !  root directory for the location of the bulk files
-   rn_pref     = 101000._wp !  reference atmospheric pressure   [N/m2]/
+   cn_dir      = './fluxes/' !  root directory for the location of the bulk files
+   rn_pref     = 101000.     !  reference atmospheric pressure   [N/m2]/
    ln_ref_apr  = .false.    !  ref. pressure: global mean Patm (T) or a constant (F)
-   ln_apr_obc  = .false.    !  inverse barometer added to OBC ssh data
+   ln_apr_obc  = .true.    !  inverse barometer added to OBC ssh data
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist_cfg
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist_cfg	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist_cfg	(revision 5539)
@@ -94,13 +94,12 @@
    sn_snow     = 'ncar_precip.15JUNE2009_fill' ,        -1         , 'SNOW'    ,   .false.    , .true. , 'yearly'  , 'weights_core2_agrif_bilinear.nc'     , ''       , ''
    sn_tdif     = 'taudif_core'                 ,        24         , 'taudif'  ,   .false.    , .true. , 'yearly'  , 'weights_core2_agrif_bilinear.nc'     , ''       , ''
+
    cn_dir      = './'      !  root directory for the location of the bulk files
-   ln_2m       = .false.   !  air temperature and humidity referenced at 2m (T) instead 10m (F)
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
-   ln_bulk2z   = .false.   !  Air temperature/humidity and wind vectors are referenced at heights rn_zqt and rn_zu
-   rn_zqt      = 3.        !  Air temperature and humidity reference height (m) (ln_bulk2z)
-   rn_zu       = 4.        !  Wind vector reference height (m)                  (ln_bulk2z)
+   rn_zqt      = 10.        !  Air temperature and humidity reference height (m)
+   rn_zu       = 10.        !  Wind vector reference height (m)
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
    rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
-   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
+   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity
                            !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds)
 /
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_demo.xml
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_demo.xml	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_demo.xml	(revision 5539)
@@ -64,5 +64,5 @@
 	  <field field_ref="sst"    name="sst_1d_inst" operation="instant" />     <!-- instant value -->
 	  <field field_ref="sst"    name="sst_1d_max"  operation="maximum" />     <!-- max -->
-	  <field field_ref="suoce"  />   <!-- include a U-grid variable in the list -->
+	  <field field_ref="ssu"  />   <!-- include a U-grid variable in the list -->
 	</file>
 	
@@ -78,6 +78,8 @@
     
    <axis_definition>  
-      <axis id="deptht"  long_name="Vertical T levels"  unit="m" positive="down" />
-      <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10   " axis_ref="deptht" />
+      <axis_group id="deptht" long_name="Vertical T levels" unit="m" positive="down" > 
+         <axis id="deptht" /> 
+         <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10" /> 
+      </axis_group>
       <axis id="depthu"  long_name="Vertical U levels"  unit="m" positive="down" />
       <axis id="depthv"  long_name="Vertical V levels"  unit="m" positive="down" />
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/1_namelist_cfg
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/1_namelist_cfg	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/1_namelist_cfg	(revision 5539)
@@ -94,13 +94,12 @@
    sn_snow     = 'ncar_precip.15JUNE2009_fill' ,        -1         , 'SNOW'    ,   .false.    , .true. , 'yearly'  , 'weights_core2_agrif_bilinear.nc'     , ''       , ''
    sn_tdif     = 'taudif_core'                 ,        24         , 'taudif'  ,   .false.    , .true. , 'yearly'  , 'weights_core2_agrif_bilinear.nc'     , ''       , ''
+
    cn_dir      = './'      !  root directory for the location of the bulk files
-   ln_2m       = .false.   !  air temperature and humidity referenced at 2m (T) instead 10m (F)
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
-   ln_bulk2z   = .false.   !  Air temperature/humidity and wind vectors are referenced at heights rn_zqt and rn_zu
-   rn_zqt      = 3.        !  Air temperature and humidity reference height (m) (ln_bulk2z)
-   rn_zu       = 4.        !  Wind vector reference height (m)                  (ln_bulk2z)
+   rn_zqt      = 10.        !  Air temperature and humidity reference height (m)
+   rn_zu       = 10.        !  Wind vector reference height (m)
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
    rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
-   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
+   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity
                            !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds)
 /
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 5539)
@@ -190,4 +190,6 @@
           <field field_ref="salinity_cat"    name="salincat" />
           <field field_ref="brinevol_cat"    name="sibricat" />
+	  <field field_ref="icetemp_cat"     name="sitemcat" />
+	  <field field_ref="snwtemp_cat"     name="sntemcat" />
 
 	</file>
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5539)
@@ -42,4 +42,6 @@
          <field id="sssmin"       long_name="min of sea surface salinity"   field_ref="sss"   operation="minimum"                 />
          <field id="sbs"          long_name="sea bottom salinity"                                                     unit="1e-3" />
+
+         <field id="taubot"       long_name="bottom stress module"                                                    unit="N/m2" /> 
 
          <field id="ssh"          long_name="sea surface height"             standard_name="sea_surface_height_above_geoid"             unit="m" />
@@ -279,4 +281,6 @@
          <field id="salinity_cat" long_name="Sea-Ice Bulk salinity for categories"                         unit="g/kg"   axis_ref="ncatice" />
          <field id="brinevol_cat" long_name="Brine volume for categories"                                  unit="%"      axis_ref="ncatice" />
+         <field id="icetemp_cat"  long_name="Ice temperature for categories"                               unit="degC"   axis_ref="ncatice" />
+         <field id="snwtemp_cat"  long_name="Snow temperature for categories"                              unit="degC"   axis_ref="ncatice" />
 
          <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     />
@@ -366,4 +370,5 @@
          <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             />
          <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             />
+         <field id="ubar"         long_name="ocean barotropic current along i-axis"                                                              unit="m/s"                             />
          <field id="uocetr_eff"   long_name="Effective ocean transport along i-axis"                 standard_name="ocean_volume_x_transport"    unit="m3/s"       grid_ref="grid_U_3D" />
          <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_U_3D" />
@@ -400,4 +405,5 @@
          <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             />
          <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             />
+         <field id="vbar"         long_name="ocean barotropic current along j-axis"                                                              unit="m/s"                             />
          <field id="vocetr_eff"   long_name="Effective ocean transport along j-axis"                 standard_name="ocean_volume_y_transport"    unit="m3/s"       grid_ref="grid_V_3D" />
          <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_V_3D" />
Index: /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5539)
@@ -476,5 +476,5 @@
 
    cn_dir      = './'       !  root directory for the location of the bulk files
-   rn_pref     = 101000._wp !  reference atmospheric pressure   [N/m2]/
+   rn_pref     = 101000.    !  reference atmospheric pressure   [N/m2]/
    ln_ref_apr  = .false.    !  ref. pressure: global mean Patm (T) or a constant (F)
    ln_apr_obc  = .false.    !  inverse barometer added to OBC ssh data
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5539)
@@ -30,5 +30,5 @@
 
    PUBLIC   lim_thd_dh      ! called by lim_thd
-   PUBLIC   lim_thd_snwblow ! called in sbcblk/sbccpl and here
+   PUBLIC   lim_thd_snwblow ! called in sbcblk/sbcclio/sbccpl and here
 
    INTERFACE lim_thd_snwblow
@@ -122,24 +122,20 @@
       END SELECT
 
-      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
-      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw )
+      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
+      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
       CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_alloc( jpij, nlay_i, icount )
-      
+       
       dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
       dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
- 
-      zqprec (:) = 0._wp ; zq_su  (:) = 0._wp ; zq_bo  (:) = 0._wp ; zf_tt  (:) = 0._wp
-      zq_rema(:) = 0._wp
-
-      zdh_s_pre(:) = 0._wp
-      zdh_s_mel(:) = 0._wp
-      zdh_s_sub(:) = 0._wp
-      zqh_s    (:) = 0._wp      
-      zqh_i    (:) = 0._wp   
-
-      zh_i      (:,:) = 0._wp       
-      zdeltah   (:,:) = 0._wp       
-      icount    (:,:) = 0
+
+      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
+      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp
+      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
+      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp     
+
+      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
+      icount (:,:) = 0
+
 
       ! Initialize enthalpy at nlay_i+1
@@ -223,6 +219,7 @@
       ! Martin Vancoppenolle, December 2006
 
+      CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing
+
       zdeltah(:,:) = 0._wp
-      CALL lim_thd_snwblow( 1. - at_i_1d, zsnw ) ! snow distribution over ice after wind blowing
       DO ji = kideb, kiut
          !-----------
@@ -689,6 +686,6 @@
       WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
       
-      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
-      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw )
+      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
+      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
       CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_dealloc( jpij, nlay_i, icount )
@@ -703,12 +700,12 @@
    !!--------------------------------------------------------------------------
    SUBROUTINE lim_thd_snwblow_2d( pin, pout )
-      REAL(wp), DIMENSION(:,:), INTENT(in)  :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
-      REAL(wp), DIMENSION(:,:), INTENT(out) :: pout
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
+      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
       pout = ( 1._wp - ( pin )**rn_betas )
    END SUBROUTINE lim_thd_snwblow_2d
 
    SUBROUTINE lim_thd_snwblow_1d( pin, pout )
-      REAL(wp), DIMENSION(:), INTENT(in)  :: pin
-      REAL(wp), DIMENSION(:), INTENT(out) :: pout
+      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
+      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
       pout = ( 1._wp - ( pin )**rn_betas )
    END SUBROUTINE lim_thd_snwblow_1d
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5539)
@@ -756,5 +756,9 @@
 
       ! --- diagnose the change in non-solar flux due to surface temperature change --- !
-      IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:)  - zqns_ice_b(:) ) * a_i_1d(:) 
+      IF ( ln_it_qnsice ) THEN
+         DO ji = kideb, kiut
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
+         END DO
+      END IF
 
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
@@ -790,5 +794,5 @@
       CALL wrk_dealloc( jpij, numeqmin, numeqmax )
       CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
-      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
+      CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
       CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
       CALL wrk_dealloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5539)
@@ -60,5 +60,5 @@
       REAL(wp) ::  z1_365
       REAL(wp) ::  ztmp
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s
       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace
       !!-------------------------------------------------------------------
@@ -66,5 +66,5 @@
       IF( nn_timing == 1 )  CALL timing_start('limwri')
 
-      CALL wrk_alloc( jpi, jpj, jpl, zoi, zei )
+      CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
       CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi )
 
@@ -243,4 +243,16 @@
       CALL iom_put( "snowthic_cat"     , ht_s        )        ! snow depth for categories
       CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories
+
+      ! ice temperature
+      IF ( iom_use( "icetemp_cat" ) ) THEN 
+         zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i
+         CALL iom_put( "icetemp_cat"   , zt_i - rt0  )
+      ENDIF
+      
+      ! snow temperature
+      IF ( iom_use( "snwtemp_cat" ) ) THEN 
+         zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s
+         CALL iom_put( "snwtemp_cat"   , zt_s - rt0  )
+      ENDIF
 
       ! Compute ice age
@@ -280,5 +292,5 @@
       !     not yet implemented
       
-      CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei )
+      CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
       CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db )
 
Index: anches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/domain.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 5538)
+++ 	(revision )
@@ -1,431 +1,0 @@
-MODULE domain
-   !!==============================================================================
-   !!                       ***  MODULE domain   ***
-   !! Ocean initialization : domain initialization
-   !!==============================================================================
-
-   !!----------------------------------------------------------------------
-   !!   dom_init       : initialize the space and time domain
-   !!   dom_nam        : read and contral domain namelists
-   !!   dom_ctl        : control print for the ocean domain
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! 
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distributed memory computing library
-
-   USE domstp          ! domain: set the time-step
-   USE domrea          ! domain: write the meshmask file
-   USE dommsk          ! domain : mask
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC dom_init       ! called by opa.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE dom_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_init  ***
-      !!                    
-      !! ** Purpose :   Domain initialization. Call the routines that are 
-      !!      required to create the arrays which define the space and time
-      !!      domain of the ocean model.
-      !!
-      !! ** Method  :
-      !!      - dom_stp: defined the model time step
-      !!      - dom_rea: read the meshmask file if nmsh=1
-      !!
-      !! History :
-      !!        !  90-10  (C. Levy - G. Madec)  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  92-01  (M. Imbard) insert time step initialization
-      !!        !  96-06  (G. Madec) generalized vertical coordinate 
-      !!        !  97-02  (G. Madec) creation of domwri.F
-      !!        !  01-05  (E.Durand - G. Madec) insert closed sea
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   jk                ! dummy loop argument
-      INTEGER ::   iconf = 0         ! temporary integers
-      !!----------------------------------------------------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_init : domain initialization'
-         WRITE(numout,*) '~~~~~~~~'
-      ENDIF
-
-      CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
-      CALL dom_zgr      ! Vertical mesh and bathymetry option
-      CALL dom_rea      ! Create a domain file
-
-     !
-      ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines
-      !        but could be usefull in many other routines
-      e12t    (:,:) = e1t(:,:) * e2t(:,:)
-      e1e2t   (:,:) = e1t(:,:) * e2t(:,:)
-      e12u    (:,:) = e1u(:,:) * e2u(:,:)
-      e12v    (:,:) = e1v(:,:) * e2v(:,:)
-      e12f    (:,:) = e1f(:,:) * e2f(:,:)
-      r1_e12t (:,:) = 1._wp    / e12t(:,:)
-      r1_e12u (:,:) = 1._wp    / e12u(:,:)
-      r1_e12v (:,:) = 1._wp    / e12v(:,:)
-      r1_e12f (:,:) = 1._wp    / e12f(:,:)
-      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
-      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
-      !
-      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
-      hv(:,:) = 0._wp
-      DO jk = 1, jpk
-         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
-
-      CALL dom_stp      ! Time step
-      CALL dom_msk      ! Masks
-      CALL dom_ctl      ! Domain control
-
-   END SUBROUTINE dom_init
-
-   SUBROUTINE dom_nam
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_nam  ***
-      !!                    
-      !! ** Purpose :   read domaine namelists and print the variables.
-      !!
-      !! ** input   : - namrun namelist
-      !!              - namdom namelist
-      !!              - namcla namelist
-      !!----------------------------------------------------------------------
-      USE ioipsl
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
-         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
-         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
-         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
-      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
-         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
-         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs, &
-         &             jphgr_msh, &
-         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
-         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
-         &             ppa2, ppkth2, ppacr2
-      NAMELIST/namcla/ nn_cla
-#if defined key_netcdf4
-      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
-#endif
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
-      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
-      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namrun )
-      !
-      IF(lwp) THEN                  ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
-         WRITE(numout,*) '~~~~~~~ '
-         WRITE(numout,*) '   Namelist namrun'  
-         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
-         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
-         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
-         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
-         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
-         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
-         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
-         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
-         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
-         WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
-         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
-         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
-         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
-         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
-         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
-         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
-      ENDIF
-      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
-      cexper = cn_exp
-      nrstdt = nn_rstctl
-      nit000 = nn_it000
-      nitend = nn_itend
-      ndate0 = nn_date0
-      nleapy = nn_leapy
-      ninist = nn_istate
-      nstock = nn_stock
-      nstocklist = nn_stocklist
-      nwrite = nn_write
-
-
-      !                             ! control of output frequency
-      IF ( nstock == 0 .OR. nstock > nitend ) THEN
-         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
-         CALL ctl_warn( ctmp1 )
-         nstock = nitend
-      ENDIF
-      IF ( nwrite == 0 ) THEN
-         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
-         CALL ctl_warn( ctmp1 )
-         nwrite = nitend
-      ENDIF
-
-      ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day)
-      ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
-      adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
-
-#if defined key_agrif
-      IF( Agrif_Root() ) THEN
-#endif
-      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
-      CASE (  1 ) 
-         CALL ioconf_calendar('gregorian')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
-      CASE (  0 )
-         CALL ioconf_calendar('noleap')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
-      CASE ( 30 )
-         CALL ioconf_calendar('360d')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
-      END SELECT
-#if defined key_agrif
-      ENDIF
-#endif
-
-      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
-      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
-903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
-      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
-904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namdom )
-
-      IF(lwp) THEN
-         WRITE(numout,*) 
-         WRITE(numout,*) '   Namelist namdom : space & time domain'
-         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
-         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
-         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
-         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
-         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
-         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
-         WRITE(numout,*) '           = 0   no file created                 '
-         WRITE(numout,*) '           = 1   mesh_mask                       '
-         WRITE(numout,*) '           = 2   mesh and mask                   '
-         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask      '
-         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
-         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
-         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
-         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
-         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin
-         WRITE(numout,*) '                  bottom  tracer rdt       rdtmax    = ', rn_rdtmax
-         WRITE(numout,*) '                  depth of transition      rn_rdth   = ', rn_rdth
-         WRITE(numout,*) '      suppression of closed seas (=0)      nn_closea = ', nn_closea
-         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
-         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
-         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
-         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
-         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
-         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
-         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
-         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
-         WRITE(numout,*) '                                        ppa0            = ', ppa0
-         WRITE(numout,*) '                                        ppa1            = ', ppa1
-         WRITE(numout,*) '                                        ppkth           = ', ppkth
-         WRITE(numout,*) '                                        ppacr           = ', ppacr
-         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
-         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
-         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
-         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
-         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
-         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
-      ENDIF
-
-      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
-      e3zps_min = rn_e3zps_min
-      e3zps_rat = rn_e3zps_rat
-      nmsh      = nn_msh
-      nacc      = nn_acc
-      atfp      = rn_atfp
-      rdt       = rn_rdt
-      rdtmin    = rn_rdtmin
-      rdtmax    = rn_rdtmin
-      rdth      = rn_rdth
-
-      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
-      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
-905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
-      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
-906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namcla )
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist namcla'
-         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
-      ENDIF
-
-#if defined key_netcdf4
-      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
-      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
-      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
-907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
-      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
-908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namnc4 )
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
-         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
-         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
-         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
-         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
-      ENDIF
-
-      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
-      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
-      snc4set%ni   = nn_nchunks_i
-      snc4set%nj   = nn_nchunks_j
-      snc4set%nk   = nn_nchunks_k
-      snc4set%luse = ln_nc4zip
-#else
-      snc4set%luse = .FALSE.        ! No NetCDF 4 case
-#endif
-      !
-   END SUBROUTINE dom_nam
-
-   SUBROUTINE dom_zgr
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_zgr  ***
-      !!                   
-      !! ** Purpose :  set the depth of model levels and the resulting 
-      !!      vertical scale factors.
-      !!
-      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
-      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
-      !!              - vertical coordinate (gdep., e3.) depending on the 
-      !!                coordinate chosen :
-      !!                   ln_zco=T   z-coordinate  
-      !!                   ln_zps=T   z-coordinate with partial steps
-      !!                   ln_zco=T   s-coordinate 
-      !!
-      !! ** Action  :   define gdep., e3., mbathy and bathy
-      !!----------------------------------------------------------------------
-      INTEGER ::   ioptio = 0   ! temporary integer
-      INTEGER ::   ios
-      !!
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
-      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
-      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzgr )
-
-      IF(lwp) THEN                     ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_zgr : vertical coordinate'
-         WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
-         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
-         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
-         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
-         WRITE(numout,*) '             ice shelf cavity               ln_isfcav = ', ln_isfcav
-      ENDIF
-
-      ioptio = 0                       ! Check Vertical coordinate options
-      IF( ln_zco ) ioptio = ioptio + 1
-      IF( ln_zps ) ioptio = ioptio + 1
-      IF( ln_sco ) ioptio = ioptio + 1
-      IF( ln_isfcav ) ioptio = 33
-      IF ( ioptio /= 1  )   CALL ctl_stop( ' none or several vertical coordinate options used' )
-      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
-
-   END SUBROUTINE dom_zgr
-
-   SUBROUTINE dom_ctl
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_ctl  ***
-      !!
-      !! ** Purpose :   Domain control.
-      !!
-      !! ** Method  :   compute and print extrema of masked scale factors
-      !!
-      !! History :
-      !!   8.5  !  02-08  (G. Madec)    Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
-      INTEGER, DIMENSION(2) ::   iloc      ! 
-      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
-      !!----------------------------------------------------------------------
-
-      ! Extrema of the scale factors
-
-      IF(lwp)WRITE(numout,*)
-      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
-      IF(lwp)WRITE(numout,*) '~~~~~~~'
-
-      IF (lk_mpp) THEN
-         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
-         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
-         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
-         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
-      ELSE
-         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-
-         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iimi1 = iloc(1) + nimpp - 1
-         ijmi1 = iloc(2) + njmpp - 1
-         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iimi2 = iloc(1) + nimpp - 1
-         ijmi2 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iima1 = iloc(1) + nimpp - 1
-         ijma1 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iima2 = iloc(1) + nimpp - 1
-         ijma2 = iloc(2) + njmpp - 1
-      ENDIF
-
-      IF(lwp) THEN
-         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
-         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
-         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
-         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
-      ENDIF
-
-   END SUBROUTINE dom_ctl
-
-   !!======================================================================
-END MODULE domain
Index: anches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/dommsk.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 5538)
+++ 	(revision )
@@ -1,132 +1,0 @@
-MODULE dommsk
-   !!======================================================================
-   !!                       ***  MODULE dommsk   ***
-   !! Ocean initialization : domain land/sea masks, off-line case 
-   !!======================================================================
-   !! History :  3.3  ! 2010-10  (C. Ethe)  adapted from OPA_SRC/DOM/dommsk
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dom_msk        : compute land/ocean mask
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE lib_mpp         ! MPP library
-   USE in_out_manager  ! I/O manager
-   USE wrk_nemo  
-   USE lbclnk  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dom_msk    ! routine called by inidom.F90
-
-   REAL(wp)        :: rn_shlat   = 2.   ! type of lateral boundary condition on velocity
-   LOGICAL, PUBLIC :: ln_vorlat  = .false.   !  consistency of vorticity boundary condition 
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dom_msk
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE dom_msk  ***
-      !!
-      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
-      !!
-      !! ** Method  :   The interior ocean/land mask is computed from tmask
-      !!              setting to zero the duplicated row and lines due to
-      !!              MPP exchange halos, est-west cyclic and north fold
-      !!              boundary conditions.
-      !!
-      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
-      !!               tpol     : ???
-      !!----------------------------------------------------------------------
-      !
-      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
-      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
-      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
-      !
-      !!---------------------------------------------------------------------
-      
-      CALL wrk_alloc( jpi, jpj, imsk )
-      !
-      ! Interior domain mask (used for global sum)
-      ! --------------------
-      ssmask(:,:)  = tmask(:,:,1)
-      tmask_i(:,:) = tmask(:,:,1)
-      iif = jpreci                        ! thickness of exchange halos in i-axis
-      iil = nlci - jpreci + 1
-      ijf = jprecj                        ! thickness of exchange halos in j-axis
-      ijl = nlcj - jprecj + 1
-      !
-      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
-      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
-      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
-      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
-      !
-      !                                   ! north fold mask
-      tpol(1:jpiglo) = 1._wp
-      !                                
-      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
-      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
-      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
-         IF( mjg(ijl-1) == jpjglo-1 ) THEN
-            DO ji = iif+1, iil-1
-               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
-            END DO
-         ENDIF
-      ENDIF 
-      !
-      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
-      ! least 1 wet u point
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector loop
-            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
-            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
-         END DO
-         DO ji = 1, jpim1      ! NO vector opt.
-            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
-               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
-         END DO
-      END DO
-      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask_i, 'V', 1._wp )
-      CALL lbc_lnk( fmask_i, 'F', 1._wp )
-
-      ! 3. Ocean/land mask at wu-, wv- and w points 
-      !----------------------------------------------
-      wmask (:,:,1) = tmask(:,:,1) ! ????????
-      wumask(:,:,1) = umask(:,:,1) ! ????????
-      wvmask(:,:,1) = vmask(:,:,1) ! ????????
-      DO jk=2,jpk
-         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
-         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
-         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
-      END DO
-      !
-      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
-         imsk(:,:) = INT( tmask_i(:,:) )
-         WRITE(numout,*) ' tmask_i : '
-         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
-         WRITE (numout,*)
-         WRITE (numout,*) ' dommsk: tmask for each level'
-         WRITE (numout,*) ' ----------------------------'
-         DO jk = 1, jpk
-            imsk(:,:) = INT( tmask(:,:,jk) )
-            WRITE(numout,*)
-            WRITE(numout,*) ' level = ',jk
-            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
-         END DO
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, imsk )
-      !
-   END SUBROUTINE dom_msk
-   !!======================================================================
-END MODULE dommsk
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 5539)
@@ -1,35 +1,41 @@
 MODULE domrea
-   !!======================================================================
-   !!                       ***  MODULE domrea  ***
-   !! Ocean initialization : read the ocean domain meshmask file(s)
-   !!======================================================================
-   !! History :  3.3  ! 2010-05  (C. Ethe)  Full reorganization of the off-line
+   !!==============================================================================
+   !!                       ***  MODULE domrea   ***
+   !! Ocean initialization : domain initialization
+   !!==============================================================================
+
    !!----------------------------------------------------------------------
-
+   !!   dom_init       : initialize the space and time domain
+   !!   dom_nam        : read and contral domain namelists
+   !!   dom_ctl        : control print for the ocean domain
    !!----------------------------------------------------------------------
-   !!   dom_rea        : read mesh and mask file(s)
-   !!                    nmsh = 1  :   mesh_mask file
-   !!                         = 2  :   mesh and mask file
-   !!                         = 3  :   mesh_hgr, mesh_zgr and mask
-   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! 
    USE dom_oce         ! ocean space and time domain
-   USE dommsk          ! domain: masks
+   USE phycst          ! physical constants
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distributed memory computing library
+
+   USE domstp          ! domain: set the time-step
+
    USE lbclnk          ! lateral boundary condition - MPP exchanges
    USE trc_oce         ! shared ocean/biogeochemical variables
-   USE lib_mpp 
-   USE in_out_manager
    USE wrk_nemo  
-
+   
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   dom_rea    ! routine called by inidom.F90
-  !! * Substitutions
+   !! * Routine accessibility
+   PUBLIC dom_rea       ! called by opa.F90
+
+   !! * Substitutions
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OFF 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
+
 CONTAINS
 
@@ -37,4 +43,395 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dom_rea  ***
+      !!                    
+      !! ** Purpose :   Domain initialization. Call the routines that are 
+      !!      required to create the arrays which define the space and time
+      !!      domain of the ocean model.
+      !!
+      !! ** Method  :
+      !!      - dom_stp: defined the model time step
+      !!      - dom_rea: read the meshmask file if nmsh=1
+      !!
+      !! History :
+      !!        !  90-10  (C. Levy - G. Madec)  Original code
+      !!        !  91-11  (G. Madec)
+      !!        !  92-01  (M. Imbard) insert time step initialization
+      !!        !  96-06  (G. Madec) generalized vertical coordinate 
+      !!        !  97-02  (G. Madec) creation of domwri.F
+      !!        !  01-05  (E.Durand - G. Madec) insert closed sea
+      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   jk                ! dummy loop argument
+      INTEGER ::   iconf = 0         ! temporary integers
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_init : domain initialization'
+         WRITE(numout,*) '~~~~~~~~'
+      ENDIF
+
+      CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
+      CALL dom_zgr      ! Vertical mesh and bathymetry option
+      CALL dom_grd      ! Create a domain file
+
+     !
+      ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines
+      !        but could be usefull in many other routines
+      e12t    (:,:) = e1t(:,:) * e2t(:,:)
+      e1e2t   (:,:) = e1t(:,:) * e2t(:,:)
+      e12u    (:,:) = e1u(:,:) * e2u(:,:)
+      e12v    (:,:) = e1v(:,:) * e2v(:,:)
+      e12f    (:,:) = e1f(:,:) * e2f(:,:)
+      r1_e12t (:,:) = 1._wp    / e12t(:,:)
+      r1_e12u (:,:) = 1._wp    / e12u(:,:)
+      r1_e12v (:,:) = 1._wp    / e12v(:,:)
+      r1_e12f (:,:) = 1._wp    / e12f(:,:)
+      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
+      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
+      !
+      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
+      hv(:,:) = 0._wp
+      DO jk = 1, jpk
+         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !                                        ! Inverse of the local depth
+      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
+      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
+
+      CALL dom_stp      ! Time step
+      CALL dom_msk      ! Masks
+      CALL dom_ctl      ! Domain control
+
+   END SUBROUTINE dom_rea
+
+   SUBROUTINE dom_nam
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read domaine namelists and print the variables.
+      !!
+      !! ** input   : - namrun namelist
+      !!              - namdom namelist
+      !!              - namcla namelist
+      !!----------------------------------------------------------------------
+      USE ioipsl
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
+         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
+         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
+      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
+         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
+         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs, &
+         &             jphgr_msh, &
+         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
+         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
+         &             ppa2, ppkth2, ppacr2
+      NAMELIST/namcla/ nn_cla
+#if defined key_netcdf4
+      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
+#endif
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
+      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namrun )
+      !
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namrun'  
+         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
+         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
+         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
+         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
+         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
+         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
+         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
+         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
+         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
+         WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
+         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
+         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
+         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
+         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
+         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
+         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
+      ENDIF
+      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
+      cexper = cn_exp
+      nrstdt = nn_rstctl
+      nit000 = nn_it000
+      nitend = nn_itend
+      ndate0 = nn_date0
+      nleapy = nn_leapy
+      ninist = nn_istate
+      nstock = nn_stock
+      nstocklist = nn_stocklist
+      nwrite = nn_write
+
+
+      !                             ! control of output frequency
+      IF ( nstock == 0 .OR. nstock > nitend ) THEN
+         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
+         CALL ctl_warn( ctmp1 )
+         nstock = nitend
+      ENDIF
+      IF ( nwrite == 0 ) THEN
+         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
+         CALL ctl_warn( ctmp1 )
+         nwrite = nitend
+      ENDIF
+
+      ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day)
+      ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
+      adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
+
+#if defined key_agrif
+      IF( Agrif_Root() ) THEN
+#endif
+      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
+      CASE (  1 ) 
+         CALL ioconf_calendar('gregorian')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
+      CASE (  0 )
+         CALL ioconf_calendar('noleap')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
+      CASE ( 30 )
+         CALL ioconf_calendar('360d')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
+      END SELECT
+#if defined key_agrif
+      ENDIF
+#endif
+
+      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namdom )
+
+      IF(lwp) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) '   Namelist namdom : space & time domain'
+         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
+         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
+         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
+         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
+         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
+         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
+         WRITE(numout,*) '           = 0   no file created                 '
+         WRITE(numout,*) '           = 1   mesh_mask                       '
+         WRITE(numout,*) '           = 2   mesh and mask                   '
+         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask      '
+         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
+         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
+         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
+         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
+         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin
+         WRITE(numout,*) '                  bottom  tracer rdt       rdtmax    = ', rn_rdtmax
+         WRITE(numout,*) '                  depth of transition      rn_rdth   = ', rn_rdth
+         WRITE(numout,*) '      suppression of closed seas (=0)      nn_closea = ', nn_closea
+         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
+         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
+         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
+         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
+         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
+         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
+         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
+         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
+         WRITE(numout,*) '                                        ppa0            = ', ppa0
+         WRITE(numout,*) '                                        ppa1            = ', ppa1
+         WRITE(numout,*) '                                        ppkth           = ', ppkth
+         WRITE(numout,*) '                                        ppacr           = ', ppacr
+         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
+         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
+         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
+         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
+         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
+         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
+      ENDIF
+
+      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
+      e3zps_min = rn_e3zps_min
+      e3zps_rat = rn_e3zps_rat
+      nmsh      = nn_msh
+      nacc      = nn_acc
+      atfp      = rn_atfp
+      rdt       = rn_rdt
+      rdtmin    = rn_rdtmin
+      rdtmax    = rn_rdtmin
+      rdth      = rn_rdth
+
+      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
+      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
+      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namcla )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namcla'
+         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
+      ENDIF
+
+#if defined key_netcdf4
+      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
+      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
+      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
+907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
+      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
+908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namnc4 )
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
+         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
+         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
+         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
+         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
+      ENDIF
+
+      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
+      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
+      snc4set%ni   = nn_nchunks_i
+      snc4set%nj   = nn_nchunks_j
+      snc4set%nk   = nn_nchunks_k
+      snc4set%luse = ln_nc4zip
+#else
+      snc4set%luse = .FALSE.        ! No NetCDF 4 case
+#endif
+      !
+   END SUBROUTINE dom_nam
+
+   SUBROUTINE dom_zgr
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_zgr  ***
+      !!                   
+      !! ** Purpose :  set the depth of model levels and the resulting 
+      !!      vertical scale factors.
+      !!
+      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
+      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
+      !!              - vertical coordinate (gdep., e3.) depending on the 
+      !!                coordinate chosen :
+      !!                   ln_zco=T   z-coordinate  
+      !!                   ln_zps=T   z-coordinate with partial steps
+      !!                   ln_zco=T   s-coordinate 
+      !!
+      !! ** Action  :   define gdep., e3., mbathy and bathy
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio = 0   ! temporary integer
+      INTEGER ::   ios
+      !!
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
+      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
+      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzgr )
+
+      IF(lwp) THEN                     ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_zgr : vertical coordinate'
+         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
+         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
+         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
+         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
+         WRITE(numout,*) '             ice shelf cavity               ln_isfcav = ', ln_isfcav
+      ENDIF
+
+      ioptio = 0                       ! Check Vertical coordinate options
+      IF( ln_zco ) ioptio = ioptio + 1
+      IF( ln_zps ) ioptio = ioptio + 1
+      IF( ln_sco ) ioptio = ioptio + 1
+      IF( ln_isfcav ) ioptio = 33
+      IF ( ioptio /= 1  )   CALL ctl_stop( ' none or several vertical coordinate options used' )
+      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
+
+   END SUBROUTINE dom_zgr
+
+   SUBROUTINE dom_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_ctl  ***
+      !!
+      !! ** Purpose :   Domain control.
+      !!
+      !! ** Method  :   compute and print extrema of masked scale factors
+      !!
+      !! History :
+      !!   8.5  !  02-08  (G. Madec)    Original code
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
+      INTEGER, DIMENSION(2) ::   iloc      ! 
+      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
+      !!----------------------------------------------------------------------
+
+      ! Extrema of the scale factors
+
+      IF(lwp)WRITE(numout,*)
+      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
+      IF(lwp)WRITE(numout,*) '~~~~~~~'
+
+      IF (lk_mpp) THEN
+         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
+         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
+         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
+         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+      ELSE
+         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+
+         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi1 = iloc(1) + nimpp - 1
+         ijmi1 = iloc(2) + njmpp - 1
+         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi2 = iloc(1) + nimpp - 1
+         ijmi2 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima1 = iloc(1) + nimpp - 1
+         ijma1 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima2 = iloc(1) + nimpp - 1
+         ijma2 = iloc(2) + njmpp - 1
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
+         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
+         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
+         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
+      ENDIF
+
+   END SUBROUTINE dom_ctl
+
+   SUBROUTINE dom_grd
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_grd  ***
       !!                   
       !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
@@ -344,5 +741,5 @@
       CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw )
       !
-   END SUBROUTINE dom_rea
+   END SUBROUTINE dom_grd
 
 
@@ -388,4 +785,101 @@
    END SUBROUTINE zgr_bot_level
 
+   SUBROUTINE dom_msk
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk  ***
+      !!
+      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
+      !!
+      !! ** Method  :   The interior ocean/land mask is computed from tmask
+      !!              setting to zero the duplicated row and lines due to
+      !!              MPP exchange halos, est-west cyclic and north fold
+      !!              boundary conditions.
+      !!
+      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
+      !!               tpol     : ???
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
+      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
+      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
+      !
+      !!---------------------------------------------------------------------
+      
+      CALL wrk_alloc( jpi, jpj, imsk )
+      !
+      ! Interior domain mask (used for global sum)
+      ! --------------------
+      ssmask(:,:)  = tmask(:,:,1)
+      tmask_i(:,:) = tmask(:,:,1)
+      iif = jpreci                        ! thickness of exchange halos in i-axis
+      iil = nlci - jpreci + 1
+      ijf = jprecj                        ! thickness of exchange halos in j-axis
+      ijl = nlcj - jprecj + 1
+      !
+      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
+      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
+      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
+      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
+      !
+      !                                   ! north fold mask
+      tpol(1:jpiglo) = 1._wp
+      !                                
+      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
+      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
+      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
+         IF( mjg(ijl-1) == jpjglo-1 ) THEN
+            DO ji = iif+1, iil-1
+               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
+            END DO
+         ENDIF
+      ENDIF 
+      !
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
+      ! least 1 wet u point
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector loop
+            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+         END DO
+         DO ji = 1, jpim1      ! NO vector opt.
+            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
+               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
+         END DO
+      END DO
+      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( vmask_i, 'V', 1._wp )
+      CALL lbc_lnk( fmask_i, 'F', 1._wp )
+
+      ! 3. Ocean/land mask at wu-, wv- and w points 
+      !----------------------------------------------
+      wmask (:,:,1) = tmask(:,:,1) ! ????????
+      wumask(:,:,1) = umask(:,:,1) ! ????????
+      wvmask(:,:,1) = vmask(:,:,1) ! ????????
+      DO jk=2,jpk
+         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
+      END DO
+      !
+      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
+         imsk(:,:) = INT( tmask_i(:,:) )
+         WRITE(numout,*) ' tmask_i : '
+         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
+         WRITE (numout,*)
+         WRITE (numout,*) ' dommsk: tmask for each level'
+         WRITE (numout,*) ' ----------------------------'
+         DO jk = 1, jpk
+            imsk(:,:) = INT( tmask(:,:,jk) )
+            WRITE(numout,*)
+            WRITE(numout,*) ' level = ',jk
+            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
+         END DO
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, imsk )
+      !
+   END SUBROUTINE dom_msk
+
    !!======================================================================
 END MODULE domrea
+
Index: anches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/istate.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/istate.F90	(revision 5538)
+++ 	(revision )
@@ -1,51 +1,0 @@
-MODULE istate
-   !!======================================================================
-   !!                     ***  MODULE  istate  ***
-   !! Ocean state   :  initial state setting, off-line case
-   !!=====================================================================
-   !! History :  3.3  ! 2010-10  (C. Ethe)  original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   istate_init   : initial state set to zero
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers 
-   USE dom_oce         ! ocean space and time domain 
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   istate_init   ! routine called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!---------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE istate_init
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE istate_init  ***
-      !! 
-      !! ** Purpose :   Initialization to zero of the dynamics and tracers.
-      !!----------------------------------------------------------------------
-      !
-      !     now fields         !     after fields      !
-      un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   !
-      vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   !
-      wn   (:,:,:)   = 0._wp   !                       !
-      hdivn(:,:,:)   = 0._wp   !                       !
-      tsn  (:,:,:,:) = 0._wp   !                       !
-      !
-      rhd  (:,:,:) = 0.e0
-      rhop (:,:,:) = 0.e0
-      rn2  (:,:,:) = 0.e0 
-      !
-   END SUBROUTINE istate_init
-
-   !!=====================================================================
-END MODULE istate
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 5539)
@@ -18,6 +18,6 @@
    USE c1d             ! 1D configuration
    USE domcfg          ! domain configuration               (dom_cfg routine)
-   USE domain          ! domain initialization             (dom_init routine)
-   USE istate          ! initial state setting          (istate_init routine)
+   USE domain          ! domain initialization from coordinate & bathymetry (dom_init routine)
+   USE domrea          ! domain initialization from mesh_mask            (dom_init routine)
    USE eosbn2          ! equation of state            (eos bn2 routine)
    !              ! ocean physics
@@ -34,5 +34,4 @@
    USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
    USE dtadyn          ! Lecture and interpolation of the dynamical fields
-   USE stpctl          ! time stepping control            (stp_ctl routine)
    !              ! I/O & MPP
    USE iom             ! I/O library
@@ -95,5 +94,5 @@
       istp = nit000
       ! 
-      CALL iom_init( "nemo" )            ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
+      CALL iom_init( cxios_context )            ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
       ! 
       DO WHILE ( istp <= nitend .AND. nstop == 0 )    ! time stepping
@@ -108,5 +107,5 @@
       END DO
 #if defined key_iomput
-      CALL iom_context_finalize( "nemo" ) ! needed for XIOS+AGRIF
+      CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
 #endif
 
@@ -143,4 +142,5 @@
       INTEGER ::   ilocal_comm   ! local integer
       INTEGER ::   ios
+      LOGICAL ::   llexist
       CHARACTER(len=80), DIMENSION(16) ::   cltxt
       !!
@@ -182,5 +182,5 @@
       !                             !--------------------------------------------!
 #if defined key_iomput
-      CALL  xios_initialize( "nemo",return_comm=ilocal_comm )
+      CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )
       narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )   ! Nodes selection
 #else
@@ -269,5 +269,10 @@
       IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration
                             CALL     dom_cfg    ! Domain configuration
-                            CALL     dom_init   ! Domain
+      !
+      INQUIRE( FILE='coordinates.nc', EXIST = llexist )   ! Check if coordinate file exist
+      !
+      IF( llexist )  THEN  ;  CALL  dom_init   !  compute the grid from coordinates and bathymetry
+      ELSE                 ;  CALL  dom_rea    !  read grid from the meskmask
+      ENDIF
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
@@ -276,6 +281,6 @@
       IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
 
-      !                                     ! Ocean physics
                             CALL     sbc_init   ! Forcings : surface module
+
 #if ! defined key_degrad
                             CALL ldf_tra_init   ! Lateral ocean tracer physics
@@ -283,27 +288,15 @@
       IF( lk_ldfslp )       CALL ldf_slp_init   ! slope of lateral mixing
 
-      !                                     ! Active tracers
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
       IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
 
-                            CALL trc_nam_run  ! Needed to get restart parameters for passive tracers
-      IF( ln_rsttr ) THEN
-        neuler = 1   ! Set time-step indicator at nit000 (leap-frog)
-        CALL trc_rst_cal( nit000, 'READ' )   ! calendar
-      ELSE
-        neuler = 0                  ! Set time-step indicator at nit000 (euler)
-        CALL day_init               ! set calendar
-      ENDIF
-      !                                     ! Dynamics
+                            CALL trc_nam_run    ! Needed to get restart parameters for passive tracers
+                            CALL trc_rst_cal( nit000, 'READ' )   ! calendar
                             CALL dta_dyn_init   ! Initialization for the dynamics
 
-      !                                     ! Passive tracers
                             CALL     trc_init   ! Passive tracers initialization
-      !
-      ! Initialise diaptr as some variables are used in if statements later (in
-      ! various advection and diffusion routines.
-                            CALL dia_ptr_init
-      !
-      IF(lwp) WRITE(numout,cform_aaa)       ! Flag AAAAAAA
+                            CALL dia_ptr_init   ! Initialise diaptr as some variables are used 
+      !                                         ! in various advection and diffusion routines
+      IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA
       !
       IF( nn_timing == 1 )  CALL timing_stop( 'nemo_init')
@@ -660,4 +653,51 @@
    END SUBROUTINE nemo_northcomms
 #endif
+
+   SUBROUTINE istate_init
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_init  ***
+      !!
+      !! ** Purpose :   Initialization to zero of the dynamics and tracers.
+      !!----------------------------------------------------------------------
+      !
+      !     now fields         !     after fields      !
+      un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   !
+      vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   !
+      wn   (:,:,:)   = 0._wp   !                       !
+      hdivn(:,:,:)   = 0._wp   !                       !
+      tsn  (:,:,:,:) = 0._wp   !                       !
+      !
+      rhd  (:,:,:) = 0.e0
+      rhop (:,:,:) = 0.e0
+      rn2  (:,:,:) = 0.e0
+      !
+   END SUBROUTINE istate_init
+
+   SUBROUTINE stp_ctl( kt, kindic )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!
+      !! ** Actions :   'time.step' file containing the last ocean time-step
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT(inout) ::   kindic  ! indicator of solver convergence
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'stp_ctl : time-stepping control'
+         WRITE(numout,*) '~~~~~~~'
+         ! open time.step file
+         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+      ENDIF
+      !
+      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
+      IF(lwp) REWIND( numstp )                       ! --------------------------
+      !
+   END SUBROUTINE stp_ctl
    !!======================================================================
 END MODULE nemogcm
Index: anches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/stpctl.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OFF_SRC/stpctl.F90	(revision 5538)
+++ 	(revision )
@@ -1,62 +1,0 @@
-MODULE stpctl
-   !!======================================================================
-   !!                       ***  MODULE  stpctl  ***
-   !! Ocean run control :  Off-line case, only save the time step in numstp
-   !!======================================================================
-   !! History :  OPA  ! 1991-03  (G. Madec) Original code
-   !!            6.0  ! 1992-06  (M. Imbard)
-   !!            8.0  ! 1997-06  (A.M. Treguier)
-   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
-   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   stp_ctl      : Control the run
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! distributed memory computing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   stp_ctl    ! routine called by opa.F90
-   
-   !!----------------------------------------------------------------------
-   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE stp_ctl( kt, kindic )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE stp_ctl  ***
-      !!                     
-      !! ** Purpose :   Control the run
-      !!
-      !! ** Method  : - Save the time step in numstp
-      !!
-      !! ** Actions :   'time.step' file containing the last ocean time-step
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in   ) ::   kt      ! ocean time-step index
-      INTEGER, INTENT(inout) ::   kindic  ! indicator of solver convergence
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'stp_ctl : time-stepping control'
-         WRITE(numout,*) '~~~~~~~'
-         ! open time.step file
-         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-      ENDIF
-      !
-      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
-      IF(lwp) REWIND( numstp )                       !  --------------------------
-      !
-   END SUBROUTINE stp_ctl
-
-   !!======================================================================
-END MODULE stpctl
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 5539)
@@ -176,5 +176,5 @@
 
      !open output file
-     IF( lwp ) THEN
+     IF( lwm ) THEN
         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
@@ -283,5 +283,5 @@
            DO jsec=1,nb_sec
 
-              IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec))
+              IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
             
               !nullify transports values after writing
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 5539)
@@ -166,16 +166,12 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process: all this will come back in input files
-            ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
-            ! which had a grid-size of 362x292.
-            ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-            ! The original domain can still be used optionally if the isf code is not activated. 
-            ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
-            ! 
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-            ENDIF                                       ! Adjust j­indices to account for more southerly starting latitude
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
+            !
             ii0 = 283           ;   ii1 = 283
-            ij0 = 201 + isrow   ;   ij1 = 201 + isrow
+            ij0 = 241 - isrow   ;   ij1 = 241 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
@@ -223,11 +219,11 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            !  This dirty section will be suppressed by simplification process:
-            !  all this will come back in input files
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 39
-            ENDIF                                      
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
             ii0 = 282           ;   ii1 = 282
-            ij0 = 201 + isrow   ;   ij1 = 201 + isrow
+            ij0 = 240 - isrow   ;   ij1 = 240 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
@@ -275,11 +271,11 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            !  This dirty section will be suppressed by simplification process:
-            !  all this will come back in input files
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 39
-            ENDIF                                      
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
             ii0 = 331           ;   ii1 = 331
-            ij0 = 176 + isrow   ;   ij1 = 176 + isrow
+            ij0 = 215 - isrow   ;   ij1 = 215 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
@@ -327,11 +323,11 @@
          CASE ( 1 )                                  !  ORCA_R1 configurations
             !                                        ! =======================
-            !  This dirty section will be suppressed by simplification process:
-            !  all this will come back in input files
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 39
-            ENDIF                                      
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
             ii0 = 297           ;   ii1 = 297
-            ij0 = 230 + isrow   ;   ij1 = 230 + isrow
+            ij0 = 269 - isrow   ;   ij1 = 269 - isrow
             !                                        ! =======================
          CASE DEFAULT                                !    ORCA R05 or R025
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5539)
@@ -46,4 +46,6 @@
    USE iom
    USE ioipsl
+   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
+
 #if defined key_lim2
    USE limwri_2 
@@ -125,4 +127,5 @@
       !!
       INTEGER                      ::   ji, jj, jk              ! dummy loop indices
+      INTEGER                      ::   jkbot                   !
       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
       !!
@@ -148,4 +151,7 @@
          CALL iom_put( "e3w" , fse3w_n(:,:,:) )
       ENDIF
+
+      CALL iom_put( "ssh" , sshn )                 ! sea surface height
+      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       
       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
@@ -154,5 +160,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_tem)
+               jkbot = mbkt(ji,jj)
+               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
             END DO
          END DO
@@ -165,8 +172,25 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_sal)
+               jkbot = mbkt(ji,jj)
+               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
             END DO
          END DO
          CALL iom_put( "sbs", z2d )                ! bottom salinity
+      ENDIF
+
+      IF ( iom_use("taubot") ) THEN                ! bottom stress
+         z2d(:,:) = 0._wp
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
+                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )      
+               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
+                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
+               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
+               !
+            ENDDO
+         ENDDO
+         CALL lbc_lnk( z2d, 'T', 1. )
+         CALL iom_put( "taubot", z2d )           
       ENDIF
          
@@ -176,9 +200,15 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = un(ji,jj,MAX(mbathy(ji,jj),1))
+               jkbot = mbku(ji,jj)
+               z2d(ji,jj) = un(ji,jj,jkbot)
             END DO
          END DO
          CALL iom_put( "sbu", z2d )                ! bottom i-current
       ENDIF
+#if defined key_dynspg_ts
+      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
+#else
+      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
+#endif
       
       CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
@@ -187,8 +217,25 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj) = vn(ji,jj,MAX(mbathy(ji,jj),1))
+               jkbot = mbkv(ji,jj)
+               z2d(ji,jj) = vn(ji,jj,jkbot)
             END DO
          END DO
          CALL iom_put( "sbv", z2d )                ! bottom j-current
+      ENDIF
+#if defined key_dynspg_ts
+      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
+#else
+      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
+#endif
+
+      CALL iom_put( "woce", wn )                   ! vertical velocity
+      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
+         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
+         z2d(:,:) = rau0 * e12t(:,:)
+         DO jk = 1, jpk
+            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
+         END DO
+         CALL iom_put( "w_masstr" , z3d )  
+         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
       ENDIF
 
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 5539)
@@ -72,4 +72,5 @@
       !!----------------------------------------------------------------------
       INTEGER ::   jc            ! dummy loop indices
+      INTEGER :: isrow           ! local index
       !!----------------------------------------------------------------------
       
@@ -91,7 +92,13 @@
          CASE ( 1 )                                  ! ORCA_R1 configuration
             !                                        ! =======================
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
+            !
             ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea
-            ncsi1(1)   = 332  ; ncsj1(1)   = 203
-            ncsi2(1)   = 344  ; ncsj2(1)   = 235
+            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow
+            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow
             ncsir(1,1) = 1    ; ncsjr(1,1) = 1
             !                                        
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5539)
@@ -162,51 +162,47 @@
             !                                             ! =====================
             ! This dirty section will be suppressed by simplification process: all this will come back in input files
-            ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
             ! which had a grid-size of 362x292.
-            ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-            ! The original domain can still be used optionally if the isf code is not activated. 
-            ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
             ! 
-            IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-            ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-            ENDIF     
+            isrow = 332 - jpjglo
             !
             ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km)
-            ij0 = 201 + isrow   ;   ij1 = 201 + isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
+            ij0 = 201 + isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km'
 
             ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km)
-            ij0 = 208 + isrow   ;   ij1 = 208 + isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            ij0 = 208 + isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km'
 
             ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km)
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
+            ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km'
 
             ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
+            ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km'
 
             ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km)
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
+            ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km'
 
             ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km)
-            ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
+            ij0 = 124 + isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km'
 
             ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km)
-            ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
+            ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km'
 
             ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km)
-            ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
+            ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km'
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5539)
@@ -402,48 +402,44 @@
       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
          !                                                 ! Increased lateral friction near of some straits
-         ! This dirty section will be suppressed by simplification process: all this will come back in input files
-         ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
-         ! which had a grid-size of 362x292.
-         ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-         ! The original domain can still be used optionally if the isf code is not activated. 
-         ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
-         ! 
-         IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-         ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-         ENDIF     
-
+         ! This dirty section will be suppressed by simplification process:
+         ! all this will come back in input files
+         ! Currently these hard-wired indices relate to configuration with
+         ! extend grid (jpjglo=332)
+         !
+         isrow = 332 - jpjglo
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
          IF(lwp) WRITE(numout,*) '      Gibraltar '
          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait 
-         ij0 = 201 + isrow   ;   ij1 = 201 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 201 + isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Bhosporus '
          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait 
-         ij0 = 208 + isrow   ;   ij1 = 208 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 208 + isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Makassar (Top) '
          ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top) 
-         ij0 = 149 + isrow   ;   ij1 = 150 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+         ij0 = 149 + isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
 
          IF(lwp) WRITE(numout,*) '      Lombok '
          ii0 =  44           ;   ii1 =  44        ! Lombok Strait 
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Ombai '
          ii0 =  53           ;   ii1 =  53        ! Ombai Strait 
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      Timor Passage '
          ii0 =  56           ;   ii1 =  56        ! Timor Passage 
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
 
          IF(lwp) WRITE(numout,*) '      West Halmahera '
          ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait 
-         ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+         ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
 
          IF(lwp) WRITE(numout,*) '      East Halmahera '
          ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait 
-         ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+         ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
          !
       ENDIF
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5539)
@@ -1123,17 +1123,13 @@
       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
          !                                             ! =====================
-         ! This dirty section will be suppressed by simplification process: all this will come back in input files
-         ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration 
+         ! This dirty section will be suppressed by simplification process:
+         ! all this will come back in input files
+         ! Currently these hard-wired indices relate to configuration with
+         ! extend grid (jpjglo=332)
          ! which had a grid-size of 362x292.
-         ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6. 
-         ! The original domain can still be used optionally if the isf code is not activated. 
-         ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
-         ! 
-         IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
-         ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain. 
-         ENDIF     
+         isrow = 332 - jpjglo
          !
          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified)
-         ij0 = 201 + isrow   ;   ij1 = 201 + isrow
+         ij0 = 241 - isrow   ;   ij1 = 241 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1156,5 +1152,5 @@
          !
          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
-         ij0 = 208 + isrow   ;   ij1 = 208 + isrow
+         ij0 = 248 - isrow   ;   ij1 = 248 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1177,5 +1173,5 @@
          !
          ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified)
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow
+         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1193,5 +1189,5 @@
          !
          ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
-         ij0 = 124 + isrow   ;   ij1 = 125 + isrow
+         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1209,5 +1205,5 @@
          !
          ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified)
-         ij0 = 124 + isrow  ;   ij1 = 125  + isrow  
+         ij0 = 164 - isrow  ;   ij1 = 165  - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1224,6 +1220,6 @@
          END DO
          !
-         ii0 =  56             ;   ii1 =  56        ! Timor Passage (e1v was modified)
-         ij0 = 124  + isrow    ;   ij1 = 125  + isrow  
+         ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified)
+         ij0 = 164 - isrow    ;   ij1 = 165  - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1240,6 +1236,6 @@
          END DO
          !
-         ii0 =  55             ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
-         ij0 = 141  + isrow    ;   ij1 = 142  + isrow  
+         ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
+         ij0 = 181 - isrow    ;   ij1 = 182 - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
@@ -1256,6 +1252,6 @@
          END DO
          !
-         ii0 =  58             ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
-         ij0 = 141  + isrow    ;   ij1 = 142  + isrow  
+         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
+         ij0 = 181 - isrow    ;   ij1 = 182 - isrow  
          DO jk = 1, jpkm1
             DO jj = mj0(ij0), mj1(ij1)
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5539)
@@ -484,19 +484,4 @@
             CALL iom_close( inum )
             mbathy(:,:) = INT( bathy(:,:) )
-            !
-            ! CL :  add Amazon deeper
-            IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
-                 ii0 = 230   ;   ii1 = 245                   ! Amazon area 
-                 ij0 = 140   ;   ij1 = 155                   ! no ocean shallower than 30 meters
-                 DO ji = mi0(ii0), mi1(ii1)
-                    DO jj = mj0(ij0), mj1(ij1)
-                       IF( bathy(ji,jj) .LE. 30. .AND. bathy(ji,jj) .GT. 0.0 ) bathy(ji,jj) = 30._wp
-                    END DO
-                 END DO
-                 IF(lwp) WRITE(numout,*)
-                 IF(lwp) WRITE(numout,*) '             orca_r1: Amazon area not shallower than 30 meters for: '
-                 IF(lwp) WRITE(numout,*) '                    Longitude index ',ii0, ii0
-                 IF(lwp) WRITE(numout,*) '                    Latitude index  ',ij0, ij0
-            ENDIF
             !                                                ! =====================
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5539)
@@ -237,5 +237,5 @@
       !
       CALL wrk_dealloc( jpi  , jpj+2, zwu               )
-      CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 )
+      CALL wrk_dealloc( jpi+4, jpj  , zwv, kistart = -1 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('div_cur')
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5539)
@@ -266,6 +266,6 @@
                ! Add volume filter correction: compatibility with tracer advection scheme
                ! => time filter + conservation correction (only at the first level)
-               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-            !
+               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
+                              &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
             ENDIF
             !
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5539)
@@ -21,5 +21,4 @@
    USE domvvl          ! Variable volume
    USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE iom             ! I/O library
    USE restart         ! only for lrst_oce
    USE in_out_manager  ! I/O manager
@@ -31,5 +30,4 @@
    USE bdy_par         
    USE bdydyn2d        ! bdy_ssh routine
-   USE iom
 #if defined key_agrif
    USE agrif_opa_update
@@ -137,6 +135,4 @@
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "ssh" , sshn )   ! sea surface height
-      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       !
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
@@ -228,22 +224,4 @@
 #endif
       !
-      !                                           !------------------------------!
-      !                                           !           outputs            !
-      !                                           !------------------------------!
-      CALL iom_put( "woce", wn )   ! vertical velocity
-      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
-         CALL wrk_alloc( jpi, jpj, z2d ) 
-         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
-         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e12t(:,:)
-         DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
-         END DO
-         CALL iom_put( "w_masstr" , z3d )  
-         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-         CALL wrk_dealloc( jpi, jpj, z2d  ) 
-         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
-      ENDIF
-      !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
 
@@ -290,5 +268,5 @@
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
-         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * ssmask(:,:)
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) - rnf_b(:,:) + rnf(:,:) ) * ssmask(:,:)
          sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5539)
@@ -91,5 +91,6 @@
    INTEGER ::   nitrst                !: time step at which restart file should be written
    LOGICAL ::   lrst_oce              !: logical to control the oce restart write 
-   INTEGER ::   numror, numrow        !: logical unit for cean restart (read and write)
+   INTEGER ::   numror = 0            !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90)
+   INTEGER ::   numrow                !: logical unit for ocean restart (write)
    INTEGER ::   nrst_lst              !: number of restart to output next
 
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 5539)
@@ -663,4 +663,5 @@
 
       ! --- evaporation minus precipitation --- !
+      zsnw(:,:) = 0._wp
       CALL lim_thd_snwblow( pfrld, zsnw )          ! snow redistribution by wind
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * ( 1._wp - zsnw )
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 5539)
@@ -613,4 +613,5 @@
 
       ! --- evaporation minus precipitation --- !
+      zsnw(:,:) = 0._wp
       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5539)
@@ -780,5 +780,5 @@
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
       !
-      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'S_QsrOce' ) + cpl_freq( 'S_QsrMix' )
+      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
       IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
          &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
@@ -1379,5 +1379,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ! for LIM3
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
       !!----------------------------------------------------------------------
       !
@@ -1524,4 +1524,5 @@
 
       ! --- heat flux associated with emp --- !
+      zsnw(:,:) = 0._wp
       CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
       zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
@@ -1556,5 +1557,4 @@
 
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
-
 #else
 
@@ -1624,4 +1624,17 @@
          ENDDO
       ENDIF
+
+#if defined key_lim3
+      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
+      ! --- solar flux over ocean --- !
+      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
+      zqsr_oce = 0._wp
+      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
+
+      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
+      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
+
+      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
+#endif
 
       IF( ln_mixcpl ) THEN
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 5539)
@@ -509,5 +509,5 @@
       CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
 
-      CALL wrk_dealloc( jpi,jpj, ztmp )
+      CALL wrk_dealloc( jpi,jpj, ztmp, zpice )
       CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
       !
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5539)
@@ -216,8 +216,8 @@
             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
                                  CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
+            ! clem: evap_ice is forced to 0 in coupled mode for now 
+            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models
+            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp
             IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
-            ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
-            evap_ice  (:,:,:) = 0._wp
-            devap_ice (:,:,:) = 0._wp
          END SELECT
          CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5539)
@@ -53,4 +53,5 @@
    USE timing           ! Timing
    USE sbcwave          ! Wave module
+   USE bdy_par          ! Require lk_bdy
 
    IMPLICIT NONE
@@ -343,5 +344,7 @@
       !                                            ! ---------------------------------------- !
       !
-      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
+      IF ( .NOT. lk_bdy ) then
+         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
+      ENDIF
                                                          ! (caution called before sbc_ssm)
       !
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5539)
@@ -209,5 +209,5 @@
       zfact = 0.5_wp
       !
-      IF( ln_rnf_depth ) THEN      !==   runoff distributed over several levels   ==!
+      IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==!
          IF( lk_vvl ) THEN             ! variable volume case
             DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90	(revision 5539)
@@ -112,7 +112,4 @@
    INTEGER         :: nn_trc_ord = 1          ! order of autoregressive processes
 
-   ! Public array with density correction
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: drho_ran
-
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -263,5 +260,4 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
 
-      IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF')
       ! Read namsto namelist : stochastic parameterization
       REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
@@ -341,4 +337,5 @@
       ENDIF
       IF( ln_sto_eos ) THEN
+         IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF')
          IF(lwp) WRITE(numout,*) '       - stochastic equation of state'
          ALLOCATE(jsto_eosi(nn_sto_eos))
@@ -668,7 +665,4 @@
       ! ------------------------------------------
       IF( ln_rststo ) CALL sto_rst_read
-
-      ! Allocate drho_ran
-      ALLOCATE(drho_ran(jpi,jpj,jpk))
 
    END SUBROUTINE sto_par_init
@@ -781,6 +775,4 @@
                CALL iom_rstput( kt, nitrst, numstow, clsto3d , sto3d(:,:,:,jsto) )
             END DO
-            ! Save drho_ran in restart file
-            CALL iom_rstput( kt, nitrst, numstow, 'drho' , drho_ran(:,:,:) )
             ! close the restart file
             CALL iom_close( numstow )
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 5539)
@@ -27,4 +27,5 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
    USE zdf_oce         ! ocean vertical mixing
    USE domvvl          ! variable volume
@@ -278,5 +279,5 @@
 
       !!     
-      LOGICAL  ::   ll_tra_hpg, ll_traqsr   ! local logical
+      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
       REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
@@ -293,7 +294,9 @@
          ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
          ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
+         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
       ELSE                          
          ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
          ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
+         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
       ENDIF
       !
@@ -319,18 +322,23 @@
                   !
                   IF( jk == 1 ) THEN           ! first level 
-                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
+                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) )
                      ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
                   ENDIF
+
                   IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
                      &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
 
-                   ze3t_f = 1.e0 / ze3t_f
-                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
-                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
-                   !
-                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
-                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
-                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
-                   ENDIF
+                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   & 		      ! river runoffs
+                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
+                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
+
+                  ze3t_f = 1.e0 / ze3t_f
+                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
+                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
+                  !
+                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
+                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
+                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
+                  ENDIF
                END DO
             END DO
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5539)
@@ -101,9 +101,10 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-
+      IF( lk_bdy     )  THEN
+         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
+                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+      ENDIF
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5539)
@@ -27,4 +27,5 @@
    USE sbc_oce          ! surface boundary condition: ocean
    USE sbctide          ! Tide initialisation
+   USE sbcapr           ! surface boundary condition: ssh_ib required by bdydta 
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90	(revision 5539)
@@ -121,4 +121,5 @@
 
    LOGICAL ::   linit = .FALSE.
+   LOGICAL ::   ldebug = .FALSE.
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -486,4 +487,5 @@
       
       IF( SUM( tree(ii)%ishape ) == 0 ) THEN                    ! create a new branch 
+         IF(ldebug) PRINT *, 'create new branch ', ii,ishape, isrt, itype
          tree(ii)%itype = itype                                        ! define the type of this branch 
          tree(ii)%ishape(:) = ishape(:)                                ! define the shape of this branch 
@@ -515,4 +517,5 @@
          tree(ii)%current%next%in_use = .FALSE.                        ! this leaf is not yet used
          tree(ii)%current%next%indic = tree(ii)%current%indic + 1      ! number of this leaf
+         IF(ldebug) PRINT *, 'add a leaf ', ii, tree(ii)%current%indic
          tree(ii)%current%next%prev => tree(ii)%current                ! previous leaf of the new leaf is the current leaf
          tree(ii)%current%next%next => NULL()                          ! next leaf is not yet defined
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/daymod.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/daymod.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/daymod.F90	(revision 5539)
@@ -80,15 +80,8 @@
       ndt05   = NINT(0.5 * rdttra(1))
 
-      ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)
-      ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
-      adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
-      IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error
-      !
-      IF(lwp) THEN
-         WRITE(numout,*) ' *** Info used values : '
-         WRITE(numout,*) '   date ndastp                                      : ', ndastp
-         WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
-         WRITE(numout,*)
-      ENDIF
+      ! ==> clem: here we read the ocean restart for the date (only if it exists)
+      !           It is not clean and another solution should be found
+      CALL day_rst( nit000, 'READ' )
+      ! ==>
 
       ! set the calendar from ndastp (read in restart file and namelist)
@@ -285,4 +278,98 @@
       !
    END SUBROUTINE day
+
+
+   SUBROUTINE day_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ts_rst  ***
+      !!
+      !!  ** Purpose : Read or write calendar in restart file:
+      !!
+      !!  WRITE(READ) mode:
+      !!       kt        : number of time step since the begining of the experiment at the
+      !!                   end of the current(previous) run
+      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
+      !!                   end of the current(previous) run (REAL -> keep fractions of day)
+      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
+      !!
+      !!   According to namelist parameter nrstdt,
+      !!       nrstdt = 0  no control on the date (nit000 is  arbitrary).
+      !!       nrstdt = 1  we verify that nit000 is equal to the last
+      !!                   time step of previous run + 1.
+      !!       In both those options, the  exact duration of the experiment
+      !!       since the beginning (cumulated duration of all previous restart runs)
+      !!       is not stored in the restart and is assumed to be (nit000-1)*rdt.
+      !!       This is valid is the time step has remained constant.
+      !!
+      !!       nrstdt = 2  the duration of the experiment in days (adatrj)
+      !!                    has been stored in the restart file.
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      !
+      REAL(wp) ::   zkt, zndastp
+      !!----------------------------------------------------------------------
+
+      IF( TRIM(cdrw) == 'READ' ) THEN
+
+         IF( iom_varid( numror, 'kt', ldstop = .FALSE. ) > 0 ) THEN
+            ! Get Calendar informations
+            CALL iom_get( numror, 'kt', zkt )   ! last time-step of previous run
+            IF(lwp) THEN
+               WRITE(numout,*) ' *** Info read in restart : '
+               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
+               WRITE(numout,*) ' *** restart option'
+               SELECT CASE ( nrstdt )
+               CASE ( 0 )   ;   WRITE(numout,*) ' nrstdt = 0 : no control of nit000'
+               CASE ( 1 )   ;   WRITE(numout,*) ' nrstdt = 1 : no control the date at nit000 (use ndate0 read in the namelist)'
+               CASE ( 2 )   ;   WRITE(numout,*) ' nrstdt = 2 : calendar parameters read in restart'
+               END SELECT
+               WRITE(numout,*)
+            ENDIF
+            ! Control of date
+            IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 )                                         &
+                 &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart',                 &
+                 &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' )
+            ! define ndastp and adatrj
+            IF ( nrstdt == 2 ) THEN
+               ! read the parameters correspondting to nit000 - 1 (last time step of previous run)
+               CALL iom_get( numror, 'ndastp', zndastp )
+               ndastp = NINT( zndastp )
+               CALL iom_get( numror, 'adatrj', adatrj  )
+            ELSE
+               ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day)
+               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
+               adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
+               ! note this is wrong if time step has changed during run
+            ENDIF
+         ELSE
+            ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day)
+            ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
+            adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
+         ENDIF
+         IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) ' *** Info used values : '
+            WRITE(numout,*) '   date ndastp                                      : ', ndastp
+            WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
+            WRITE(numout,*)
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+         !
+         IF( kt == nitrst ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'rst_write : write oce restart file  kt =', kt
+            IF(lwp) WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         ! calendar control
+         CALL iom_rstput( kt, nitrst, numrow, 'kt'     , REAL( kt    , wp) )   ! time-step
+         CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) )   ! date
+         CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj            )   ! number of elapsed days since
+         !                                                                     ! the begining of the run [s]
+      ENDIF
+      !
+   END SUBROUTINE day_rst
    !!======================================================================
 END MODULE daymod
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90	(revision 5539)
@@ -59,4 +59,9 @@
    USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges
    USE icbstp          ! handle bergs, calving, themodynamics and transport
+#if defined key_bdy
+   USE bdyini          ! open boundary cond. setting       (bdy_init routine). clem: mandatory for LIM3
+   USE bdydta          ! open boundary cond. setting   (bdy_dta_init routine). clem: mandatory for LIM3
+#endif
+   USE bdy_par
 
    IMPLICIT NONE
@@ -354,4 +359,11 @@
 
                             CALL sbc_init   ! Forcings : surface module 
+                            
+      ! ==> clem: open boundaries init. is mandatory for LIM3 because ice BDY is not decoupled from  
+      !           the environment of ocean BDY. Therefore bdy is called in both OPA and SAS modules. 
+      !           This is not clean and should be changed in the future. 
+      IF( lk_bdy        )   CALL     bdy_init
+      IF( lk_bdy        )   CALL bdy_dta_init
+      ! ==>
       
       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
@@ -502,5 +514,10 @@
       USE diawri    , ONLY: dia_wri_alloc
       USE dom_oce   , ONLY: dom_oce_alloc
-      USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
+#if defined key_bdy   
+      USE bdy_oce   , ONLY: bdy_oce_alloc
+      USE oce         ! clem: mandatory for LIM3 because needed for bdy arrays
+#else
+      USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
+#endif
       !
       INTEGER :: ierr,ierr1,ierr2,ierr3,ierr4,ierr5,ierr6
@@ -510,11 +527,16 @@
       ierr =        dia_wri_alloc   ()
       ierr = ierr + dom_oce_alloc   ()          ! ocean domain
-      ALLOCATE( snwice_mass(jpi,jpj)  , snwice_mass_b(jpi,jpj),             &
-         &      snwice_fmass(jpi,jpj), STAT= ierr1 )
+#if defined key_bdy
+      ierr = ierr + bdy_oce_alloc   ()          ! bdy masks (incl. initialization)
+      ierr = ierr + oce_alloc       ()          ! (tsn...)
+#endif
+
+#if ! defined key_bdy
+       ALLOCATE( snwice_mass(jpi,jpj)  , snwice_mass_b(jpi,jpj),             &
+         &      snwice_fmass(jpi,jpj)  , STAT= ierr1 )
       !
       ! lim code currently uses surface temperature and salinity in tsn array for initialisation
-      ! and ub, vb arrays in ice dynamics
-      ! so allocate enough of arrays to use
-      !
+      ! and ub, vb arrays in ice dynamics, so allocate enough of arrays to use
+      ! clem: should not be needed. To be checked out
       jpm = MAX(jp_tem, jp_sal)
       ALLOCATE( tsn(jpi,jpj,1,jpm)  , STAT=ierr2 )
@@ -523,6 +545,6 @@
       ALLOCATE( tsb(jpi,jpj,1,jpm)  , STAT=ierr5 )
       ALLOCATE( sshn(jpi,jpj)       , STAT=ierr6 )
-
       ierr = ierr + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 + ierr6 
+#endif
       !
       IF( lk_mpp    )   CALL mpp_sum( ierr )
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/step.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/step.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/SAS_SRC/step.F90	(revision 5539)
@@ -37,4 +37,9 @@
 
    USE timing           ! Timing            
+
+   USE bdy_par          ! clem: mandatory for LIM3
+#if defined key_bdy
+   USE bdydta           ! clem: mandatory for LIM3
+#endif
 
    IMPLICIT NONE
@@ -81,4 +86,12 @@
                              CALL iom_setkt( kstp - nit000 + 1, cxios_context )   ! tell iom we are at time step kstp
 
+      ! ==> clem: open boundaries is mandatory for LIM3 because ice BDY is not decoupled from  
+      !           the environment of ocean BDY. Therefore bdy is called in both OPA and SAS modules.
+      !           From SAS: ocean bdy data are wrong  (but we do not care) and ice bdy data are OK.  
+      !           This is not clean and should be changed in the future. 
+#if defined key_bdy
+      IF( lk_bdy     )       CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+#endif
+      ! ==>
                              CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
 
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5539)
@@ -117,5 +117,9 @@
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
             CALL fld_read( kt, 1, sf_dust )
-            dust(:,:) = sf_dust(1)%fnow(:,:,1)
+            IF( nn_ice_tr == -1 .AND. .NOT. ln_ironice ) THEN
+               dust(:,:) = sf_dust(1)%fnow(:,:,1)
+            ELSE
+               dust(:,:) = sf_dust(1)%fnow(:,:,1) * ( 1.0 - fr_i(:,:) )
+            ENDIF
          ENDIF
       ENDIF
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 5539)
@@ -185,5 +185,6 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !
-      INTEGER :: ji, jj, jk, jn, jl, jc                     ! dummy loop indicesa
+      INTEGER :: ji , jj, jk, jn, jl, jc                     ! dummy loop indicesa
+      INTEGER :: isrow                                      ! local index
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta       ! 3D  workspace
 
@@ -201,4 +202,13 @@
             !
             SELECT CASE ( jp_cfg )
+            !                                           ! =======================
+            CASE ( 1 )                                  ! eORCA_R1 configuration
+            !                                           ! =======================
+            isrow = 332 - jpjglo
+            !
+                                                        ! Caspian Sea
+            nctsi1(1)   = 332  ; nctsj1(1)   = 243 - isrow
+            nctsi2(1)   = 344  ; nctsj2(1)   = 275 - isrow
+            !                                        
             !                                           ! =======================
             CASE ( 2 )                                  !  ORCA_R2 configuration
Index: /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 5539)
@@ -207,35 +207,32 @@
          ENDIF
 
-         CALL iom_open( TRIM(cn_trcrst_indir)//'/'//cn_trcrst_in, numrtr, kiolib = jlibalt )
-
-         CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
-
-         IF(lwp) THEN
-            WRITE(numout,*) ' *** Info read in restart : '
-            WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
-            WRITE(numout,*) ' *** restart option'
-            SELECT CASE ( nn_rsttr )
-            CASE ( 0 )   ;   WRITE(numout,*) ' nn_rsttr = 0 : no control of nittrc000'
-            CASE ( 1 )   ;   WRITE(numout,*) ' nn_rsttr = 1 : no control the date at nittrc000 (use ndate0 read in the namelist)'
-            CASE ( 2 )   ;   WRITE(numout,*) ' nn_rsttr = 2 : calendar parameters read in restart'
-            END SELECT
-            WRITE(numout,*)
-         ENDIF
-         ! Control of date 
-         IF( nittrc000  - NINT( zkt ) /= nn_dttrc .AND.  nn_rsttr /= 0 )                                  &
-            &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 &
-            &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' )
-         IF( lk_offline ) THEN      ! set the date in offline mode
-            ! Check dynamics and tracer time-step consistency and force Euler restart if changed
-            IF( iom_varid( numrtr, 'rdttrc1', ldstop = .FALSE. ) > 0 )   THEN
-               CALL iom_get( numrtr, 'rdttrc1', zrdttrc1 )
-               IF( zrdttrc1 /= rdt * nn_dttrc )   neuler = 0
-            ENDIF
-            !                          ! define ndastp and adatrj
-            IF( nn_rsttr == 2 ) THEN
+         IF( ln_rsttr ) THEN
+            CALL iom_open( TRIM(cn_trcrst_indir)//'/'//cn_trcrst_in, numrtr, kiolib = jlibalt )
+            CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
+
+            IF(lwp) THEN
+               WRITE(numout,*) ' *** Info read in restart : '
+               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
+               WRITE(numout,*) ' *** restart option'
+               SELECT CASE ( nn_rsttr )
+               CASE ( 0 )   ;   WRITE(numout,*) ' nn_rsttr = 0 : no control of nittrc000'
+               CASE ( 1 )   ;   WRITE(numout,*) ' nn_rsttr = 1 : no control the date at nittrc000 (use ndate0 read in the namelist)'
+               CASE ( 2 )   ;   WRITE(numout,*) ' nn_rsttr = 2 : calendar parameters read in restart'
+               END SELECT
+               WRITE(numout,*)
+            ENDIF
+            ! Control of date 
+            IF( nittrc000  - NINT( zkt ) /= nn_dttrc .AND.  nn_rsttr /= 0 )                                  &
+               &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 &
+               &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' )
+         ENDIF
+         !
+         IF( lk_offline ) THEN    
+            !                                          ! set the date in offline mode
+            IF( ln_rsttr .AND. nn_rsttr == 2 ) THEN
                CALL iom_get( numrtr, 'ndastp', zndastp ) 
                ndastp = NINT( zndastp )
                CALL iom_get( numrtr, 'adatrj', adatrj  )
-            ELSE
+             ELSE
                ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam
                adatrj = ( REAL( nittrc000-1, wp ) * rdttra(1) ) / rday
@@ -248,4 +245,8 @@
               WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
               WRITE(numout,*)
+            ENDIF
+            !
+            IF( ln_rsttr )  THEN   ;    neuler = 1
+            ELSE                   ;    neuler = 0
             ENDIF
             !
Index: /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-XC40_METO
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-XC40_METO	(revision 5539)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-XC40_METO	(revision 5539)
@@ -0,0 +1,95 @@
+#/bin/bash 
+#!
+#PBS -N nemo_sette
+#PBS -l walltime=00:30:00
+#PBS -j oe
+#PBS -q QUEUE
+#PBS -l SELECT
+
+  export PBS_O_WORKDIR=$(readlink -f $PBS_O_WORKDIR)
+  export OMP_NUM_THREADS=1
+  cd $PBS_O_WORKDIR
+  export XIO_HOME=/projects/ocean/nemo/xios/xios_1.0_r618_20150619/XIOS
+#
+  echo " ";
+  export OMP_NUM_THREADS=1
+  O_PER_NODE=32
+  X_PER_NODE=4
+  OCORES=NPROCS
+  XCORES=NXIOPROCS
+  if [ $OCORES -le 32 ] ; then O_PER_NODE=$OCORES; fi
+  export SETTE_DIR=DEF_SETTE_DIR
+
+###############################################################
+#
+#
+# load sette functions (only post_test_tidyup needed)
+#
+  . ${SETTE_DIR}/all_functions.sh
+###############################################################
+#
+# Don't remove neither change the following line
+# BODY
+#
+# Test specific settings. Do not hand edit these lines; the fcm_job.sh script will set these
+# (via sed operating on this template job file). Note that the number of compute nodes required
+# is also set by the fcm_job.sh on the PBS select header line above.
+#
+# These variables are needed by post_test_tidyup function in all_functions.sh
+#
+  export INPUT_DIR=DEF_INPUT_DIR
+  export CONFIG_DIR=DEF_CONFIG_DIR
+  export NEMO_VALIDATION_DIR=DEF_NEMO_VALIDATION
+  export NEW_CONF=DEF_NEW_CONF
+  export CMP_NAM=DEF_CMP_NAM
+  export TEST_NAME=DEF_TEST_NAME
+  export EXE_DIR=DEF_EXE_DIR
+  ulimit -c unlimited
+  ulimit -s unlimited
+  export FORT_FMT_RECL=132
+#
+# end of set up
+###############################################################
+#
+# change to the working directory 
+#
+  cd $EXE_DIR
+  echo Directory is `pwd`
+  
+  if [ $XCORES -gt 0 ]; then
+#
+#  Run MPMD case
+#
+     #XIOS will run on a separate node so will run in parallel queue
+     if [ ! -f ./xios_server.exe ] && [ -f ${XIO_HOME}/bin/xios_server.exe ]; then
+        cp ${XIO_HOME}/bin/xios_server.exe .
+     fi
+     if [ ! -f ./xios_server.exe ]; then
+        echo "./xios_server.exe not found"
+        echo "run aborted"
+        exit
+     fi
+       echo time aprun -b -n $XCORES -N $X_PER_NODE ./xios_server.exe : -n $OCORES -N $O_PER_NODE ./opa
+            time aprun -b -n $XCORES -N $X_PER_NODE ./xios_server.exe : -n $OCORES -N $O_PER_NODE ./opa
+#
+  else
+#
+# Run SPMD case
+#
+    if [ $OCORES -gt 32 ] ; then
+        #Still more than one node so run in parallel queue
+        echo time aprun -b  -n $OCORES -N $O_PER_NODE ./opa
+             time aprun -b  -n $OCORES -N $O_PER_NODE ./opa
+    else
+       echo time mpiexec -n $OCORES ./opa
+            time mpiexec -n $OCORES ./opa
+    fi
+  fi
+#
+
+#
+  post_test_tidyup
+# END_BODY
+# Don't remove neither change the previous line
+  exit
+
Index: /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/prepare_job.sh
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/prepare_job.sh	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/prepare_job.sh	(revision 5539)
@@ -226,5 +226,56 @@
                                         # round up the number of nodes required.
                                         #
-                                        NB_NODES=$( echo $NB_NODES $NXIO_PROC | awk '{print ($1 + ( $2 / 4 ) + 1)}')                    
+                                        NB_NODES=$( echo $NB_NODES $NXIO_PROC | awk '{print ($1 + ( $2 / 4 ) + 1)}') 
+                                fi
+                                ;;
+                        XC40_METO*) #Setup for Met Office XC40 with any compiler
+                                # ocean cores are packed 32 to a node
+                                # If we need more than one node then have to use parallel queue and XIOS must have a node to itself
+                                NB_REM=$( echo $NB_PROC | awk '{print ( $1 % 32 ) }')
+                                if [ ${NB_REM} == 0 ] ; then
+                                        # number of processes required is an integer multiple of 32
+                                        #
+                                        NB_NODES=$( echo $NB_PROC $NXIO_PROC | awk '{print ($1) / 32}')
+                                else
+                                        #
+                                        # number of processes required is not an integer multiple of 32
+                                        # round up the number of nodes required.
+                                        #
+                                        NB_NODES=$( echo $NB_PROC $NXIO_PROC | awk '{printf("%d",($1) / 32 + 1 )}')
+                                fi
+                                # xios cores are sparsely packed at 4 to a node
+                                if [ $NXIO_PROC == 0 ] ; then
+                                    NB_XNODES=0
+                                else
+                                    NB_REM=$( echo $NXIO_PROC | awk '{print ( $1 % 4 ) }')
+                                    if [ ${NB_REM} == 0 ] ; then
+                                            # number of processes required is an integer multiple of 4                           
+                                            #
+                                            NB_XNODES=$( echo $NXIO_PROC | awk '{print (( $1 / 4 ) + 1)}') 
+                                    else
+                                            #
+                                            # number of processes required is not an integer multiple of 4                             
+                                            # round up the number of nodes required.
+                                            #
+                                            NB_XNODES=$( echo $NXIO_PROC | awk '{printf("%d",($1) / 4 + 1) }')                    
+                                    fi
+                                fi
+                                if [ ${NB_XNODES} -ge 1 ] ; then
+                                   NB_NODES=$((NB_NODES+NB_XNODES))
+                                fi
+                                echo NB_XNODES=${NB_XNODES} 
+                                echo Total NB_NODES=${NB_NODES}
+                                if [ ${NB_NODES} -eq 1 ] ; then
+                                   QUEUE=shared
+                                   #Not using XIOS in detatched mode and using less than one node so should be ok on shared node
+                                   #Load snplauncher module to allow use of mpiexec
+                                   SELECT="select=1:ncpus=$((NXIO_PROC + NB_PROC))":mem=15GB
+                                   module load cray-snplauncher
+                                   echo 'Shared Queue'
+                                else
+                                   QUEUE=normal
+                                   SELECT="select=$NB_NODES"
+                                   module unload cray-snplauncher #Make sure snplauncher module is not loaded
+                                   echo 'Normal Queue'
                                 fi
                                 ;;
@@ -279,4 +330,9 @@
                     cat run_sette_test.job | sed -e"s/NPROC_NODE/${NB_PROC_NODE}/" \
                                                  -e"s:QUEUE:${QUEUE}:" > run_sette_test1.job
+                    mv run_sette_test1.job run_sette_test.job
+                    ;;
+              XC40_METO*)
+                    cat run_sette_test.job | sed -e"s/QUEUE/${QUEUE}/" \
+                                                 -e"s/SELECT/${SELECT}/" > run_sette_test1.job
                     mv run_sette_test1.job run_sette_test.job
                     ;;
Index: /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/sette.sh
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/sette.sh	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/SETTE/sette.sh	(revision 5539)
@@ -857,4 +857,6 @@
     export TEST_NAME="REPRO_4_8"
     . ./prepare_exe_dir.sh
+    JOB_FILE=${EXE_DIR}/run_job.sh
+    if [ -f ${JOB_FILE} ] ; then \rm ${JOB_FILE} ; fi
     cd ${EXE_DIR}
     set_namelist namelist_cfg nn_it000 1
Index: /branches/UKMO/restart_datestamp/NEMOGCM/TOOLS/MISCELLANEOUS/chk_wrk_alloc.sh
===================================================================
--- /branches/UKMO/restart_datestamp/NEMOGCM/TOOLS/MISCELLANEOUS/chk_wrk_alloc.sh	(revision 5538)
+++ /branches/UKMO/restart_datestamp/NEMOGCM/TOOLS/MISCELLANEOUS/chk_wrk_alloc.sh	(revision 5539)
@@ -18,4 +18,6 @@
 for ff in $( grep -il "^ *use  *wrk_nemo" $( find . -name "*90" )  $( find . -name "*h90" ) )
 do
+    ierr=0
+    
     # number of lines with wrk_alloc
     n1=$( grep -ic "call *wrk_alloc *(" $ff )
@@ -28,15 +30,48 @@
 	n2=$( sed -e "s/wrk_alloc/wrk_dealloc/" $ff | grep -ic "call *wrk_dealloc *(" )
 	# we should get n2 = 2 * n1...
-	[ $(( 2 * $n1 )) -ne $n2 ] && echo "problem with wrk_alloc in $ff" 
-	
+	if [ $(( 2 * $n1 )) -ne $n2 ]
+	then
+	    ierr=1
+	    echo "problem with wrk_alloc in $ff" 
+	fi
 	# same story but for wrk_dealloc
 	nn2=$( sed -e "s/wrk_dealloc/wrk_alloc/" $ff | grep -ic "call *wrk_alloc *(" )
 	if [ $(( 2 * $nn1 )) -ne $nn2 ]
 	then
-	    echo "problem with wrk_dealloc in $ff"
+	    ierr=1
+	    echo "problem with wrk_dealloc in $ff" 
+	fi
+
+	if [ $ierr -eq 0 ] # check that wrk_alloc block is the same as wrk_dealloc block
+	then
+	    grep -i "call *wrk_alloc *("   $ff | sed -e "s/ //g" | sed -e "s/!.*//g" > txt1$$
+	    grep -i "call *wrk_dealloc *(" $ff | sed -e "s/wrk_dealloc/wrk_alloc/"  | sed -e "s/ //g" | sed -e "s/!.*//g" > txt2$$
+	    cmp txt1$$ txt2$$
+	    if [ $? -ne 0 ]
+	    then
+		echo "different syntax in wrk_alloc and wrk_dealloc in $ff"
+		echo
+		for ll in $( seq 1 $n1 )  # compare each line
+		do
+		    sed -n ${ll}p txt1$$ > ll1$$
+		    sed -n ${ll}p txt2$$ > ll2$$
+		    cmp ll1$$ ll2$$ > /dev/null
+		    if [ $? -ne 0 ]
+		    then
+			grep -i "call *wrk_alloc *("   $ff | sed -n ${ll}p
+			grep -i "call *wrk_dealloc *(" $ff | sed -n ${ll}p
+			echo
+		    fi
+		    rm -f ll1$$ ll2$$
+		done
+	    fi
+	    rm -f txt1$$ txt2$$
+	else
 	    grep -i "call *wrk_alloc *(" $ff
-	    grep -i "call *wrk_dealloc *(" $ff
+	    echo
+	    grep -i "call *wrk_dealloc *(" $ff	
 	    echo
 	fi
+	
     fi
     
