Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 7351)
@@ -2,5 +2,5 @@
    !!==============================================================================
    !!                       ***  MODULE domzgr   ***
-   !! Ocean initialization : domain initialization
+   !! Ocean domain : definition of the vertical coordinate system
    !!==============================================================================
    !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
@@ -18,4 +18,5 @@
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case  
    !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye  
+   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying
    !!----------------------------------------------------------------------
 
@@ -36,6 +37,8 @@
    USE oce               ! ocean variables
    USE dom_oce           ! ocean domain
+   USE wet_dry           ! wetting and drying
    USE closea            ! closed seas
    USE c1d               ! 1D vertical configuration
+   !
    USE in_out_manager    ! I/O manager
    USE iom               ! I/O library
@@ -73,5 +76,4 @@
 
   !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -102,5 +104,5 @@
       INTEGER ::   ios
       !
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
       !!----------------------------------------------------------------------
       !
@@ -120,10 +122,13 @@
          WRITE(numout,*) 'dom_zgr : vertical coordinate'
          WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
-         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
-         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
-         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
-         WRITE(numout,*) '             ice shelf cavities             ln_isfcav = ', ln_isfcav
+         WRITE(numout,*) '   Namelist namzgr : set vertical coordinate'
+         WRITE(numout,*) '      z-coordinate - full steps      ln_zco    = ', ln_zco
+         WRITE(numout,*) '      z-coordinate - partial steps   ln_zps    = ', ln_zps
+         WRITE(numout,*) '      s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
+         WRITE(numout,*) '      ice shelf cavities             ln_isfcav = ', ln_isfcav
+         WRITE(numout,*) '      linear free surface            ln_linssh = ', ln_linssh
       ENDIF
+
+      IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time'
 
       ioptio = 0                       ! Check Vertical coordinate options
@@ -132,4 +137,9 @@
       IF( ln_sco      )   ioptio = ioptio + 1
       IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
+      !
+      ioptio = 0
+      IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
+      IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
+      IF( ioptio > 0 )   CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
       !
       ! Build the vertical coordinate system
@@ -155,18 +165,18 @@
       !
       IF( nprint == 1 .AND. lwp )   THEN
-         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+         WRITE(numout,*) ' MIN val mbathy  ', MINVAL(  mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
          WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
-            &                   ' w ',   MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ),  &
-            &                   ' u ',   MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ),  &
-            &                   ' uw',   MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)),   &
-            &                   ' w ',   MINVAL( e3w_0(:,:,:) )
+            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  &
+            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  &
+            &                          ' uw', MINVAL(  e3uw_0(:,:,:) ), ' vw', MINVAL(  e3vw_0(:,:,:)),   &
+            &                          ' w ', MINVAL(   e3w_0(:,:,:) )
 
          WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
-            &                   ' w ',   MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ),  &
-            &                   ' u ',   MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ),  &
-            &                   ' uw',   MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),   &
-            &                   ' w ',   MAXVAL( e3w_0(:,:,:) )
+            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  &
+            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  &
+            &                          ' uw', MAXVAL(  e3uw_0(:,:,:) ), ' vw', MAXVAL(  e3vw_0(:,:,:) ),  &
+            &                          ' w ', MAXVAL(   e3w_0(:,:,:) )
       ENDIF
       !
@@ -379,5 +389,5 @@
       !!              - bathy : meter bathymetry (in meters)
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
       INTEGER  ::   inum                      ! temporary logical unit
       INTEGER  ::   ierror                    ! error flag
@@ -498,4 +508,6 @@
             CALL iom_close( inum )
             mbathy(:,:) = INT( bathy(:,:) )
+            ! initialisation isf variables
+            risfdep(:,:)=0._wp ; misfdep(:,:)=1             
             !                                                ! =====================
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
@@ -534,6 +546,7 @@
             CALL iom_close( inum )
             !                                                
-            risfdep(:,:)=0._wp         
-            misfdep(:,:)=1             
+            ! initialisation isf variables
+            risfdep(:,:)=0._wp ; misfdep(:,:)=1             
+            !
             IF ( ln_isfcav ) THEN
                CALL iom_open ( 'isf_draft_meter.nc', inum ) 
@@ -541,4 +554,11 @@
                CALL iom_close( inum )
                WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
+
+               ! set grounded point to 0 
+               ! (a treshold could be set here if needed, or set it offline based on the grounded fraction)
+               WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin )
+                  misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
+                  mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
+               END WHERE
             END IF
             !       
@@ -578,11 +598,4 @@
       !                       
       IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
-         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
-         IF ( ln_isfcav ) THEN
-            WHERE (bathy == risfdep)
-               bathy   = 0.0_wp ; risfdep = 0.0_wp
-            END WHERE
-         END IF
-         ! end patch
          IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
          ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
@@ -674,9 +687,7 @@
       !!              - update bathy : meter bathymetry (in meters)
       !!----------------------------------------------------------------------
-      !!
       INTEGER ::   ji, jj, jl                    ! dummy loop indices
       INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
       REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy
-
       !!----------------------------------------------------------------------
       !
@@ -775,12 +786,4 @@
          IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
       ENDIF
-
-      IF( lwp .AND. nprint == 1 ) THEN      ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
-         WRITE(numout,*) ' ------------------'
-         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
-         WRITE(numout,*)
-      ENDIF
       !
       CALL wrk_dealloc( jpi, jpj, zbathy )
@@ -803,5 +806,4 @@
       !!                                     (min value = 1 over land)
       !!----------------------------------------------------------------------
-      !!
       INTEGER ::   ji, jj   ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk
@@ -835,7 +837,8 @@
    END SUBROUTINE zgr_bot_level
 
-      SUBROUTINE zgr_top_level
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_bot_level  ***
+
+   SUBROUTINE zgr_top_level
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_top_level  ***
       !!
       !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
