Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90	(revision 6060)
@@ -26,13 +26,13 @@
    !!
    !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE iom             !
-   USE ioipsl, ONLY :   ymds2ju   ! for calendar
-   USE prtctl          ! Print control
-   USE trc_oce, ONLY : lk_offline ! offline flag
-   USE timing          ! Timing
-   USE restart         ! restart
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE iom            !
+   USE ioipsl  , ONLY :   ymds2ju   ! for calendar
+   USE prtctl         ! Print control
+   USE trc_oce , ONLY : lk_offline ! offline flag
+   USE timing         ! Timing
+   USE restart        ! restart
 
    IMPLICIT NONE
@@ -43,5 +43,5 @@
    PUBLIC   day_mth    ! Needed by TAM
 
-   INTEGER, PUBLIC ::   nsecd, nsecd05, ndt, ndt05 ! (PUBLIC for TAM)
+   INTEGER, PUBLIC ::   nsecd, nsecd05, ndt, ndt05   !: (PUBLIC for TAM)
 
    !!----------------------------------------------------------------------
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6060)
@@ -11,4 +11,5 @@
    !!                             to the optimization of BDY communications
    !!            3.7  ! 2015-11  (G. Madec) introduce surface and scale factor ratio
+   !!             -   ! 2015-11  (G. Madec, A. Coward)  time varying zgr by default
    !!----------------------------------------------------------------------
 
@@ -47,6 +48,6 @@
    !! Free surface parameters
    !! =======================
-   LOGICAL, PUBLIC :: ln_dynspg_exp      !: Explicit free surface flag
-   LOGICAL, PUBLIC :: ln_dynspg_ts       !: Split-Explicit free surface flag
+   LOGICAL , PUBLIC :: ln_dynspg_exp     !: Explicit free surface flag
+   LOGICAL , PUBLIC :: ln_dynspg_ts      !: Split-Explicit free surface flag
 
    !! Time splitting parameters
@@ -61,5 +62,5 @@
    !! Horizontal grid parameters for domhgr
    !! =====================================
-   INTEGER       ::   jphgr_msh        !: type of horizontal mesh
+   INTEGER       ::   jphgr_msh          !: type of horizontal mesh
    !                                       !  = 0 curvilinear coordinate on the sphere read in coordinate.nc
    !                                       !  = 1 geographical mesh on the sphere with regular grid-spacing
@@ -68,29 +69,29 @@
    !                                       !  = 4 Mercator grid with T/U point at the equator
 
-   REAL(wp)      ::   ppglam0              !: longitude of first raw and column T-point (jphgr_msh = 1)
-   REAL(wp)      ::   ppgphi0              !: latitude  of first raw and column T-point (jphgr_msh = 1)
+   REAL(wp)      ::   ppglam0            !: longitude of first raw and column T-point (jphgr_msh = 1)
+   REAL(wp)      ::   ppgphi0            !: latitude  of first raw and column T-point (jphgr_msh = 1)
    !                                                        !  used for Coriolis & Beta parameters (jphgr_msh = 2 or 3)
-   REAL(wp)      ::   ppe1_deg             !: zonal      grid-spacing (degrees)
-   REAL(wp)      ::   ppe2_deg             !: meridional grid-spacing (degrees)
-   REAL(wp)      ::   ppe1_m               !: zonal      grid-spacing (degrees)
-   REAL(wp)      ::   ppe2_m               !: meridional grid-spacing (degrees)
+   REAL(wp)      ::   ppe1_deg           !: zonal      grid-spacing (degrees)
+   REAL(wp)      ::   ppe2_deg           !: meridional grid-spacing (degrees)
+   REAL(wp)      ::   ppe1_m             !: zonal      grid-spacing (degrees)
+   REAL(wp)      ::   ppe2_m             !: meridional grid-spacing (degrees)
 
    !! Vertical grid parameter for domzgr
    !! ==================================
-   REAL(wp)      ::   ppsur                !: ORCA r4, r2 and r05 coefficients
-   REAL(wp)      ::   ppa0                 !: (default coefficients)
-   REAL(wp)      ::   ppa1                 !:
-   REAL(wp)      ::   ppkth                !:
-   REAL(wp)      ::   ppacr                !:
+   REAL(wp)      ::   ppsur              !: ORCA r4, r2 and r05 coefficients
+   REAL(wp)      ::   ppa0               !: (default coefficients)
+   REAL(wp)      ::   ppa1               !:
+   REAL(wp)      ::   ppkth              !:
+   REAL(wp)      ::   ppacr              !:
    !
    !  If both ppa0 ppa1 and ppsur are specified to 0, then
    !  they are computed from ppdzmin, pphmax , ppkth, ppacr in dom_zgr
-   REAL(wp)      ::   ppdzmin              !: Minimum vertical spacing
-   REAL(wp)      ::   pphmax               !: Maximum depth
+   REAL(wp)      ::   ppdzmin            !: Minimum vertical spacing
+   REAL(wp)      ::   pphmax             !: Maximum depth
    !
-   LOGICAL       ::   ldbletanh            !: Use/do not use double tanf function for vertical coordinates
-   REAL(wp)      ::   ppa2                 !: Double tanh function parameters
-   REAL(wp)      ::   ppkth2               !:
-   REAL(wp)      ::   ppacr2               !:
+   LOGICAL       ::   ldbletanh          !: Use/do not use double tanf function for vertical coordinates
+   REAL(wp)      ::   ppa2               !: Double tanh function parameters
+   REAL(wp)      ::   ppkth2             !:
+   REAL(wp)      ::   ppacr2             !:
 
    !                                    !! old non-DOCTOR names still used in the model
@@ -106,7 +107,8 @@
    REAL(wp), PUBLIC ::   rdth            !: depth variation of tracer step
 
-   !                                                  !!! associated variables
-   INTEGER , PUBLIC                 ::   neuler        !: restart euler forward option (0=Euler)
-   REAL(wp), PUBLIC                 ::   atfp1         !: asselin time filter coeff. (atfp1= 1-2*atfp)
+   !                                   !!! associated variables
+   INTEGER , PUBLIC ::   neuler          !: restart euler forward option (0=Euler)
+   REAL(wp), PUBLIC ::   atfp1           !: asselin time filter coeff. (atfp1= 1-2*atfp)
+   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdttra  !: vertical profile of tracer time step
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dtra  !: = 2*rdttra except at nit000 (=rdttra) if neuler=0
@@ -177,53 +179,37 @@
    !! vertical coordinate and scale factors
    !! ---------------------------------------------------------------------
-   !                                 !!* Namelist namzgr : vertical coordinate *
-   LOGICAL, PUBLIC ::   ln_zco        !: z-coordinate - full step
-   LOGICAL, PUBLIC ::   ln_zps        !: z-coordinate - partial step
-   LOGICAL, PUBLIC ::   ln_sco        !: s-coordinate or hybrid z-s coordinate
-   LOGICAL, PUBLIC ::   ln_isfcav     !: presence of ISF 
-
-   !! All coordinates
-   !! ---------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
-#if defined key_vvl
-   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .TRUE.    !: variable grid flag
-
-   !! All coordinates
-   !! ---------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_n           !: now depth of T-points (sum of e3w) (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_n, gdepw_n   !: now depth at T-W  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_b, gdepw_b   !: before depth at T-W  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_n              !: now    vertical scale factors at  t       point  (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_n  , e3v_n     !:            -      -      -    -   u --v   points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_n  , e3f_n     !:            -      -      -    -   w --f   points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_n , e3vw_n    !:            -      -      -    -   uw--vw  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_b              !: before     -      -      -    -   t       points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_b              !: before     -      -      -    -   t       points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_b  , e3v_b     !:   -        -      -      -    -   u --v   points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_b , e3vw_b    !:   -        -      -      -    -   uw--vw  points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_a              !: after      -      -      -    -   t       point  (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_a  , e3v_a     !:   -        -      -      -    -   u --v   points (m)
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
-#endif
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur  , hvr     !: Now    inverse of u and v-points ocean depth (1/m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu   , hv      !:        depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht             !:        depth at t-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehur_a, ehvr_a !: After  inverse of u and v-points ocean depth (1/m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehu_a , ehv_a  !:        depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehur_b, ehvr_b !: Before inverse of u and v-points ocean depth (1/m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehu_b , ehv_b  !:        depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0           !: reference depth at t-       points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters)
+   !                                !!* Namelist namzgr : vertical coordinate *
+   LOGICAL, PUBLIC ::   ln_zco       !: z-coordinate - full step
+   LOGICAL, PUBLIC ::   ln_zps       !: z-coordinate - partial step
+   LOGICAL, PUBLIC ::   ln_sco       !: s-coordinate or hybrid z-s coordinate
+   LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF 
+   LOGICAL, PUBLIC ::   ln_linssh    !: variable grid flag
+
+   !                                                        !  ref.   ! before  !   now   ! after  !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+
+   !                                                        !  ref.   ! before  !   now   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+   
+   !                                                      !  ref. ! before  !   now   !  after  !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0            ,    ht_n             !: t-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: u-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
 
-   !! z-coordinate with full steps (also used in the other cases as reference z-coordinate)
+   !! 1D reference  vertical coordinate
    !! =-----------------====------
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
@@ -231,4 +217,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3tp    , e3wp     !: ocean bottom level thickness at T and W points
 
+!!gm  This should be removed from here....  ==>>> only used in domzgr at initialization phase
    !! s-coordinate and hybrid z-s-coordinate
    !! =----------------======---------------
@@ -244,4 +231,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu       !: and quasi-uniform spacing             t--u points (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1                !: Maximum grid stiffness ratio
+!!gm end
 
    !!----------------------------------------------------------------------
@@ -249,9 +237,7 @@
    !! ---------------------------------------------------------------------
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters)
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt, mbku, mbkv   !: vertical index of the bottom last T-, U- & V ocean level
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                              !: land/ocean mask of barotropic stream function
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level                (ISF)
@@ -352,24 +338,22 @@
          &        ff (jpi,jpj)                                                         , STAT=ierr(3) )
          !
-      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         &
-         &      gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         &
-         &      gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) )
-         !
-#if defined key_vvl
-      ALLOCATE( gdep3w_n(jpi,jpj,jpk) , e3t_n (jpi,jpj,jpk) , e3u_n (jpi,jpj,jpk) ,                           &
-         &      gdept_n (jpi,jpj,jpk) , e3v_n (jpi,jpj,jpk) , e3w_n (jpi,jpj,jpk) ,                           &
-         &      gdepw_n (jpi,jpj,jpk) , e3f_n (jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , e3uw_n(jpi,jpj,jpk) ,     &
-         &      e3t_b   (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) ,                           &
-         &      e3uw_b  (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,                                                 &
-         &      gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) ,                           &
-         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) ,                           &
-         &      ehu_a   (jpi,jpj)     , ehv_a (jpi,jpj),                                                     &
-         &      ehur_a  (jpi,jpj)     , ehvr_a(jpi,jpj),                                                     &
-         &      ehu_b   (jpi,jpj)     , ehv_b (jpi,jpj),                                                     &
-         &      ehur_b  (jpi,jpj)     , ehvr_b(jpi,jpj),                                  STAT=ierr(5) )                          
-#endif
-         !
-      ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , ht_0(jpi,jpj) ,     &
-         &      hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , ht  (jpi,jpj) , STAT=ierr(6)  )
+      ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,     &
+         &      gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) ,                             &
+         &      gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) )
+         !
+      ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) ,   &
+         &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   &                         &
+         &      e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) ,   &                         &
+         &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             &
+         !                                                          !
+         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         &
+         &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &               
+         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
+         !
+      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           &
+         &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     &
+         &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     &
+         &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  )
+         !
          !
       ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) ,                                     &
