Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_cfg
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_cfg	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_cfg	(revision 15278)
@@ -15,5 +15,5 @@
 !-----------------------------------------------------------------------
    !
-   ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens. 
+   ln_e3_dep   = .false.   ! =T : e3=dk[depth] in discret sens. 
    !                       !      ===>>> will become the only possibility in v4.0
    !                       ! =F : e3 analytical derivative of depth function
@@ -34,31 +34,6 @@
 &namzgr        !   vertical coordinate
 !-----------------------------------------------------------------------
-   ln_mes      = .true.    !  Multi-Envelope s-coordinate
+   ln_zps      = .true.    !  z-coordinate with partial steps
    ln_linssh   = .true.    !  linear free surface
-/
-!-----------------------------------------------------------------------
-&namzgr_mes    !   MEs-coordinate                     
-!-----------------------------------------------------------------------
-   ln_envl     =   .TRUE. , .TRUE. , .TRUE. , .TRUE., .FALSE.  ! (T/F) If the envelope is used
-   nn_strt     =     1    ,   1    ,    1   ,   1    ,   1     ! Stretch. funct.: Madec 1996 (0) or
-                                                               ! Song & Haidvogel 1994 (1) or                    
-                                                               ! Siddorn & Furner 2012 (2)
-   nn_slev     =     3    ,   37   ,   17   ,   18   ,    0    ! number of s-lev between env(n-1)
-                                                               ! and env(n)
-   rn_e_hc     =     0.0  ,    0.0 ,   0.0  ,   0.0  ,   0.0   ! critical depth for transition to        
-                                                               ! stretch. coord.
-   rn_e_th     =     1.0  ,    1.2 ,   2.0  ,   0.0  ,   0.0   ! surf. control param.:                  
-                                                               ! SH94 or MD96: 0<=th<=20
-                                                               ! SF12: thickness surf. cell
-   rn_e_bb     =     0.9  ,    0.9 ,   0.85 ,   0.0  ,   0.0   ! bot. control param.:
-                                                               ! SH94 or MD96: 0<=bb<=1
-                                                               ! SF12: offset for calculating Zb
-   rn_e_al     =     0.0  ,    0.0 ,   0.0  ,   0.0  ,   0.0   ! alpha stretching param with SF12
-   rn_e_ba     =     0.0  ,    0.0 ,   0.0  ,   0.0  ,   0.0   ! SF12 bathymetry scaling factor for 
-                                                               ! calculating Zb
-   rn_bot_min  = 10.0       ! minimum depth of the ocean bottom (>0) (m)
-   rn_bot_max  = 5900.0     ! maximum depth of the ocean bottom (= ocean depth) (>0) (m)
-
-   ln_loc_mes  = .TRUE.
 /
 !-----------------------------------------------------------------------
@@ -67,22 +42,4 @@
    nn_bathy    =  1
    nn_msh      =  1
-   jphgr_msh   =  0
-   ldbletanh   = .TRUE.
-   ppglam0     =  999999.d0
-   ppgphi0     =  999999.d0
-   ppe1_deg    =  999999.d0
-   ppe2_deg    =  999999.d0
-   ppe1_m      =  999999.d0
-   ppe2_m      =  999999.d0                 
-   ppa0        =  103.9530096000000
-   ppa1        =  2.415951269000000
-   ppa2        =  100.7609285000000
-   ppacr       =  7.0
-   ppacr2      =  13.0
-   ppdzmin     =  999999.0
-   pphmax      =  999999.0
-   ppkth       =  15.35101370000000
-   ppkth2      =  48.02989372000000
-   ppsur       = -3958.951371276829
    rn_atfp     =  0.1
    rn_e3zps_min=  25.0
@@ -90,4 +47,22 @@
    rn_hmin     = -8.0
    rn_rdt      =  1350.0
+   jphgr_msh   =       0               !  type of horizontal mesh
+   ppglam0     =  999999.0             !  longitude of first raw and column T-point (jphgr_msh = 1)
+   ppgphi0     =  999999.0             ! latitude  of first raw and column T-point (jphgr_msh = 1)
+   ppe1_deg    =  999999.0             !  zonal      grid-spacing (degrees)
+   ppe2_deg    =  999999.0             !  meridional grid-spacing (degrees)
+   ppe1_m      =  999999.0             !  zonal      grid-spacing (degrees)
+   ppe2_m      =  999999.0             !  meridional grid-spacing (degrees)
+   ppsur       =   -3958.951371276829  !  ORCA r4, r2 and r05 coefficients
+   ppa0        =    103.9530096000000  ! (default coefficients)
+   ppa1        =      2.415951269000000  !
+   ppkth       =     15.35101370000000 !
+   ppacr       =       7.0             !
+   ppdzmin     =  999999.0             !  Minimum vertical spacing
+   pphmax      =  999999.0             !  Maximum depth
+   ldbletanh   =   .TRUE.              !  Use/do not use double tanf function for vertical coordinates
+   ppa2        =     100.7609285000000 !  Double tanh function parameters
+   ppkth2      =      48.02989372000000  !
+   ppacr2      =      13.              !
 /
 !-----------------------------------------------------------------------
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_ref
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_ref	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/namelist_ref	(revision 15278)
@@ -107,4 +107,12 @@
    ln_isfcav   = .false.   !  ice shelf cavity
    ln_linssh   = .false.   !  linear free surface
+   ln_loc_zgr  = .false.   !  to localise the chosen vert. coord. system
+                           !  If TRUE, the user must provide:
+                           !     1) a bathy_meter.nc file including
+                           !        *) s2z_msk: 2D mask for s(=2), s2z(=1) and z(=0) areas
+                           !        *) s2z_wgt: 2D field with 1 in s-areas, distance-based
+                           !                    weights (0<wgt<1) in s2z-areas and 0 elsewhere.
+                           !     2) a domain_cfg_global.nc file including all the geometrical
+                           !        information describing the vertical grid used globally.
 /
 !-----------------------------------------------------------------------
@@ -160,12 +168,9 @@
    rn_bot_max  = 2201.5    ! maximum depth of the ocean bottom (= ocean depth) (>0) (m)
 