@@ -847,5 +850,4 @@
       !!                                     (min value = 1)
       !!----------------------------------------------------------------------
-      !!
       INTEGER ::   ji, jj   ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:) ::  zmik
@@ -881,9 +883,10 @@
    END SUBROUTINE zgr_top_level
 
+
    SUBROUTINE zgr_zco
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zgr_zco  ***
       !!
-      !! ** Purpose :   define the z-coordinate system
+      !! ** Purpose :   define the reference z-coordinate system
       !!
       !! ** Method  :   set 3D coord. arrays to reference 1D array 
@@ -895,14 +898,14 @@
       !
       DO jk = 1, jpk
-         gdept_0 (:,:,jk) = gdept_1d(jk)
-         gdepw_0 (:,:,jk) = gdepw_1d(jk)
-         gdep3w_0(:,:,jk) = gdepw_1d(jk)
-         e3t_0   (:,:,jk) = e3t_1d  (jk)
-         e3u_0   (:,:,jk) = e3t_1d  (jk)
-         e3v_0   (:,:,jk) = e3t_1d  (jk)
-         e3f_0   (:,:,jk) = e3t_1d  (jk)
-         e3w_0   (:,:,jk) = e3w_1d  (jk)
-         e3uw_0  (:,:,jk) = e3w_1d  (jk)
-         e3vw_0  (:,:,jk) = e3w_1d  (jk)
+         gdept_0(:,:,jk) = gdept_1d(jk)
+         gdepw_0(:,:,jk) = gdepw_1d(jk)
+         gde3w_0(:,:,jk) = gdepw_1d(jk)
+         e3t_0  (:,:,jk) = e3t_1d  (jk)
+         e3u_0  (:,:,jk) = e3t_1d  (jk)
+         e3v_0  (:,:,jk) = e3t_1d  (jk)
+         e3f_0  (:,:,jk) = e3t_1d  (jk)
+         e3w_0  (:,:,jk) = e3w_1d  (jk)
+         e3uw_0 (:,:,jk) = e3w_1d  (jk)
+         e3vw_0 (:,:,jk) = e3w_1d  (jk)
       END DO
       !
@@ -917,5 +920,5 @@
       !!                     
       !! ** Purpose :   the depth and vertical scale factor in partial step
-      !!      z-coordinate case
+      !!              reference z-coordinate case
       !!
       !! ** Method  :   Partial steps : computes the 3D vertical scale factors
@@ -957,13 +960,10 @@
       !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
       !!----------------------------------------------------------------------
-      !!
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       INTEGER  ::   ik, it, ikb, ikt ! temporary integers
-      LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
-      REAL(wp) ::   zmax             ! Maximum depth
       REAL(wp) ::   zdiff            ! temporary scalar
-      REAL(wp) ::   zrefdep          ! temporary scalar
+      REAL(wp) ::   zmax             ! temporary scalar
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
       !!---------------------------------------------------------------------
@@ -971,5 +971,5 @@
       IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zprt )
+      CALL wrk_alloc( jpi,jpj,jpk,   zprt )
       !
       IF(lwp) WRITE(numout,*)
@@ -977,13 +977,4 @@
       IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
-
-      ll_print = .FALSE.                   ! Local variable for debugging
-      
-      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
-         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
-      ENDIF
-
 
       ! bathymetry in level (from bathy_meter)
@@ -1004,6 +995,4 @@
       END DO
 
-      IF ( ln_isfcav ) CALL zgr_isf
-
       ! Scale factors and depth at T- and W-points
       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
@@ -1013,5 +1002,736 @@
          e3w_0  (:,:,jk) = e3w_1d  (jk)
       END DO
