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 15125)
+++ NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90	(revision 15126)
@@ -4,16 +4,12 @@
    !! Ocean initialization : Multiple Enveloped s coordinate (MES)
    !!==============================================================================
-   !!  NEMO      3.6  ! 2017-05  (D. Bruciaferri)   
+   !!  NEMO      4.0  ! 2021-07  (D. Bruciaferri)   
    !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dom_mes          : define the MES vertical coordinate system
-   !!       fs_ma96      : tanh stretch function Madec 1996
-   !!       fs_sh94      : Song and Haidvogel 1994 stretch function
-   !!---------------------------------------------------------------------
+   !
    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
@@ -28,26 +24,15 @@
 
    PUBLIC   zgr_mes        ! called by domzgr.F90
-   ! Bruciaferri, Shapiro & Wobus (2018) envelopes and stretching parameters
-   CHARACTER(lc)      :: ctlmes                  ! control message error (lc=256)
-   INTEGER, PARAMETER :: max_nn_env = 5          ! Maximum allowed number of envelopes.
-   REAL(wp)           :: rn_bot_min              ! minimum depth of ocean bottom (>0) (m)
-   REAL(wp)           :: rn_bot_max              ! maximum depth of ocean bottom (= ocean depth) (>0) (m)
-   LOGICAL            :: ln_envl(max_nn_env)     ! Array with flags (T/F) specifing which envelope is used
-   INTEGER            :: nn_strt(max_nn_env)     ! Array specifing the stretching function for each envelope:
-                                                 ! Madec 1996 (0), Song and Haidvogel 1994 (1) 
-   REAL(wp)           :: rn_ebot_min(max_nn_env) ! minimum depth of envelopes (>0) (m)
-   REAL(wp)           :: rn_ebot_max(max_nn_env) ! maximum depth of envelopes (= envelopes depths) (>0) (m)
-   INTEGER            :: nn_slev(max_nn_env)     ! Array specifing number of levels of each enveloped vertical zone
-   REAL(wp)           :: rn_e_hc(max_nn_env)     ! Array specifing critical depth for transition to stretched 
-                                                 ! coordinates of each envelope
-   REAL(wp)           :: rn_e_th(max_nn_env)     ! Array specifing surface control parameter (0<=th<=20) of SH94
-                                                 ! or rn_zs of SF12 for each vertical sub-zone
-   REAL(wp)           :: rn_e_bb(max_nn_env)     ! Array specifing bottom control parameter (0<=bb<=1) of SH94
-                                                 ! or rn_zb_b of SF12 for each vertical sub-zone
-   REAL(wp)           :: rn_e_ba(max_nn_env)     ! Array specifing bottom control parameter rn_zb_a of SF12 for
-                                                 ! each vertical sub-zone
-   REAL(wp)           :: rn_e_al(max_nn_env)     ! Array specifing stretching parameter rn_alpha of SF12 for
-                                                 ! each vertical sub-zone
-   LOGICAL            :: ln_loc_mes              ! To use localised MEs (.TRUE.) or not (.FALSE.
+   ! 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
@@ -68,11 +53,4 @@
       ! 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(:,:,:) :: 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]
       !-----------------------------------------------------------------------
       !
@@ -80,16 +58,22 @@
       CALL zgr_mes_build
 
-      if( ln_loc_mes ) THEN
-         !   Local MEs
+      ! 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(ln_isfcav,                                      &  ! ice-cavities
-                 &         gdept_1d_in, gdepw_1d_in, e3t_1d_in, e3w_1d_in, &  ! 1D ref. vert. coord.
-                 &         gdept_in   , gdepw_in   ,                       &  ! 3D t & w-points depth
-                 &         e3t_in     , e3u_in     , e3v_in   , e3f_in   , &  ! vertical scale factors
-                 &         e3w_in     , e3uw_in    , e3vw_in               )   
+         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)
+         !
       END IF
 
@@ -98,955 +82,5 @@
 ! =====================================================================================================
 