-   ln_loc_mes  = .FALSE.   ! To use MEs-coordinates only in localised regions. 
-                           ! If TRUE, the user must provide:
-                           !     1) a bathy_meter.nc file including
-                           !        *) s2z_msk: 2D mask for s(=2), s2z(=1) and z(=0) areas
-                           !        *) s2z_wgt: 2D field with 1 in s-areas, distance-based 
-                           !                    weights (0<wgt<1) in s2z-areas and 0 elsewhere.
-                           !     2) a domain_cfg_global.nc file including all the geometrical
-                           !        information describing the vertical grid used globally.  
+   ln_pst_mes   = .FALSE.  ! To use partial-steps algorithm where MEs-coordinates are used.
+   ln_pst_l2g   = .FALSE.  ! To use partial-steps algorithm in the transition area 
+                           ! (requires ln_loc_zgr = .TRUE.)
+   rn_e3pst_min = 20.      ! partial step thickness is set larger than the minimum of
+   rn_e3pst_rat = 0.1      ! rn_e3pst_min and rn_e3pst_rat*e3t, with 0<rn_e3pst_rat<1
 /
 !-----------------------------------------------------------------------
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/dom_oce.f90
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/dom_oce.f90	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/dom_oce.f90	(revision 15278)
@@ -178,4 +178,6 @@
    LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF 
    LOGICAL, PUBLIC ::   ln_linssh    !: variable grid flag
+   LOGICAL, PUBLIC ::   ln_loc_zgr   !: To localise (.TRUE.) or not (.FALSE.) the 
+                                     !  chosen vertical coordinate system
 
    !                                                        !  ref.   ! before  !   now   ! after  !
@@ -245,8 +247,10 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mes_msk      ! mask for local-ME coordinates (ln_loc_mes=TRUE)  
-                                                                       !    MEs-area (=1),
-                                                                       !    z-area (=0)
-
+   ! LOCAL VERTICAL COORDINATE SYSTEM 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: l2g_msk ! mask for identifying areas using
+                                                                  ! -) the local coord. system (=2),
+                                                                  ! -) transitioning from local to global (=1),
+                                                                  ! -) using global coord. system (=0)
+                                                                  ! If ln_loc_zgr=.FALSE. then l2g_msk(:,:)=2
 
    !!----------------------------------------------------------------------
@@ -370,7 +374,8 @@
 
       ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       &
-         &     tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                       & 
-         &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
-         &     mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , mes_msk(jpi,jpj) , STAT=ierr(9) )
+         &      tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                      &
+         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj), &
+         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) ,                   &
+         &      l2g_msk(jpi,jpj) , STAT=ierr(9) )
 
 ! (ISF) Allocation of basic array   
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domwri.f90
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domwri.f90	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domwri.f90	(revision 15278)
@@ -16,5 +16,4 @@
    !!----------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
-   USE mes, ONLY : ln_loc_mes
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O library
@@ -596,5 +595,5 @@
          END DO
       END DO
-      IF( ln_loc_mes ) zx1 = zx1 * mes_msk 
+      IF( ln_loc_zgr .AND. ln_mes ) WHERE (l2g_msk(:,:) == 0.0)  zx1(:,:) = 0.0 
       CALL lbc_lnk( zx1, 'T', 1. )
       !
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domzgr.f90
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domzgr.f90	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domzgr.f90	(revision 15278)
@@ -113,5 +113,5 @@
       INTEGER ::   ios
       !
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_mes, ln_isfcav, ln_linssh
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_mes, ln_isfcav, ln_linssh, ln_loc_zgr
       !!----------------------------------------------------------------------
       !
@@ -141,4 +141,5 @@
 
       IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time'
+      IF( ln_loc_zgr.AND. lwp) WRITE(numout,*) '   the chosen vert. coord. system is implemented LOCALLY'
 
       ioptio = 0                       ! Check Vertical coordinate options
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90	(revision 15278)
@@ -44,5 +44,9 @@
    REAL(wp)           :: rn_e_al(max_nn_env)     ! Array specifing stretching parameter rn_alpha of SF12 for
                                                  ! each vertical sub-zone
-   LOGICAL, PUBLIC    :: ln_loc_mes              ! To use localised MEs (.TRUE.) or not (.FALSE.
+   LOGICAL, PUBLIC    :: ln_pst_mes              ! To apply partial steps to MEs (.TRUE.) or not (.FALSE.).
+   LOGICAL, PUBLIC    :: ln_pst_l2g              ! To apply partial steps to the transition zone (.TRUE.) 
+                                                 ! or not (.FALSE.) when ln_loc_zgr = .TRUE. 
+   REAL, PUBLIC       :: rn_e3pst_min            ! partial steps parameters (require ln_pst_mes = .TRUE.)
+   REAL, PUBLIC       :: rn_e3pst_rat            ! partial steps parameters (require ln_pst_mes = .TRUE.)
    !
    REAL(wp), POINTER, DIMENSION(:,:,:) :: envlt  ! array for the envelopes
@@ -669,8 +673,9 @@
       REAL(wp)                     :: rn_bot_max ! max depth of ocean bottom (= ocean depth) (>0) (m)
 
-      NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl, &
-                          nn_strt    , nn_slev    , rn_e_hc, &
-                          rn_e_th    , rn_e_bb    , rn_e_ba, &
-                          rn_e_al    , ln_loc_mes
+      NAMELIST/namzgr_mes/rn_bot_min  , rn_bot_max  , ln_envl   , &
+          &               nn_strt     , nn_slev     , rn_e_hc   , &
+          &               rn_e_th     , rn_e_bb     , rn_e_ba   , &
+          &               rn_e_al     , ln_pst_mes  , ln_pst_l2g, &
+          &               rn_e3pst_min, rn_e3pst_rat
       !!----------------------------------------------------------------------
       !
@@ -720,4 +725,6 @@
       WRITE(ctlmes,*) 'TOT number of levels (jpk) IS DIFFERENT from sum over nn_slev(:)'
       IF ( SUM(nn_slev(:)) /= jpk ) CALL ctl_stop( ctlmes )
+      !
+      IF( .NOT.(ln_loc_zgr) .AND. ln_pst_l2g) ln_pst_l2g = .FALSE.
       !
       ! 3) Reading Bathymetry and envelopes
@@ -810,4 +817,25 @@
             WRITE(numout,*) '          ------------------------------------------------------------------'
          ENDDO