+      
+      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
+      IF ( ln_isfcav ) CALL zgr_isf
+
+      ! Scale factors and depth at T- and W-points
+      IF ( .NOT. ln_isfcav ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ik = mbathy(ji,jj)
+               IF( ik > 0 ) THEN               ! ocean point only
+                  ! max ocean level case
+                  IF( ik == jpkm1 ) THEN
+                     zdepwp = bathy(ji,jj)
+                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
+                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
+                     e3t_0(ji,jj,ik  ) = ze3tp
+                     e3t_0(ji,jj,ik+1) = ze3tp
+                     e3w_0(ji,jj,ik  ) = ze3wp
+                     e3w_0(ji,jj,ik+1) = ze3tp
+                     gdepw_0(ji,jj,ik+1) = zdepwp
+                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
+                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
+                     !
+                  ELSE                         ! standard case
+                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
+                     ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+                     ENDIF
+   !gm Bug?  check the gdepw_1d
+                     !       ... on ik
+                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
+                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
+                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
+                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
+                     e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
+                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+                     !       ... on ik+1
+                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+         !
+         it = 0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ik = mbathy(ji,jj)
+               IF( ik > 0 ) THEN               ! ocean point only
+                  e3tp (ji,jj) = e3t_0(ji,jj,ik)
+                  e3wp (ji,jj) = e3w_0(ji,jj,ik)
+                  ! test
+                  zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
+                  IF( zdiff <= 0._wp .AND. lwp ) THEN 
+                     it = it + 1
+                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
+                     WRITE(numout,*) ' bathy = ', bathy(ji,jj)
+                     WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+                     WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      END IF
+      !
+      ! Scale factors and depth at U-, V-, UW and VW-points
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3u_0 (:,:,jk) = e3t_1d(jk)
+         e3v_0 (:,:,jk) = e3t_1d(jk)
+         e3uw_0(:,:,jk) = e3w_1d(jk)
+         e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+
+      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
+               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
+               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
+               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      IF ( ln_isfcav ) THEN
+      ! (ISF) define e3uw (adapted for 2 cells in the water column)
+         DO jj = 2, jpjm1 
+            DO ji = 2, fs_jpim1   ! vector opt. 
+               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
+               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
+               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
+               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
+            END DO
+         END DO
+      END IF
+
+      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
+      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
+      !
+
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
+         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+      
+      ! Scale factor at F-point
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
+      !
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+!!gm  bug ? :  must be a do loop with mj0,mj1
       ! 
+      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_zps :   e r r o r   e3t_0 <= 0' )
+      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
+      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
+      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
+     
+      ! Compute gde3w_0 (vertical sum of e3w)
+      IF ( ln_isfcav ) THEN ! if cavity
+         WHERE( misfdep == 0 )   misfdep = 1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
+               DO jk = 2, misfdep(ji,jj)
+                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+               END DO
+               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
+               DO jk = misfdep(ji,jj) + 1, jpk
+                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+      ELSE ! no cavity
+         gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
+         DO jk = 2, jpk
+            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
+         END DO
+      END IF
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   zprt )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
+      !
+   END SUBROUTINE zgr_zps
+
+
+   SUBROUTINE zgr_isf
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_isf  ***
+      !!   
+      !! ** Purpose :   check the bathymetry in levels
+      !!   
+      !! ** Method  :   THe water column have to contained at least 2 cells
+      !!                Bathymetry and isfdraft are modified (dig/close) to respect
+      !!                this criterion.
+      !!   
+      !! ** Action  : - test compatibility between isfdraft and bathy 
+      !!              - bathy and isfdraft are modified
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices
+      INTEGER  ::   ik, it               ! temporary integers
+      INTEGER  ::   icompt, ibtest       ! (ISF)
+      INTEGER  ::   ibtestim1, ibtestip1 ! (ISF)
+      INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF)
+      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t
+      REAL(wp) ::   zmax             ! Maximum and minimum depth
+      REAL(wp) ::   zbathydiff       ! isf temporary scalar
+      REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar
+      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
+      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t
+      REAL(wp) ::   zdiff            ! temporary scalar
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
+      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('zgr_isf')
+      !
+      CALL wrk_alloc( jpi,jpj,   zbathy, zmask, zrisfdep)
+      CALL wrk_alloc( jpi,jpj,   zmisfdep, zmbathy )
+
+      ! (ISF) compute misfdep
+      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
+      ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
+      END WHERE  
+
+      ! Compute misfdep for ocean points (i.e. first wet level) 
+      ! find the first ocean level such that the first level thickness 
+      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 
+      ! e3t_0 is the reference level thickness 
+      DO jk = 2, jpkm1 
+         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
+         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
+      END DO 
+      WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) )
+         risfdep(:,:) = 0. ; misfdep(:,:) = 1
+      END WHERE
+
+      ! remove very shallow ice shelf (less than ~ 10m if 75L)
+      WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1)
+         misfdep = 0; risfdep = 0.0_wp;
+         mbathy  = 0; bathy   = 0.0_wp;
+      END WHERE
+      WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp)
+         misfdep = 0; risfdep = 0.0_wp;
+         mbathy  = 0; bathy   = 0.0_wp;
+      END WHERE
+ 
+! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
+      icompt = 0 
+! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
+      DO jl = 1, 10     
+         ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)
+         WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin)
+            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
+            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
+         END WHERE
+         WHERE (mbathy(:,:) <= 0) 
+            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
+            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
+         END WHERE
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
+            misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west 
+            misfdep(jpi,:) = misfdep(  2  ,:) 
+            mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west
+            mbathy(jpi,:)  = mbathy(  2  ,:)
+         ENDIF
+
+         ! split last cell if possible (only where water column is 2 cell or less)
+         ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
+         IF ( .NOT. ln_iscpl) THEN
+            DO jk = jpkm1, 1, -1
+               zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+               WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
+                  mbathy(:,:) = jk
+                  bathy(:,:)  = zmax
+               END WHERE
+            END DO
+         END IF
+ 
+         ! split top cell if possible (only where water column is 2 cell or less)
+         DO jk = 2, jpkm1
+            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
+               misfdep(:,:) = jk
+               risfdep(:,:) = zmax
+            END WHERE
+         END DO
+
+ 
+ ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! find the minimum change option:
+               ! test bathy
+               IF (risfdep(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN
+                     zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
+                         &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
+                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
+                        IF (zbathydiff <= zrisfdepdiff) THEN
+                           bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
+                           mbathy(ji,jj)= mbathy(ji,jj) + 1
+                        ELSE
+                           risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
+                           misfdep(ji,jj) = misfdep(ji,jj) - 1
+                        END IF
+                     ENDIF
+                  ELSE
+                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
+                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
+                        misfdep(ji,jj) = misfdep(ji,jj) - 1
+                     END IF
+                  END IF
+               END IF
+            END DO
+         END DO
+ 
+         ! At least 2 levels for water thickness at T, U, and V point.
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! find the minimum change option:
+               ! test bathy
+               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) &
+                         &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) & 
+                         &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy(ji,jj) = mbathy(ji,jj) + 1
+                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)= misfdep(ji,jj) - 1
+                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj)= misfdep(ji,jj) - 1
+                     risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+ 
+ ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) &
+                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
+                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy(ji,jj) = mbathy(ji,jj) + 1
+                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
+                        risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
+                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+ 
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+ ! point V misdep(ji,jj) == mbathy(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
+                           &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) &
+                           &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
+                        bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)   = misfdep(ji,jj) - 1
+                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
+                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+ 
+ 
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF 
+ 
+ ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                  zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
+                       &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
+                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
+                       &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
+                  IF (zbathydiff <= zrisfdepdiff) THEN
+                     mbathy(ji,jj) = mbathy(ji,jj) + 1
+                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+                  ELSE
+                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
+                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
+                  END IF
+                  ELSE
+                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
+                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+ 
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF 
+ 
+ ! point U misfdep(ji,jj) == bathy(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
+                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) &
+                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
+                        bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
+                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
+                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+ 
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+      END DO
+      ! end dig bathy/ice shelf to be compatible
+      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
+      DO jl = 1,20
+ 
+ ! remove single point "bay" on isf coast line in the ice shelf draft'
+         DO jk = 2, jpk
+            WHERE (misfdep==0) misfdep=jpk
+            zmask=0._wp
+            WHERE (misfdep <= jk) zmask=1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF (misfdep(ji,jj) == jk) THEN
+                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
+                     IF (ibtest <= 1) THEN
+                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
+                        IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk
+                     END IF
+                  END IF
+               END DO
+            END DO
+         END DO
+         WHERE (misfdep==jpk)
+             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
+         END WHERE
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+ 
+ ! remove single point "bay" on bathy coast line beneath an ice shelf'
+         DO jk = jpk,1,-1
+            zmask=0._wp
+            WHERE (mbathy >= jk ) zmask=1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN
+                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
+                     IF (ibtest <= 1) THEN
+                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
+                        IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0
+                     END IF
+                  END IF
+               END DO
+            END DO
+         END DO
+         WHERE (mbathy==0)
+             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
+         END WHERE
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+ 
+ ! fill hole in ice shelf
+         zmisfdep = misfdep
+         zrisfdep = risfdep
+         WHERE (zmisfdep <= 1._wp) zmisfdep=jpk
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
+               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
+               IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk
+               IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk
+               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk
+               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk
+               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
+               IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN
+                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
+               END IF
+               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN
+                  misfdep(ji,jj) = ibtest
+                  risfdep(ji,jj) = gdepw_1d(ibtest)
+               ENDIF
+            ENDDO
+         ENDDO
+ 
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep, 'T', 1. ) 
+            CALL lbc_lnk( bathy,   'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy,  'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ !
+ !! fill hole in bathymetry
+         zmbathy (:,:)=mbathy (:,:)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
+               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
+               IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0
+               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
+               IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN
+                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
+               END IF
+               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN
+                  mbathy(ji,jj) = ibtest
+                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
+               ENDIF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep, 'T', 1. ) 
+            CALL lbc_lnk( bathy,   'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy,  'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie U point water column less than 2 cells), mask U
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
+                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep, 'T', 1. ) 
+            CALL lbc_lnk( bathy,   'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy,  'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie U point water column less than 2 cells), mask U
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
+                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie V point water column less than 2 cells), mask V
+         DO jj = 1, jpjm1
+            DO ji = 1, jpi
+               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
+                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie V point water column less than 2 cells), mask V
+         DO jj = 1, jpjm1
+            DO ji = 1, jpi
+               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
+                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check, mask T
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
+                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
+               END IF
+            END DO
+         END DO
+ 
+         WHERE (mbathy(:,:) == 1)
+            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
+         END WHERE
+      END DO 
+! end check compatibility ice shelf/bathy
+      ! remove very shallow ice shelf (less than ~ 10m if 75L)
+      WHERE (risfdep(:,:) <= 10._wp)
+         misfdep = 1; risfdep = 0.0_wp;
+      END WHERE
+
+      IF( icompt == 0 ) THEN 
+         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
+      ELSE 
+         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
+      ENDIF 
+
+      ! compute scale factor and depth at T- and W- points
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1035,4 +1755,5 @@
                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
                   ENDIF