@@ -386,5 +370,4 @@
       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      &
          &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), &
-         &     bmask  (jpi,jpj) ,                                                       &
          &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )
 
@@ -405,3 +388,2 @@
    !!======================================================================
 END MODULE dom_oce
-
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6060)
@@ -13,4 +13,5 @@
    !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
    !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
+   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
    !!----------------------------------------------------------------------
    
@@ -36,4 +37,5 @@
    !
    USE in_out_manager  ! I/O manager
+   USE wrk_nemo        ! Memory Allocation
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
@@ -45,6 +47,4 @@
    PUBLIC   dom_init   ! called by opa.F90
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!-------------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -70,6 +70,7 @@
       !!              - 1D configuration, move Coriolis, u and v at T-point
       !!----------------------------------------------------------------------
-      INTEGER ::   jk          ! dummy loop argument
+      INTEGER ::   jk          ! dummy loop indices
       INTEGER ::   iconf = 0   ! local integers
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0
       !!----------------------------------------------------------------------
       !
@@ -82,15 +83,17 @@
       ENDIF
       !
-                             CALL dom_nam      ! read namelist ( namrun, namdom )
-                             CALL dom_clo      ! Closed seas and lake
-                             CALL dom_hgr      ! Horizontal mesh
-                             CALL dom_zgr      ! Vertical mesh and bathymetry
-                             CALL dom_msk      ! Masks
-      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
-      !
-      ht_0(:,:) = 0._wp                        ! Reference ocean depth at T-points
-      hu_0(:,:) = 0._wp                        ! Reference ocean depth at U-points
-      hv_0(:,:) = 0._wp                        ! Reference ocean depth at V-points
-      DO jk = 1, jpk
+      !                       !==  Reference coordinate system  ==!
+      !
+                     CALL dom_nam               ! read namelist ( namrun, namdom )
+                     CALL dom_clo               ! Closed seas and lake
+                     CALL dom_hgr               ! Horizontal mesh
+                     CALL dom_zgr               ! Vertical mesh and bathymetry
+                     CALL dom_msk               ! Masks
+      IF( ln_sco )   CALL dom_stiff             ! Maximum stiffness ratio/hydrostatic consistency
+      !
+      ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
+      hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
+      hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpk
          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
          hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
@@ -98,24 +101,45 @@
       END DO
       !
-      IF( lk_vvl         )   CALL dom_vvl_init ! Vertical variable mesh
-      !
-      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
-      !
-      !
-      hu(:,:) = 0._wp                          ! Ocean depth at U-points
-      hv(:,:) = 0._wp                          ! Ocean depth at V-points
-      ht(:,:) = 0._wp                          ! Ocean depth at T-points
-      DO jk = 1, jpkm1
-         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
-
-                             CALL dom_stp      ! time step
-      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
-      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
+      !              !==  time varying part of coordinate system  ==!
+      !
+      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all
+         !       before        !          now          !       after         !
+         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
+         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
+         ;                     ;   gde3w_n = gde3w_0   !        ---          !
+         !                                                                  
+         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
+         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
+         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
+         ;                     ;     e3f_n =   e3f_0   !        ---          !
+         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
+         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
+         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
+         !
+         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
+         !
+         z1_hu_0(:,:) = umask_i(:,:) / ( hu_0(:,:) + 1._wp - umask_i(:,:) )     ! _i mask due to ISF
+         z1_hv_0(:,:) = vmask_i(:,:) / ( hv_0(:,:) + 1._wp - vmask_i(:,:) )
+         !
+         !        before       !          now          !       after         !
+         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
+         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   ! 
+         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
+         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
+         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
+         !
+         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
+         !
+      ELSE                         ! time varying : initialize before/now/after variables
+         !
+         CALL dom_vvl_init 
+         !
+      ENDIF
+      !
+      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
+      !
+                             CALL dom_stp       ! time step
+      IF( nmsh /= 0      )   CALL dom_wri       ! Create a domain file
+      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
       !
       IF( nn_timing == 1 )   CALL timing_stop('dom_init')
@@ -403,5 +427,5 @@
       INTEGER  ::   ji, jj, jk 
       REAL(wp) ::   zrxmax
-      REAL(wp), DIMENSION(4) :: zr1
+      REAL(wp), DIMENSION(4) ::   zr1
       !!----------------------------------------------------------------------
       rx1(:,:) = 0._wp
@@ -412,22 +436,22 @@
          DO jj = 2, jpjm1
             DO jk = 1, jpkm1
-               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
-                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
-                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
-                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
-               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
-                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
-                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
-                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
-               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
-                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
-                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
-                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
-               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
-                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
-                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
-                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
-               zrxmax = MAXVAL(zr1(1:4))
-               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
+               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
+                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
+                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
+                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
+               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
+                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
+                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
+                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
+               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
+                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
+                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
+                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
+               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
+                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
+                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
+                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
+               zrxmax = MAXVAL( zr1(1:4) )
+               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
             END DO
          END DO
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90	(revision 6060)
@@ -118,22 +118,24 @@
          WRITE(numout,*) '          south-west indices    jpizoom = ', jpizoom,   &
             &                                           ' jpjzoom = ', jpjzoom
-         WRITE(numout,*)
-         WRITE(numout,*) '          conversion local  ==> data i-index domain'
-         WRITE(numout,25)              (mig(ji),ji = 1,jpi)
-         WRITE(numout,*)
-         WRITE(numout,*) '          conversion data   ==> local  i-index domain'
-         WRITE(numout,*) '             starting index'
-         WRITE(numout,25)              (mi0(ji),ji = 1,jpidta)
-         WRITE(numout,*) '             ending index'
-         WRITE(numout,25)              (mi1(ji),ji = 1,jpidta)
-         WRITE(numout,*)
-         WRITE(numout,*) '          conversion local  ==> data j-index domain'
-         WRITE(numout,25)              (mjg(jj),jj = 1,jpj)
-         WRITE(numout,*)
-         WRITE(numout,*) '          conversion data  ==> local j-index domain'
-         WRITE(numout,*) '             starting index'
-         WRITE(numout,25)              (mj0(jj),jj = 1,jpjdta)
-         WRITE(numout,*) '             ending index'
-         WRITE(numout,25)              (mj1(jj),jj = 1,jpjdta)
+         IF( nn_print >= 1 ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) '          conversion local  ==> data i-index domain'
+            WRITE(numout,25)              (mig(ji),ji = 1,jpi)
+            WRITE(numout,*)
+            WRITE(numout,*) '          conversion data   ==> local  i-index domain'
+            WRITE(numout,*) '             starting index'
+            WRITE(numout,25)              (mi0(ji),ji = 1,jpidta)
+            WRITE(numout,*) '             ending index'
+            WRITE(numout,25)              (mi1(ji),ji = 1,jpidta)
+            WRITE(numout,*)
+            WRITE(numout,*) '          conversion local  ==> data j-index domain'
+            WRITE(numout,25)              (mjg(jj),jj = 1,jpj)
+            WRITE(numout,*)
+            WRITE(numout,*) '          conversion data  ==> local j-index domain'
+            WRITE(numout,*) '             starting index'
+            WRITE(numout,25)              (mj0(jj),jj = 1,jpjdta)
+            WRITE(numout,*) '             ending index'
+            WRITE(numout,25)              (mj1(jj),jj = 1,jpjdta)
+         ENDIF
       ENDIF
  25   FORMAT( 100(10x,19i4,/) )
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6060)
@@ -348,5 +348,5 @@
       e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
 