+         IF( ln_loc_zgr) THEN
+           WRITE(numout,*) ''
+           WRITE(numout,*) '          LOCALISED MEs '
+           WRITE(numout,*) '          ============= '
+         ENDIF
+         IF( ln_pst_mes .OR. ln_pst_l2g) THEN
+           WRITE(numout,*) ''
+           WRITE(numout,*) '          APPLYING PARTIAL-STEPS to '
+           WRITE(numout,*) '          ========================= '
+           WRITE(numout,*) ''
+           WRITE(numout,*) '            a) MEs area       : ', ln_pst_mes
+           IF( ln_loc_zgr) WRITE(numout,*) '            b) Transition area: ', ln_pst_l2g
+           WRITE(numout,*) ''
+           WRITE(numout,*) '            Partial step thickness is set larger than the minimum of'
+           WRITE(numout,*) ''
+           WRITE(numout,*) '                  rn_e3pst_min = ', rn_e3pst_min
+           WRITE(numout,*) '            and'
+           WRITE(numout,*) '                  rn_e3pst_rat = ', rn_e3pst_rat
+           WRITE(numout,*) ''
+           WRITE(numout,*) '            with 0 < rn_e3pst_rat < 1'
+         ENDIF
       ENDIF
 
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrloc.F90
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrloc.F90	(revision 15278)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrloc.F90	(revision 15278)
@@ -0,0 +1,262 @@
+MODULE zgrloc
+   !!==============================================================================
+   !!                       ***  MODULE zgrmes   ***
+   !! Ocean initialization : Multiple Enveloped s coordinate (MES)
+   !!==============================================================================
+   !!  NEMO      4.0  ! 2021-07  (D. Bruciaferri)   
+   !!----------------------------------------------------------------------
+   !
+   USE oce               ! ocean variables
+   USE dom_oce           ! ocean domain
+   USE depth_e3          ! depth <=> e3
+   USE closea            ! closed seas
+   USE mes               ! MEs-coordinates
+   !
+   USE in_out_manager    ! I/O manager
+   USE iom               ! I/O library
+   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp           ! distributed memory computing library
+   USE wrk_nemo          ! Memory allocation
+   USE timing            ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zgr_loc        ! called by zgrmes.F90
+   !
+   ! INPUT GLOBAL VERTICAL COORDINATE SYSTEM
+   !
+   REAL(wp), POINTER, DIMENSION(:)     :: e3t_1d_in  , e3w_1d_in   !: ref. scale factors (m)
+   REAL(wp), POINTER, DIMENSION(:)     :: gdept_1d_in, gdepw_1d_in !: ref. depth (m)
+   !
+   REAL(wp), POINTER, DIMENSION(:,:)   :: mbathy_in
+   !
+   REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m]
+   REAL(wp), POINTER, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in         !:   -      -
+   ! 
+   REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_in, gdepw_in               !: depths [m]
+   !
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+
+CONTAINS
+
+! =====================================================================================================
+
+   SUBROUTINE zgr_loc
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE zgr_loc  ***
+      !!
+      !! ** Purpose : Wrap the steps to generate a localised vertical 
+      !!              coordinate systems
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk     ! dummy loop index
+      !-----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in)
+      CALL wrk_alloc( jpi, jpj, mbathy_in)
+      CALL wrk_alloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in)
+      CALL wrk_alloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in)
+      !
+      ! Reading the input vertical grid that will be used globally
+      CALL zgr_dom_read
+      ! Creating the local cooridnate system within the global input vertical grid
+      IF ( ln_mes ) CALL zgr_mes_local
+      !
+      CALL wrk_dealloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in)
+      CALL wrk_dealloc( jpi, jpj, mbathy_in)
+      CALL wrk_dealloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in)
+      CALL wrk_dealloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in)
+
+   END SUBROUTINE zgr_loc
+
+! =====================================================================================================
+
+   SUBROUTINE zgr_dom_read
+     !!---------------------------------------------------------------------
+     !!              ***  ROUTINE zgr_dom_read  ***
+     !!
+     !! ** Purpose :   Read the vertical information in the domain configuration file
+     !!
+     !!----------------------------------------------------------------------
+     INTEGER  ::   jk     ! dummy loop index
+     INTEGER  ::   inum   ! local logical unit
+     REAL(WP) ::   z_cav
+     REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
+     !!----------------------------------------------------------------------
+     !
+     IF(lwp) THEN
+       WRITE(numout,*)
+       WRITE(numout,*) '   zgr_dom_read : read the vertical coordinates in domain_cfg_global.nc file'
+       WRITE(numout,*) '   ~~~~~~~~'
+     ENDIF
+     !
+     CALL iom_open( 'domain_cfg_global.nc', inum )
+     !
+     !                          !* ocean cavities under iceshelves
+     CALL iom_get( inum, 'ln_isfcav', z_cav )
+     IF( z_cav == 0._wp ) THEN   ;   ln_isfcav = .false.   ;   ELSE   ;   ln_isfcav = .true.   ;   ENDIF
+     !
+     !                          !* vertical scale factors
+     CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d_in ) ! 1D reference coordinate
+     CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d_in )
+     !
+     CALL iom_get( inum, jpdom_global, 'bottom_level' , mbathy_in )     ! 2D mbathy
+     !
+     CALL iom_get( inum, jpdom_global, 'e3t_0'  , e3t_in  )     ! 3D coordinate
+     CALL iom_get( inum, jpdom_global, 'e3u_0'  , e3u_in  )
+     CALL iom_get( inum, jpdom_global, 'e3v_0'  , e3v_in  )
+     CALL iom_get( inum, jpdom_global, 'e3f_0'  , e3f_in  )
+     CALL iom_get( inum, jpdom_global, 'e3w_0'  , e3w_in  )
+     CALL iom_get( inum, jpdom_global, 'e3uw_0' , e3uw_in )
+     CALL iom_get( inum, jpdom_global, 'e3vw_0' , e3vw_in )
+     !
+     !                          !* depths
+     !                                   !- old depth definition (obsolescent feature)
+     IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
+        & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
+        & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
+        & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
+        CALL ctl_warn( 'zgr_dom_read : old definition of depths and scale factors used ', & 
+           &           '           depths at t- and w-points read in the domain configuration file')
+        CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d_in )   
+        CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d_in )
+        CALL iom_get( inum, jpdom_global , 'gdept_0' , gdept_in )
+        CALL iom_get( inum, jpdom_global , 'gdepw_0' , gdepw_in )
+        !
+     ELSE                                !- depths computed from e3. scale factors
+        CALL e3_to_depth( e3t_1d_in, e3w_1d_in, gdept_1d_in, gdepw_1d_in )    ! 1D reference depth
+        CALL e3_to_depth( e3t_in   , e3w_in   , gdept_in   , gdepw_in    )    ! 3D depths
+        IF(lwp) THEN
+          WRITE(numout,*)
+          WRITE(numout,*) '              GLOBAL reference 1D z-coordinate depth and scale factors:'
+          WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
+          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d_in(jk), gdepw_1d_in(jk), e3t_1d_in(jk), e3w_1d_in(jk), jk = 1, jpk )
+        ENDIF
+     ENDIF
+     !
+     CALL iom_close( inum )
+     !
+   END SUBROUTINE zgr_dom_read
+
+! =====================================================================================================
+
+   SUBROUTINE zgr_mes_local
+     !!---------------------------------------------------------------------
+     !!              ***  ROUTINE zgr_dom_merge  ***
+     !!
+     !! ** Purpose :   Create a local MEs grid within a gloabal grid 
+     !!                using different vertical coordinates.
+     !!
+     !!----------------------------------------------------------------------
+     INTEGER                      ::   ji, jj, jk ! dummy loop index
+     INTEGER                      ::   inum       ! local logical unit
+     REAL(wp), DIMENSION(jpi,jpj) ::   l2g_wgt    ! weigths for computing model
+                                                  ! levels in transition area
+     !!----------------------------------------------------------------------
+
+     IF(lwp) THEN
+       WRITE(numout,*)
+       WRITE(numout,*) '   zgr_mes_merge : read the vertical coordinates in domain_cfg_global.nc file'
+       WRITE(numout,*) '   ~~~~~~~~'
+     ENDIF
+     !
+     CALL iom_open( 'bathy_meter.nc', inum )
+     CALL iom_get( inum, jpdom_data, 's2z_msk', l2g_msk)
+     CALL iom_get( inum, jpdom_data, 's2z_wgt', l2g_wgt)
+
+     DO jj = 1,jpj
+        DO ji = 1,jpi
+           SELECT CASE (INT(l2g_msk(ji,jj)))
+             CASE (0) ! global zps area
+                  mbathy (ji,jj  ) = mbathy_in(ji,jj)
+                  gdept_0(ji,jj,:) = gdept_in(ji,jj,:)
+                  gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:)
+                  e3t_0  (ji,jj,:) = e3t_in (ji,jj,:)
+                  e3w_0  (ji,jj,:) = e3w_in (ji,jj,:)
+                  e3u_0  (ji,jj,:) = e3u_in (ji,jj,:)
+                  e3v_0  (ji,jj,:) = e3v_in (ji,jj,:)
+                  e3f_0  (ji,jj,:) = e3f_in (ji,jj,:)
+                  e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:)
+                  e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:)
+             CASE (1) ! MEs to zps transition area 
+                  gdept_0(ji,jj,:) =           l2g_wgt(ji,jj)   * gdept_0(ji,jj,:) + &
+                    &                ( 1._wp - l2g_wgt(ji,jj) ) * gdept_in(ji,jj,:)
+                  gdepw_0(ji,jj,:) =           l2g_wgt(ji,jj)   * gdepw_0(ji,jj,:) + &
+                    &                ( 1._wp - l2g_wgt(ji,jj) ) * gdepw_in(ji,jj,:)
+             CASE (2) ! MEs area
+                  CYCLE     
+           END SELECT
+        END DO
+     END DO
+     !
+     ! e3t, e3w for transition zone
+     ! as finite differences
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+              DO jk = 1,jpkm1
+                 e3t_0(ji,jj,jk)   = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk)
+                 e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk)
+              ENDDO
+              ! Surface
+              jk = 1
+              e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1))
+              !
+              ! Bottom
+              jk = jpk
+              e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk))
+           END IF
+        END DO
+     END DO
+     !
+     ! MBATHY transition zone
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+              DO jk = 1, jpkm1
+                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
+                 IF( scobot(ji,jj) == 0.e0              ) mbathy(ji,jj) = 0
+                 IF( scobot(ji,jj) < 0.e0               ) mbathy(ji,jj) = INT( scobot(ji,jj)) ! do we need it?
+              END DO
+           END IF
+        END DO
+     END DO
+     !
+     ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0 
+     ! for transition zone
+     !
+     DO jj = 1, jpjm1
+        DO ji = 1, jpim1
+           IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+              DO jk = 1, jpk
+                 e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
+                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) /         &
+                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
+
+                 e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
+                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) /         &
+                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
+
+                 e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
+                                   MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) /         &
+                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
+
+                 e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
+                                   MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) /         &
+                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
+
+                 e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)         +       &
+                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk)      +       &
+                                  MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) +       &
+                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) )    /       &
+                                  MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  &
+                                +   MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))  )
+              END DO
+           END IF
+        END DO
+     END DO
+
+   END SUBROUTINE zgr_mes_local
+
+END MODULE zgrloc
Index: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90
===================================================================
--- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90	(revision 15276)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90	(revision 15278)
@@ -12,4 +12,5 @@
    USE closea            ! closed seas
    USE mes               ! MEs-coordinates