+      !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
 !gm Bug?  check the gdepw_1d
                   !       ... on ik
@@ -1040,12 +1761,9 @@
                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
-                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
-                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
-                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
-                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+                  e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  )
+                  e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1)
                   !       ... on ik+1
                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
-                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
                ENDIF
             ENDIF
@@ -1073,687 +1791,58 @@
       END DO
       !
-      IF ( ln_isfcav ) THEN
       ! (ISF) Definition of e3t, u, v, w for ISF case
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               ik = misfdep(ji,jj) 
-               IF( ik > 1 ) THEN               ! ice shelf point only 
-                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
-                  gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+            ik = misfdep(ji,jj) 
+            IF( ik > 1 ) THEN               ! ice shelf point only 
+               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
+               gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
 !gm Bug?  check the gdepw_0 
-               !       ... on ik 
-                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
-                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
-                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
-                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
-                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
-
-                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
-                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
-                  ENDIF 
-               !       ... on ik / ik-1 
-                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
-                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
+            !       ... on ik 
+               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
+                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
+                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
+               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
+               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
+
+               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
+                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
+               ENDIF 
+            !       ... on ik / ik-1 
+               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
+               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
-                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+            ENDIF 
+         END DO 
+      END DO 
+   
+      it = 0 
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+            ik = misfdep(ji,jj) 
+            IF( ik > 1 ) THEN               ! ice shelf point only 
+               e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
+               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
+            ! test 
+               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
+               IF( zdiff <= 0. .AND. lwp ) THEN  
+                  it = it + 1 
+                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
+                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
+                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
+                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
                ENDIF 
-            END DO 
+            ENDIF 
          END DO 