-      IF( lwp .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
+      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
          WRITE(numout,*)
          WRITE(numout,*) '          longitude and e1 scale factors'
@@ -391,6 +391,6 @@
          !
 #if defined key_agrif
-         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
-            IF( .NOT. Agrif_Root() ) THEN
+         IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN       ! for EEL6 configuration only
+            IF( .NOT.Agrif_Root() ) THEN
               zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
             ENDIF
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6060)
@@ -7,6 +7,5 @@
    !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
    !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
-   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
-   !!                 !                      pression of the double computation of bmask
+   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask
    !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
    !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
@@ -25,7 +24,8 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp
+   USE lib_mpp         !
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
@@ -34,5 +34,5 @@
    PRIVATE
 
-   PUBLIC   dom_msk         ! routine called by inidom.F90
+   PUBLIC   dom_msk    ! routine called by inidom.F90
 
    !                            !!* Namelist namlbc : lateral boundary condition *
@@ -89,16 +89,9 @@
       !!
       !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
-      !!      are defined with the proper value at lateral domain boundaries,
-      !!      but bmask. indeed, bmask defined the domain over which the
-      !!      barotropic stream function is computed. this domain cannot
-      !!      contain identical columns because the matrix associated with
-      !!      the barotropic stream function equation is then no more inverti-
-      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
-      !!      even IF nperio is not zero.
+      !!      are defined with the proper value at lateral domain boundaries.
       !!
       !!      In case of open boundaries (lk_bdy=T):
       !!        - tmask is set to 1 on the points to be computed bay the open
       !!          boundaries routines.
-      !!        - bmask is  set to 0 on the open boundaries.
       !!
       !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
@@ -107,6 +100,4 @@
       !!               fmask    : land/ocean mask at f-point (=0. or 1.)
       !!                          =rn_shlat along lateral boundaries
-      !!               bmask    : land/ocean mask at barotropic stream
-      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
       !!               tmask_i  : interior ocean mask
       !!----------------------------------------------------------------------
@@ -254,50 +245,4 @@
       END DO
 
-      ! 4. ocean/land mask for the elliptic equation
-      ! --------------------------------------------
-      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
-      !
-      !                               ! Boundary conditions
-      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
-      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
-         bmask( 1 ,:) = 0._wp
-         bmask(jpi,:) = 0._wp
-      ENDIF
-      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
-         bmask(:, 1 ) = 0._wp
-      ENDIF
-      !                                    ! north fold : 
-      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
-         DO ji = 1, jpi                      
-            ii = ji + nimpp - 1
-            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
-            bmask(ji,jpj  ) = 0._wp
-         END DO
-      ENDIF
-      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
-         bmask(:,jpj) = 0._wp
-      ENDIF
-      !
-      IF( lk_mpp ) THEN                    ! mpp specificities
-         !                                      ! bmask is set to zero on the overlap region
-         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
-         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
-         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
-         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
-         !
-         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
-            DO ji = 1, nlci
-               ii = ji + nimpp - 1
-               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
-               bmask(ji,nlcj  ) = 0._wp
-            END DO
-         ENDIF
-         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
-            DO ji = 1, nlci
-               bmask(ji,nlcj  ) = 0._wp
-            END DO
-         ENDIF
-      ENDIF
-
       ! Lateral boundary conditions on velocity (modify fmask)
       ! ---------------------------------------     
@@ -399,51 +344,6 @@
       !
       CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
-
+      !
       ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
-            
-      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
-         WRITE(numout,*)
-         WRITE(numout,*) ' dom_msk: vmask for each level'
-         WRITE(numout,*) ' -----------------------------'
-         DO jk = 1, jpk
-            imsk(:,:) = INT( vmask(:,:,jk) )
-            WRITE(numout,*)
-            WRITE(numout,*) ' level = ',jk
-            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
-               &                              1, jpj, 1, 1, numout)
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) ' dom_msk: fmask for each level'
-         WRITE(numout,*) ' -----------------------------'
-         DO jk = 1, jpk
-            imsk(:,:) = INT( fmask(:,:,jk) )
-            WRITE(numout,*)
-            WRITE(numout,*) ' level = ',jk
-            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
-               &                              1, jpj, 1, 1, numout )
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) ' dom_msk: bmask '
-         WRITE(numout,*) ' ---------------'
-         WRITE(numout,*)
-         imsk(:,:) = INT( bmask(:,:) )
-         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
-            &                              1, jpj, 1, 1, numout )
-      ENDIF
       !
       CALL wrk_dealloc( jpi, jpj, imsk )
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90	(revision 6060)
@@ -22,6 +22,4 @@
    PUBLIC   dom_stp   ! routine called by inidom.F90
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -84,5 +82,5 @@
             IF(lwp) WRITE(numout,*)'               accelerating the convergence'
             IF(lwp) WRITE(numout,*)'               dynamics time step = ', rdt/3600., ' hours'
-            IF( ln_sco .AND. rdtmin /= rdtmax .AND. lk_vvl )   &
+            IF( ln_sco .AND. rdtmin /= rdtmax .AND. .NOT.ln_linssh )   &
                  & CALL ctl_stop ( ' depth dependent acceleration of convergence not implemented in s-coordinates &
                  &                   nor in variable volume' )
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6060)
@@ -20,9 +20,11 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
+   USE phycst          ! physical constant
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce         ! ocean surface boundary condition
+   USE restart         ! ocean restart
+   !
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O manager library
-   USE restart         ! ocean restart
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -60,12 +62,10 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 
+   !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -74,5 +74,5 @@
       !!                ***  FUNCTION dom_vvl_alloc  ***
       !!----------------------------------------------------------------------
-      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
+      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
@@ -81,6 +81,6 @@
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
-         un_td = 0.0_wp
-         vn_td = 0.0_wp
+         un_td = 0._wp
+         vn_td = 0._wp
       ENDIF
       IF( ln_vvl_ztilde ) THEN
@@ -89,5 +89,5 @@
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
       ENDIF
-
+      !
    END FUNCTION dom_vvl_alloc
 
@@ -103,11 +103,11 @@
       !!               - interpolate scale factors
       !!
-      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
-      !!              - Regrid: fse3(u/v)_n
-      !!                        fse3(u/v)_b       
-      !!                        fse3w_n           
-      !!                        fse3(u/v)w_b      
-      !!                        fse3(u/v)w_n      
-      !!                        fsdept_n, fsdepw_n and fsde3w_n
+      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
+      !!              - Regrid: e3(u/v)_n
+      !!                        e3(u/v)_b       
+      !!                        e3w_n           
+      !!                        e3(u/v)w_b      
+      !!                        e3(u/v)w_n      
+      !!                        gdept_n, gdepw_n and gde3w_n
       !!              - h(t/u/v)_0
       !!              - frq_rst_e3t and frq_rst_hdv
@@ -115,108 +115,110 @@
       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      USE phycst,  ONLY : rpi, rsmall, rad
-      !! * Local declarations
-      INTEGER ::   ji,jj,jk
+      INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
       REAL(wp)::   zcoef
       !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_init')
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-
-      ! choose vertical coordinate (z_star, z_tilde or layer)
-      ! ==========================
-      CALL dom_vvl_ctl
-
-      ! Allocate module arrays
-      ! ======================
+      !
+      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
+      !
+      !                    ! Allocate module arrays
       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
-
-      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
-      ! ============================================================================
+      !
+      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
       CALL dom_vvl_rst( nit000, 'READ' )
-      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
-
-      ! Reconstruction of all vertical scale factors at now and before time steps
-      ! =============================================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
-      ! Vertical scale factor interpolations
-      ! ------------------------------------
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
-      ! t- and w- points depth
-      ! ----------------------
-      ! set the isf depth as it is in the initial step
-      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
-      fsdepw_n(:,:,1) = 0.0_wp
-      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
-      fsdepw_b(:,:,1) = 0.0_wp
-
-      DO jk = 2, jpk
+      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      !
+      !                    !== Set of all other vertical scale factors  ==!  (now and before)
+      !                                ! Horizontal interpolation of e3t
+      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
+      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
+      !                                ! Vertical interpolation of e3t,u,v 
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
+      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
+      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
+      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      !
+      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
+      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
+      gdepw_n(:,:,1) = 0.0_wp
+      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
+      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
+      gdepw_b(:,:,1) = 0.0_wp
+      DO jk = 2, jpk                               ! vertical sum
          DO jj = 1,jpj
             DO ji = 1,jpi
-              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
-                                                     ! 0.5 where jk = mikt  
-               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
-               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
-               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
-               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
-               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk)) 
+               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
+               !                             ! 0.5 where jk = mikt     
+!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
+               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
+               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
+                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
+               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
+               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
+               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
+                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
             END DO
          END DO
       END DO