+   USE zgrloc            ! localised vert. coor. system 
    !
    USE in_out_manager    ! I/O manager
@@ -24,15 +25,4 @@
 
    PUBLIC   zgr_mes        ! called by domzgr.F90
-   ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.)
-   !
-   REAL(wp), POINTER, DIMENSION(:)     :: e3t_1d_in  , e3w_1d_in   !: ref. scale factors (m)
-   REAL(wp), POINTER, DIMENSION(:)     :: gdept_1d_in, gdepw_1d_in !: ref. depth (m)
-   !
-   REAL(wp), POINTER, DIMENSION(:,:)   :: mbathy_in 
-   !
-   REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m]
-   REAL(wp), POINTER, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in         !:   -      -
-   ! 
-   REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_in, gdepw_in               !: depths [m]
 
    !! * Substitutions
@@ -50,30 +40,62 @@
       !!              MEs-coordinates systems
       !!----------------------------------------------------------------------
-      !
-      ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.)
-      !
+      INTEGER  ::   ji, jj, jk     ! dummy loop index
       !-----------------------------------------------------------------------
       !
+      ! Initialise the mask to use MEs-coord. globally
+      l2g_msk(:,:) = 2.0
+     
       ! Generating a global MEs vertical grid
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       CALL mes_build
-
+      
       ! Local MEs