-      ! 
-         it = 0 
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               ik = misfdep(ji,jj) 
-               IF( ik > 1 ) THEN               ! ice shelf point only 
-                  e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
-                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
-               ! test 
-                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
-                  IF( zdiff <= 0. .AND. lwp ) THEN  
-                     it = it + 1 
-                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
-                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
-                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
-                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
-                  ENDIF 
-               ENDIF 
-            END DO 
-         END DO 
-      END IF
-      ! END (ISF)
-
-      ! Scale factors and depth at U-, V-, UW and VW-points
-      DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3u_0 (:,:,jk) = e3t_1d(jk)
-         e3v_0 (:,:,jk) = e3t_1d(jk)
-         e3uw_0(:,:,jk) = e3w_1d(jk)
-         e3vw_0(:,:,jk) = e3w_1d(jk)
-      END DO
-      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
-               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
-               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
-               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
-            END DO
-         END DO
-      END DO
-      IF ( ln_isfcav ) THEN
-      ! (ISF) define e3uw (adapted for 2 cells in the water column)
-         DO jj = 2, jpjm1 
-            DO ji = 2, fs_jpim1   ! vector opt. 
-               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
-               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
-               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
-                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
-               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
-               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
-               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
-                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
-            END DO
-         END DO
-      END IF
-
-      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
-      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
-      !
-      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
-         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
-         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
-         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
-      END DO
-      
-      ! Scale factor at F-point
-      DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3f_0(:,:,jk) = e3t_1d(jk)
-      END DO
-      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
-      !
-      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
-      END DO
-!!gm  bug ? :  must be a do loop with mj0,mj1
-      ! 
-      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_zps :   e r r o r   e3t_0 <= 0' )
-      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
-      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
-      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
-     
-      ! Compute gdep3w_0 (vertical sum of e3w)
-      IF ( ln_isfcav ) THEN ! if cavity
-         WHERE (misfdep == 0) misfdep = 1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
-               DO jk = 2, misfdep(ji,jj)
-                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
-               END DO
-               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
-               DO jk = misfdep(ji,jj) + 1, jpk
-                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
-               END DO
-            END DO
-         END DO
-      ELSE ! no cavity
-         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
-         DO jk = 2, jpk
-            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
-         END DO
-      END IF
-      !                                               ! ================= !
-      IF(lwp .AND. ll_print) THEN                     !   Control print   !
-         !                                            ! ================= !
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX( mbathy(ji,jj), 1 )
-               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
-               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
-               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
-               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
-               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
-               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-      ENDIF  
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
-      !
-   END SUBROUTINE zgr_zps
-
-   SUBROUTINE zgr_isf
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_isf  ***
-      !!   
-      !! ** Purpose :   check the bathymetry in levels
-      !!   
-      !! ** Method  :   THe water column have to contained at least 2 cells
-      !!                Bathymetry and isfdraft are modified (dig/close) to respect
-      !!                this criterion.
-      !!                 
-      !!   
-      !! ** Action  : - test compatibility between isfdraft and bathy 
-      !!              - bathy and isfdraft are modified
-      !!----------------------------------------------------------------------
-      !!   
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  ::   ik, it           ! temporary integers
-      INTEGER  ::   id, jd, nprocd
-      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
-      LOGICAL  ::   ll_print         ! Allow  control print for debugging
-      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
-      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
-      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
-      REAL(wp) ::   zdiff            ! temporary scalar
-      REAL(wp) ::   zrefdep          ! temporary scalar
-      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
-      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
-      !
-      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
-      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
-
-
-      ! (ISF) compute misfdep
-      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
-      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
-      END WHERE  
-
-      ! Compute misfdep for ocean points (i.e. first wet level) 
-      ! find the first ocean level such that the first level thickness 
-      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 
-      ! e3t_0 is the reference level thickness 
-      DO jk = 2, jpkm1 
-         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
-         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
       END DO 