-   SUBROUTINE zgr_mes_build
-      !!-----------------------------------------------------------------------------
-      !!                  ***  ROUTINE zgr_mes  ***
-      !!                     
-      !! ** Purpose :   define the Multi Enveloped S-coordinate (MES) system
-      !!
-      !! ** Method  :   s-coordinate defined respect to multiple
-      !!                externally defined envelope
-      !!
-      !!      The depth of model levels is defined as the product of an
-      !!      analytical function by the externally defined envelopes.
-      !! 
-      !!      - Read envelopes (in meters) at t-point and compute the
-      !!        envelopes at u-, v-, and f-points.
-      !!            hbatu = mi( hbatt )
-      !!            hbatv = mj( hbatt )
-      !!            hbatf = mi( mj( hbatt ) )
-      !!
-      !!      - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
-      !!        function and its finite difference derivative.
-      !!            z_gsigt(k) = fssig (k    )
-      !!            z_gsigw(k) = fssig (k-0.5)
-      !!            z_esigt(k) = z_gsigw(k+1) - z_gsigw(k)
-      !!            z_esigw(k) = z_gsigt(k+1) - z_gsigt(k)
-      !!
-      !!      Three options for stretching are given:
-      !! 
-      !!          *) nn_strt = 0 : Madec et al 1996 cosh/tanh function
-      !!          *) nn_strt = 1 : Song and Haidvogel 1994 sinh/tanh function  
-      !!          *) nn_strt = 2 : Siddorn and Furner gamma function
-      !!
-      !!----------------------------------------------------------------------
-      !!        SKETCH of the GEOMETRY OF THE S-S-Z COORDINATE SYSTEM
-      !!
-      !!        Lines represent W-levels:
-      !!
-      !!        0: First s-level part. The total number
-      !!           of levels in this part is controlled
-      !!           by the nn_slev(1) namelist parameter.
-      !!
-      !!           Depth first W-lev: 0 m (surfcace)
-      !!           Depth last  W-lev: depth of 1st envelope
-      !!
-      !!        o: Second  s-level part. The total number
-      !!           of levels in this part is controlled
-      !!           by the nn_slev(2) namelist parameter.
-      !!        
-      !!           Depth last  W-lev: depth of 2nd envelope
-      !!
-      !!        @: Third s-level part. The total number
-      !!           of levels in this part is controlled
-      !!           by the nn_slev(3) namelist parameter.
-      !!        
-      !!           Depth last  W-lev: depth of 3rd envelope  
-      !!    
-      !! z |~~~~~~~~~~~~~~~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~~~~          SURFACE                 
-      !!   | 
-      !!   |-----------------------0-----------------------   ln_envl(1) = true, nn_slev(1) = 3
-      !!   |
-      !!   |=======================0=======================          ENVELOPE 1
-      !!   | 
-      !!   |-----------------------o-----------------------
-      !!   |                                                 ln_envl(2) = true, nn_slev(2) = 3
-      !!   |-----------------------o-----------------------
-      !!   |
-      !!   |¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬o¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬          ENVELOPE 2
-      !!   |
-      !!   |-----------------------@-----------------------  ln_envl(3) = true, nn_slev(3) = 2
-      !!   |
-      !!   |_______________________@_______________________          ENVELOPE 3
-      !!  \|/
-      !!
-      !!
-      !!-----------------------------------------------------------------------------
-      !! Author: Diego Bruciaferri (University of Plymouth), September 2017
-      !!-----------------------------------------------------------------------------
-
-      !
-      CHARACTER(lc)                       ::   env_name                 ! name of the externally defined envelope
-      INTEGER                             ::   ji, jj, jk, jl, je       ! dummy loop argument
-      INTEGER                             ::   ios                      ! Local integer output status for namelist read
-      INTEGER                             ::   inum                     ! temporary logical unit
-      INTEGER                             ::   tot_env, cor_env
-      INTEGER                             ::   num_s, s_1st, ind        ! for loops over envelopes
-      INTEGER                             ::   num_s_up, num_s_dw       ! for loops over envelopes
-      REAL(wp)                            ::   D1s_up, D1s_dw           ! for loops over envelopes
-      INTEGER                             ::   ibcbeg, ibcend           ! boundary condition for cubic splines
-      INTEGER, PARAMETER                  ::   n = 2                    ! number of considered points for each interval
-      REAL(wp)                            ::   c(4,n)                   ! matrix for cubic splines coefficents
-      REAL(wp)                            ::   d(2,2)                   ! matrix for depth values and depth derivative 
-                                                                        ! at intervals' boundaries
-      REAL(wp)                            ::   tau(n)                   ! abscissas of intervals' boundaries
-      REAL(wp)                            ::   coefa, coefb             ! for loops over envelopes
-      REAL(wp)                            ::   alpha, env_hc            ! for loops over envelopes
-      REAL(wp)                            ::   dep_up, dep_dw           ! for loops over envelopes
-      REAL(wp)                            ::   dep0, dep1, dep2, dep3   ! for loops over envelopes 
-                                                                        ! for cubic splines
-      REAL(wp)                            ::   zcoeft, zcoefw, sw, st   ! temporary scalars
-      REAL(wp)                            ::   max_env_up, min_env_up   ! to identify geopotential levels
-      REAL(wp)                            ::   max_env_dw, min_env_dw   ! to identify geopotential levels
-      REAL(wp), POINTER, DIMENSION(:)     ::   z_gsigw1, z_gsigt1       ! MES-coordinates analytical trans. 
-                                                                        ! (0.<= s <= 1.) for T- and W-points
-      REAL(wp), POINTER, DIMENSION(:)     ::   z_esigw1, z_esigt1       ! MES analytical scale factors  
-                                                                        ! of T- and W-points
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: env_up, env_dw             ! for loops over envelopes
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: env0, env1, env2, env3     ! for loops over envelopes
-                                                                        ! for cubic splines
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: envlt                      ! array for the envelopes
-      INTEGER                             :: gst_envl(max_nn_env)       ! Array to deal with a ghost last envelope 
-
-      NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl, nn_strt, &
-                          rn_ebot_min, rn_ebot_max, nn_slev, rn_e_hc, &
-                          rn_e_th, rn_e_bb, rn_e_ba, rn_e_al, ln_loc_mes
-
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_mes')
-      !
-      CALL wrk_alloc( jpi, jpj, max_nn_env , envlt )
-      CALL wrk_alloc( jpi, jpj, env_up, env_dw, env0, env1, env2, env3 )
-      CALL wrk_alloc( jpk, z_gsigw1, z_gsigt1, z_esigw1, z_esigt1 )
-      !
-      ! Namelist namzgr_mes in reference namelist : envelopes and
-      ! sigma-stretching parameters
-      REWIND( numnam_ref )         
-      READ  ( numnam_ref, namzgr_mes, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in reference namelist', lwp )
-      ! Namelist namzgr_mes in configuration namelist : envelopes and
-      ! sigma-stretching parameters
-      REWIND( numnam_cfg )
-      READ  ( numnam_cfg, namzgr_mes, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzgr_mes )
-
-      IF(lwp) THEN                           ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr_mes : Multi Enveloped S-coordinate (Bruciaferri, Shapiro and Wobus 2017)'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namzgr_mes'
-         WRITE(numout,*) ''
-         WRITE(numout,*) '   Minimum depth of the ocean   rn_bot_min, ', rn_bot_min
-         WRITE(numout,*) '   Maximum depth of the ocean   rn_bot_max, ', rn_bot_max
-         WRITE(numout,*) ''
-         WRITE(numout,*) '-------------------------------------------------------------------------------'
-         DO je = 1, max_nn_env
-            WRITE(numout,*) 'SUBDOMAIN ', je,':'
-            IF ( je == 1) THEN
-               WRITE(numout,*) '   Envelope up  : envelope 0, free surface'
-               WRITE(numout,*) '   Envelope down: envelope ',je,',ln_envl(',je,') = ',ln_envl(je)
-            ELSE
-               WRITE(numout,*) '   Envelope up  : envelope ',je-1,',ln_envl(',je,') = ',ln_envl(je-1)
-               WRITE(numout,*) '   Envelope down: envelope ',je  ,',ln_envl(',je,') = ',ln_envl(je)
-            END IF
-            WRITE(numout,*) '   min dep of envlp down rn_ebot_min(',je,') = ',rn_ebot_min(je)
-            WRITE(numout,*) '   max dep of envlp down rn_ebot_max(',je,') =',rn_ebot_max(je)
-            WRITE(numout,*) '   num. of MEs-lev.          nn_slev(',je,') = ',nn_slev(je)
-            IF ( isodd(je) ) THEN
-               WRITE(numout,*) '   Stretched s-coordinates: '
-            ELSE
-               WRITE(numout,*) '   Stretched CUBIC SPLINES: '
-            END IF
-            IF (nn_strt(je) == 0) WRITE(numout,*) '     M96  stretching function'
-            IF (nn_strt(je) == 1) WRITE(numout,*) '     SH94 stretching function'
-            IF (nn_strt(je) == 2) WRITE(numout,*) '     SF12 stretching function'
-            WRITE(numout,*) '     critical depth        rn_e_hc(',je,') = ',rn_e_hc(je)
-            WRITE(numout,*) '     surface stretc. coef. rn_e_th(',je,') = ',rn_e_th(je)
-            IF (nn_strt(je) == 2) THEN
-               WRITE(numout,*) '     bottom  stretc. coef. rn_e_ba(',je,') = ',rn_e_ba(je)
-            END IF
-            WRITE(numout,*) '     bottom  stretc. coef. rn_e_bb(',je,') = ',rn_e_bb(je)
-            IF (nn_strt(je) == 2) THEN
-               WRITE(numout,*) '     bottom  stretc. coef. rn_e_al(',je,') = ',rn_e_al(je)
-            END IF
-            WRITE(numout,*) '-------------------------------------------------------------------------------'
-         ENDDO
-      ENDIF
-
-      ! Check if namelist is defined correctly.
-      ! Not strictly needed but we force the user
-      ! to define the namelist correctly.
-
-      tot_env = 0 ! total number of requested envelopes
-      cor_env = 0 ! total number of correctly defined envelopes
-      DO je = 1, max_nn_env
-         IF ( ln_envl(je) )                        tot_env = tot_env + 1
-         IF ( ln_envl(je) .AND. cor_env == (je-1)) cor_env = cor_env + 1
-      ENDDO
-      WRITE(ctlmes,*) 'number of REQUESTED envelopes and number of CORRECTLY defined envelopes are DIFFERENT'
-      IF ( tot_env /= cor_env ) CALL ctl_stop( ctlmes )
-
-      !Checking consistency of user defined parameters
-      WRITE(ctlmes,*) 'TOT number of levels (jpk) IS DIFFERENT from sum over nn_slev(:)'
-      IF ( SUM(nn_slev(:)) /= jpk ) CALL ctl_stop( ctlmes )
-
-      ! Determining if there is a "ghost" envelope: 
-      gst_envl(:) = 0
-      IF ( nn_slev(tot_env) == 0 ) gst_envl(tot_env) = 1
-
-      ! Reading bathymetry and envelopes.
-      ! In future should be included in zgr_bat.
-
-      IF( ntopo == 1 ) THEN
-
-        CALL iom_open ( 'bathy_meter.nc', inum )
-        CALL iom_get  ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr  )
-
-        DO je = 1, tot_env
-           WRITE (env_name, '(A6,I0)') 'hbatt_', je
-           CALL iom_get  ( inum, jpdom_data, TRIM(env_name) , envlt(:,:,je), lrowattr=ln_use_jattr  )
-        ENDDO
-
-        CALL iom_close( inum )
-
-      ELSE
-
-        WRITE(ctlmes,*) 'parameter , ntopo = ', ntopo
-        CALL ctl_stop( ctlmes )
-
-      ENDIF
-      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy ) !==  NO closed seas or lakes  ==!  
-
-      ! Checking consistency of envelopes
-
-      DO je = 1, tot_env-1
-         WRITE(ctlmes,*) 'Envelope ', je+1, ' is shallower that Envelope ', je
-         IF (MAXVAL(envlt(:,:,je+1)) < MAXVAL(envlt(:,:,je))) CALL ctl_stop( ctlmes )
-      ENDDO
-
-      ! Set maximum and minimum ocean depth
-      bathy(:,:) = MIN( rn_bot_max, bathy(:,:) )
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_bot_min, bathy(ji,jj) )
-         END DO
-      END DO
-
-      ! Initializing to 0.0 some arrays:
-      !
-      ! *) gdept_1d, gdepw_1d : depth of T- and W-point (m)
-      ! *) e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
-      ! *) z_gsigt1, z_gsigw1 : MES-coordinates (0.<= s <= 1.) of T- and W-points
-
-      gdept_1d = 0._wp   ;   gdepw_1d = 0._wp   ;
-      e3t_1d   = 0._wp   ;   e3w_1d   = 0._wp   ;
-      z_gsigw1 = 0._wp   ;   z_gsigt1 = 0._wp   ;
-      z_esigw1 = 0._wp   ;   z_esigt1 = 0._wp   ;
-
-      ! Initialize mbathy to the maximum ocean level available
-
-      mbathy(:,:) = jpkm1
-
-      ! Defining scosrf and scobot
-
-      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
-      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
-
-      !==================================================================
-      ! 1. Defining 1D MEs-levels and computing their depths 
-      !    (they will be used in diawri.F90)
-      !==================================================================
-
-      ! Computing gdept_1d and gdepw_1d
-
-      DO je = 1, tot_env ! LOOP over all the requested envelopes
-         IF (je == 1) THEN
-            s_1st  = 1
-            num_s  = nn_slev(je)
-            dep_up = 0._wp      ! surface
-         ELSE
-            s_1st = s_1st + num_s - 1
-            num_s  = nn_slev(je)  + 1
-            dep_up = rn_ebot_max(je-1)
-         ENDIF
-
-         dep_dw  = rn_ebot_max(je)
-         env_hc  = rn_e_hc(je)
-         coefa   = rn_e_th(je)
-         coefb   = rn_e_bb(je)
-         alpha   = rn_e_al(je)
-
-         IF (nn_strt(je) == 2) THEN
-            coefa = coefa / (dep_dw - dep_up)
-            coefb = coefb + ((dep_dw - dep_up)*rn_e_ba(je))
-            coefb = 1.0_wp - (coefb/(dep_dw - dep_up))
-         ENDIF
-
-         IF ( isodd(je) ) THEN
-            IF ( gst_envl(je) == 0 ) THEN
-               DO jk = 1, num_s
-                  ind = s_1st+jk-1
-                  ! W-GRID
-                  z_gsigw1(ind) = -s_coord(sigma(1,jk,num_s), coefa, &
-                                           coefb, alpha, num_s, nn_strt(je))
-                  ! T-GRID
-                  z_gsigt1(ind) = -s_coord(sigma(0,jk,num_s), coefa, &
-                                           coefb, alpha, num_s, nn_strt(je))
-                  zcoefw = z_gsigw1(ind)
-                  zcoeft = z_gsigt1(ind)
-                              
-                  IF (nn_strt(je) /= 2) THEN
-                     gdept_1d(ind) = z_mes(dep_up, dep_dw, zcoeft, &
-                                           -sigma(0,jk,num_s), env_hc)
-                     gdepw_1d(ind) = z_mes(dep_up, dep_dw, zcoefw, &
-                                           -sigma(1,jk,num_s), env_hc)
-                  ELSE
-                     gdept_1d(ind) = z_mes(dep_up, dep_dw, zcoeft, &
-                                           -sigma(0,jk,num_s), 0._wp)
-                     gdepw_1d(ind) = z_mes(dep_up, dep_dw, zcoefw, &
-                                           -sigma(1,jk,num_s), 0._wp)
-                  END IF
-               ENDDO
-            ENDIF
-
-         ELSE
-
-            ! Ensure that for splines we do not use SF12
-            !IF (nn_strt(je) > 1) nn_strt(je) = 1
-
-            ! Interval's endpoints TAU are 0 and 1.
-            tau(1) = 0._wp
-            tau(n) = 1._wp
-
-            IF ( je == 2 ) THEN
-               dep0 = 0._wp
-            ELSE 
-               dep0 = rn_ebot_max(je-2)
-            END IF 
-            dep1   = dep_up
-            dep2   = dep_dw
-            dep3   = rn_ebot_max(je+1)
-
-            ! Computing derivative analytically
-            !
-            ! 1. Derivative for upper envelope
-            ibcbeg = 1
-            IF (nn_strt(je-1) == 2) THEN
-               alpha    = rn_e_al(je-1)
-               num_s_up = nn_slev(je-1)
-               coefa    = rn_e_th(je-1) / (dep1-dep0)
-               coefb    = rn_e_bb(je-1) + ((dep1-dep0)*rn_e_ba(je-1))
-               coefb    = 1.0_wp-(coefb/(dep1-dep0))
-            ELSE
-               alpha    = 0._wp
-               num_s_up = 1
-               coefa    = rn_e_th(je-1)
-               coefb    = rn_e_bb(je-1)
-            END IF
-            D1s_up = D1s_coord(-1._wp, coefa, coefb, &
-                               alpha, num_s_up, nn_strt(je-1))
-            d(1,1) = dep_up
-            IF (nn_strt(je-1) == 2) THEN
-               d(2,1) = D1z_mes(dep0, dep1, D1s_up, 0._wp, 1)
-            ELSE
-               d(2,1) = D1z_mes(dep0, dep1, D1s_up, rn_e_hc(je-1), 1)
-            END IF
-            !
-            ! 2. Derivative for deeper envelope
-            IF ( gst_envl(je+1) == 0 ) THEN
-               ibcend = 1 ! Bound. cond for the 1st derivative          
-               IF (nn_strt(je+1) == 2) THEN
-                  alpha    = rn_e_al(je+1)
-                  num_s_dw = nn_slev(je+1)
-                  coefa    = rn_e_th(je+1) / (dep3-dep2)
-                  coefb    = rn_e_bb(je+1) + ((dep3-dep2)*rn_e_ba(je+1))
-                  coefb    = 1.0_wp-(coefb/(dep3-dep2))
-               ELSE
-                  alpha    = 0._wp
-                  num_s_dw = 1
-                  coefa    = rn_e_th(je+1)
-                  coefb    = rn_e_bb(je+1)
-               END IF
-               D1s_dw = D1s_coord(0._wp, coefa, coefb, &
-                                  alpha, num_s_dw, nn_strt(je+1))
-               d(1,2) = dep_dw
-               IF (nn_strt(je+1) == 2) THEN
-                  d(2,2) = D1z_mes(dep2, dep3, D1s_dw, 0._wp, 1)
-               ELSE
-                  d(2,2) = D1z_mes(dep2, dep3, D1s_dw, rn_e_hc(je+1), 1)
-               END IF
-            ELSE
-               ibcend = 2     ! Bound. cond for the 2nd derivative
-               d(2,2) = 0._wp ! 2nd der = 0
-            ENDIF
-            !
-            ! Computing spline's coefficients
-            IF ( lwp ) THEN
-               WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:"
-               WRITE(numout,*) "ibcbeg: ", ibcbeg
-               WRITE(numout,*) "ibcend: ", ibcend
-            ENDIF
-            c = cub_spl(tau, d, n, ibcbeg, ibcend )
-
-            IF ( lwp ) THEN
-               WRITE(numout,*) "Subzone between depth ", dep1, " and depth ", dep2 
-               WRITE(numout,*) ""
-               WRITE(numout,*) "D1s_up = ", D1s_up
-               WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep1, " is ", d(2,1)
-               WRITE(numout,*) "D1s_dw = ", D1s_dw
-               WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep2, " is ", d(2,2)
-               WRITE(numout,*) ""
-               WRITE(numout,*) "CUBIC SPLINE COEFFICENTS"
-               WRITE(numout,*) "c1 = ", c(1,1)
-               WRITE(numout,*) "c2 = ", c(2,1)
-               WRITE(numout,*) "c3 = ", c(3,1)
-               WRITE(numout,*) "c4 = ", c(4,1)
-               WRITE(numout,*) ""
-               WRITE(numout,*) "CUB. SPL. at s = ", -sigma(1, 1, num_s),      &
-                               " is ", z_cub_spl(-sigma(1, 1, num_s), tau, c)
-               WRITE(numout,*) "1st deriv. of CUB. SPL. at s = ", -sigma(1, 1, num_s) ,     & 
-                               " is ", D1z_cub_spl(-sigma(1, 1, num_s), tau, c, 1)
-               WRITE(numout,*) "CUB. SPL. at s = ", -sigma(1, num_s, num_s) , &
-                               " is ", z_cub_spl(-sigma(1, num_s, num_s), tau, c)
-               WRITE(numout,*) "1st deriv. of CUB. SPL. at s = ", -sigma(1, num_s, num_s) , &
-                               " is ", D1z_cub_spl(-sigma(1, num_s, num_s), tau, c, 1)
-            END IF
-
-            env_hc  = rn_e_hc(je)
-            coefa   = rn_e_th(je)
-            coefb   = rn_e_bb(je)
-            alpha   = rn_e_al(je)
-
-            DO jk = 1, num_s
-               ind = s_1st+jk-1
-               ! W-GRID
-               z_gsigw1(s_1st+jk-1) = -s_coord( sigma(1,jk,num_s), coefa, &
-                                                coefb, alpha, num_s, nn_strt(je) )
-               ! T-GRID
-               z_gsigt1(s_1st+jk-1) = -s_coord( sigma(0,jk,num_s), coefa, &
-                                                coefb, alpha, num_s, nn_strt(je) )
-               st  = z_gsigt1(ind)
-               sw  = z_gsigw1(ind)
-               gdept_1d(ind) = z_cub_spl(st, tau, c)
-               gdepw_1d(ind) = z_cub_spl(sw, tau, c)
-            END DO
- 
-         END IF
-
-      END DO
-
-      !==================================================================
-      ! 2. Computing derivative of MEs-coordinate and scale factors
-      !    as finite differences
-      !==================================================================
-
-      ! Computing derivative of MEs-coordinate as finite differences
-      DO jk = 1, jpkm1
-         z_esigt1(jk)   = z_gsigw1(jk+1) - z_gsigw1(jk)
-         z_esigw1(jk+1) = z_gsigt1(jk+1) - z_gsigt1(jk)
-      END DO
-      ! Surface
-      z_esigw1( 1 ) = 2._wp * ( z_gsigt1( 1 ) - z_gsigw1( 1 ) )
-      ! Bottom
-      z_esigt1(jpk) = 2._wp * ( z_gsigt1(jpk) - z_gsigw1(jpk) )
-
-      ! Computing scale factors as finite differences
-      DO jk=1,jpkm1
-         e3t_1d(jk)   = gdepw_1d(jk+1) - gdepw_1d(jk)
-         e3w_1d(jk+1) = gdept_1d(jk+1) - gdept_1d(jk)
-      ENDDO
-      ! Bottom:
-      jk = jpk
-      e3t_1d(jk) = 2.0_wp * (gdept_1d(jk) - gdepw_1d(jk))
-      ! and Surface  (H. Liu)
-      jk = 1
-      e3w_1d(jk) = 2.0_wp * (gdept_1d(1) - gdepw_1d(1))
-
- 
-      ! Control print
-      IF ( lwp ) THEN
-        WRITE(numout,*) ""
-        WRITE(numout,*) "1a. MEs-levels and first derivative"
-        WRITE(numout,*) ""
-        WRITE(numout,*) "      k     z_gsigw1      z_esigw1     z_gsigt1      z_esigt1"
-        WRITE(numout,*) ""
-        DO jk = 1, jpk
-           WRITE(numout,*) jk, z_gsigw1(jk), z_esigw1(jk), z_gsigt1(jk), z_esigt1(jk)
-        ENDDO
-        WRITE(numout,*) "-----------------------------------------------------------------"
-        
-        WRITE(numout,*) ""
-        WRITE(numout,*) "2a. Depth of MEs-levels and scale factors"
-        WRITE(numout,*) ""
-        WRITE(numout,*) "      k     gdepw_1d      e3w_1d       gdepw_1d      e3w_1d"
-        WRITE(numout,*) ""
-        DO jk = 1, jpk
-           WRITE(numout,*) jk, gdepw_1d(jk), e3w_1d(jk), gdept_1d(jk), e3t_1d(jk)
-        ENDDO
-        WRITE(numout,*) "-----------------------------------------------------------------"
-      ENDIF
-
-      ! Checking monotonicity for cubic splines
-      DO jk = 1, jpk-1
-         IF ( gdept_1d(jk+1) < gdept_1d(jk) ) THEN
-            WRITE(ctlmes,*) 'NOT MONOTONIC gdept_1d: change envelopes'
-            CALL ctl_stop( ctlmes )
-         END IF
-         IF ( gdepw_1d(jk+1) < gdepw_1d(jk) ) THEN
-            WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_1d: change envelopes'
-            CALL ctl_stop( ctlmes )
-         END IF
-      END DO
-
-      !========================
-      ! 3. Compute 3D MEs mesh
-      !========================
-
-      DO je = 1, tot_env ! LOOP over all the requested envelopes
-
-         IF (je == 1) THEN
-            s_1st  = 1
-            num_s  = nn_slev(je)
-            env_up = scosrf       ! surface
-         ELSE
-            s_1st  = s_1st + num_s - 1
-            num_s  = nn_slev(je) + 1
-            env_up = envlt(:,:,je-1)
-         ENDIF
-
-         env_dw  = envlt(:,:,je)
-         env_hc  = rn_e_hc(je)
-         alpha   = rn_e_al(je)
-
-         IF ( isodd(je) ) THEN
-
-            IF ( gst_envl(je) == 0 ) THEN
-               DO jk = 1, num_s
-                  ind = s_1st+jk-1
-                  DO jj = 1,jpj
-                     DO ji = 1,jpi
-                        IF (env_dw(ji,jj) < env_hc) THEN   
-                           ! shallow water, uniform sigma
-                           zcoefw = -sigma(1, jk, num_s) 
-                           zcoeft = -sigma(0, jk, num_s)
-                           !IF (nn_strt(je) == 2) THEN
-                           !   ! shallow water, z coordinates
-                           !   ! SF12 in AMM* configurations uses 
-                           !   ! this (ln_sigcrit=.false.)
-                           !   zcoefw = zcoefw*(env_hc/(env_dw(ji,jj)-env_up(ji,jj)))
-                           !   zcoeft = zcoeft*(env_hc/(env_dw(ji,jj)-env_up(ji,jj)))
-                           !ENDIF
-                        ELSE                               ! deep water, stretched s
-                           IF (nn_strt(je) == 2) THEN
-                              coefa = rn_e_th(je) /  (env_dw(ji,jj)-env_up(ji,jj))
-                              coefb = rn_e_bb(je) + ((env_dw(ji,jj)-env_up(ji,jj))*rn_e_ba(je))
-                              coefb = 1.0_wp-(coefb/(env_dw(ji,jj)-env_up(ji,jj)))
-                              ! W-GRID
-                              zcoefw = -s_coord( sigma(1,jk,num_s), coefa, &
-                                                 coefb, alpha, num_s, nn_strt(je) )
-                              ! T-GRID
-                              zcoeft = -s_coord( sigma(0,jk,num_s), coefa, &
-                                                 coefb, alpha, num_s, nn_strt(je) )
-                           ELSE
-                              zcoefw = z_gsigw1(ind)
-                              zcoeft = z_gsigt1(ind)
-                           END IF
-                        ENDIF
-
-                        IF (nn_strt(je) /= 2) THEN
-                           gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, &
-                                                      -sigma(0, jk, num_s), env_hc)
-                                                
-                           gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, &
-                                                      -sigma(1, jk, num_s), env_hc)
-                        ELSE
-                           gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, &
-                                                      -sigma(0, jk, num_s), 0._wp)
-
-                           gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, &
-                                                      -sigma(1, jk, num_s), 0._wp)
-                        END IF
-                     END DO
-                  END DO
-               END DO
-            END IF
-
-         ELSE
-
-            ! Interval's endpoints TAU are 0 and 1.
-            tau(1) = 0._wp
-            tau(n) = 1._wp
-
-            IF ( je == 2 ) THEN
-               env0 = scosrf
-            ELSE
-               env0 = envlt(:,:,je-2) 
-            END IF
-            env1   = env_up
-            env2   = env_dw
-            env3   = envlt(:,:,je+1)
-
-            DO jj = 1,jpj
-               DO ji = 1,jpi
-                  d(1,1) = env_up(ji,jj)
-                  d(1,2) = env_dw(ji,jj)
-                  !
-                  ! Computing derivative analytically
-                  !
-                  ! 1. Derivative for upper envelope
-                  ibcbeg = 1
-                  IF (nn_strt(je-1) == 2) THEN
-                     alpha    = rn_e_al(je-1)
-                     num_s_up = nn_slev(je-1)
-                     coefa    = rn_e_th(je-1) / &
-                                (env1(ji,jj)-env0(ji,jj))
-                     coefb    = rn_e_bb(je-1) + &
-                               ((env1(ji,jj)-env0(ji,jj))*rn_e_ba(je-1))
-                     coefb    = 1.0_wp-(coefb/(env1(ji,jj)-env0(ji,jj)))                   
-                  ELSE
-                     alpha    = 0._wp
-                     num_s_up = 1
-                     coefa    = rn_e_th(je-1)
-                     coefb    = rn_e_bb(je-1)
-                  END IF
-                  D1s_up = D1s_coord(-1._wp, coefa, coefb, &
-                                     alpha, num_s_up, nn_strt(je-1))
-                  IF (nn_strt(je-1) == 2) THEN
-                     d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, &
-                                      0._wp, 1)
-                  ELSE
-                     d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, &
-                                      rn_e_hc(je-1), 1)
-                  END IF
-                  !
-                  ! 2. Derivative for lower envelope
-                  IF ( gst_envl(je+1) == 0 ) THEN
-                     ibcend = 1 ! Bound. cond for the 1st derivative
-                     IF (nn_strt(je+1) == 2) THEN
-                        alpha    = rn_e_al(je+1)
-                        num_s_dw = nn_slev(je+1)
-                        coefa    = rn_e_th(je+1) / &
-                                   (env3(ji,jj)-env2(ji,jj))
-                        coefb    = rn_e_bb(je+1) + &
-                                   ((env3(ji,jj)-env2(ji,jj))*rn_e_ba(je+1))
-                        coefb    = 1.0_wp-(coefb/(env3(ji,jj)-env2(ji,jj)))
-                     ELSE
-                        alpha    = 0._wp
-                        num_s_dw = 1
-                        coefa    = rn_e_th(je+1)
-                        coefb    = rn_e_bb(je+1)
-                     END IF
-                     D1s_dw = D1s_coord(0._wp, coefa, coefb, &
-                                        alpha, num_s_dw, nn_strt(je+1))
-                  
-                     IF (nn_strt(je+1) == 2) THEN
-                        d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, &
-                                         0._wp, 1)
-                     ELSE
-                        d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, &
-                                         rn_e_hc(je+1), 1)
-                     END IF
-                  ELSE
-                     ibcend = 2     ! Bound. cond for the 2nd derivative
-                     d(2,2) = 0._wp ! 2nd der = 0
-                  ENDIF
-
-                  ! Computing spline's coefficients
-                  !IF ( lwp ) THEN
-                  !   WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:"
-                  !   WRITE(numout,*) "ibcbeg: ", ibcbeg
-                  !   WRITE(numout,*) "ibcend: ", ibcend
-                  !ENDIF
-                  c = cub_spl(tau, d, n, ibcbeg, ibcend )
-
-                  ! Computing model level depths
-                  DO jk = 1, num_s
-                     ind = s_1st+jk-1
-                     st  = z_gsigt1(ind)
-                     sw  = z_gsigw1(ind)
-                     gdept_0(ji,jj,ind) = z_cub_spl(st, tau, c)
-                     gdepw_0(ji,jj,ind) = z_cub_spl(sw, tau, c)
-                  END DO
-
-                  ! Control print
-                  IF ( lwp ) THEN
-                     IF ( jj == 25 .AND. ji == 7 ) THEN
-                        WRITE(numout,*) ""
-                  !   WRITE(numout,*) "env_up = ", env1(ji,jj)
-                        WRITE(numout,*) "1st deriv. at TAU = 0._wp is ", d(2,1)
-                  !   WRITE(numout,*) "env_dw = ", env2(ji,jj)
-                        WRITE(numout,*) "1st deriv. at TAU = 1._wp is ", d(2,2)
-                  !   WRITE(numout,*) "CUBIC SPLINE at s = ", -sigma(1, 1, num_s) , " is ", &
-                  !                   z_cub_spl(-sigma(1, 1, num_s), tau, c)
-                        WRITE(numout,*) "1st deriv. at s = ", -sigma(1, 1, num_s) , " is ", &
-                                        D1z_cub_spl(-sigma(1, 1, num_s), tau, c, 1)
-                  !   WRITE(numout,*) "CUBIC SPLINE at s = ", -sigma(1, num_s, num_s) , " is ", &
-                  !                   z_cub_spl(-sigma(1, num_s, num_s), tau, c)
-                        WRITE(numout,*) "1st deriv. at s = ", -sigma(1, num_s, num_s) , " is ", &
-                                        D1z_cub_spl(-sigma(1, num_s, num_s), tau, c, 1)
-                     END IF
-                  END IF
- 
-               END DO
-            END DO
-
-         END IF
-
-      END DO
-
-      ! =============================================
-      ! 4. COMPUTE e3t, e3w for ALL general levels
-      !    as finite differences
-      ! =============================================
-      !
-      DO jk=1,jpkm1
-         e3t_0(:,:,jk)   = gdepw_0(:,:,jk+1) - gdepw_0(:,:,jk)
-         e3w_0(:,:,jk+1) = gdept_0(:,:,jk+1) - gdept_0(:,:,jk)
-      ENDDO
-      ! Surface
-      jk = 1
-      e3w_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,1) - gdepw_0(:,:,1))
-      !
-      ! Bottom
-      jk = jpk
-      e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk))
-
-      IF ( lwp ) THEN
-         jj = 25 
-         ji = 7
-         WRITE(numout,*) ""
-         WRITE(numout,*) "      k    gdepw_0(k)    gdept_0(k)    e3w_0(k)    e3t_0(k)"
-         WRITE(numout,*) ""
-         DO jk = 1, jpk
-            WRITE(numout,*) jk, gdepw_0(ji,jj,jk), gdept_0(ji,jj,jk), e3w_0(ji,jj,jk), e3t_0(ji,jj,jk)
-         ENDDO
-      END IF
-
-      !==============================
-      ! Computing gdept3w
-      !==============================
-
-      gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
-      DO jk = 2, jpk
-         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
-      END DO
-      e3t_0(:,:,jpk) = 2.0_wp * (gdept_0(:,:,jpk) - gdepw_0(:,:,jpk))
-
-      !==============================
-      ! Computing mbathy
-      !==============================
-
-      IF( lwp ) WRITE(numout,*) ''
-      IF( lwp ) WRITE(numout,*) 'MES mbathy:'
-      IF( lwp ) WRITE(numout,*) ''
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            DO jk = 1, jpkm1
-               !IF( lwp ) WRITE(numout,*) 'scobot: ', scobot(ji,jj), 'gdept_0: ', gdept_0(ji,jj,jk)
-               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
-            !IF( lwp ) WRITE(numout,*) 'mbathy(',ji,',',jj,') = ', mbathy(ji,jj)
-         END DO
-      END DO
-
-      !==============================================
-      ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0
-      !==============================================
-
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-            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))  )
-            ENDDO
-         ENDDO
-      ENDDO     
-
-      ! Adjusting for geopotential levels, if any
-
-      DO je = 1, tot_env ! LOOP over all the requested envelopes
-
-         IF (je == 1) THEN
-            s_1st  = 1
-            num_s  = nn_slev(je)
-            max_env_up = 0._wp      ! surface
-            min_env_up = 0._wp      ! surface
-         ELSE
-            s_1st  = s_1st + num_s - 1
-            num_s  = nn_slev(je) + 1
-            max_env_up = rn_ebot_max(je-1)
-            min_env_up = rn_ebot_min(je-1)
-         ENDIF
-
-         max_env_dw  = rn_ebot_max(je)
-         min_env_dw  = rn_ebot_min(je)
-
-         IF ( max_env_up == min_env_up .AND. max_env_dw == min_env_dw ) THEN
-
-            DO jj = 1,jpjm1
-               DO ji = 1,jpim1
-                  DO jk = 1, num_s
-
-                     ind = s_1st+jk-1
-                     e3u_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind))
-                     e3uw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji+1,jj,ind))
-                     e3v_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji,jj+1,ind))
-                     e3vw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji,jj+1,ind))
-                     e3f_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind))
-
-                  ENDDO
-               ENDDO
-            ENDDO
-
-         ENDIF
-
-      ENDDO
-
-      !==============================
-      ! Computing gdept3w
-      !==============================
-
-      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.
-
-      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
-            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 DO
-      END DO
-
-      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
-         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
-
-      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
-         WRITE(numout,*) 'zgr_sco : mbathy'
-         WRITE(numout,*) '~~~~~~~'
-         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( 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( 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
-
-      !================================================================================
-      ! Check coordinates makes sense
-      !================================================================================
-
-      DO ji = 1, jpi
-         DO jj = 1, jpj
-            DO jk = 1, mbathy(ji,jj)
-               ! check coordinate is monotonically increasing
-               IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
-                  WRITE(ctmp1,*) 'ERROR zgr_mes :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                  WRITE(numout,*) 'ERROR zgr_mes :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                  WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
-                  WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
-                  CALL ctl_stop( ctmp1 )
-               ENDIF
-               ! and check it has never gone negative
-               IF ( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
-                  WRITE(ctmp1,*) 'ERROR zgr_mes :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
-                  WRITE(numout,*) 'ERROR zgr_mes :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
-                  WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
-                  WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
-                  CALL ctl_stop( ctmp1 )
-               ENDIF
-!              ! and check it never exceeds the total depth
-!               IF ( fsdepw(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,:)
-!                  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
-!                  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,:)
-!                  CALL ctl_stop( ctmp1 )
-!               ENDIF
-!            END DO
-         END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi, jpj, max_nn_env , envlt )
-      CALL wrk_dealloc( jpi, jpj, env_up, env_dw )
-      CALL wrk_dealloc( jpk, z_gsigw1, z_gsigt1, z_esigw1, z_esigt1 )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_mes')
-
-   END SUBROUTINE zgr_mes_build
-
-! =====================================================================================================
-
-   SUBROUTINE zgr_dom_read(ld_isfcav,                             &   ! ice-cavities
-      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d, &   ! 1D reference vertical coordinate
-      &                    pdept , pdepw ,                        &   ! 3D t & w-points depth
-      &                    pe3t  , pe3u  , pe3v   , pe3f ,        &   ! vertical scale factors
-      &                    pe3w  , pe3uw , pe3vw                  )   !     -      -      -
+   SUBROUTINE zgr_dom_read
      !!---------------------------------------------------------------------
      !!              ***  ROUTINE zgr_dom_read  ***
@@ -1055,11 +89,4 @@
      !!
      !!----------------------------------------------------------------------
-     LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
-     REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
-     REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
-     REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m]
-     REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
-     REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
-     !
      INTEGER  ::   jk     ! dummy loop index
      INTEGER  ::   inum   ! local logical unit
@@ -1078,17 +105,19 @@
      !                          !* ocean cavities under iceshelves
      CALL iom_get( inum, 'ln_isfcav', z_cav )
-     IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF
+     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'  , pe3t_1d  )                     ! 1D reference coordinate
-     CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  )
-     !
-     CALL iom_get( inum, jpdom_global, 'e3t_0'  , pe3t  )    ! 3D coordinate
-     CALL iom_get( inum, jpdom_global, 'e3u_0'  , pe3u  )
-     CALL iom_get( inum, jpdom_global, 'e3v_0'  , pe3v  )
-     CALL iom_get( inum, jpdom_global, 'e3f_0'  , pe3f  )
-     CALL iom_get( inum, jpdom_global, 'e3w_0'  , pe3w  )
-     CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw )
-     CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw )
+     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
@@ -1100,17 +129,17 @@
         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', pdept_1d )   
