Index: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/dom_oce.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_domcfg/src/dom_oce.F90	(revision 11201)
+++ NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/dom_oce.F90	(revision 11308)
@@ -36,4 +36,9 @@
    REAL(wp), PUBLIC ::   rn_e3zps_rat    !: minimum thickness ration for partial steps
    INTEGER , PUBLIC ::   nn_msh          !: = 1 create a mesh-mask file
+
+   INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file
+   INTEGER, PUBLIC  ::   nperio          !: type of lateral boundary condition
+   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters)
+   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps
 
    INTEGER, PUBLIC :: nn_interp
Index: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domain.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_domcfg/src/domain.F90	(revision 11201)
+++ NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domain.F90	(revision 11308)
@@ -105,6 +105,6 @@
          &             ln_cfmeta, ln_iscpl
       NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp,                        &
-         &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin,           &
-         &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                  &
+         &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,                       &
+         &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                                &
          &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
          &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
@@ -204,5 +204,4 @@
          WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
          WRITE(numout,*) '      min number of ocean level (<0)       '
-         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
          WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
          WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
@@ -415,5 +414,15 @@
          END DO
       END DO
-      CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , z2d , ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum, 'bathy_metry_e3'   , z2d , ktype = jp_r4 )
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            z2d (ji,jj) = SUM ( e3t_0(ji,jj, 1:mikt(ji,jj)-1 ) ) * ssmask(ji,jj) 
+         END DO
+      END DO
+      CALL iom_rstput( 0, 0, inum, 'isf_draft_e3'   , z2d , ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum, 'isf_draft'   , risfdep , ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , bathy , ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum, 'hw',bathy-risfdep, ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum, 'mhw',mbkt*ssmask-mikt*ssmask, ktype = jp_i4 )
       !
       !                              !== closed sea ==!