-      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
-         risfdep(:,:) = 0. ; misfdep(:,:) = 1
-      END WHERE
- 
-! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
-      icompt = 0 
-! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
-      DO jl = 1, 10     
-         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
-            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
-            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
-         END WHERE
-         WHERE (mbathy(:,:) <= 0) 
-            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
-            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
-         ENDWHERE
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
-         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
-            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west 
-            misfdep(jpi,:) = misfdep(  2  ,:) 
-         ENDIF
-
-         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
-            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
-            mbathy(jpi,:) = mbathy(  2  ,:)
-         ENDIF
-
-         ! split last cell if possible (only where water column is 2 cell or less)
-         DO jk = jpkm1, 1, -1
-            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
-               mbathy(:,:) = jk
-               bathy(:,:)  = zmax
-            END WHERE
-         END DO
- 
-         ! split top cell if possible (only where water column is 2 cell or less)
-         DO jk = 2, jpkm1
-            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
-               misfdep(:,:) = jk
-               risfdep(:,:) = zmax
-            END WHERE
-         END DO
-
- 
- ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! find the minimum change option:
-               ! test bathy
-               IF (risfdep(ji,jj) .GT. 1) THEN
-               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
-                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
-               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
-                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
- 
-                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
-                     IF (zbathydiff .LE. zrisfdepdiff) THEN
-                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
-                        mbathy(ji,jj)= mbathy(ji,jj) + 1
-                     ELSE
-                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
-                        misfdep(ji,jj) = misfdep(ji,jj) - 1
-                     END IF
-                  END IF
-               END IF
-            END DO
-         END DO
- 
-          ! At least 2 levels for water thickness at T, U, and V point.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! find the minimum change option:
-               ! test bathy
-               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)&
-                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  &
-                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj)= misfdep(ji,jj) - 1
-                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
-                  END IF
-               ENDIF
-            END DO
-         END DO
- 
- ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   &
-                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  &
-                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) &
-                   &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
-                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &
-                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
-                  END IF
-               ENDIF
-            END DO
-         END DO
- 
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
- ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
-                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  &
-                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
-                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
-                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
-                  END IF
-               ENDIF
-            END DO
-         END DO
- 
- 
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- 
- ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
-                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
-                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
-                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
-                  END IF
-               ENDIF
-            ENDDO
-         ENDDO
- 
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- 
- ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
-                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  &
-                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
-                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  &
-                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
-                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) &
-                      &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
-                  END IF
-               ENDIF
-            ENDDO
-         ENDDO
- 
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
-      END DO
-      ! end dig bathy/ice shelf to be compatible
-      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
-      DO jl = 1,20
- 
- ! remove single point "bay" on isf coast line in the ice shelf draft'
-         DO jk = 2, jpk
-            WHERE (misfdep==0) misfdep=jpk
-            zmask=0
-            WHERE (misfdep .LE. jk) zmask=1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF (misfdep(ji,jj) .EQ. jk) THEN
-                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
-                     IF (ibtest .LE. 1) THEN
-                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
-                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
-                     END IF
-                  END IF
-               END DO
-            END DO
-         END DO
-         WHERE (misfdep==jpk)
-             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
-         END WHERE
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
- 
- ! remove single point "bay" on bathy coast line beneath an ice shelf'
-         DO jk = jpk,1,-1
-            zmask=0
-            WHERE (mbathy .GE. jk ) zmask=1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
-                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
-                     IF (ibtest .LE. 1) THEN
-                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
-                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
-                     END IF
-                  END IF
-               END DO
-            END DO
-         END DO
-         WHERE (mbathy==0)
-             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
-         END WHERE
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
- 
- ! fill hole in ice shelf
-         zmisfdep = misfdep
-         zrisfdep = risfdep
-         WHERE (zmisfdep .LE. 1) zmisfdep=jpk
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
-               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
-               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
-                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
-               END IF
-               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
-                  misfdep(ji,jj) = ibtest
-                  risfdep(ji,jj) = gdepw_1d(ibtest)
-               ENDIF
-            ENDDO
-         ENDDO
- 
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- !
- !! fill hole in bathymetry
-         zmbathy (:,:)=mbathy (:,:)
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
-               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
-               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
-               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
-               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
-                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
-               END IF
-               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
-                  mbathy(ji,jj) = ibtest
-                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
-               ENDIF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie U point water column less than 2 cells), mask U
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
-                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie U point water column less than 2 cells), mask U
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
-                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie V point water column less than 2 cells), mask V
-         DO jj = 1, jpjm1
-            DO ji = 1, jpi
-               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
-                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie V point water column less than 2 cells), mask V
-         DO jj = 1, jpjm1
-            DO ji = 1, jpi
-               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
-                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check, mask T
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
-                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
-               END IF
-            END DO
-         END DO
- 
-         WHERE (mbathy(:,:) == 1)
-            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
-         END WHERE
-      END DO 
-! end check compatibility ice shelf/bathy
-      ! remove very shallow ice shelf (less than ~ 10m if 75L)
-      WHERE (misfdep(:,:) <= 5)
-         misfdep = 1; risfdep = 0.0_wp;
-      END WHERE
-
-      IF( icompt == 0 ) THEN 
-         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
-      ELSE 
-         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
-      ENDIF 
 
       CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
       CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
-
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
-      
-   END SUBROUTINE
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zgr_isf')
+      !      
+   END SUBROUTINE zgr_isf
+
 
    SUBROUTINE zgr_sco
@@ -1801,5 +1890,4 @@
       !!
       !!----------------------------------------------------------------------
-      !
       INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
       INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
@@ -1810,12 +1898,12 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
-
+      !!
       NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
-                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
+         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
      !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
       !
-      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
+      CALL wrk_alloc( jpi,jpj,   zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
       !
       REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
@@ -1862,9 +1950,11 @@
       bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
-         END DO
-      END DO
+      IF( .NOT.ln_wd ) THEN                      
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
+            END DO
+         END DO
+      END IF
       !                                        ! =============================
       !                                        ! Define the envelop bathymetry   (hbatt)
@@ -1873,19 +1963,29 @@
       zenv(:,:) = bathy(:,:)
       !
+      IF( .NOT.ln_wd ) THEN    
       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           IF( bathy(ji,jj) == 0._wp ) THEN
-             iip1 = MIN( ji+1, jpi )
-             ijp1 = MIN( jj+1, jpj )
-             iim1 = MAX( ji-1, 1 )
-             ijm1 = MAX( jj-1, 1 )
-             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
-        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
-               zenv(ji,jj) = rn_sbot_min
-             ENDIF
-           ENDIF
-         END DO
-      END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( bathy(ji,jj) == 0._wp ) THEN
+                  iip1 = MIN( ji+1, jpi )
+                  ijp1 = MIN( jj+1, jpj )
+                  iim1 = MAX( ji-1, 1 )
+                  ijm1 = MAX( jj-1, 1 )
+!!gm BUG fix see ticket #1617
+                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
+                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
+                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
+                     &    zenv(ji,jj) = rn_sbot_min
+!!gm
+!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
+!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
+!!gm               zenv(ji,jj) = rn_sbot_min
+!!gm             ENDIF
+!!gm end
+              ENDIF
+            END DO
+         END DO
+      END IF
+
       ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
       CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
@@ -1975,15 +2075,4 @@
       ENDIF
       !
-      IF(lwp) THEN                             ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
-         WRITE(numout,*)
-         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
-         IF( nprint == 1 )   THEN        
-            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
-            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
-         ENDIF
-      ENDIF
-
       !                                        ! ==============================
       !                                        !   hbatu, hbatv, hbatf fields
@@ -1991,5 +2080,10 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+         IF( .NOT.ln_wd ) THEN
+           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+         ELSE
+           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min
+           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld
+         ENDIF
       ENDIF
       hbatu(:,:) = rn_sbot_min
@@ -2004,4 +2098,19 @@
         END DO
       END DO
+
+      IF( ln_wd ) THEN               !avoid the zero depth on T- (U-,V-,F-) points
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) &
+              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) &
+              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) &
+              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) &
+              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1
+          END DO
+        END DO
+      END IF
       ! 
       ! Apply lateral boundary condition