-      IF( ln_loc_mes ) THEN
-         !
-         CALL wrk_alloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in)
-         CALL wrk_alloc( jpi, jpj, mbathy_in)
-         CALL wrk_alloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in)
-         CALL wrk_alloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in)
-         !
-         ! Reading the input vertical grid that will be used globally
-         CALL zgr_dom_read
-         ! Creating the local MEs within the global input vertical grid
-         CALL zgr_mes_local
-         !
-         CALL wrk_dealloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in  , e3w_1d_in)
-         CALL wrk_dealloc( jpi, jpj, mbathy_in)
-         CALL wrk_dealloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in)
-         CALL wrk_dealloc( jpi, jpj, jpk, e3u_in  , e3v_in  , e3f_in, e3uw_in, e3vw_in)
-         !
+      ! ~~~~~~~~~
+      IF( ln_loc_zgr )  CALL zgr_loc
+      
+      ! MEs with partial steps
+      ! ~~~~~~~~~~~~~~~~~~~~~~
+      IF(ln_pst_mes .OR. ln_pst_l2g)  CALL zgr_pst_mes
+      !
+      ! Computing gdep3w
+      gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
+      DO jk = 2, jpk
+         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
+      END DO
+      !
+      ! From here equal to sco code - domzgr.F90 line 2183
+      !
+      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
+      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
+      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
+      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
+      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
+      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
+      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
+      !
+      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
+      !
+      IF ( ln_loc_zgr ) THEN
+         IF ( .NOT. ln_pst_l2g ) THEN
+            ! Only in the transition zone since MEs
+            ! zone has been already taken care of
+            IF ( lwp ) WRITE(numout,*) 'Refine mbathy'
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+                     DO jk = 1, jpkm1
+                        IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
+                     END DO
+                     IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
+                  END IF
+               END DO
+            END DO
+         END IF
       END IF
 
@@ -82,233 +104,264 @@
 ! =====================================================================================================
 
-   SUBROUTINE zgr_dom_read
+   SUBROUTINE zgr_pst_mes
      !!---------------------------------------------------------------------
-     !!              ***  ROUTINE zgr_dom_read  ***
+     !!              ***  ROUTINE zgr_pst_mes  ***
      !!
-     !! ** Purpose :   Read the vertical information in the domain configuration file
+     !! ** Purpose : Apply partial steps to MEs coordinates 
+     !!              (and transition zone if using local MEs)
      !!
      !!----------------------------------------------------------------------
-     INTEGER  ::   jk     ! dummy loop index
-     INTEGER  ::   inum   ! local logical unit
-     REAL(WP) ::   z_cav
-     REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
+     INTEGER                             ::   ji, jj, jk, jp, jl   ! dummy loop index
+     INTEGER                             ::   jbku, jbkv, jbkf     ! dummy loop index
+     REAL(wp)                            ::   zdepth, e3psm, e3psr ! local variables
+     REAL(wp)                            ::   zgdept, zgdepw, zgdepwp1
+     REAL(wp)                            ::   ze3t, ze3w, e3u_j, e3u_jp1
+     REAL(wp)                            ::   e3t_i, e3t_ip1, e3w_i, e3w_ip1
+     REAL(wp)                            ::   e3t_j, e3t_jp1, e3w_j, e3w_jp1
+     REAL(wp)                            ::   ps_i, ps_ip1, ps_j, ps_jp1
+     REAL(wp)                            ::   ps_uj, ps_ujp1, ps_p, ps_s
+     !
+     REAL(wp), POINTER, DIMENSION(:,:)   :: pst_msk              !: mask for areas where to
+                                                                 !: apply partial steps
+     REAL(wp), POINTER, DIMENSION(:,:)   :: apst                 !: tracking points where
+                                                                 !: we apply partial steps
+     ! ARRAYS TO BACKUP ORIGINAL MEs COORDINATE SYSTEM 
+     REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_mes, e3w_mes     !: scale factors [m]
+     REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_mes, gdepw_mes !: depths [m]
      !!----------------------------------------------------------------------