-
-      ! Before depth and Inverse of the local depth of the water column at u- and v- points
-      ! -----------------------------------------------------------------------------------
-      hu_b(:,:) = 0.
-      hv_b(:,:) = 0.
-      DO jk = 1, jpkm1
-         hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
-         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+      !
+      !                    !==  thickness of the water column  !!   (ocean portion only)
+      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
+      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
+      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
+      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
+      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
+         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
+         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
       END DO
-      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
-      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
-
-      ! Restoring frequencies for z_tilde coordinate
-      ! ============================================
+      !
+      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
+      r1_hu_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )    ! _i mask due to ISF
+      r1_hu_n(:,:) = umask_i(:,:) / ( hu_n(:,:) + 1._wp - umask_i(:,:) )
+      r1_hv_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
+      r1_hv_n(:,:) = vmask_i(:,:) / ( hv_n(:,:) + 1._wp - vmask_i(:,:) )
+
+      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
       IF( ln_vvl_ztilde ) THEN
-         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
-         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
-         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
-         IF( ln_vvl_ztilde_as_zstar ) THEN
-            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
-            frq_rst_e3t(:,:) = 0.0_wp 
-            frq_rst_hdv(:,:) = 1.0_wp / rdt
-         ENDIF
-         IF ( ln_vvl_zstar_at_eqtor ) THEN
+!!gm : idea: add here a READ in a file of custumized restoring frequency
+         !                                   ! Values in days provided via the namelist
+         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
+         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
+         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
+         !
+         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
+            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
+            frq_rst_hdv(:,:) = 1._wp / rdt
+         ENDIF
+         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
             DO jj = 1, jpj
                DO ji = 1, jpi
+!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
                   IF( ABS(gphit(ji,jj)) >= 6.) THEN
                      ! values outside the equatorial band and transition zone (ztilde)
                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
-                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
+                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
                      ! values inside the equatorial band (ztilde as zstar)
                      frq_rst_e3t(ji,jj) =  0.0_wp
                      frq_rst_hdv(ji,jj) =  1.0_wp / rdt
-                  ELSE
-                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
+                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
+                     !                                      ! (linearly transition from z-tilde to z-star)
                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
@@ -229,6 +231,6 @@
                END DO
             END DO
-            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
-               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
+            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
+               ii0 = 103   ;   ii1 = 111       
                ij0 = 128   ;   ij1 = 135   ;   
                frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
@@ -237,7 +239,7 @@
          ENDIF
       ENDIF
-
+      !
       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
-
+      !
    END SUBROUTINE dom_vvl_init
 
@@ -261,25 +263,27 @@
       !!               - tilde_e3t_a: after increment of vertical scale factor 
       !!                              in z_tilde case
-      !!               - fse3(t/u/v)_a
+      !!               - e3(t/u/v)_a
       !!
       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
-      !! * Arguments
-      INTEGER, INTENT( in )                  :: kt                    ! time step
-      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
-      !! * Local declarations
-      INTEGER                                :: ji, jj, jk            ! dummy loop indices
-      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
-      REAL(wp)                               :: z2dt                  ! temporary scalars
-      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
-      LOGICAL                                :: ll_do_bclinic         ! temporary logical
-      !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
-      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
-
-      IF(kt == nit000)   THEN
+      INTEGER, INTENT( in )           ::   kt      ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
+      !
+      INTEGER                ::   ji, jj, jk            ! dummy loop indices
+      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
+      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
+      LOGICAL                ::   ll_do_bclinic         ! local logical
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze3t
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zht, z_scale, zwu, zwv, zhdiv
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_sf_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,zht,   z_scale, zwu, zwv, zhdiv )
+      CALL wrk_alloc( jpi,jpj,jpk,   ze3t )
+
+      IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
@@ -289,5 +293,5 @@
       ll_do_bclinic = .TRUE.
       IF( PRESENT(kcall) ) THEN
-         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
+         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
       ENDIF
 
@@ -295,18 +299,16 @@
       ! After acale factors at t-points !
       ! ******************************* !
-
       !                                                ! --------------------------------------------- !
-                                                       ! z_star coordinate and barotropic z-tilde part !
+      !                                                ! z_star coordinate and barotropic z-tilde part !
       !                                                ! --------------------------------------------- !
-
+      !
       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
       DO jk = 1, jpkm1
-         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
-         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
+         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
       END DO
-
+      !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
          !                                                            ! ------baroclinic part------ !
-
          ! I - initialization
          ! ==================
@@ -314,9 +316,9 @@
          ! 1 - barotropic divergence
          ! -------------------------
-         zhdiv(:,:) = 0.
-         zht(:,:)   = 0.
+         zhdiv(:,:) = 0._wp
+         zht(:,:)   = 0._wp
          DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
+            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
@@ -325,15 +327,15 @@
          ! --------------------------------------------------
          IF( ln_vvl_ztilde ) THEN
-            IF( kt .GT. nit000 ) THEN
+            IF( kt > nit000 ) THEN
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
                END DO
             ENDIF
-         END IF
+         ENDIF
 
          ! II - after z_tilde increments of vertical scale factors
          ! =======================================================
-         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
+         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
 
          ! 1 - High frequency divergence term
@@ -341,11 +343,11 @@
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
-            END DO
-         END IF
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+            END DO
+         ENDIF
 
          ! 2 - Restoring term (z-tilde case only)
@@ -355,12 +357,11 @@
                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
             END DO
-         END IF
+         ENDIF
 
          ! 3 - Thickness diffusion term
          ! ----------------------------
-         zwu(:,:) = 0.0_wp
-         zwv(:,:) = 0.0_wp
-         ! a - first derivative: diffusive fluxes
-         DO jk = 1, jpkm1
+         zwu(:,:) = 0._wp
+         zwv(:,:) = 0._wp
+         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -374,6 +375,5 @@
             END DO
          END DO
-         ! b - correction for last oceanic u-v points
-         DO jj = 1, jpj
+         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
             DO ji = 1, jpi
                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
@@ -381,6 +381,5 @@
             END DO
          END DO
-         ! c - second derivative: divergence of diffusive fluxes
-         DO jk = 1, jpkm1
+         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -391,6 +390,6 @@
             END DO
          END DO
-         ! d - thickness diffusion transport: boundary conditions
-         !     (stored for tracer advction and continuity equation)
+         !                       ! d - thickness diffusion transport: boundary conditions
+         !                             (stored for tracer advction and continuity equation)
          CALL lbc_lnk( un_td , 'U' , -1._wp)
          CALL lbc_lnk( vn_td , 'V' , -1._wp)
@@ -410,5 +409,5 @@
          ! Maximum deformation control
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ze3t(:,:,jpk) = 0.0_wp
+         ze3t(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
@@ -462,5 +461,5 @@
          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
          END DO
 
@@ -470,5 +469,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -485,25 +484,25 @@
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
+            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
+            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
          END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
          !
          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
@@ -524,6 +523,6 @@
       ! *********************************** !
 
-      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
+      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
 
       ! *********************************** !
@@ -531,19 +530,20 @@
       ! *********************************** !
 
-      hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
-      hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
-      DO jk = 1, jpkm1
-         hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
+      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
+      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
+         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
-      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
-
-      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
-
+!!gm BUG ?  don't understand the use of umask_i here .....
+      r1_hu_a(:,:) = umask_i(:,:) / ( hu_a(:,:) + 1._wp - umask_i(:,:) )
+      r1_hv_a(:,:) = vmask_i(:,:) / ( hv_a(:,:) + 1._wp - vmask_i(:,:) )
+      !
+      CALL wrk_dealloc( jpi,jpj,       zht, z_scale, zwu, zwv, zhdiv )
+      CALL wrk_dealloc( jpi,jpj,jpk,   ze3t )
+      !
       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
-
+      !
    END SUBROUTINE dom_vvl_sf_nxt
 
@@ -561,11 +561,11 @@
       !!               - recompute depths and water height fields
       !!
-      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
+      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
       !!               - Recompute:
-      !!                    fse3(u/v)_b       
-      !!                    fse3w_n           
-      !!                    fse3(u/v)w_b      
-      !!                    fse3(u/v)w_n      
-      !!                    fsdept_n, fsdepw_n  and fsde3w_n
+      !!                    e3(u/v)_b       
+      !!                    e3w_n           
+      !!                    e3(u/v)w_b      
+      !!                    e3(u/v)w_n      
+      !!                    gdept_n, gdepw_n  and gde3w_n
       !!                    h(u/v) and h(u/v)r
       !!
@@ -573,11 +573,12 @@
       !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in )               :: kt       ! time step
-      !! * Local declarations
-      INTEGER                             :: ji,jj,jk       ! dummy loop indices
-      REAL(wp)                            :: zcoef
-      !!----------------------------------------------------------------------
-
+      INTEGER, INTENT( in ) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zcoef        ! local scalar
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
       !
@@ -590,5 +591,5 @@
       ! Time filter and swap of scale factors
       ! =====================================
