Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4810)
+++ branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5034)
@@ -35,4 +35,5 @@
    USE oce               ! ocean variables
    USE dom_oce           ! ocean domain
+   USE sbc_oce           ! surface variable (isf)
    USE closea            ! closed seas
    USE c1d               ! 1D vertical configuration
@@ -101,5 +102,5 @@
       INTEGER ::   ios
       !
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
       !!----------------------------------------------------------------------
       !
@@ -120,7 +121,8 @@
          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,*) '             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
       ENDIF
 
@@ -145,4 +147,5 @@
       IF( .NOT.lk_c1d )   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
       !
       IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
@@ -294,4 +297,16 @@
          gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
       ENDIF
+
+! 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) 
+      END DO
+      e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
+
+      DO jk = 2, jpk
+         e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
+      END DO
+      e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
 
 !!gm BUG in s-coordinate this does not work!
@@ -450,4 +465,28 @@
             END DO
          END DO
+         risfdep(:,:)=0.e0
+         misfdep(:,:)=1
+         !
+         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
+         IF( cp_cfg == "isomip" ) THEN 
+	        ! 
+	        risfdep(:,:)=200.e0 
+	        misfdep(:,:)=1 
+	        ij0 = 1 ; ij1 = 40 
+	        DO jj = mj0(ij0), mj1(ij1) 
+	           risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
+                END DO 
+            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
+	        ! 
+         ELSEIF ( cp_cfg == "isomip2" ) THEN
+         ! 
+            risfdep(:,:)=0.e0
+            misfdep(:,:)=1
+            ij0 = 1 ; ij1 = 40
+            DO jj = mj0(ij0), mj1(ij1)
+               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
+            END DO
+            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
+         END IF
          !
          !                                            ! ================ !
@@ -492,5 +531,14 @@
             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close( inum )
-            !                                                
+            !  
+            risfdep(:,:)=0._wp         
+            misfdep(:,:)=1             
+            IF ( ln_isfcav ) THEN
+               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
+               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
+               CALL iom_close( inum )
+               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
+            END IF
+            !       
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
                !
@@ -530,4 +578,9 @@
       !                       
       IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
+         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
+         WHERE (bathy == risfdep)
+            bathy   = 0.0_wp ; risfdep = 0.0_wp
+         END WHERE
+         ! 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
@@ -783,4 +836,49 @@
    END SUBROUTINE zgr_bot_level
 
+      SUBROUTINE zgr_top_level
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_bot_level  ***
+      !!
+      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
+      !!
+      !! ** Method  :   computes from misfdep with a minimum value of 1
+      !!
+      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest 
+      !!                                     ocean level at t-, u- & v-points
+      !!                                     (min value = 1)
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER ::   ji, jj   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level')
+      !
+      CALL wrk_alloc( jpi, jpj, zmik )
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
+      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
+      !
+      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1)
+      !                                      ! top k-index of W-level (=mikt)
+      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level
+         DO ji = 1, jpim1
+            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
+            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
+            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  )
+         END DO
+      END DO
+
+      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
+      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 )
+      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 )
+      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 )
+      !
+      CALL wrk_dealloc( jpi, jpj, zmik )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level')
+      !
+   END SUBROUTINE zgr_top_level
 
    SUBROUTINE zgr_zco
@@ -861,12 +959,17 @@
       !!----------------------------------------------------------------------
       !!
-      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      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             ! Maximum depth
+      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
       REAL(wp) ::   zdiff            ! temporary scalar
       REAL(wp) ::   zrefdep          ! temporary scalar
+      REAL(wp) ::   zbathydiff, zrisfdepdiff 
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH)
+      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
       !!---------------------------------------------------------------------
@@ -875,4 +978,6 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zprt )
+      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
+      CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep)
       !
       IF(lwp) WRITE(numout,*)
@@ -889,5 +994,4 @@
       ENDIF
 
-
       ! bathymetry in level (from bathy_meter)
       ! ===================
@@ -906,4 +1010,466 @@
          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
       END DO
+      ! (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 = 1, 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 ) 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 
 
       ! Scale factors and depth at T- and W-points
@@ -938,7 +1504,7 @@
 !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) ))
+                  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) ) 
@@ -973,4 +1539,52 @@
          END DO
       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  ) = 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)
+            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 (ISF)
 
       ! Scale factors and depth at U-, V-, UW and VW-points
@@ -991,4 +1605,17 @@
          END DO
       END DO
+      ! (ISF) define e3uw
+      DO jk = 2,jpk                          
+         DO jj = 1, jpjm1 
+            DO ji = 1, fs_jpim1   ! vector opt. 
+               e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) &
+                 &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) )
+               e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) &
+                 &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) )
+            END DO 
+         END DO 
+      END DO
+      !End (ISF)
+      
       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 )
@@ -1032,9 +1659,17 @@
      
       ! Compute gdep3w_0 (vertical sum of e3w)
-      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
-        
+      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
       !                                               ! ================= !
       IF(lwp .AND. ll_print) THEN                     !   Control print   !
@@ -1066,4 +1701,6 @@
       !
       CALL wrk_dealloc( jpi, jpj, jpk, zprt )
+      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
+      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
@@ -1445,6 +2082,6 @@
             DO jk = 1, jpkm1
                IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
-               IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0
-            END DO
+            END DO
+            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
          END DO
       END DO