-     !
-     IF(lwp) THEN
-       WRITE(numout,*)
-       WRITE(numout,*) '   zgr_dom_read : read the vertical coordinates in domain_cfg_global.nc file'
-       WRITE(numout,*) '   ~~~~~~~~'
-     ENDIF
-     !
-     CALL iom_open( 'domain_cfg_global.nc', inum )
-     !
-     !                          !* ocean cavities under iceshelves
-     CALL iom_get( inum, 'ln_isfcav', z_cav )
-     IF( z_cav == 0._wp ) THEN   ;   ln_isfcav = .false.   ;   ELSE   ;   ln_isfcav = .true.   ;   ENDIF
-     !
-     !                          !* vertical scale factors
-     CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d_in ) ! 1D reference coordinate
-     CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d_in )
-     !
-     CALL iom_get( inum, jpdom_global, 'bottom_level' , mbathy_in )     ! 2D mbathy
-     !
-     CALL iom_get( inum, jpdom_global, 'e3t_0'  , e3t_in  )     ! 3D coordinate
-     CALL iom_get( inum, jpdom_global, 'e3u_0'  , e3u_in  )
-     CALL iom_get( inum, jpdom_global, 'e3v_0'  , e3v_in  )
-     CALL iom_get( inum, jpdom_global, 'e3f_0'  , e3f_in  )
-     CALL iom_get( inum, jpdom_global, 'e3w_0'  , e3w_in  )
-     CALL iom_get( inum, jpdom_global, 'e3uw_0' , e3uw_in )
-     CALL iom_get( inum, jpdom_global, 'e3vw_0' , e3vw_in )
-     !
-     !                          !* depths
-     !                                   !- old depth definition (obsolescent feature)
-     IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
-        & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
-        & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
-        & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
-        CALL ctl_warn( 'zgr_dom_read : old definition of depths and scale factors used ', & 
-           &           '           depths at t- and w-points read in the domain configuration file')
-        CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d_in )   
-        CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d_in )
-        CALL iom_get( inum, jpdom_global , 'gdept_0' , gdept_in )
-        CALL iom_get( inum, jpdom_global , 'gdepw_0' , gdepw_in )
-        !
-     ELSE                                !- depths computed from e3. scale factors
-        CALL e3_to_depth( e3t_1d_in, e3w_1d_in, gdept_1d_in, gdepw_1d_in )    ! 1D reference depth
-        CALL e3_to_depth( e3t_in   , e3w_in   , gdept_in   , gdepw_in    )    ! 3D depths
-        IF(lwp) THEN
-          WRITE(numout,*)
-          WRITE(numout,*) '              GLOBAL reference 1D z-coordinate depth and scale factors:'
-          WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
-          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d_in(jk), gdepw_1d_in(jk), e3t_1d_in(jk), e3w_1d_in(jk), jk = 1, jpk )
-        ENDIF
-     ENDIF
-     !
-     CALL iom_close( inum )
-     !
-   END SUBROUTINE zgr_dom_read
-
-! =====================================================================================================
-
-   SUBROUTINE zgr_mes_local
-     !!---------------------------------------------------------------------
-     !!              ***  ROUTINE zgr_dom_merge  ***
-     !!
-     !! ** Purpose :   Create a local MEs grid within a gloabal grid 
-     !!                using different vertical coordinates.
-     !!
-     !!----------------------------------------------------------------------
-     INTEGER                      ::   ji, jj, jk ! dummy loop index
-     INTEGER                      ::   inum       ! local logical unit
-     REAL(wp), DIMENSION(jpi,jpj) ::   s2z_msk    ! mask for MEs-area (=2),
-                                                  !          transition area (=1),
-                                                  !          z-area (=0)
-     REAL(wp), DIMENSION(jpi,jpj) ::   s2z_wgt    ! weigths for computing model
-                                                  ! levels in transition area
-     !!----------------------------------------------------------------------
-
-     IF(lwp) THEN
-       WRITE(numout,*)
-       WRITE(numout,*) '   zgr_mes_merge : read the vertical coordinates in domain_cfg_global.nc file'
-       WRITE(numout,*) '   ~~~~~~~~'
-     ENDIF
-     !
-     CALL iom_open( 'bathy_meter.nc', inum )
-     CALL iom_get( inum, jpdom_data, 's2z_msk', s2z_msk)
-     CALL iom_get( inum, jpdom_data, 's2z_wgt', s2z_wgt)
-
-     mes_msk(:,:) = 1.
-     WHERE (s2z_msk(:,:).eq.0.0)  mes_msk(:,:) = 0.0
-
-     DO jj = 1,jpj
-        DO ji = 1,jpi
-           SELECT CASE (INT(s2z_msk(ji,jj)))
-             CASE (0) ! global zps area
-                  mbathy (ji,jj  ) = mbathy_in(ji,jj)
-                  gdept_0(ji,jj,:) = gdept_in(ji,jj,:)
-                  gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:)
-                  e3t_0  (ji,jj,:) = e3t_in (ji,jj,:)
-                  e3w_0  (ji,jj,:) = e3w_in (ji,jj,:)
-                  e3u_0  (ji,jj,:) = e3u_in (ji,jj,:)
-                  e3v_0  (ji,jj,:) = e3v_in (ji,jj,:)
-                  e3f_0  (ji,jj,:) = e3f_in (ji,jj,:)
-                  e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:)
-                  e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:)
-             CASE (1) ! MEs to zps transition area 
-                  gdept_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdept_0(ji,jj,:) + &
-                    &                ( 1._wp - s2z_wgt(ji,jj) ) * gdept_in(ji,jj,:)
-                  gdepw_0(ji,jj,:) =           s2z_wgt(ji,jj)   * gdepw_0(ji,jj,:) + &
-                    &                ( 1._wp - s2z_wgt(ji,jj) ) * gdepw_in(ji,jj,:)
-             CASE (2) ! MEs area
-                  CYCLE     
-           END SELECT
-        END DO
-     END DO
-     !
-     ! e3t, e3w for transition zone
-     ! as finite differences
+
+     e3psm = rn_e3pst_min
+     e3psr = rn_e3pst_rat
+
+     CALL wrk_alloc( jpi, jpj, pst_msk, apst)
+     CALL wrk_alloc( jpi, jpj, jpk, gdept_mes, gdepw_mes, e3t_mes, e3w_mes)
+
+     ! Initialise the mask to identify where to apply partial steps
+     pst_msk(:,:) = 0.0
+     IF(ln_pst_mes)  WHERE (l2g_msk(:,:) == 2.0)  pst_msk(:,:) = 1.0
+     IF(ln_pst_l2g)  WHERE (l2g_msk(:,:) == 1.0)  pst_msk(:,:) = 1.0
+
+
+     ! Compute mbathy for ocean points (i.e. the number of ocean levels)
+     ! find the number of ocean levels such that the last level thickness
+     ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0
      DO jj = 1, jpj
         DO ji = 1, jpi