-      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
+      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          IF( neuler == 0 .AND. kt == nit000 ) THEN
@@ -600,10 +601,10 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
-      fsdept_b(:,:,:) = fsdept_n(:,:,:)
-      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
-
-      fse3t_n(:,:,:) = fse3t_a(:,:,:)
-      fse3u_n(:,:,:) = fse3u_a(:,:,:)
-      fse3v_n(:,:,:) = fse3v_a(:,:,:)
+      gdept_b(:,:,:) = gdept_n(:,:,:)
+      gdepw_b(:,:,:) = gdepw_n(:,:,:)
+
+      e3t_n(:,:,:) = e3t_a(:,:,:)
+      e3u_n(:,:,:) = e3u_a(:,:,:)
+      e3v_n(:,:,:) = e3v_a(:,:,:)
 
       ! Compute all missing vertical scale factor and depths
@@ -611,22 +612,21 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
+      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
+      
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      
       ! Vertical scale factor interpolations
-      ! ------------------------------------
-      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
-      ! t- and w- points depth
-      ! ----------------------
-      ! set the isf depth as it is in the initial step
-      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
-      fsdepw_n(:,:,1) = 0.0_wp
-      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-
+      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
+      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+
+      ! t- and w- points depth (set the isf depth as it is in the initial step)
+      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
+      gdepw_n(:,:,1) = 0.0_wp
+      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
       DO jk = 2, jpk
          DO jj = 1,jpj
@@ -635,44 +635,38 @@
                                                                  ! 1 for jk = mikt
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
-               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
-               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
+               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
+               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
+                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
+               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
             END DO
          END DO
       END DO
 
-      ! Local depth and Inverse of the local depth of the water column at u- and v- points
-      ! ----------------------------------------------------------------------------------
-      hu (:,:) = hu_a (:,:)
-      hv (:,:) = hv_a (:,:)
-
-      ! Inverse of the local depth
-      hur(:,:) = hur_a(:,:)
-      hvr(:,:) = hvr_a(:,:)
-
-      ! Local depth of the water column at t- points
-      ! --------------------------------------------
-      ht(:,:) = 0.
-      DO jk = 1, jpkm1
-         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+      ! Local depth and Inverse of the local depth of the water
+      ! -------------------------------------------------------
+      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
+      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
+      !
+      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
+      DO jk = 2, jpkm1
+         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
 
       ! Write outputs
       ! =============
-      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) )
-      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) )
-      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) )
-      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) )
-      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
+      CALL iom_put(     "e3t",   e3t_n(:,:,:) )
+      CALL iom_put(     "e3u",   e3u_n(:,:,:) )
+      CALL iom_put(     "e3v",   e3v_n(:,:,:) )
+      CALL iom_put(     "e3w",   e3w_n(:,:,:) )
+      CALL iom_put( "tpt_dep", gde3w_n(:,:,:) )
       IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
+         CALL iom_put( "e3tdef", ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100. * tmask(:,:,:) ) ** 2 )
 
       ! write restart file
       ! ==================
-      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
-
+      IF( lrst_oce )   CALL dom_vvl_rst( kt, 'WRITE' )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_sf_swp')
+      !
    END SUBROUTINE dom_vvl_sf_swp
 
@@ -696,5 +690,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
+      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_interpol')
       !
       SELECT CASE ( pout )    !==  type of interpolation  ==!
@@ -770,5 +764,5 @@
       END SELECT
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
+      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_interpol')
       !
    END SUBROUTINE dom_vvl_interpol
@@ -801,6 +795,6 @@
             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
             !
-            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
+            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
@@ -810,39 +804,39 @@
             !                             ! --------- !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
                ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
                WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  fse3t_n(:,:,:) = e3t_0(:,:,:)
-                  fse3t_b(:,:,:) = e3t_0(:,:,:)
+                  e3t_n(:,:,:) = e3t_0(:,:,:)
+                  e3t_b(:,:,:) = e3t_0(:,:,:)
                END WHERE
                IF( neuler == 0 ) THEN
-                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
+                  e3t_b(:,:,:) = e3t_n(:,:,:)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
-               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
-               fse3t_n(:,:,:) = fse3t_b(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
+               e3t_n(:,:,:) = e3t_b(:,:,:)
                neuler = 0
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
-               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
+               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
-               fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
+               e3t_b(:,:,:) = e3t_n(:,:,:)
                neuler = 0
             ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               DO jk=1,jpk
-                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
-                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
-                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
+               DO jk = 1, jpk
+                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
+                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
                END DO
-               fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               e3t_b(:,:,:) = e3t_n(:,:,:)
                neuler = 0
             ENDIF
@@ -875,6 +869,6 @@
             !
          ELSE                                   !* Initialize at "rest"
-            fse3t_b(:,:,:) = e3t_0(:,:,:)
-            fse3t_n(:,:,:) = e3t_0(:,:,:)
+            e3t_b(:,:,:) = e3t_0(:,:,:)
+            e3t_n(:,:,:) = e3t_0(:,:,:)
             sshn(:,:) = 0.0_wp
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
@@ -891,6 +885,6 @@
          !                                           ! all cases !
          !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
-         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) )
          !                                           ! ----------------------- !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
@@ -975,5 +969,5 @@
       !
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
-      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
+      IF( .NOT. ln_vvl_zstar .AND. nn_isf /= 0)   CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
       IF(lwp) THEN                   ! Print the choice
@@ -989,5 +983,5 @@
       !
 #if defined key_agrif
-      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
+      IF(.NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' )
 #endif
       !
@@ -996,5 +990,2 @@
    !!======================================================================
 END MODULE domvvl
-
-
-
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6060)
@@ -2,5 +2,5 @@
    !!==============================================================================
    !!                       ***  MODULE domzgr   ***
-   !! Ocean initialization : domain initialization
+   !! Ocean domain : definition of the vertical coordinate system
    !!==============================================================================
    !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
@@ -38,4 +38,5 @@
    USE closea            ! closed seas
    USE c1d               ! 1D vertical configuration
+   !
    USE in_out_manager    ! I/O manager
    USE iom               ! I/O library
@@ -73,5 +74,4 @@
 
   !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -102,5 +102,5 @@
       INTEGER ::   ios
       !
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
       !!----------------------------------------------------------------------
       !
@@ -120,10 +120,13 @@
          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 cavities             ln_isfcav = ', ln_isfcav
+         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 cavities             ln_isfcav = ', ln_isfcav
+         WRITE(numout,*) '      linear free surface            ln_linssh = ', ln_linssh
       ENDIF
+
+      IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time'
 
       ioptio = 0                       ! Check Vertical coordinate options
@@ -155,18 +158,18 @@
       !
       IF( nprint == 1 .AND. lwp )   THEN
-         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+         WRITE(numout,*) ' MIN val mbathy  ', MINVAL(  mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
          WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
-            &                   ' w ',   MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ),  &
-            &                   ' u ',   MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ),  &
-            &                   ' uw',   MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)),   &
-            &                   ' w ',   MINVAL( e3w_0(:,:,:) )
+            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  &
+            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  &
+            &                          ' uw', MINVAL(  e3uw_0(:,:,:) ), ' vw', MINVAL(  e3vw_0(:,:,:)),   &
+            &                          ' w ', MINVAL(   e3w_0(:,:,:) )
 
          WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
-            &                   ' w ',   MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ),  &
-            &                   ' u ',   MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ),  &
-            &                   ' uw',   MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),   &
-            &                   ' w ',   MAXVAL( e3w_0(:,:,:) )
+            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  &
+            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  &
+            &                          ' uw', MAXVAL(  e3uw_0(:,:,:) ), ' vw', MAXVAL(  e3vw_0(:,:,:) ),  &
+            &                          ' w ', MAXVAL(   e3w_0(:,:,:) )
       ENDIF
       !
@@ -674,9 +677,7 @@
       !!              - update bathy : meter bathymetry (in meters)
       !!----------------------------------------------------------------------
-      !!
       INTEGER ::   ji, jj, jl                    ! dummy loop indices
       INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
       REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy
-
       !!----------------------------------------------------------------------
       !
@@ -775,12 +776,4 @@
          IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
       ENDIF
-
-      IF( lwp .AND. nprint == 1 ) THEN      ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
-         WRITE(numout,*) ' ------------------'
-         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
-         WRITE(numout,*)
-      ENDIF
       !
       CALL wrk_dealloc( jpi, jpj, zbathy )
@@ -803,5 +796,4 @@
       !!                                     (min value = 1 over land)
       !!----------------------------------------------------------------------
-      !!
       INTEGER ::   ji, jj   ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk
@@ -835,5 +827,6 @@
    END SUBROUTINE zgr_bot_level
 
-      SUBROUTINE zgr_top_level
+
+   SUBROUTINE zgr_top_level
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE zgr_bot_level  ***
@@ -847,5 +840,4 @@
       !!                                     (min value = 1)
       !!----------------------------------------------------------------------
-      !!
       INTEGER ::   ji, jj   ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:) ::  zmik
@@ -881,9 +873,10 @@
    END SUBROUTINE zgr_top_level
 
+
    SUBROUTINE zgr_zco
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zgr_zco  ***
       !!
-      !! ** Purpose :   define the z-coordinate system
+      !! ** Purpose :   define the reference z-coordinate system
       !!
       !! ** Method  :   set 3D coord. arrays to reference 1D array 
@@ -895,14 +888,14 @@
       !
       DO jk = 1, jpk