@@ -2011,4 +2120,5 @@
          DO ji = 1, jpi
             IF( hbatu(ji,jj) == 0._wp ) THEN
+               !No worries about the following line when ln_wd == .true.
                IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
                IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
@@ -2036,8 +2146,10 @@
 
 !!bug:  key_helsinki a verifer
-      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
-      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
-      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
-      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+      IF( .NOT.ln_wd ) THEN
+        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
+        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
+        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
+        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+      END IF
 
       IF( nprint == 1 .AND. lwp )   THEN
@@ -2080,49 +2192,39 @@
       CALL lbc_lnk( e3uw_0, 'U', 1._wp )
       CALL lbc_lnk( e3vw_0, 'V', 1._wp )
-
-      fsdepw(:,:,:) = gdepw_0 (:,:,:)
-      fsde3w(:,:,:) = gdep3w_0(:,:,:)
-      !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      !
+      IF( .NOT.ln_wd ) THEN
+        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
+        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
+        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
+        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
+        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
+        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
+        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
+      END IF
 
 #if defined key_agrif
-      ! Ensure meaningful vertical scale factors in ghost lines/columns
-      IF( .NOT. Agrif_Root() ) THEN
-         !  
-         IF((nbondi == -1).OR.(nbondi == 2)) THEN
-            e3u_0(1,:,:) = e3u_0(2,:,:)
-         ENDIF
-         !
-         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
-            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
-         ENDIF
-         !
-         IF((nbondj == -1).OR.(nbondj == 2)) THEN
-            e3v_0(:,1,:) = e3v_0(:,2,:)
-         ENDIF
-         !
-         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
-            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
-         ENDIF
-         !
-      ENDIF
+      IF( .NOT. Agrif_Root() ) THEN    ! Ensure meaningful vertical scale factors in ghost lines/columns
+         IF( nbondi == -1 .OR. nbondi == 2 )   e3u_0(  1   ,  :   ,:) = e3u_0(  2   ,  :   ,:)
+         IF( nbondi ==  1 .OR. nbondi == 2 )   e3u_0(nlci-1,  :   ,:) = e3u_0(nlci-2,  :   ,:)
+         IF( nbondj == -1 .OR. nbondj == 2 )   e3v_0(  :   ,  1   ,:) = e3v_0(  :   ,  2   ,:)
+         IF( nbondj ==  1 .OR. nbondj == 2 )   e3v_0(  :   ,nlcj-1,:) = e3v_0(  :   ,nlcj-2,:)
+       ENDIF
 #endif
 
-      fsdept(:,:,:) = gdept_0 (:,:,:)
-      fsdepw(:,:,:) = gdepw_0 (:,:,:)
-      fsde3w(:,:,:) = gdep3w_0(:,:,:)
-      fse3t (:,:,:) = e3t_0   (:,:,:)
-      fse3u (:,:,:) = e3u_0   (:,:,:)
-      fse3v (:,:,:) = e3v_0   (:,:,:)
-      fse3f (:,:,:) = e3f_0   (:,:,:)
-      fse3w (:,:,:) = e3w_0   (:,:,:)
-      fse3uw(:,:,:) = e3uw_0  (:,:,:)
-      fse3vw(:,:,:) = e3vw_0  (:,:,:)
+!!gm   I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays)
+!!gm   and only that !!!!!
+!!gm   THIS should be removed from here !
+      gdept_n(:,:,:) = gdept_0(:,:,:)
+      gdepw_n(:,:,:) = gdepw_0(:,:,:)
+      gde3w_n(:,:,:) = gde3w_0(:,:,:)
+      e3t_n  (:,:,:) = e3t_0  (:,:,:)
+      e3u_n  (:,:,:) = e3u_0  (:,:,:)
+      e3v_n  (:,:,:) = e3v_0  (:,:,:)
+      e3f_n  (:,:,:) = e3f_0  (:,:,:)
+      e3w_n  (:,:,:) = e3w_0  (:,:,:)
+      e3uw_n (:,:,:) = e3uw_0 (:,:,:)
+      e3vw_n (:,:,:) = e3vw_0 (:,:,:)
+!!gm and obviously in the following, use the _0 arrays until the end of this subroutine
+!! gm end
 !!
       ! HYBRID : 
@@ -2130,7 +2232,15 @@
          DO ji = 1, jpi
             DO jk = 1, jpkm1
-               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
-            END DO
-            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
+               IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
+            END DO
+            IF( ln_wd ) THEN
+              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
+                mbathy(ji,jj) = 0
+              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN
+                mbathy(ji,jj) = 2
+              ENDIF
+            ELSE
+              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0
+            ENDIF
          END DO
       END DO
@@ -2141,15 +2251,15 @@
          WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
          WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
-            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
-            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
-            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
+            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
+            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
+            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
             &                          ' w ', MINVAL( e3w_0  (:,:,:) )
 
          WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