Index: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domisf.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domisf.F90	(revision 11308)
+++ NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domisf.F90	(revision 11308)
@@ -0,0 +1,563 @@
+MODULE domisf
+   !!==============================================================================
+   !!                       ***  MODULE domisf   ***
+   !! Ocean domain : 
+   !!==============================================================================
+   !! History :  4.1  ! 2019-07  (P. Mathiot) re-write of the geometry definition under ice shelf
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   zgr_isf_nam    : read ice shelf namelist; print and compatibility option check
+   !!   zgr_isf_zspace : check compatibility between bathymetry and ice shelf draft and adjust ice shelf draft if needed (z space) 
+   !!   zgr_isf_kspace : adjust ice shelf draft and compute top level to fit 2 cell in the water column (k space) 
+   !!   zgr_isf        : compute top partial cell scale factor
+   !!   zgr_isf_e3uv_w : correct e3uw and e3vw in case of 2 cell in water column under an ice shelf
+   !!---------------------------------------------------------------------
+   USE dom_oce
+   USE domutil           ! flood filling algorithm
+   USE domngb            ! find nearest neighbourg
+   USE in_out_manager    ! I/O manager
+   USE lbclnk            ! lbclnk and lbclnk_multi
+   USE lib_mpp           ! 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zgr_isf_nam    !: read domisf namelist
+   PUBLIC   zgr_isf_zspace !: check compatibility ice shelf draft/bathymetry and dig ice shelf draft if needed
+   PUBLIC   zgr_isf_kspace !: compute misfdep and dig ice shelf draft if needed
+   PUBLIC   zps_isf        !: compute zps scale factor
+   PUBLIC   zps_isf_e3uv_w !: compute e3uw and e3vw in case of 2 cell in water column
+   PUBLIC   zgr_isf_subgl
+
+   INTEGER          :: nn_kisfmax    = 999              !: maximal number of level change allowed by ln_isfconnect option 
+   REAL(wp), PUBLIC :: rn_isfdep_min = 10.0_wp          !: ice shelf minimal thickness 
+   REAL(wp), PUBLIC :: rn_glhw_min   = 1.0e-3           !: threshold on hw to define grounding line water
+   REAL(wp), PUBLIC :: rn_isfhw_min  = 1.0e-3           !: threshold on hw to define isf draft into the cavity
+   REAL(wp), PUBLIC :: rn_isfshallow = 0.0_wp           !: threshold to define shallow ice shelf cavity
+   REAL(wp), PUBLIC :: rn_zisfmax    = 6000.0_wp        !: maximun meter of ice we are allowed to dig to assure connectivity
+   REAL(wp), PUBLIC :: rn_isfsubgllon, rn_isfsubgllat   !: seed for the detection of the 3d open ocean
+   LOGICAL , PUBLIC :: ln_isfcheminey= .FALSE.          !: remove cheminey
+   LOGICAL , PUBLIC :: ln_isfconnect = .FALSE.          !: assure connectivity 
+   LOGICAL , PUBLIC :: ln_isfchannel = .FALSE.          !: remove channel in water column (on z space not level space)
+   LOGICAL , PUBLIC :: ln_isfsubgl   = .FALSE.          !: remove subglacial lakes
+
+CONTAINS
+
+   SUBROUTINE zgr_isf_nam
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zps_isf  ***
+      !!   
+      !! ** Purpose :   Read namzgr_isf namelist and 
+      !!                check compatibility with other option
+      !!   
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER  :: ios, ierr 
+      !!---------------------------------------------------------------------
+      NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, rn_isfdep_min, rn_glhw_min, rn_isfhw_min, &
+         &                ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
+      !
+      ! 0.0 read namelist
+      REWIND( numnam_ref )              ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition
+      READ  ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition
+      READ  ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzgr_isf )
+      !
+      IF (lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' dom_isf : '
+         WRITE(numout,*) ' ~~~~~~~   '
+         WRITE(numout,*) ' nn_kisfmax       = ',nn_kisfmax           !
+         WRITE(numout,*) ' rn_zisfmax       = ',rn_zisfmax           !
+         WRITE(numout,*) ' rn_isfdep_min    = ',rn_isfdep_min        !
+         WRITE(numout,*) ' rn_isfhw_min     = ',rn_isfhw_min         !
+         WRITE(numout,*) ' rn_glhw_min      = ',rn_glhw_min          !
+         WRITE(numout,*) ' ln_isfcheminey   = ',ln_isfcheminey       !
+         WRITE(numout,*) ' ln_isfconnect    = ',ln_isfconnect        !
+         WRITE(numout,*) ' ln_isfchannel    = ',ln_isfchannel        !
+      END IF
+      !
+      ! 0.1 compatibility option
+      ierr = 0
+      IF ( ln_zco .OR. ln_sco ) ierr = ierr + 1
+      IF ( ierr > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
+
+   END SUBROUTINE zgr_isf_nam
+
+   SUBROUTINE zgr_isf_zspace
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dom_isf  ***
+      !!   
+      !! ** Purpose :   compute the grounding line position
+      !!   
+      !! ** Method  :   bathy set to 0 and isfdraft are modified (dig/close) to respect
+      !!                these criterions.
+      !!                digging and filling base on depth criterion only
+      !!                          1.0 = set iceshelf to the minimum depth allowed
+      !!                          1.1 = ground ice shelf if water column less than X m
+      !!                          1.2 = ensure a minimum thickness for iceshelf cavity in shallow water
+      !!                          1.3 = remove channels and single point 'bay'
+      !!   
+      !! ** Action  : - test compatibility between isfdraft and bathy 
+      !!              - bathy and isfdraft are modified
+      !!----------------------------------------------------------------------
+      !!  
+      INTEGER  :: ji, jj                             ! loop indexes
+      INTEGER  :: imskjp1, imskjm1, imskip1, imskim1 ! local variable
+      INTEGER, DIMENSION(jpi,jpj) :: imask           ! isf mask
+      !
+      REAL(wp) ::   zisfdep_min  ! minimal ice shelf draft allowed
+      !!---------------------------------------------------------------------
+      !
+      ! 1.0 set iceshelf to the minimum depth allowed
+      zisfdep_min = MAX(rn_isfdep_min,e3t_1d(1))
+      WHERE(risfdep(:,:) > 0.0_wp .AND. risfdep(:,:) < zisfdep_min)
+         risfdep(:,:) = zisfdep_min
+      END WHERE
+      !  
+      ! 1.1 ground ice shelf if water column less than rn_glhw_min m 
+      ! => set the grounding line position
+      WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min ) 
+         bathy  (:,:) = 0._wp ; risfdep(:,:) = 0._wp
+      END WHERE
+      !
+      ! 1.2 ensure a minimum thickness for iceshelf cavity 
+      ! => avoid to negative e3t if ssh + sum(e3t_0*tmask) < 0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( bathy(ji,jj) - risfdep(ji,jj) < rn_isfhw_min ) THEN
+               risfdep(ji,jj) = bathy(ji,jj) - rn_isfhw_min
+               ! sanity check on risfdep (if < zisfdep_min)
+               ! => we ground it as it failed to respect condition 1.0 and 1.1
+               IF ( risfdep(ji,jj) < zisfdep_min ) THEN 
+                  bathy(ji,jj)=0._wp ; risfdep(ji,jj)=0._wp
+               END IF
+            END IF
+         END DO
+      END DO
+      !
+      ! 1.3 Remove channels and single point 'bay' using bathy mask.
+      ! => channel could be created if connectivity is enforced later.
+      IF (ln_isfchannel) THEN
+         imask(:,:) = 0
+         WHERE ( bathy(:,:) > 0._wp .AND. risfdep(:,:) > 0._wp )
+            imask(:,:) = 1
+         END WHERE
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF( imask(ji,jj) == 1 ) THEN
+                  ! define connexion
+                  imskip1=imask(ji,jj)*imask(ji+1,jj  )  ! 1 = connexion
+                  imskjp1=imask(ji,jj)*imask(ji  ,jj+1)  ! 1 = connexion
+                  imskim1=imask(ji,jj)*imask(ji-1,jj  )  ! 1 = connexion
+                  imskjm1=imask(ji,jj)*imask(ji  ,jj-1)  ! 1 = connexion
+                  ! zonal channel and single bay
+                  IF ((imskip1+imskim1>=1) .AND. (imskjp1+imskjm1==0)) THEN
+                     bathy(ji,jj)=0._wp ; risfdep(ji,jj)=0._wp
+                  END IF
+                  ! meridionnal channel and single bay
+                  IF ((imskjp1+imskjm1>=1) .AND. (imskip1+imskim1==0)) THEN
+                     bathy(ji,jj)=0._wp ; risfdep(ji,jj)=0._wp
+                  END IF
+               END IF
+            END DO
+         END DO
+         ! ensure halo correct 
+         CALL lbc_lnk_multi( 'domisf', risfdep, 'T', 1._wp, bathy  , 'T', 1._wp )
+      END IF
+      !
+   END SUBROUTINE zgr_isf_zspace
+
+   SUBROUTINE zgr_isf_kspace
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dom_isf  ***
+      !!   
+      !! ** Purpose :   compute misfdep and update isf draft if needed
+      !!   
+      !! ** Method  :   The water column has to contain at least 2 wet cells in the water column under an ice shelf
+      !!                isf draft is modified (dig/close) to respect this criterion.
+      !!                compute level
+      !!                          2.0 = compute level
+      !!                          2.1 = ensure misfdep is not below bathymetry after step 2.0
+      !!                digging and filling base on level criterion only                
+      !!                          3.0 = dig to fit the 2 water level criterion (closed pool possible after this step)
+      !!                          3.1 = dig enough to ensure connectivity of all the cell beneath an ice shelf 
+      !!                                (most of the close pool should be remove after this step)
+      !!                          3.2 = fill chimney
+      !!   
+      !! ** Action  : - compute misfdep
+      !!              - isf draft is modified if needed
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jn
+      INTEGER  ::   icompt
+      INTEGER  ::   ibtest,           ibtestim1, ibtestip1, ibtestjm1, ibtestjp1
+      INTEGER  ::   ibathy, ibathyij, ibathyim1, ibathyip1, ibathyjm1, ibathyjp1
+      INTEGER  ::   imskjp1  , imskjm1  , imskip1  , imskim1
+      INTEGER  ::   imskip1_r, imskim1_r, imskjp1_r, imskjm1_r
+      INTEGER , DIMENSION(jpi,jpj) :: imbathy, imisfdep, imask
+      !
+      REAL(wp) ::   zdepth
+      REAL(wp), DIMENSION(jpi,jpj) :: zrisfdep, zdummy   ! 2D workspace (ISH)
+      !!---------------------------------------------------------------------
+      !
+      ! 2 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 
+      !
+      ! 2.0 compute misfdep (taking into account the minimal cell thickness allowed)
+      ! set misfdep to 1 on open water and initialisation beneath ice shelf
+      WHERE( risfdep(:,:) == 0._wp ) ;   misfdep(:,:) = 1   ! open water or grounded : set misfdep to 1  
+      ELSEWHERE                      ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
+      END WHERE  
+      !
+      DO jk = 2, jpkm1 
+         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
+         WHERE( risfdep(:,:) > 0._wp .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
+      END DO 
+      !
+      ! 2.1 fill isolated grid point in the bathymetry
+      ! will be done again later on in zgr_bat_ctl, but need to be done here to adjust misfdep respectively
+      icompt = 0
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
+               &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
+               IF( ibtest < mbathy(ji,jj) ) THEN
+                  mbathy(ji,jj) = ibtest
+                  icompt = icompt + 1
+               END IF
+         END DO
+      END DO
+      !
+      ! ensure halo correct 
+      zdummy(:,:) = FLOAT( mbathy(:,:) ) ; CALL lbc_lnk('domisf', zdummy, 'T', 1._wp ) ; mbathy(:,:) = INT( zdummy(:,:) )
+      !
+      IF( lk_mpp )   CALL mpp_sum('domisf', icompt )
+      IF( icompt == 0 ) THEN
+         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
+      ELSE
+         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
+      ENDIF
+      !
+      ! 2.2 be sure misfdep not below mbathy 
+      ! warning because of condition 4 we could have 'wet cell with misfdep below mbathy
+      ! risfdep of these cells will be fix later on (see 3)
+      WHERE( misfdep > mbathy ) misfdep(:,:) = MAX( 1, mbathy(:,:) )
+      !
+      ! 3.0 Assure 2 wet cells in the water column at T point and along the edge.
+      ! find the deepest isfdep level that fit the 2 wet cell on the water column
+      ! on all the sides (still need 4 pass)
+      ! It need 4 iterations: if un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1
+      DO jn = 1, 4
+         imisfdep = misfdep
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ! ISF cell only
+               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
+                  ibathyij  = mbathy(ji  ,jj)
+                  !
+                  ! set ground edge value to jpk to skip it later on
+                  ibathyip1 = mbathy(ji+1,jj) ; IF ( ibathyip1 < misfdep(ji,jj) ) ibathyip1 = jpk ! no wet cell in common on this edge
+                  ibathyim1 = mbathy(ji-1,jj) ; IF ( ibathyim1 < misfdep(ji,jj) ) ibathyim1 = jpk ! no wet cell in common on this edge
+                  ibathyjp1 = mbathy(ji,jj+1) ; IF ( ibathyjp1 < misfdep(ji,jj) ) ibathyjp1 = jpk ! no wet cell in common on this edge
+                  ibathyjm1 = mbathy(ji,jj-1) ; IF ( ibathyjm1 < misfdep(ji,jj) ) ibathyjm1 = jpk ! no wet cell in common on this edge
+                  !
+                  ! find shallowest bathy level among the current cell and the neigbourging cells
+                  ibathy = MIN(ibathyij,ibathyip1,ibathyim1,ibathyjp1,ibathyjm1)
+                  !
+                  ! update misfdep and risfdep if needed
+                  ! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the water column
+                  jk = MIN(misfdep(ji,jj),ibathy-1)
+                  IF ( jk < misfdep(ji,jj) ) THEN
+                     imisfdep(ji,jj) = jk
+                     risfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+         misfdep=imisfdep
+         !
+         ! ensure halo correct before new pass 
+         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
+         CALL lbc_lnk('domisf', risfdep, 'T', 1. )
+      END DO ! jn
+      !
+      ! 3.1 condition block to inssure connectivity everywhere beneath an ice shelf
+      IF (ln_isfconnect) THEN
+         imask(:,:) = 1
+         imbathy  = mbathy
+         imisfdep = misfdep
+         zrisfdep = risfdep
+         WHERE ( mbathy(:,:) == 0 )
+            imask(:,:) = 0
+            imbathy(:,:) = jpk
+         END WHERE
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
+                  !
+                  ! what it should be
+                  imskip1 = imask(ji,jj) * imask(ji+1,jj  )  ! 1 = should be connected
+                  imskim1 = imask(ji,jj) * imask(ji-1,jj  )  ! 1 = should be connected
+                  imskjp1 = imask(ji,jj) * imask(ji  ,jj+1)  ! 1 = should be connected
+                  imskjm1 = imask(ji,jj) * imask(ji  ,jj-1)  ! 1 = should be connected
+                  !
+                  ! what it is
+                  imskip1_r=jpk ; imskim1_r=jpk; imskjp1_r=jpk; imskjm1_r=jpk
+                  IF (misfdep(ji,jj) > imbathy(ji+1,jj  )) imskip1_r=1.0 ! 1 = no effective connection
+                  IF (misfdep(ji,jj) > imbathy(ji-1,jj  )) imskim1_r=1.0 ! 1 = no effective connection
+                  IF (misfdep(ji,jj) > imbathy(ji  ,jj+1)) imskjp1_r=1.0 ! 1 = no effective connection
+                  IF (misfdep(ji,jj) > imbathy(ji  ,jj-1)) imskjm1_r=1.0 ! 1 = no effective connection
+                  !
+                  ! defining level needed for connectivity
+                  ! imskip1 * imskip1_r == 1 means connections need to be enforce
+                  jk=MIN(imbathy(ji+1,jj  ) * imskip1_r * imskip1, &
+                     &   imbathy(ji-1,jj  ) * imskim1_r * imskim1, &
+                     &   imbathy(ji  ,jj+1) * imskjp1_r * imskjp1, &
+                     &   imbathy(ji  ,jj-1) * imskjm1_r * imskjm1, &
+                     &   jpk+1 ) ! add jpk in the MIN to avoid out of boundary later on
+                  !
+                  ! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep
+                  imisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1)
+                  !
+                  ! update isf draft if needed (need to be done here because next test is in meter and level)
+                  IF ( imisfdep(ji,jj) < misfdep(ji,jj) ) THEN
+                     jk = imisfdep(ji,jj)
+                     zrisfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+                  END IF
+   
+                  ! sanity check
+                  ! check if we dig more than nn_kisfmax level or reach the surface
+                  ! check if we dig more than rn_zisfmax meter
+                  ! => if this is the case, undo what has been done before
+                  IF (      (misfdep(ji,jj)-imisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) &
+                     & .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj)  )) ) THEN
+                     imisfdep(ji,jj)=misfdep(ji,jj) 
+                     zrisfdep(ji,jj)=risfdep(ji,jj)
+                  END IF
+               END IF
+            END DO
+         END DO
+         misfdep=imisfdep
+         risfdep=zrisfdep
+         !
+         ! ensure halo correct 
+         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
+         CALL lbc_lnk('domisf', risfdep, 'T', 1. )
+      END IF
+      !
+      ! 3.2 fill hole in ice shelf (ie cell with no velocity point)
+      !      => misfdep = MIN(misfdep at U, U-1, V, V-1)
+      !         risfdep = gdepw(misfdep) (ie full cell)
+      IF (ln_isfcheminey) THEN
+         imisfdep = misfdep
+         WHERE (mbathy == 0) imisfdep=jpk   ! grounded
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ibtest    = imisfdep(ji  ,jj  )
+               ibtestim1 = imisfdep(ji-1,jj  ) ; ibtestip1 = imisfdep(ji+1,jj  )
+               ibtestjm1 = imisfdep(ji  ,jj-1) ; ibtestjp1 = imisfdep(ji  ,jj+1)
+               !
+               ! correction in case isfdraft(ii,jj) deeper than bathy U-1/U+1/V-1/V+1
+               IF( imisfdep(ji,jj) > mbathy(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1
+               IF( imisfdep(ji,jj) > mbathy(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1
+               IF( imisfdep(ji,jj) > mbathy(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
+               IF( imisfdep(ji,jj) > mbathy(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
+               !
+               ! correction in case bathy(ii,jj) shallower than isfdraft U-1/U+1/V-1/V+1
+               IF( mbathy(ji,jj) < imisfdep(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1
+               IF( mbathy(ji,jj) < imisfdep(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1
+               IF( mbathy(ji,jj) < imisfdep(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
+               IF( mbathy(ji,jj) < imisfdep(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
+               ! 
+               ! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job) 
+               ! if misfdep is not changed, nothing is done
+               ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1))
+               IF (misfdep(ji,jj) < ibtest) THEN
+                  misfdep(ji,jj) = ibtest
+                  risfdep(ji,jj) = gdepw_1d(ibtest)
+               END IF
+            ENDDO
+         ENDDO
+         !
+         ! if surround is fully masked, we mask it
+         WHERE( misfdep(:,:) == jpk)   ! ground case (1)
+            mbathy(:,:) = 0 ; bathy(:,:) = 0.0_wp ; misfdep(:,:) = 1 ; risfdep(:,:) = 0.0_wp
+         END WHERE
+         !
+         ! ensure halo correct
+         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
+         zdummy(:,:) = FLOAT( mbathy (:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); mbathy (:,:) = INT( zdummy(:,:) )
+         CALL lbc_lnk('domisf', risfdep, 'T', 1. )
+         CALL lbc_lnk('domisf', bathy  , 'T', 1. )
+      END IF
+      !
+   END SUBROUTINE zgr_isf_kspace
+
+   SUBROUTINE zgr_isf_subgl
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zps_isf_subg  ***
+      !!   
+      !! ** Purpose :   remove subglacial lakes
+      !!   
+      !! ** Method  :   Use a flood filling algorithm to detect the open ocean
+      !!                and reverse the mask
+      !!   
+      !! ** Action  : - ckeck if seed on land
+      !!              - detect main ocean
+      !!              - mask subglacial lake (closed sea under an ice shelf)
+      !!              - mask properly mbathy, misfdep, bathy, risfdep
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER  :: jk                       ! loop index
+      INTEGER  :: iiseed, ijseed           ! seed for flood filling algorithm
+      INTEGER, DIMENSION(jpi,jpj) :: imask ! 2d mask
+      !
+      REAL(wp)               :: zdist      ! distance between the seed and the nearest neighbourg                
+      REAL(wp), DIMENSION(1) :: zchk       ! sanity check variable
+      !!----------------------------------------------------------------------
+      !
+      ! check closed wet pool
+      CALL dom_ngb(rn_isfsubgllon, rn_isfsubgllat, iiseed, ijseed, zdist, 'T')
+      !
+      ! fill open ocean
+      CALL fill_pool( iiseed, ijseed, 1, tmask, -1._wp )
+      ! at this point tmask:,:,:) can have 3 different values :
+      !              0 where there where already 0
+      !             -1 where the ocean points are connected
+      !              1 where ocean points in tmask are not connected
+      IF (lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'dom_isf_subg : removal of subglacial lakes '
+         WRITE(numout,*)'~~~~~~~~~~~~'
+         WRITE(numout,*)'Number of disconected points : ', COUNT(  (tmask(:,:,:) == 1) )
+         WRITE(numout,*)'lon/lat seed to detect main ocean is: ', rn_isfsubgllon, rn_isfsubgllat
+         WRITE(numout,*)'i/j     seed to detect main ocean is: ', iiseed, ijseed
+         WRITE(numout,*)
+      END IF
+      !
+      DO jk = 1, jpk
+         ! remove only subglacial lake (ie similar to close sea only below an ice shelf)
+         WHERE (tmask(:,:,jk) > 0 .AND. misfdep(:,:) > 1) tmask(:,:,jk) = 0
+         WHERE (tmask(:,:,jk) < 0) tmask(:,:,jk) = 1
+         !
+      END DO
+      !
+      ! surface mask
+      imask(:,:)   = MAXVAL( tmask(:,:,:), DIM=3 )
+      !
+      ! mask mbathy, misfdep, bathy and risfdep 
+      bathy(:,:)   = bathy(:,:)   * imask(:,:) 
+      mbathy(:,:)  = mbathy(:,:)  * imask(:,:)
+      risfdep(:,:) = risfdep(:,:) * imask(:,:) 
+      misfdep(:,:) = misfdep(:,:) * imask(:,:)
+      !
+      ! sanity check on seed location (land or not)
+      zchk = 0._wp
+      IF (mi0(iiseed) == mi1(iiseed) .AND. mj0(ijseed) == mj1(ijseed)) zchk = tmask(mi0(iiseed),mj0(ijseed),1)
+      CALL mpp_max('domisf',zchk)
+      IF (zchk(1) == 0._wp) CALL ctl_stop( 'STOP', 'open sea seed is on land, please update namelist (rn_lon_opnsea,rn_lat_opnsea)' )
+      !
+   END SUBROUTINE zgr_isf_subgl
+
+
+   SUBROUTINE zps_isf
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zps_isf  ***
+      !!   
+      !! ** Purpose : compute top partial cell
+      !!   
+      !! ** Method  : same method as for the bottom adapted for the top
+      !!   
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER :: jj, ji   ! loop indices
+      INTEGER :: ik       ! local top level
+      INTEGER :: it       ! counter index
+      !
+      REAL(wp) :: zdiff
+      !!----------------------------------------------------------------------
+      !
+      ! compute top partial cell
+      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 Bu 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  ) = 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)
+               gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik)
+               gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,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
+            ENDIF
+         END DO
+      END DO
+      !
+   END SUBROUTINE zps_isf
+
+   SUBROUTINE zps_isf_e3uv_w
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zps_isf  ***
+      !!   
+      !! ** Purpose :   correct e3uw and e3vw for case with 2 wet cell in the water column under an ice shelf
+      !!   
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER :: ji , jj   ! loop indices
+      INTEGER :: ikt, ikb  ! local top/bottom level
+      !!----------------------------------------------------------------------
+      !
+      ! define e3uw (adapted for 2 cells in the water column)
+      DO jj = 2, jpjm1
+         DO ji = 2, 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 SUBROUTINE zps_isf_e3uv_w
+
+END MODULE domisf
Index: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/dommsk.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_domcfg/src/dommsk.F90	(revision 11201)
+++ NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/dommsk.F90	(revision 11308)
@@ -24,4 +24,5 @@
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
+   USE domisf         ! domain: ice shelf
    USE domwri         ! domain: write the meshmask file
    USE bdy_oce        ! open boundary
@@ -142,6 +143,8 @@
             END DO
          END DO
-      END DO
-     
+      END DO    
+ 
+      IF ( ln_isfsubgl ) CALL zgr_isf_subgl
+
 !SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read !
 !!gm I don't understand why...  
Index: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domzgr.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90	(revision 11201)
+++ NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domzgr.F90	(revision 11308)
@@ -17,5 +17,5 @@
    !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
    !!            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.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capability
    !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying
    !!----------------------------------------------------------------------
@@ -43,4 +43,5 @@
    USE lib_fortran
    USE dombat
+   USE domisf
 
    IMPLICIT NONE
@@ -48,5 +49,5 @@
 
    PUBLIC   dom_zgr        ! called by dom_init.F90
-
+   !
    !                              !!* Namelist namzgr_sco *
    LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
@@ -57,8 +58,4 @@
    REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1)
    REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates
-   INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file
-   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters)
-   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps
-   INTEGER, PUBLIC  ::   nperio          !: type of lateral boundary condition
 
    ! Song and Haidvogel 1994 stretching parameters
@@ -148,4 +145,5 @@
       IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
       !
+      IF ( ln_isfcav ) CALL zgr_isf_nam
       ioptio = 0
       IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
@@ -160,8 +158,8 @@
       IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
       IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
+                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
       !
       ! final adjustment of mbathy & check 
       ! -----------------------------------
-                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
                           CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
                           CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
@@ -313,8 +311,4 @@
       IF ( ln_isfcav .OR. ln_e3_dep ) THEN      ! e3. = dk[gdep]   
          !
-!==>>>   need to be like this to compute the pressure gradient with ISF. 
-!        If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
-!        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
-!
          DO jk = 1, jpkm1
             e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
@@ -578,12 +572,4 @@
                CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
                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
             !       
@@ -624,5 +610,5 @@
          ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
          ENDIF
-         zhmin = gdepw_1d(ik+1)                                                         ! minimum depth = ik+1 w-levels 
+         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels 
          WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
          ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
@@ -931,4 +917,8 @@
       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
 
+      ! compute position of the ice shelf grounding line
+      ! set bathy and isfdraft to 0 where grounded
+      IF ( ln_isfcav ) CALL zgr_isf_zspace
+
       ! bathymetry in level (from bathy_meter)
       ! ===================
@@ -948,4 +938,9 @@
       END DO
 
+      ! Check compatibility between bathy and iceshelf draft
+      ! insure at least 2 wet level on the vertical under an ice shelf
+      ! compute misfdep and adjust isf draft if needed
+      IF ( ln_isfcav ) CALL zgr_isf_kspace
+
       ! Scale factors and depth at T- and W-points
       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
@@ -956,731 +951,5 @@
       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, 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, 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('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
-      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', 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, 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('domzgr', 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
-      !
-      DEALLOCATE( zprt )
-      !
-   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), ALLOCATABLE, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
-      INTEGER , ALLOCATABLE, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
-      !!---------------------------------------------------------------------
-      !
-      ALLOCATE( zbathy(jpi,jpj), zmask(jpi,jpj), zrisfdep(jpi,jpj) )
-      ALLOCATE( zmisfdep(jpi,jpj), zmbathy(jpi,jpj) )
-      !
-      ! (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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk('domzgr', bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:)  = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy,  'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-
-            CALL lbc_lnk( 'domzgr',risfdep, 'T', 1. ) 
-            CALL lbc_lnk( 'domzgr',bathy,   'T', 1. )
-
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy,  'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-
-            CALL lbc_lnk( 'domzgr',risfdep, 'T', 1. ) 
-            CALL lbc_lnk( 'domzgr',bathy,   'T', 1. )
-
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy,  'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-
-            CALL lbc_lnk('domzgr', risfdep, 'T', 1. ) 
-            CALL lbc_lnk('domzgr', bathy,   'T', 1. )
-
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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('domzgr', zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-
-            CALL lbc_lnk('domzgr', risfdep,'T', 1. ) 
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
-            CALL lbc_lnk('domzgr', bathy,  'T', 1. )
-
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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( 'domzgr',zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-
-            CALL lbc_lnk( 'domzgr',risfdep,'T', 1. ) 
-            CALL lbc_lnk( 'domzgr',bathy,  'T', 1. )
-
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( 'domzgr',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
@@ -1704,5 +973,4 @@
                   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
@@ -1710,9 +978,12 @@
                      &                             * ((gdept_1d(     ik  ) - gdepw_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)
+                  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
@@ -1740,57 +1011,94 @@
       END DO
       !
-      ! (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) 
-!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  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
-               gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik)
-               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
-               e3w_0  (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_1d(ik-2)
-               gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,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 
-            ENDIF 
-         END DO 
-      END DO 
-
-      DEALLOCATE( zbathy, zmask, zrisfdep )
-      DEALLOCATE( zmisfdep, zmbathy )
-      !
-   END SUBROUTINE zgr_isf
-
+      ! compute top scale factor if ice shelf
+      IF (ln_isfcav) CALL zps_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, 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
+
+      ! update e3uw in case only 2 cells in the water column
+      IF ( ln_isfcav ) CALL zps_isf_e3uv_w
+      !
+      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
+      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', 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, 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('domzgr', 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   ! to check if it can be simplify
+         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
+      !
+      DEALLOCATE( zprt )
+      !
+   END SUBROUTINE zgr_zps
 
    SUBROUTINE zgr_sco