-        CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
-        CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept)
-        CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw)
+        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( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth
-        CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths
+        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, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
+          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
@@ -1122,594 +151,161 @@
 ! =====================================================================================================
 
-   FUNCTION isodd( a ) RESULT(odd)
-      ! Notice that this method uses a btest intrinsic function. If you look at any number in binary 
-      ! notation, you'll note that the Least Significant Bit (LSB) will always be set for odd numbers 
-      ! and cleared for even numbers. Hence, all that is necessary to determine if a number is even or 
-      ! odd is to check the LSB. Since bitwise operations like AND, OR, XOR etc. are usually much faster 
-      ! than the arithmetic operations, this method will run a lot faster than the one with modulo 
-      ! operation.
-      ! http://forums.devshed.com/software-design-43/quick-algorithm-determine-odd-29843.html
-         
-      INTEGER, INTENT (in) :: a
-      LOGICAL              :: odd
-      odd = btest(a, 0)
-   END FUNCTION isodd
-
-! =====================================================================================================
- 
-   FUNCTION iseven( a ) RESULT(even)
-      INTEGER, INTENT (in) :: a
-      LOGICAL              :: even
-      even = .NOT. isodd(a)
-   END FUNCTION iseven
-
-! =====================================================================================================
-
-   FUNCTION sech( x ) RESULT( seh )
-
-      REAL(wp), INTENT(in   ) :: x
-      REAL(wp)                :: seh
-
-      seh = 1._wp / COSH(x)
-      
-   END FUNCTION sech
-
-! =====================================================================================================
-
-   FUNCTION sigma( vgrid, pk, kmax ) RESULT( ps1 )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE sigma  ***
-      !!       
-      !! ** Purpose :   provide the analytical function for sigma-coordinate
-      !!                (not stretched s-coordinate).
-      !!          
-      !! ** Method  :   the function provide the non-dimensional position of
-      !!                T and W points (i.e. between 0 and 1).
-      !!                T-points at integer values (between 1 and jpk)
-      !!                W-points at integer values - 1/2 (between 0.5 and
-      !!                jpk-0.5)
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT (in) ::   pk    ! continuous k coordinate
-      INTEGER, INTENT (in) ::   kmax  ! number of levels
-      INTEGER, INTENT (in) ::   vgrid ! type of vertical grid: 0 = T, 1 = W
-      REAL(wp)             ::   kindx ! index of T or W points
-      REAL(wp)             ::   ps1   ! value of sigma coordinate (-1 <= ps1 <=0)
-      !
-      SELECT CASE (vgrid)
-
-        CASE (0) ! T-points at integer values (between 1 and jpk)
-             kindx = REAL(pk,wp)
-        CASE (1) ! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
-             kindx = REAL(pk,wp) - 0.5_wp
-        CASE DEFAULT
-             WRITE(ctlmes,*) 'No valid vertical grid option selected'
-
-      END SELECT
-
-
-       ps1 = -(kindx - 0.5) / REAL(kmax-1) ! The surface is at the first W level
-                                           ! while the bottom coincides with the
-                                           ! last W level
-   END FUNCTION sigma
-
-! =====================================================================================================
-
-   FUNCTION s_coord( s, ca, cb, alpha, kmax, ftype ) RESULT ( pf1 )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE s_coord ***
-      !!
-      !! ** Purpose :   provide the analytical stretching function 
-      !!                for s-coordinate.
-      !!
-      !! ** Method  :   if ftype = 2: Siddorn and Furner 2012
-      !!                              analytical function is used
-      !!                if ftype = 1: Song and Haidvogel 1994 
-      !!                              analytical function is used
-      !!                if ftype = 0: Madec et al. 1996
-      !!                              analytical function is used
-      !!                              (pag 65 of NEMO Manual)
-      !!                              Reference  : Madec, Lott, Delecluse and
-      !!                              Crepon, 1996. JPO, 26, 1393-1408
-      !!
-      !!                s MUST be NEGATIVE: -1 <= s <= 0
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in) ::   s              ! continuous not stretched 
-                                               ! "s" coordinate
-      REAL(wp), INTENT(in) ::   ca             ! surface stretch. coeff
-      REAL(wp), INTENT(in) ::   cb             ! bottom stretch. coeff
-      REAL(wp), INTENT(in) ::   alpha          ! alpha stretch. coeff for SF12
-      INTEGER, INTENT (in) ::   kmax           ! number of levels
-      INTEGER,  INTENT(in) ::   ftype          ! type of stretching function
-      REAL(wp)             ::   pf1            ! "s" stretched value
-      ! SF12 stretch. funct. parameters
-      REAL(wp)             ::   psmth          ! smoothing parameter for transition
-                                               ! between shallow and deep areas
-      REAL(wp)             ::   za1,za2,za3    ! local variables
-      REAL(wp)             ::   zn1,zn2,ps     ! local variables
-      REAL(wp)             ::   za,zb,zx       ! local variables
-      !!----------------------------------------------------------------------
-      SELECT CASE (ftype)
-
-        CASE (0) ! M96  stretching function
-           pf1 =   (   TANH( ca * ( s + cb )  ) - TANH( cb * ca ) ) &
-            &    * (   COSH( ca )                                   &
-            &        + COSH( ca * ( 2.e0 * cb - 1.e0 ) )          ) &
-            &    / ( 2._wp * SINH( ca ) )
-        CASE (1) ! SH94 stretching function
-           IF ( ca == 0 ) THEN      ! uniform sigma
-              pf1 = s
-           ELSE                     ! stretched sigma
-              pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) + &
-            &       cb * ( ( TANH(ca*(s + 0.5_wp)) - TANH(0.5_wp*ca) ) / &
-            &       (2._wp*TANH(0.5_wp*ca)) )
+   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)
+
+     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
+     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))
+              !
+              ! Bottom
+              jk = jpk
+              e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk))
            END IF