-           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
-              DO jk = 1,jpkm1
-                 e3t_0(ji,jj,jk)   = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk)
-                 e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk)
-              ENDDO
-              ! Surface
-              jk = 1
-              e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1))
+           IF ( pst_msk(ji,jj) == 1.0 .AND. bathy(ji,jj) > 0._wp) THEN
+              DO jk = jpkm1, 1, -1
+                 zdepth = gdepw_0(ji,jj,jk) + MIN( e3psm, e3t_0(ji,jj,jk)*e3psr )
+                 IF ( bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1
+              END DO
+            END IF
+        END DO
+     END DO
+
+     ! ===================================================== 
+     ! Scale factors and depth at T- and W- points
+     ! =====================================================
+     !
+     ! Backup to the reference MEs-coordinate
+     e3t_mes  (:,:,:) = e3t_0  (:,:,:)
+     e3w_mes  (:,:,:) = e3w_0  (:,:,:)
+     gdept_mes(:,:,:) = gdept_0(:,:,:)
+     gdepw_mes(:,:,:) = gdepw_0(:,:,:)
+     
+     ! Mask identifying points where we apply partial steps:
+     !  bottom cell is     a partial cell = 0 
+     !  bottom cell is not a partial cell = 1
+     apst(:,:) = 1._wp
+     
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           jk = mbathy(ji,jj)
+           !
+           ! We use partial steps only at ocean points where 
+           ! envelopes differ from the bottom topography 
+           IF ( pst_msk(ji,jj) == 1.0 .AND. jk > 0 ) THEN    
               !
-              ! Bottom
-              jk = jpk
-              e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk))
+              IF ( gdepw_0(ji,jj,jk+1) /= bathy(ji,jj) ) apst(ji,jj) = 0     
+              !
+              IF ( jk == jpkm1 ) THEN          ! max ocean level case
+                 zgdepw = bathy(ji,jj)
+                 ze3t   = bathy(ji,jj) - gdepw_mes(ji,jj,jk)
+                 ze3w   = 0.5_wp * e3w_mes(ji,jj,jk) * ( 1._wp + ( ze3t/e3t_mes(ji,jj,jk) ) )
+                 e3t_0(ji,jj,jk  ) = ze3t
+                 e3w_0(ji,jj,jk  ) = ze3w
+                 gdepw_0(ji,jj,jk+1) = zgdepw
+                 gdept_0(ji,jj,jk  ) = gdept_mes(ji,jj,jk-1) + ze3w
+                 gdept_0(ji,jj,jk+1) = gdept_0(ji,jj,jk) + ze3t
+                 !
+              ELSE                             ! standard case
+                 zgdepw   = gdepw_mes(ji,jj,jk  )
+                 zgdepwp1 = gdepw_mes(ji,jj,jk+1)
+                 zgdept   = gdept_mes(ji,jj,jk  )
+                 ze3t     = e3t_mes  (ji,jj,jk  )
+                 ze3w     = e3w_mes  (ji,jj,jk  )
+                 IF ( bathy(ji,jj) <= zgdepwp1 ) gdepw_0(ji,jj,jk+1) = bathy(ji,jj)
+                 !gm Bug?  check the gdepw_1d
+                 !       ... on jk
+                 gdept_0(ji,jj,jk) = zgdepw + ( gdepw_0(ji,jj,jk+1) - zgdepw      ) &
+                    &                       * ((zgdept-zgdepw) / (zgdepwp1-zgdepw))
+                 e3t_0  (ji,jj,jk) = ze3t   * (gdepw_0(ji,jj,jk+1) - zgdepw) / (zgdepwp1 - zgdepw)
+                 e3w_0  (ji,jj,jk) = 0.5_wp * (gdepw_0(ji,jj,jk+1) + zgdepwp1 - 2._wp * zgdepw) &
+                    &                       * (ze3w / (zgdepwp1 - zgdepw))
+                 !       ... on jk+1
+                 gdept_0(ji,jj,jk+1) = gdept_0(ji,jj,jk) + e3t_0(ji,jj,jk)
+              ENDIF
+           ENDIF
+        ENDDO
+     ENDDO     
+     !
+     ! set value at mbathy+1
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           jk = mbathy(ji,jj)
+           IF ( pst_msk(ji,jj) == 1.0 .AND. jk > 0) THEN
+               e3t_0  (ji,jj,jk+1) = e3t_0  (ji,jj,jk)
+               e3w_0  (ji,jj,jk+1) = e3t_0  (ji,jj,jk)
            END IF
         END DO
      END DO
      !