-         gdept_0 (:,:,jk) = gdept_1d(jk)
-         gdepw_0 (:,:,jk) = gdepw_1d(jk)
-         gdep3w_0(:,:,jk) = gdepw_1d(jk)
-         e3t_0   (:,:,jk) = e3t_1d  (jk)
-         e3u_0   (:,:,jk) = e3t_1d  (jk)
-         e3v_0   (:,:,jk) = e3t_1d  (jk)
-         e3f_0   (:,:,jk) = e3t_1d  (jk)
-         e3w_0   (:,:,jk) = e3w_1d  (jk)
-         e3uw_0  (:,:,jk) = e3w_1d  (jk)
-         e3vw_0  (:,:,jk) = e3w_1d  (jk)
+         gdept_0(:,:,jk) = gdept_1d(jk)
+         gdepw_0(:,:,jk) = gdepw_1d(jk)
+         gde3w_0(:,:,jk) = gdepw_1d(jk)
+         e3t_0  (:,:,jk) = e3t_1d  (jk)
+         e3u_0  (:,:,jk) = e3t_1d  (jk)
+         e3v_0  (:,:,jk) = e3t_1d  (jk)
+         e3f_0  (:,:,jk) = e3t_1d  (jk)
+         e3w_0  (:,:,jk) = e3w_1d  (jk)
+         e3uw_0 (:,:,jk) = e3w_1d  (jk)
+         e3vw_0 (:,:,jk) = e3w_1d  (jk)
       END DO
       !
@@ -917,5 +910,5 @@
       !!                     
       !! ** Purpose :   the depth and vertical scale factor in partial step
-      !!      z-coordinate case
+      !!              reference z-coordinate case
       !!
       !! ** Method  :   Partial steps : computes the 3D vertical scale factors
@@ -957,8 +950,6 @@
       !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
       !!----------------------------------------------------------------------
-      !!
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       INTEGER  ::   ik, it, ikb, ikt ! temporary integers
-      LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
@@ -971,5 +962,5 @@
       IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zprt )
+      CALL wrk_alloc( jpi,jpj,jpk,   zprt )
       !
       IF(lwp) WRITE(numout,*)
@@ -977,13 +968,4 @@
       IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
-
-      ll_print = .FALSE.                   ! Local variable for debugging
-      
-      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
-         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
-      ENDIF
-
 
       ! bathymetry in level (from bathy_meter)
@@ -1196,58 +1178,32 @@
       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
      
-      ! Compute gdep3w_0 (vertical sum of e3w)
+      ! Compute gde3w_0 (vertical sum of e3w)
       IF ( ln_isfcav ) THEN ! if cavity
-         WHERE (misfdep == 0) misfdep = 1
+         WHERE( misfdep == 0 )   misfdep = 1
          DO jj = 1,jpj
             DO ji = 1,jpi
-               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
+               gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
                DO jk = 2, misfdep(ji,jj)
-                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
                END DO
-               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
+               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
                DO jk = misfdep(ji,jj) + 1, jpk
-                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
                END DO
             END DO
          END DO
       ELSE ! no cavity
-         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
+         gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
          DO jk = 2, jpk
-            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
+            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
          END DO
       END IF
-      !                                               ! ================= !
-      IF(lwp .AND. ll_print) THEN                     !   Control print   !
-         !                                            ! ================= !
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX( mbathy(ji,jj), 1 )
-               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
-               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
-               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
-               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
-               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
-               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-      ENDIF  
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   zprt )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
       !
    END SUBROUTINE zgr_zps
+
 
    SUBROUTINE zgr_isf
@@ -1265,10 +1221,8 @@
       !!              - bathy and isfdraft are modified
       !!----------------------------------------------------------------------
-      !!   
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       INTEGER  ::   ik, it           ! temporary integers
       INTEGER  ::   id, jd, nprocd
       INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
-      LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
@@ -1281,8 +1235,8 @@
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
-      !
-      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
-      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
+      IF( nn_timing == 1 )   CALL timing_start('zgr_isf')
+      !
+      CALL wrk_alloc( jpi,jpj,   zbathy, zmask, zrisfdep)
+      CALL wrk_alloc( jpi,jpj,   zmisfdep, zmbathy )
 
 
@@ -1300,6 +1254,6 @@
          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
       END DO 
-      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
-         risfdep(:,:) = 0. ; misfdep(:,:) = 1
+      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) > 0._wp)
+         risfdep(:,:) = 0.   ;   misfdep(:,:) = 1
       END WHERE
  
@@ -1308,5 +1262,5 @@
 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
       DO jl = 1, 10     
-         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
+         WHERE (bathy(:,:) == risfdep(:,:) )
             misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
             mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
@@ -1315,5 +1269,5 @@
             misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
             mbathy (:,:) = 0; bathy  (:,:) = 0._wp
-         ENDWHERE
+         END WHERE
          IF( lk_mpp ) THEN
             zbathy(:,:) = FLOAT( misfdep(:,:) )
@@ -1360,5 +1314,5 @@
                ! find the minimum change option:
                ! test bathy
-               IF (risfdep(ji,jj) .GT. 1) THEN
+               IF (risfdep(ji,jj) > 1) THEN
                zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
                  &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
@@ -1366,5 +1320,5 @@
                  &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
  
-                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
+                  IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN
                      IF (zbathydiff .LE. zrisfdepdiff) THEN
                         bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
@@ -1752,8 +1706,9 @@
       CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
       CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
-
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
-      
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zgr_isf')
+      !      
    END SUBROUTINE
+
 
    SUBROUTINE zgr_sco
@@ -1801,5 +1756,4 @@
       !!
       !!----------------------------------------------------------------------
-      !
       INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
       INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
@@ -1810,12 +1764,12 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
-
+      !!
       NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
-                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
+         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
      !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
       !
-      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
+      CALL wrk_alloc( jpi,jpj,   zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
       !
       REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
@@ -1876,13 +1830,19 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-           IF( bathy(ji,jj) == 0._wp ) THEN
-             iip1 = MIN( ji+1, jpi )
-             ijp1 = MIN( jj+1, jpj )
-             iim1 = MAX( ji-1, 1 )
-             ijm1 = MAX( jj-1, 1 )
-             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
-        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
-               zenv(ji,jj) = rn_sbot_min
-             ENDIF
+            IF( bathy(ji,jj) == 0._wp ) THEN
+               iip1 = MIN( ji+1, jpi )
+               ijp1 = MIN( jj+1, jpj )
+               iim1 = MAX( ji-1, 1 )
+               ijm1 = MAX( jj-1, 1 )
+!!gm BUG fix see ticket #1617
+               IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)            &
+                  &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )            &
+                  &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp )   zenv(ji,jj) = rn_sbot_min
+!!gm
+!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
+!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
+!!gm               zenv(ji,jj) = rn_sbot_min
+!!gm             ENDIF
+!!gm end
            ENDIF
          END DO
@@ -1975,15 +1935,4 @@
       ENDIF
       !
-      IF(lwp) THEN                             ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
-         WRITE(numout,*)
-         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
-         IF( nprint == 1 )   THEN        
-            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
-            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
-         ENDIF
-      ENDIF
-
       !                                        ! ==============================
       !                                        !   hbatu, hbatv, hbatf fields
@@ -2080,49 +2029,37 @@
       CALL lbc_lnk( e3uw_0, 'U', 1._wp )
       CALL lbc_lnk( e3vw_0, 'V', 1._wp )
-
-      fsdepw(:,:,:) = gdepw_0 (:,:,:)
-      fsde3w(:,:,:) = gdep3w_0(:,:,:)
-      !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      !
+      WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
+      WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
+      WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
+      WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
+      WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
+      WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
+      WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
 
 #if defined key_agrif
-      ! Ensure meaningful vertical scale factors in ghost lines/columns
-      IF( .NOT. Agrif_Root() ) THEN
-         !  
-         IF((nbondi == -1).OR.(nbondi == 2)) THEN
-            e3u_0(1,:,:) = e3u_0(2,:,:)
-         ENDIF
-         !
-         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
-            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
-         ENDIF
-         !
-         IF((nbondj == -1).OR.(nbondj == 2)) THEN
-            e3v_0(:,1,:) = e3v_0(:,2,:)
-         ENDIF
-         !
-         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
-            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
-         ENDIF
-         !
-      ENDIF
+      IF( .NOT. Agrif_Root() ) THEN    ! Ensure meaningful vertical scale factors in ghost lines/columns
+         IF( nbondi == -1 .OR. nbondi == 2 )   e3u_0(  1   ,  :   ,:) = e3u_0(  2   ,  :   ,:)
+         IF( nbondi ==  1 .OR. nbondi == 2 )   e3u_0(nlci-1,  :   ,:) = e3u_0(nlci-2,  :   ,:)
+         IF( nbondj == -1 .OR. nbondj == 2 )   e3v_0(  :   ,  1   ,:) = e3v_0(  :   ,  2   ,:)
+         IF( nbondj ==  1 .OR. nbondj == 2 )   e3v_0(  :   ,nlcj-1,:) = e3v_0(  :   ,nlcj-2,:)
+       ENDIF
 #endif
 