-        CASE (2) ! SF12 stretching function
-           psmth = 1.0_wp ! We consider only the case for efold = 0
-           ps  = -s
-           zn1 =  1. / REAL(kmax-1)
-           zn2 =  1. -  zn1
-
-           za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp)
-           za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp)
-           za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
-
-           za  = cb - za3*(ca-za1)-za2
-           za  = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) )
-           zb  = (ca - 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
-
-           pf1 = za*(ps*(1.0_wp-ps/2.0_wp))+zb*ps**3.0_wp + &
-            &    zx*( (alpha+2.0_wp)*ps**(alpha+1.0_wp) - &
-            &         (alpha+1.0_wp)*ps**(alpha+2.0_wp) )
-           pf1 = -pf1*psmth+ps*(1.0_wp-psmth)
-        CASE (3) ! stretching function for deeper part of the domain
-           IF ( ca == 0 ) THEN      ! uniform sigma
-              pf1 = s
-           ELSE                     ! stretched sigma
-              pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) 
+        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 SELECT
-       
-   END FUNCTION s_coord     
-
-! =====================================================================================================
-
-   FUNCTION D1s_coord( s, ca, cb, alpha, kmax, ftype) RESULT( pf1 )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE D1s_coord ***
-      !!
-      !! ** Purpose :   provide the 1st derivative of the analytical 
-      !!                stretching function for s-coordinate.
-      !!
-      !! ** Method  :   if ftype = 2: Siddorn and Furner 2012
-      !!                              analytical function is used
-      !!                if ftype = 1: Song and Haidvogel 1994 
-      !!                              analytical function is used
-      !!                if ftype = 0: Madec et al. 1996
-      !!                              analytical function is used
-      !!                              (pag 65 of NEMO Manual)
-      !!                              Reference  : Madec, Lott, Delecluse and
-      !!                              Crepon, 1996. JPO, 26, 1393-1408
-      !!
-      !!                s MUST be NEGATIVE: -1 <= s <= 0
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) ::   s           ! continuous not stretched 
-                                               ! "s" coordinate
-      REAL(wp), INTENT(in)    ::   ca          ! surface stretch. coeff
-      REAL(wp), INTENT(in)    ::   cb          ! bottom stretch. coeff
-      REAL(wp), INTENT(in)    ::   alpha       ! alpha stretch. coeff for SF12
-      INTEGER, INTENT (in)    ::   kmax        ! number of levels
-      INTEGER,  INTENT(in   ) ::   ftype       ! type of stretching function
-      REAL(wp)                ::   pf1         ! "s" stretched value
-      ! SF12 stretch. funct. parameters
-      REAL(wp)                ::   psmth       ! smoothing parameter for transition
-                                               ! between shallow and deep areas
-      REAL(wp)                ::   za1,za2,za3 ! local variables
-      REAL(wp)                ::   zn1,zn2,ps  ! local variables
-      REAL(wp)                ::   za,zb,zx    ! local variables
-      REAL(wp)                ::   zt1,zt2,zt3 ! local variables
-      !!----------------------------------------------------------------------
-      SELECT CASE (ftype)
-
-        CASE (0) ! M96  stretching function
-           pf1 =  ( ca * ( COSH(ca*(2._wp * cb - 1._wp)) + COSH(ca)) * &
-                     sech(ca * (s+cb))**2 ) / (2._wp * SINH(ca))
-        CASE (1) ! SH94 stretching function
-           IF ( ca == 0 ) then      ! uniform sigma
-              pf1 = 1._wp / REAL(kmax,wp)
-           ELSE                     ! stretched sigma
-              pf1 = (1._wp - cb) * ca * COSH(ca*s) / (SINH(ca) * REAL(kmax,wp)) + &
-                    cb * ca * &
-                   (sech((s+0.5_wp)*ca)**2) / (2._wp * TANH(0.5_wp*ca) * REAL(kmax,wp))
+        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 ( 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))  )
+              END DO
            END IF