-     ! MBATHY transition zone
-     DO jj = 1, jpj
-        DO ji = 1, jpi
-           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
-              DO jk = 1, jpkm1
-                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
-                 IF( scobot(ji,jj) == 0.e0              ) mbathy(ji,jj) = 0
-                 IF( scobot(ji,jj) < 0.e0               ) mbathy(ji,jj) = INT( scobot(ji,jj)) ! do we need it?
-              END DO
-           END IF
-        END DO
-     END DO
-     !
-     ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0 
-     ! for transition zone
-     !
+     ! =========================================== 
+     ! Scale factors at U-, V-, UW- and VW- points
+     ! ===========================================
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
-           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
-              DO jk = 1, jpk
-                 e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
-                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) /         &
-                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
-
-                 e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)      +         &
-                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) /         &
-                                  MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
-
-                 e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
-                                   MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) /         &
-                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))  )
-
-                 e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk)      +         &
-                                   MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) /         &
-                                   MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  )
-
-                 e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk)         +       &
-                                  MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk)      +       &
-                                  MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) +       &
-                                  MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) )    /       &
-                                  MAX(  1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))  &
-                                +   MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))  )
+           !
+           DO jp = 0,1
+              !  
+              ! e3u at the bottom
+              ! ~~~~~~~~~~~~~~~~~
+              jbku = mbathy(ji+jp,jj)
+              ! Only ocean points  
+              IF (pst_msk(ji,jj) == 1.0 .AND. jbku > 0) THEN
+                 ! work at bottom topography and level below
+                 DO jk = jbku, jbku+1
+                    ps_i   = apst(ji  ,jj  )   ! is bottom cell at T(ji  ,jj) point a partial cell or not
+                    ps_ip1 = apst(ji+1,jj  )   ! is bottom cell at T(ji+1,jj) point a partial cell or not
+                    ! correction for level jk
+                    IF ( (jk /= mbathy(ji  ,jj)) .AND. (jk /= mbathy(ji  ,jj)+1) ) ps_i   = 1
+                    IF ( (jk /= mbathy(ji+1,jj)) .AND. (jk /= mbathy(ji+1,jj)+1) ) ps_ip1 = 1
+                    !
+                    ps_p = ps_i * ps_ip1
+                    ps_s = ps_i + ps_ip1
+                    !
+                    e3t_i   = e3t_0(ji  ,jj,jk)
+                    e3t_ip1 = e3t_0(ji+1,jj,jk)
+                    e3w_i   = e3w_0(ji  ,jj,jk)
+                    e3w_ip1 = e3w_0(ji+1,jj,jk)
+                    !
+                    e3u_0(ji,jj,jk)  = MIN( e3u_0(ji,jj,jk) ,                              &
+                     &                        ps_p    * e3u_0(ji,jj,jk)                    & ! no zps cell
+                     &                     + (1-ps_p) * ( ps_ip1 * e3t_i + ps_i * e3t_ip1  & ! one of neighb is zps
+                     &                                   + (1-ps_s) * MIN(e3t_i,e3t_ip1) ) & ! both neighb are zps
+                     &                    )
+                     
+                    e3uw_0(ji,jj,jk) = MIN(e3uw_0(ji,jj,jk),                               &
+                     &                        ps_p    * e3uw_0(ji,jj,jk)                   & ! no zps 
+                     &                     + (1-ps_p) * ( ps_ip1 * e3w_i + ps_i * e3w_ip1  & ! one of neighb is zps
+                     &                                   + (1-ps_s) * MIN(e3w_i,e3w_ip1) ) & ! both side are zps
+                     &                    )
+                 END DO
+              END IF
+              ! 
+              ! e3v at the bottom
+              ! ~~~~~~~~~~~~~~~~~
+              jbkv = mbathy(ji,jj+jp)
+              ! Only ocean points  
+              IF (pst_msk(ji,jj) == 1.0 .AND. jbkv > 0) THEN
+                 ! work at bottom topography and level below
+                 DO jk = jbkv, jbkv+1
+                    ps_j   = apst(ji, jj  )   ! is bottom cell at T(ji,jj  ) point a partial cell or not
+                    ps_jp1 = apst(ji, jj+1)   ! is bottom cell at T(ji,jj+1) point a partial cell or not
+                    ! correction for level jk
+                    IF ( (jk /= mbathy(ji  ,jj)) .AND. (jk /= mbathy(ji  ,jj)+1) ) ps_j   = 1
+                    IF ( (jk /= mbathy(ji+1,jj)) .AND. (jk /= mbathy(ji+1,jj)+1) ) ps_jp1 = 1
+                    !
+                    ps_p = ps_j * ps_jp1
+                    ps_s = ps_j + ps_jp1
+                    !
+                    e3t_j   = e3t_0(ji,jj  ,jk)
+                    e3t_jp1 = e3t_0(ji,jj+1,jk)
+                    e3w_j   = e3w_0(ji,jj  ,jk)
+                    e3w_jp1 = e3w_0(ji,jj+1,jk)
+                    !
+                    e3v_0(ji,jj,jk)  = MIN( e3v_0(ji,jj,jk) ,                              &
+                     &                        ps_p    * e3v_0(ji,jj,jk)                    & ! no zps cell
+                     &                     + (1-ps_p) * ( ps_jp1 * e3t_j + ps_j * e3t_jp1  & ! one of neighb is zps
+                     &                                   + (1-ps_s) * MIN(e3t_j,e3t_jp1) ) & ! both neighb are zps
+                     &                    )
+
+                    e3vw_0(ji,jj,jk) = MIN(e3vw_0(ji,jj,jk),                               &
+                     &                        ps_p    * e3vw_0(ji,jj,jk)                   & ! no zps 
+                     &                     + (1-ps_p) * ( ps_jp1 * e3w_j + ps_j * e3w_jp1  & ! one of neighb is zps
+                     &                                   + (1-ps_s) * MIN(e3w_j,e3w_jp1) ) & ! both side are zps
+                     &                    )
+                 END DO
+              END IF
+           END DO
+        END DO
+     END DO
+
+     ! ========================= 
+     ! Scale factors at F-points
+     ! =========================   
+     DO jj = 1, jpjm1
+        DO ji = 1, jpim1
+           DO jp = 0,1
+              DO jl = 0,1
+                 jbkf = mbathy(ji+jl,jj+jp)
+                 ! Only ocean points  
+                 IF (pst_msk(ji,jj) == 1.0 .AND. jbkf > 0) THEN
+                    ! work at bottom topography and level below
+                    DO jk = jbkf, jbkf+1
+                       ps_uj   = MIN(apst(ji,jj  ),apst(ji+1,jj  )) ! is bot cell at U(ji,jj  ) pnt a part. cell or not
+                       ps_ujp1 = MIN(apst(ji,jj+1),apst(ji+1,jj+1)) ! is bot cell at U(ji,jj+1) pnt a part. cell or not
+                       ! correction for level jk
+                       IF (      (jk /= mbathy(ji  ,jj  )) .AND. (jk /= mbathy(ji  ,jj  )+1) &
+                         & .AND. (jk /= mbathy(ji+1,jj  )) .AND. (jk /= mbathy(ji+1,jj  )+1) ) ps_uj   = 1
+                       IF (      (jk /= mbathy(ji  ,jj+1)) .AND. (jk /= mbathy(ji  ,jj+1)+1) &
+                         & .AND. (jk /= mbathy(ji+1,jj+1)) .AND. (jk /= mbathy(ji+1,jj+1)+1) ) ps_ujp1 = 1
+                       !
+                       ps_p = ps_uj * ps_ujp1
+                       ps_s = ps_uj + ps_ujp1
+                       !
+                       e3u_j   = e3u_0(ji,jj  ,jk)
+                       e3u_jp1 = e3u_0(ji,jj+1,jk)
+                       !
+                       e3f_0(ji,jj,jk) = MIN( e3f_0(ji,jj,jk) ,                                &
+                        &                       ps_p    * e3f_0(ji,jj,jk)                      & ! no zps cell
+                        &                    + (1-ps_p) * ( ps_ujp1 * e3u_j + ps_uj * e3u_jp1  & ! one of neighb is zps
+                        &                                  + (1-ps_s) * MIN(e3u_j,e3u_jp1)   ) & ! both neighb are zps
+                        &                   )
+                    END DO
+                 END IF
               END DO
-           END IF
+           END DO
         END DO
      END DO
-     ! Computing gdep3w
-     gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
-     DO jk = 2, jpk
-        gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
-     END DO
-     !
-     ! From here equal to sco code - domzgr.F90 line 2079
-     ! Lateral B.C. we apply only for transition zone
-     ! since for s- and zps- areas has already been applied
-
-     CALL lbc_lnk( e3t_0 , 'T', 1._wp )
-     CALL lbc_lnk( e3u_0 , 'U', 1._wp )
-     CALL lbc_lnk( e3v_0 , 'V', 1._wp )
-     CALL lbc_lnk( e3f_0 , 'F', 1._wp )
-     CALL lbc_lnk( e3w_0 , 'W', 1._wp )
-     CALL lbc_lnk( e3uw_0, 'U', 1._wp )
-     CALL lbc_lnk( e3vw_0, 'V', 1._wp )
-
-     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
-
-     IF( lwp ) WRITE(numout,*) 'Refine mbathy'
-     DO jj = 1, jpj
-        DO ji = 1, jpi
-           IF ( s2z_msk(ji,jj) == 1._wp ) THEN
-              DO jk = 1, jpkm1
-                 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
-              END DO
-              IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
-           END IF
-        END DO
-     END DO
-
-   END SUBROUTINE zgr_mes_local
+     !
+     ! set to z-scale factor if zero (i.e. along closed boundaries)
+     DO jk = 1, jpk
+        WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_0(ji,jj,jk)
+     END DO
+     !
+     e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
+     e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
+     e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
+     e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
+     e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
+
+     ! Control of the sign
+     IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   e3t_0 <= 0' )
+     IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   e3w_0 <= 0' )
+     IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   gdept_0 < 0' )
+     IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '  zgr_pst_mes : e r r o r   gdepw_0 < 0' )
+
+     CALL wrk_dealloc( jpi, jpj, pst_msk, apst)
+     CALL wrk_dealloc( jpi, jpj, jpk, gdept_mes, gdepw_mes, e3t_mes, e3w_mes)
+
+   END SUBROUTINE zgr_pst_mes
 
 END MODULE zgrmes