-      fsdept(:,:,:) = gdept_0 (:,:,:)
-      fsdepw(:,:,:) = gdepw_0 (:,:,:)
-      fsde3w(:,:,:) = gdep3w_0(:,:,:)
-      fse3t (:,:,:) = e3t_0   (:,:,:)
-      fse3u (:,:,:) = e3u_0   (:,:,:)
-      fse3v (:,:,:) = e3v_0   (:,:,:)
-      fse3f (:,:,:) = e3f_0   (:,:,:)
-      fse3w (:,:,:) = e3w_0   (:,:,:)
-      fse3uw(:,:,:) = e3uw_0  (:,:,:)
-      fse3vw(:,:,:) = e3vw_0  (:,:,:)
+!!gm   I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays)
+!!gm   and only that !!!!!
+!!gm   THIS should be removed from here !
+      gdept_n(:,:,:) = gdept_0(:,:,:)
+      gdepw_n(:,:,:) = gdepw_0(:,:,:)
+      gde3w_n(:,:,:) = gde3w_0(:,:,:)
+      e3t_n  (:,:,:) = e3t_0  (:,:,:)
+      e3u_n  (:,:,:) = e3u_0  (:,:,:)
+      e3v_n  (:,:,:) = e3v_0  (:,:,:)
+      e3f_n  (:,:,:) = e3f_0  (:,:,:)
+      e3w_n  (:,:,:) = e3w_0  (:,:,:)
+      e3uw_n (:,:,:) = e3uw_0 (:,:,:)
+      e3vw_n (:,:,:) = e3vw_0 (:,:,:)
+!!gm and obviously in the following, use the _0 arrays until the end of this subroutine
+!! gm end
 !!
       ! HYBRID : 
@@ -2130,7 +2067,7 @@
          DO ji = 1, jpi
             DO jk = 1, jpkm1
-               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
-            END DO
-            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
+               IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
+            END DO
+            IF( scobot(ji,jj) == 0._wp                )   mbathy(ji,jj) = 0
          END DO
       END DO
@@ -2141,15 +2078,15 @@
          WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
          WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
-            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
-            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
-            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
+            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
+            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
+            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
             &                          ' w ', MINVAL( e3w_0  (:,:,:) )
 
          WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
-            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
-            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
-            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
+            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
+            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
+            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
             &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
       ENDIF
@@ -2183,50 +2120,48 @@
          END DO
       ENDIF
-
-!================================================================================
-! check the coordinate makes sense
-!================================================================================
+      !
+      !================================================================================
+      ! check the coordinate makes sense
+      !================================================================================
       DO ji = 1, jpi
          DO jj = 1, jpj
-
+            !
             IF( hbatt(ji,jj) > 0._wp) THEN
                DO jk = 1, mbathy(ji,jj)
                  ! check coordinate is monotonically increasing
-                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
+                 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
-                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
+                    WRITE(numout,*) 'e3w',e3w_n(ji,jj,:)
+                    WRITE(numout,*) 'e3t',e3t_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                  ! and check it has never gone negative
-                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
+                 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
-                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
+                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                  ! and check it never exceeds the total depth
-                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
+                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                END DO
-
+               !
                DO jk = 1, mbathy(ji,jj)-1
                  ! and check it never exceeds the total depth
-                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                END DO
-
             ENDIF
-
          END DO
       END DO
@@ -2238,7 +2173,6 @@
    END SUBROUTINE zgr_sco
 
-!!======================================================================
+
    SUBROUTINE s_sh94()
-
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE s_sh94  ***
@@ -2251,5 +2185,4 @@
       !! Reference : Song and Haidvogel 1994. 
       !!----------------------------------------------------------------------
-      !
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
@@ -2257,7 +2190,8 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
-
-      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
-      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_alloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
 
       z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
@@ -2265,8 +2199,8 @@
       z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
       z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
-
+      !
       DO ji = 1, jpi
          DO jj = 1, jpj
-
+            !
             IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
                DO jk = 1, jpk
@@ -2297,7 +2231,7 @@
                zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
                zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
-               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
-               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
+               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
             END DO
            !
@@ -2331,12 +2265,12 @@
         END DO
       END DO
-
-      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
-      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
-
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+      !
    END SUBROUTINE s_sh94
 
+
    SUBROUTINE s_sf12
-
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE s_sf12 *** 
@@ -2354,5 +2288,4 @@
       !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
       !!----------------------------------------------------------------------
-      !
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zsmth               ! smoothing around critical depth
@@ -2361,5 +2294,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
-
+      !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
@@ -2425,7 +2358,7 @@
 
           DO jk = 1, jpk
-             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
-             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
-             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
+             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
+             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
+             gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
           END DO
 
@@ -2454,5 +2387,5 @@
              e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
              !
-             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
+             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
              e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
              e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
@@ -2467,13 +2400,11 @@
       CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
       !
-      !                                               ! =============
-
-      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
-      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
-
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+      !
    END SUBROUTINE s_sf12
 
+
    SUBROUTINE s_tanh()
-
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE s_tanh*** 
@@ -2485,13 +2416,12 @@
       !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
       !!----------------------------------------------------------------------
-
-      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      INTEGER  ::   ji, jj, jk       ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
-
       REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
       REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
-
-      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
-      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
+      CALL wrk_alloc( jpk,   z_esigt, z_esigw )
 
       z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
@@ -2523,7 +2453,7 @@
          zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
          zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
-         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
-         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
+         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
+         gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
       END DO
 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
@@ -2542,9 +2472,10 @@
          END DO
       END DO
-
-      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
-      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
-
+      !
+      CALL wrk_dealloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
+      CALL wrk_dealloc( jpk,   z_esigt, z_esigw          )
+      !
    END SUBROUTINE s_tanh
+
 
    FUNCTION fssig( pk ) RESULT( pf )
@@ -2618,33 +2549,32 @@
       REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
       REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
-      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
-      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
-      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
-      REAL(wp)                ::   za1,za2,za3    ! local variables
-      REAL(wp)                ::   zn1,zn2        ! local variables
-      REAL(wp)                ::   za,zb,zx       ! local variables
-      integer                 ::   jk
-      !!----------------------------------------------------------------------
-      !
-
-      zn1  =  1./(jpk-1.)
-      zn2  =  1. -  zn1
-
+      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
+      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
+      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
+      !
+      INTEGER  ::   jk             ! dummy loop index
+      REAL(wp) ::   za1,za2,za3    ! local scalar
+      REAL(wp) ::   zn1,zn2        !   -      - 
+      REAL(wp) ::   za,zb,zx       !   -      - 
+      !!----------------------------------------------------------------------
+      !
+      zn1  =  1._wp / REAL( jpkm1, wp )
+      zn2  =  1._wp -  zn1
+      !
       za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
       za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
       za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
-     
+      !
       za = pzb - za3*(pzs-za1)-za2
       za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
       zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
       zx = 1.0_wp-za/2.0_wp-zb
- 
+      !
       DO jk = 1, jpk
-        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
-                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
-                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
+         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
+            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
+            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
         p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