-        CASE (2) ! SF12 stretching function
-           ps  = -s
-           zn1 =  1. / REAL(kmax-1)
-           zn2 =  1. -  zn1
-
-           za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp)
-           za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp)
-           za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
-
-           za  = cb - za3*(ca-za1)-za2
-           za  = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) )
-           zb  = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
-           zx  = (alpha+2.0_wp)*(alpha+1.0_wp)
-
-           zt1 = 0.5_wp*za*(ps-1._wp)*((ps**2.0_wp + 3._wp*ps + 2._wp)*ps**alpha - 2._wp)
-           zt2 = zb*(zx*ps**(alpha+1.0_wp) - zx*ps**(alpha) + 3._wp*ps**2._wp)
-           zt3 = zx*(ps-1._wp)*ps**(alpha)
-
-           pf1 = zt1 + zt2 - zt3
-
-      END SELECT
-
-   END FUNCTION D1s_coord
-
-! =====================================================================================================
-
-  FUNCTION z_mes(dep_up, dep_dw, Cs, s, hc) RESULT( z )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE z_mes ***
-      !!
-      !! ** Purpose :   provide the analytical trasformation from the 
-      !!                computational space (mes-coordinate) to the  
-      !!                physical space (depth z)
-      !!
-      !!    N.B.: z is downward positive defined as well as envelope surfaces.
-      !!          Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1.
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) :: dep_up  ! shallower envelope
-      REAL(wp), INTENT(in   ) :: dep_dw  ! deeper envelope
-      REAL(wp), INTENT(in   ) :: Cs      ! stretched s coordinate   
-      REAL(wp), INTENT(in   ) :: s       ! not stretched s coordinate
-      REAL(wp), INTENT(in   ) :: hc      ! critical depth
-      REAL                    :: z       ! downward positive depth
-      !!----------------------------------------------------------------------
-      ! 
-
-      z = dep_up + hc * s + Cs * (dep_dw - hc - dep_up)
-
-  END FUNCTION z_mes
-
-! =====================================================================================================
-
-  FUNCTION D1z_mes(dep_up, dep_dw, d1Cs, hc, kmax) RESULT( d1z )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE D1z_mes ***
-      !!
-      !! ** Purpose :   provide the 1st derivative of the analytical 
-      !!                trasformation from the computational space 
-      !!                (mes-coordinate) to the physical space (depth z)
-      !!
-      !!    N.B.: z is downward positive defined as well as envelope surfaces.
-      !!          Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1.
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) :: dep_up  ! shallower envelope
-      REAL(wp), INTENT(in   ) :: dep_dw  ! deeper envelope
-      REAL(wp), INTENT(in   ) :: d1Cs    ! 1st derivative of the 
-                                         ! stretched s coordinate   
-      REAL(wp), INTENT(in   ) :: hc      ! critical depth
-      INTEGER,  INTENT(in   ) :: kmax    ! number of levels
-      REAL(wp)                :: d1z     ! 1st derivative of downward
-                                                 ! positive depth
-      !!----------------------------------------------------------------------
-      ! 
-
-      d1z = hc / REAL(kmax,wp) + d1Cs * (dep_dw - hc - dep_up)
-      !d1z = hc + d1Cs * (dep_dw - hc - dep_up)
-
-  END FUNCTION D1z_mes
-
-! =====================================================================================================
-
-  FUNCTION cub_spl(tau, d, n, ibcbeg, ibcend) RESULT ( c )
-      !!----------------------------------------------------------------------
-      !!
-      !! CUBSPL defines an interpolatory cubic spline.
-      !!
-      !! Discussion:
-      !!
-      !!    A tridiagonal linear system for the unknown slopes S(I) of
-      !!    F at TAU(I), I=1,..., N, is generated and then solved by Gauss
-      !!    elimination, with S(I) ending up in C(2,I), for all I.
-      !!
-      !! Author: Carl de Boor
-      !!
-      !! Reference: Carl de Boor, Practical Guide to Splines,
-      !!             Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27.
-      !!
-      !! Parameters:
-      !!
-      !!    Input: 
-      !!            TAU(N) : the abscissas or X values of the data points.  
-      !!                     The entries of TAU are assumed to be strictly 
-      !!                     increasing.
-      !!
-      !!            N      : the number of data points.  N is assumed to be 
-      !!                     at least 2.
-      !! 
-      !!            IBCBEG : boundary condition indicator at TAU(1).
-      !!                     = 0 no boundary condition at TAU(1) is given.
-      !!                         In this case, the "not-a-knot condition" 
-      !!                         is used.
-      !!                     = 1 the 1st derivative at TAU(1) is equal to 
-      !!                         the input value D(2,1).
-      !!                     = 2 the 2nd derivative at TAU(1) is equal to 
-      !!                         the input value D(2,1).
-      !!
-      !!            IBCEND : boundary condition indicator at TAU(N).
-      !!                     = 0 no boundary condition at TAU(N) is given.
-      !!                         In this case, the "not-a-knot condition" 
-      !!                         is used.
-      !!                     = 1 the 1st derivative at TAU(N) is equal to 
-      !!                         the input value D(2,2).
-      !!                     = 2 the 2nd derivative at TAU(N) is equal to 
-      !!                         the input value D(2,2).                                 
-      !!
-      !!            D(1,1): value of the function at TAU(1)
-      !!            D(1,1): value of the function at TAU(N)
-      !!            D(2,1): if IBCBEG is 1 (2) it is the value of the
-      !!                    1st (2nd) derivative at TAU(1)
-      !!            D(2,2): if IBCBEG is 1 (2) it is the value of the
-      !!                    1st (2nd) derivative at TAU(N) 
-      !!
-      !!    Output: 
-      !!            C(4,N) : contains the polynomial coefficients of the 
-      !!                     cubic interpolating spline.
-      !!                     In the interval interval (TAU(I), TAU(I+1)), 
-      !!                     the spline F is given by
-      !!
-      !!                          F(X) = 
-      !!                                  C(1,I) + 
-      !!                                  C(2,I) * H +
-      !!                                  C(3,I) * H^2 / 2 + 
-      !!                                  C(4,I) * H^3 / 6.
-      !!
-      !!                     where H = X - TAU(I).  
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER,  INTENT(in   ) :: n
-      REAL(wp), INTENT(in   ) :: tau(n)
-      REAL(wp), INTENT(in   ) :: d(2,2)
-      INTEGER,  INTENT(in   ) :: ibcbeg, ibcend
-      REAL(wp)                :: c(4,n)
-      REAL(wp)                :: divdf1
-      REAL(wp)                :: divdf3
-      REAL(wp)                :: dtau
-      REAL(wp)                :: g
-      INTEGER(wp)             :: i
-      !!----------------------------------------------------------------------
-      !  Initialise c and copy d values to c
-      c(:,:) = 0.0D+00
-      c(1,1) = d(1,1)
-      c(1,n) = d(1,2)
-      c(2,1) = d(2,1)
-      c(2,n) = d(2,2)
-      ! 
-      !  C(3,*) and C(4,*) are used initially for temporary storage.
-      !  Store first differences of the TAU sequence in C(3,*).
-      !  Store first divided difference of data in C(4,*).
-      DO i = 2, n
-         c(3,i) = tau(i) - tau(i-1)
-         c(4,i) = ( c(1,i) - c(1,i-1) ) / c(3,i)
-      END DO
-
-      !  Construct the first equation from the boundary condition
-      !  at the left endpoint, of the form:
-      !
-      !    C(4,1) * S(1) + C(3,1) * S(2) = C(2,1)
-      !
-      !  IBCBEG = 0: Not-a-knot
-      IF ( ibcbeg == 0 ) THEN
-
-         IF ( n <= 2 ) THEN
-            c(4,1) = 1._wp
-            c(3,1) = 1._wp
-            c(2,1) = 2._wp * c(4,2)
-         ELSE
-            c(4,1) = c(3,3)
-            c(3,1) = c(3,2) + c(3,3)
-            c(2,1) = ( ( c(3,2) + 2._wp * c(3,1) ) * c(4,2) &
-                     * c(3,3) + c(3,2)**2 * c(4,3) ) / c(3,1)
-         END IF
-
-      !  IBCBEG = 1: derivative specified.
-      ELSE IF ( ibcbeg == 1 ) then
-
-         c(4,1) = 1._wp
-         c(3,1) = 0._wp
-
-      !  IBCBEG = 2: Second derivative prescribed at left end.
-      ELSE IF ( ibcbeg == 2 ) then
-
-         c(4,1) = 2._wp
-         c(3,1) = 1._wp
-         c(2,1) = 3._wp * c(4,2) - c(3,2) / 2._wp * c(2,1)
-
-      ELSE
-         WRITE(ctlmes,*) 'CUBSPL - Error, invalid IBCBEG input option!'
-         CALL ctl_stop( ctlmes )
-      END IF
-  
-      !  If there are interior knots, generate the corresponding
-      !  equations and carry out the forward pass of Gauss,
-      !  elimination after which the I-th equation reads:
-      !
-      !    C(4,I) * S(I) + C(3,I) * S(I+1) = C(2,I).
-
-      IF ( n > 2 ) THEN
-
-         DO i = 2, n-1
-            g = -c(3,i+1) / c(4,i-1)
-            c(2,i) = g * c(2,i-1) + 3._wp * &
-                     ( c(3,i) * c(4,i+1) + c(3,i+1) * c(4,i) )
-            c(4,i) = g * c(3,i-1) + 2._wp * ( c(3,i) + c(3,i+1))
-         END DO
-
-         !  Construct the last equation from the second boundary,
-         !  condition of the form
-         !
-         !    -G * C(4,N-1) * S(N-1) + C(4,N) * S(N) = C(2,N)
-         !
-         !  If 1st der. is prescribed at right end ( ibcend == 1 ), 
-         !  one can go directly to back-substitution, since the C
-         !  array happens to be set up just right for it at this 
-         !  point.
-
-         IF ( ibcend < 1 ) THEN
-            !  Not-a-knot and 3 <= N, and either 3 < N or also 
-            !  not-a-knot at left end point.
-            IF ( n /= 3 .OR. ibcbeg /= 0 ) THEN
-               g      = c(3,n-1) + c(3,n)
-               c(2,n) = ( ( c(3,n) + 2._wp * g ) * c(4,n) * c(3,n-1) &
-                        + c(3,n)**2 * ( c(1,n-1) - c(1,n-2) ) / &
-                        c(3,n-1) ) / g
-               g      = - g / c(4,n-1)
-               c(4,n) = c(3,n-1)
-               c(4,n) = c(4,n) + g * c(3,n-1)
-               c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
-            ELSE
-            !  N = 3 and not-a-knot also at left.
-               c(2,n) = 2._wp * c(4,n)
-               c(4,n) = 1._wp
-               g      = -1._wp / c(4,n-1)
-               c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
-               c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
-            END IF
-
-         ELSE IF ( ibcend == 2 ) THEN
-         !  IBCEND = 2: Second derivative prescribed at right endpoint.
-               c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n)
-               c(4,n) = 2._wp
-               g      = -1._wp / c(4,n-1)
-               c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
-               c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
-         END IF
-
-      ELSE
-         !  N = 2 (assumed to be at least equal to 2!).
-
-         IF ( ibcend == 2  ) THEN
-
-            c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n)
-            c(4,n) = 2._wp
-            g      = -1._wp / c(4,n-1)
-            c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
-            c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
-
-         ELSE IF ( ibcend == 0 .AND. ibcbeg /= 0 ) THEN
-
-            c(2,n) = 2._wp * c(4,n)
-            c(4,n) = 1._wp
-            g      = -1._wp / c(4,n-1)
-            c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
-            c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
-
-         ELSE IF ( ibcend == 0 .AND. ibcbeg == 0 ) THEN
-
-            c(2,n) = c(4,n)
-
-         END IF
-
-      END IF
-  
-      !  Back solve the upper triangular system 
-      !
-      !    C(4,I) * S(I) + C(3,I) * S(I+1) = B(I)
-      !
-      !  for the slopes C(2,I), given that S(N) is already known.
-      ! 
-      DO i = n-1, 1, -1
-         c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i)
-      END DO
-
-      !  Generate cubic coefficients in each interval, that is, the
-      !  derivatives at its left endpoint, from value and slope at its
-      !  endpoints.
-      DO i = 2, n
-         dtau     = c(3,i)
-         divdf1   = ( c(1,i) - c(1,i-1) ) / dtau
-         divdf3   = c(2,i-1) + c(2,i) - 2._wp * divdf1
-         c(3,i-1) = 2._wp * ( divdf1 - c(2,i-1) - divdf3 ) / dtau
-         c(4,i-1) = 6._wp * divdf3 / dtau**2
-      END DO
-
-  END FUNCTION cub_spl
-
-! =====================================================================================================
-
-  FUNCTION z_cub_spl(s, tau, c) RESULT ( z_cs )
-      !----------------------------------------------------------------------
-      !                 ***  function z_cub_spl ***
-      !
-      ! ** Purpose : evaluate the cubic spline F in the interval 
-      !              (TAU(I), TAU(I+1)). F is given by 
-      !              
-      !
-      !           F(X) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6
-      !
-      !              where H = S-TAU(I).
-      !
-      !              s MUST be positive: 0 <= s <= 1
-      !---------------------------------------------------------------------
-      INTEGER, PARAMETER      :: n = 2
-      REAL(wp), INTENT(in   ) :: s
-      REAL(wp), INTENT(in   ) :: tau(n)
-      REAL(wp), INTENT(in   ) :: c(4,n)
-      REAL(wp)                :: z_cs
-      REAL(wp)                :: c1, c2, c3, c4
-      REAL(wp)                :: H
-      !---------------------------------------------------------------------
-      c1 = c(1,1)
-      c2 = c(2,1)
-      c3 = c(3,1)
-      c4 = c(4,1)
-
-      H    = s - tau(1)
-
-      z_cs = c1 + c2 * H + (c3 * H**2._wp)/2._wp + (c4 * H**3._wp)/6._wp
-
-  END FUNCTION z_cub_spl  
-
-! =====================================================================================================
-
-  FUNCTION D1z_cub_spl(s, tau, c, kmax) RESULT ( d1z_cs )
-      !----------------------------------------------------------------------
-      !                 ***  function D1z_cub_spl ***
-      !
-      ! ** Purpose : evaluate the 1st derivative of cubic spline F  
-      !              in the interval (TAU(I), TAU(I+1)). 
-      !              The 1st derivative D1F is given by 
-      !              
-      !
-      !           D1F(S) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6
-      !
-      !              where H = S-TAU(I).
-      !
-      !              s MUST be positive: 0 <= s <= 1
-      !---------------------------------------------------------------------
-      INTEGER,  PARAMETER     :: n = 2
-      REAL(wp), INTENT(in   ) :: s
-      REAL(wp), INTENT(in   ) :: tau(n)
-      REAL(wp), INTENT(in   ) :: c(4,n)
-      INTEGER,  INTENT(in   ) :: kmax 
-      REAL(wp)                :: d1z_cs
-      REAL(wp)                :: c2, c3, c4
-      REAL(wp)                :: H
-      !---------------------------------------------------------------------
-      c2 = c(2,1) / REAL(kmax,wp)
-      c3 = c(3,1) / REAL(kmax,wp)
-      c4 = c(4,1) / REAL(kmax,wp)
-
-      H    = s - tau(1)
-
-      d1z_cs = c2 + c3 * H + (c4 * H**2._wp)/2._wp
-
-  END FUNCTION D1z_cub_spl
-
-! =====================================================================================================
-
+        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
 
 END MODULE zgrmes