-            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
-            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
-            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
+            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gde3w_0(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
+            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
+            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
             &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
       ENDIF
@@ -2183,50 +2293,48 @@
          END DO
       ENDIF
-
-!================================================================================
-! check the coordinate makes sense
-!================================================================================
+      !
+      !================================================================================
+      ! check the coordinate makes sense
+      !================================================================================
       DO ji = 1, jpi
          DO jj = 1, jpj
-
+            !
             IF( hbatt(ji,jj) > 0._wp) THEN
                DO jk = 1, mbathy(ji,jj)
                  ! check coordinate is monotonically increasing
-                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
+                 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
-                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
+                    WRITE(numout,*) 'e3w',e3w_n(ji,jj,:)
+                    WRITE(numout,*) 'e3t',e3t_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                  ! and check it has never gone negative
-                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
+                 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
-                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
+                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                  ! and check it never exceeds the total depth
-                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
+                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                END DO
-
+               !
                DO jk = 1, mbathy(ji,jj)-1
                  ! and check it never exceeds the total depth
-                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
                     WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
-                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
                     CALL ctl_stop( ctmp1 )
                  ENDIF
                END DO
-
             ENDIF
-
          END DO
       END DO
@@ -2238,7 +2346,6 @@
    END SUBROUTINE zgr_sco
 
-!!======================================================================
+
    SUBROUTINE s_sh94()
-
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE s_sh94  ***
@@ -2251,13 +2358,15 @@
       !! Reference : Song and Haidvogel 1994. 
       !!----------------------------------------------------------------------
-      !
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
+      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
-
-      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
-      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_alloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
 
       z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
@@ -2265,8 +2374,8 @@
       z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
       z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
-
+      !
       DO ji = 1, jpi
          DO jj = 1, jpj
-
+            !
             IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
                DO jk = 1, jpk
@@ -2297,7 +2406,7 @@
                zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
                zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
-               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
-               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
+               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
             END DO
            !
@@ -2307,16 +2416,43 @@
       DO ji = 1, jpim1
          DO jj = 1, jpjm1
+            ! extended for Wetting/Drying case
+            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
+            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
+            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
+            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
+            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
+                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
             DO jk = 1, jpk
-               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
-                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
+                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+               ELSE
+                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
+                        &              / ztmpu
+                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
+                        &              / ztmpu
+               END IF
+
+               IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
+                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+               ELSE
+                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
+                        &              / ztmpv
+                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
+                        &              / ztmpv
+               END IF
+
+               IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
+                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj  ,jk) + z_esigt3(ji+1,jj  ,jk)               &
+                        &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
+               ELSE
+                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
+                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
+                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
+                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
+               END IF
+
                !
                e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
@@ -2331,12 +2467,12 @@
         END DO
       END DO
-
-      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
-      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
-
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+      !
    END SUBROUTINE s_sh94
 
+
    SUBROUTINE s_sf12
-
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE s_sf12 *** 
@@ -2354,12 +2490,13 @@
       !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
       !!----------------------------------------------------------------------
-      !
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zsmth               ! smoothing around critical depth
       REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      REAL(wp) ::   ztmpu, ztmpv, ztmpf
+      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
-
+      !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
@@ -2425,7 +2562,7 @@
 
           DO jk = 1, jpk
-             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
-             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
-             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
+             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
+             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
+             gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
           END DO
 
@@ -2436,16 +2573,61 @@
         DO jj=1,jpj-1
 
-          DO jk = 1, jpk
-                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
-                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+           ! extend to suit for Wetting/Drying case
+           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
+           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
+           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
+           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
+           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
+                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
+           DO jk = 1, jpk
+              IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
+                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+              ELSE
+                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
+                       &              / ztmpu
+                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
+                       &              / ztmpu
+              END IF
+
+              IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
+                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+              ELSE
+                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
+                       &              / ztmpv
+                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
+                       &              / ztmpv
+              END IF
+
+              IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
+                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)                 &
+                       &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
+              ELSE
+                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
+                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
+                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
+                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
+              END IF
+
+             ! Code prior to wetting and drying option (for reference)
+             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
+             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
+             !
+             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
+             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
+             !
+             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
+             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
+             !
+             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
+             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
+             !
+             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
+             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
+             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
+             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
+             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
 
              e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
@@ -2454,5 +2636,5 @@
              e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
              !
-             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
+             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
              e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
              e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
@@ -2467,13 +2649,11 @@
       CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
       !
-      !                                               ! =============
-
-      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
-      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
-
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+      !
    END SUBROUTINE s_sf12
 
+
    SUBROUTINE s_tanh()
-
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE s_tanh*** 
@@ -2485,13 +2665,12 @@
       !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
       !!----------------------------------------------------------------------
-
-      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      INTEGER  ::   ji, jj, jk       ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
-
       REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
       REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
-
-      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
-      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
+      CALL wrk_alloc( jpk,   z_esigt, z_esigw )
 
       z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
@@ -2523,7 +2702,7 @@
          zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
          zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
-         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
-         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
+         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
+         gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
       END DO
 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
@@ -2542,9 +2721,10 @@
          END DO
       END DO
-
-      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
-      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
-
+      !
+      CALL wrk_dealloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
+      CALL wrk_dealloc( jpk,   z_esigt, z_esigw          )
+      !
    END SUBROUTINE s_tanh
+
 
    FUNCTION fssig( pk ) RESULT( pf )
@@ -2618,33 +2798,32 @@
       REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
       REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
-      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
-      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
-      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
-      REAL(wp)                ::   za1,za2,za3    ! local variables
-      REAL(wp)                ::   zn1,zn2        ! local variables
-      REAL(wp)                ::   za,zb,zx       ! local variables
-      integer                 ::   jk
-      !!----------------------------------------------------------------------
-      !
-
-      zn1  =  1./(jpk-1.)
-      zn2  =  1. -  zn1
-
+      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
+      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
+      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
+      !
+      INTEGER  ::   jk             ! dummy loop index
+      REAL(wp) ::   za1,za2,za3    ! local scalar
+      REAL(wp) ::   zn1,zn2        !   -      - 
+      REAL(wp) ::   za,zb,zx       !   -      - 
+      !!----------------------------------------------------------------------
+      !
+      zn1  =  1._wp / REAL( jpkm1, wp )
+      zn2  =  1._wp -  zn1
+      !
       za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
       za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
       za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
-     
+      !
       za = pzb - za3*(pzs-za1)-za2
       za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
       zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
       zx = 1.0_wp-za/2.0_wp-zb
- 
+      !
       DO jk = 1, jpk
-        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
-                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
-                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
+         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
+            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
+            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
         p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
-      ENDDO 
-
+      END DO
       !
    END FUNCTION fgamma