-      ENDDO 
-
+      END DO
       !
    END FUNCTION fgamma
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 6059)
+++ 	(revision )
@@ -1,122 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                    ***  domzgr_substitute.h90   ***
-   !!----------------------------------------------------------------------
-   !! ** purpose :   substitute fsdep. and fse.., the vert. depth and scale
-   !!      factors depending on the vertical coord. used, using CPP macro.
-   !!----------------------------------------------------------------------
-   !! History :  1.0  !  2005-10  (A. Beckmann, G. Madec) generalisation to all coord.
-   !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate
-   !!----------------------------------------------------------------------
-
-#if defined key_vvl
-! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._n)
-
-#   define  fse3t_b(i,j,k)   e3t_b(i,j,k)
-#   define  fse3u_b(i,j,k)   e3u_b(i,j,k)
-#   define  fse3v_b(i,j,k)   e3v_b(i,j,k)
-#   define  fse3uw_b(i,j,k)  e3uw_b(i,j,k)
-#   define  fse3vw_b(i,j,k)  e3vw_b(i,j,k)
-
-#   define  fsdept_n(i,j,k)  gdept_n(i,j,k)
-#   define  fsdepw_n(i,j,k)  gdepw_n(i,j,k)
-#   define  fsde3w_n(i,j,k)  gdep3w_n(i,j,k)
-#   define  fse3t_n(i,j,k)   e3t_n(i,j,k)
-#   define  fse3u_n(i,j,k)   e3u_n(i,j,k)
-#   define  fse3v_n(i,j,k)   e3v_n(i,j,k)
-#   define  fse3f_n(i,j,k)   e3f_n(i,j,k)
-#   define  fse3w_n(i,j,k)   e3w_n(i,j,k)
-#   define  fse3uw_n(i,j,k)  e3uw_n(i,j,k)
-#   define  fse3vw_n(i,j,k)  e3vw_n(i,j,k)
-
-#   define  fsdept_b(i,j,k)  gdept_b(i,j,k)
-#   define  fsdepw_b(i,j,k)  gdepw_b(i,j,k)
-#   define  fse3w_b(i,j,k)   e3w_b(i,j,k)
-
-#   define  fse3t_a(i,j,k)   e3t_a(i,j,k)
-#   define  fse3u_a(i,j,k)   e3u_a(i,j,k)
-#   define  fse3v_a(i,j,k)   e3v_a(i,j,k)
-
-#   define  fse3t_m(i,j)     e3t_m(i,j)
-
-#   define  hu_a(i,j)        ehu_a(i,j)
-#   define  hur_a(i,j)       ehur_a(i,j)
-#   define  hv_a(i,j)        ehv_a(i,j)
-#   define  hvr_a(i,j)       ehvr_a(i,j)
-#   define  hu_b(i,j)        ehu_b(i,j)
-#   define  hur_b(i,j)       ehur_b(i,j)
-#   define  hv_b(i,j)        ehv_b(i,j)
-#   define  hvr_b(i,j)       ehvr_b(i,j)
-
-! This part should be removed one day ...
-! ... In that case all occurence of the above statement functions
-!     have to be replaced in the code by xxx_n
-#   define  fsdept(i,j,k)    gdept_n(i,j,k)
-#   define  fsdepw(i,j,k)    gdepw_n(i,j,k)
-#   define  fsde3w(i,j,k)    gdep3w_n(i,j,k)
-#   define  fse3t(i,j,k)     e3t_n(i,j,k) 
-#   define  fse3u(i,j,k)     e3u_n(i,j,k) 
-#   define  fse3v(i,j,k)     e3v_n(i,j,k) 
-#   define  fse3f(i,j,k)     e3f_n(i,j,k) 
-#   define  fse3w(i,j,k)     e3w_n(i,j,k) 
-#   define  fse3uw(i,j,k)    e3uw_n(i,j,k)
-#   define  fse3vw(i,j,k)    e3vw_n(i,j,k)
-
-#else
-! z- or s-coordinate (1D or 3D + no time dependency) use reference in all cases
-
-#   define  fse3t_b(i,j,k)   e3t_0(i,j,k)
-#   define  fse3u_b(i,j,k)   e3u_0(i,j,k)
-#   define  fse3v_b(i,j,k)   e3v_0(i,j,k)
-#   define  fse3uw_b(i,j,k)  e3uw_0(i,j,k)
-#   define  fse3vw_b(i,j,k)  e3vw_0(i,j,k)
-
-#   define  fsdept_n(i,j,k)  gdept_0(i,j,k)
-#   define  fsdepw_n(i,j,k)  gdepw_0(i,j,k)
-#   define  fsde3w_n(i,j,k)  gdep3w_0(i,j,k)
-#   define  fse3t_n(i,j,k)   e3t_0(i,j,k)
-#   define  fse3u_n(i,j,k)   e3u_0(i,j,k)
-#   define  fse3v_n(i,j,k)   e3v_0(i,j,k)
-#   define  fse3f_n(i,j,k)   e3f_0(i,j,k)
-#   define  fse3w_n(i,j,k)   e3w_0(i,j,k)
-#   define  fse3uw_n(i,j,k)  e3uw_0(i,j,k)
-#   define  fse3vw_n(i,j,k)  e3vw_0(i,j,k)
-
-#   define  fsdept_b(i,j,k)  gdept_0(i,j,k)
-#   define  fsdepw_b(i,j,k)  gdepw_0(i,j,k)
-#   define  fse3w_b(i,j,k)   e3w_0(i,j,k)
-
-#   define  fse3t_a(i,j,k)   e3t_0(i,j,k)
-#   define  fse3u_a(i,j,k)   e3u_0(i,j,k)
-#   define  fse3v_a(i,j,k)   e3v_0(i,j,k)
-
-#   define  fse3t_m(i,j)     e3t_0(i,j,1)
-
-#   define  hu_a(i,j)        hu(i,j)
-#   define  hur_a(i,j)       hur(i,j)
-#   define  hv_a(i,j)        hv(i,j)
-#   define  hvr_a(i,j)       hvr(i,j)
-#   define  hu_b(i,j)        hu(i,j)
-#   define  hur_b(i,j)       hur(i,j)
-#   define  hv_b(i,j)        hv(i,j)
-#   define  hvr_b(i,j)       hvr(i,j)
-
-! This part should be removed one day ...
-! ... In that case all occurence of the above statement functions
-!     have to be replaced in the code by xxx_n
-#   define  fsdept(i,j,k)    gdept_0(i,j,k)
-#   define  fsdepw(i,j,k)    gdepw_0(i,j,k)
-#   define  fsde3w(i,j,k)    gdep3w_0(i,j,k)
-#   define  fse3t(i,j,k)     e3t_0(i,j,k)
-#   define  fse3u(i,j,k)     e3u_0(i,j,k)
-#   define  fse3v(i,j,k)     e3v_0(i,j,k)
-#   define  fse3f(i,j,k)     e3f_0(i,j,k)
-#   define  fse3w(i,j,k)     e3w_0(i,j,k)
-#   define  fse3uw(i,j,k)    e3uw_0(i,j,k)
-#   define  fse3vw(i,j,k)    e3vw_0(i,j,k)
-
-#endif
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 6060)
@@ -35,6 +35,4 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -250,25 +248,4 @@
       ENDIF
       !
-      IF( lwp .AND. kt == nit000 ) THEN
-         WRITE(numout,*) ' temperature Levitus '
-         WRITE(numout,*)
-         WRITE(numout,*)'  level = 1'
-         CALL prihre( ptsd(:,:,1    ,jp_tem), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         WRITE(numout,*)'  level = ', jpk/2
-         CALL prihre( ptsd(:,:,jpk/2,jp_tem), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         WRITE(numout,*)'  level = ', jpkm1
-         CALL prihre( ptsd(:,:,jpkm1,jp_tem), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         WRITE(numout,*)
-         WRITE(numout,*) ' salinity Levitus '
-         WRITE(numout,*)
-         WRITE(numout,*)'  level = 1'
-         CALL prihre( ptsd(:,:,1    ,jp_sal), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         WRITE(numout,*)'  level = ', jpk/2
-         CALL prihre( ptsd(:,:,jpk/2,jp_sal), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         WRITE(numout,*)'  level = ', jpkm1
-         CALL prihre( ptsd(:,:,jpkm1,jp_sal), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         WRITE(numout,*)
-      ENDIF
-      !
       IF( .NOT.ln_tsd_tradmp ) THEN                   !==   deallocate T & S structure   ==! 
          !                                              (data used only for initialisation)
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6060)
@@ -49,5 +49,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -120,41 +119,40 @@
             ENDIF
          ENDIF
-         !   
-         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
-         IF( lk_vvl ) THEN
+         !
+!!gm This is to be changed !!!!
+         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
+         IF( .NOT.ln_linssh ) THEN
             DO jk = 1, jpk
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+               e3t_b(:,:,jk) = e3t_n(:,:,jk)
             END DO
          ENDIF
+!!gm 
          ! 
       ENDIF
-      !
       ! 
       ! Initialize "now" and "before" barotropic velocities:
-      ! Do it whatever the free surface method, these arrays
-      ! being eventually used
-      !
+      ! Do it whatever the free surface method, these arrays being eventually used
       !
       un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
       ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
       !
+!!gm  the use of umsak & vmask is not necessary belox as un, vn, ub, vb are always masked
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
                !
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
+               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      un_b(:,:) = un_b(:,:) * hur  (:,:)
-      vn_b(:,:) = vn_b(:,:) * hvr  (:,:)
-      !
-      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
-      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
-      !
+      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
+      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
+      !
+      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
       !
       IF( nn_timing == 1 )   CALL timing_stop('istate_init')
@@ -184,6 +182,6 @@
       !
       DO jk = 1, jpk
-         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   &
-            &                + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
+         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) )   &
+            &                + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
          tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
       END DO
@@ -238,11 +236,7 @@
             !
             DO jk = 1, jpk
-               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
+               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
                tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
             END DO
-            !
-            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
-               &                             1     , jpi   , 5     , 1     , jpk   ,   &
-               &                             1     , 1.    , numout                  )
             !
             ! set salinity field to a constant value
@@ -314,8 +308,4 @@
             tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem)                            ! set nox temperature to tb
             !
-            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
-               &                            1     , jpi   , 5     , 1     , jpk   ,   &
-               &                            1     , 1.    , numout                  )
-            !
             ! set salinity field to a constant value
             ! --------------------------------------
@@ -363,20 +353,20 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 )         )   &
-                       &           * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2               &
-                       &       + (      15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) )       &
-                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &    
-                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   & 
-                       &           * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
+                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 )         )   &
+                       &           * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2               &
+                       &       + (      15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) )       &
+                       &                - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.)               &    
+                       &                + 7.  * (1500. - gdept_n(ji,jj,jk)) / 1500.             )   & 
+                       &           * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2
                   tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
                   tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
 
-                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 )  )  &
-                     &              * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2          &
-                     &          + (  35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000.         &
-                     &                - 1.62 * TANH( (fsdept(ji,jj,jk) - 60.  ) / 650. )       &
-                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 35.  ) / 100. )       &
-                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.)    )  &
-                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 
+                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 )  )  &
+                     &              * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2          &
+                     &          + (  35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000.         &
+                     &                - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60.  ) / 650. )       &
+                     &                + 0.2  * TANH( (gdept_n(ji,jj,jk) - 35.  ) / 100. )       &
+                     &                + 0.2  * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.)    )  &
+                     &              * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 
                   tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
                   tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
@@ -451,9 +441,9 @@
       zalfg = 0.5 * grav * rau0
       
-      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
+      zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
 
       DO jk = 2, jpkm1                                              ! Vertical integration from the surface
          zprn(:,:,jk) = zprn(:,:,jk-1)   &
-            &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
+            &         + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
       END DO  
 
