Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 4747)
@@ -193,10 +193,11 @@
          CALL iom_put( "sss2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface salinity
       END IF
+      ! multiply by umask to prevent not numerical value error in the ioserver sometimes
       IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN
-         CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport
-         CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport
+         CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport
+         CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport
       ELSE
-         CALL iom_put( "uoce" , un                         )    ! i-current
-         CALL iom_put( "voce" , vn                         )    ! j-current
+         CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:)                  )    ! i-current
+         CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:)                  )    ! j-current
          DO jj = 1, jpj
             DO ji = 1, jpi
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4747)
@@ -175,5 +175,5 @@
    LOGICAL, PUBLIC ::   ln_zps        !: z-coordinate - partial step
    LOGICAL, PUBLIC ::   ln_sco        !: s-coordinate or hybrid z-s coordinate
-   LOGICAL, PUBLIC ::   ln_isf        !: presence of ISF 
+   LOGICAL, PUBLIC ::   ln_isfcav     !: presence of ISF 
 
    !! All coordinates
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4747)
@@ -849,4 +849,9 @@
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+               ! needed to restart if land processor not computed 
+               WHERE ( tmask(:,:,:) == 0.0_wp ) 
+                  fse3t_n(:,:,:) = e3t_0(:,:,:)
+                  fse3t_b(:,:,:) = e3t_0(:,:,:)
+               END WHERE
                IF( neuler == 0 ) THEN
                   fse3t_b(:,:,:) = fse3t_n(:,:,:)
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4747)
@@ -102,7 +102,5 @@
       INTEGER ::   ios
       !
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr, nn_fwb, ln_cdgw, nn_isf
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
       !!----------------------------------------------------------------------
       !
@@ -123,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
 
@@ -533,5 +532,5 @@
             risfdep(:,:)=0._wp         
             misfdep(:,:)=1             
-            IF ( nn_isf == 1 .OR. nn_isf == 4 ) THEN
+            IF ( ln_isfcav ) THEN
                CALL iom_open ( 'isf_draft_meter.nc', inum ) 
                CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep )
@@ -577,4 +576,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
@@ -956,5 +960,5 @@
       INTEGER  ::   ik, it           ! temporary integers
       INTEGER  ::   id, jd, nprocd
-      INTEGER  ::   icompt, ibtest   ! (ISF)
+      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
@@ -963,6 +967,6 @@
       REAL(wp) ::   zdiff            ! temporary scalar
       REAL(wp) ::   zrefdep          ! temporary scalar
-      REAL(wp) ::   eps=0.99            ! small offset to avoid large pool in case bathy slightly greater than risfdep
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zbathy, zmask   ! 3D workspace (ISH)
+      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
@@ -972,5 +976,5 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zprt )
-      CALL wrk_alloc( jpi, jpj, zbathy, zmask)
+      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
       CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep)
       !
@@ -988,5 +992,4 @@
       ENDIF
 
-
       ! bathymetry in level (from bathy_meter)
       ! ===================
@@ -1006,6 +1009,6 @@
       END DO
       ! (ISF) compute misfdep
-      WHERE( risfdep(:,:) == 0._wp ) ;   misfdep(:,:) = 1   ! no ice shelf : set misfdep to 1  
-      ELSEWHERE                     ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
+      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  
 
@@ -1018,7 +1021,6 @@
          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
+      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
+         risfdep(:,:) = 0. ; misfdep(:,:) = 1
       END WHERE
  
@@ -1026,5 +1028,13 @@
       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 
+      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(:,:) )
@@ -1041,66 +1051,115 @@
             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
-
-         WHERE (mbathy == 0) 
-            risfdep = 0._wp
-            misfdep= 0
-            bathy  = 0._wp
-         ENDWHERE
-
-         WHERE (bathy(:,:) < risfdep(:,:)+eps)
-            misfdep(:,:) = 0 
-            risfdep(:,:)  = 0._wp
-            mbathy(:,:)  = 0
-            bathy(:,:)   = 0._wp
-         END WHERE
-! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
+ 
+         ! 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
-               IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
-                  bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
-                  mbathy(ji,jj)= mbathy(ji,jj) + 1
+               ! 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.
-         zmisfdep(:,:)=misfdep(:,:)
-         zmbathy (:,:)=mbathy (:,:)
  
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ! T point
-               IF( zmisfdep(ji,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN
-                  mbathy(ji,jj) = zmbathy(ji,jj) + 1 
-                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat )
+          ! 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
-               ! V point
-               IF( zmisfdep(ji,jj+1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN
-                  mbathy(ji,jj) = zmbathy(ji,jj) + 1 
-                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat )
+            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
-               ! V point -1
-               IF( zmisfdep(ji,jj-1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN
-                  mbathy(ji,jj) = zmbathy(ji,jj) + 1 
-                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat )
+            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
-               ! U point
-               IF( zmisfdep(ji+1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN
-                  mbathy(ji,jj) = zmbathy(ji,jj) + 1 
-                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat )
-               ENDIF
-               ! U point -1
-               IF( zmisfdep(ji-1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN
-                  mbathy(ji,jj) = zmbathy(ji,jj) + 1 
-                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat )
-               ENDIF
-            END DO
-         END DO
+            END DO
+         END DO
+ 
+ 
          IF( lk_mpp ) THEN 
             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
@@ -1113,19 +1172,22 @@
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF 
-
-        ! if single ocean point put as land
-         zmask=1
-         WHERE (mbathy .EQ. 0) zmask=0
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               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)=0._wp
-                  mbathy(ji,jj)=0
-                  risfdep(ji,jj)=0._wp
-                  misfdep(ji,jj)=0
-               END IF
-            END DO
-         END DO
+ 
+ ! 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(:,:) ) 
@@ -1137,7 +1199,39 @@
             CALL lbc_lnk( zbathy, 'T', 1. )
             mbathy(:,:) = INT( zbathy(:,:) )
-         END IF 
-
-         ! if single point on isf coast line
+         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
@@ -1155,36 +1249,72 @@
                END DO
             END DO
-            WHERE (misfdep==jpk) 
-               misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
-            END WHERE
-  
-            IF( lk_mpp ) THEN 
-               zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-               CALL lbc_lnk( zbathy, 'T', 1. ) 
-               misfdep(:,:) = INT( zbathy(:,:) ) 
-               CALL lbc_lnk( risfdep, 'T', 1. ) 
-               CALL lbc_lnk( bathy, 'T', 1. )
-               zbathy(:,:) = FLOAT( mbathy(:,:) )
-               CALL lbc_lnk( zbathy, 'T', 1. )
-               mbathy(:,:) = INT( zbathy(:,:) )
-            ENDIF 
-         END DO
-
-! fill hole in ice shelf
-         WHERE (misfdep==0) misfdep=jpk
-         zmisfdep(:,:)=misfdep(:,:)
-         zmbathy (:,:)=mbathy (:,:)
+         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
-               ibtest=MIN(zmisfdep(ji-1,jj), zmisfdep(ji+1,jj), zmisfdep(ji,jj-1), zmisfdep(ji,jj+1))
-               IF( ibtest > zmisfdep(ji,jj)) THEN
+               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)
+                  risfdep(ji,jj) = gdepw_1d(ibtest)
                ENDIF
             ENDDO
          ENDDO
-         WHERE (misfdep==jpk) 
-            misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0
-         END WHERE
+ 
          IF( lk_mpp ) THEN 
             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
@@ -1197,13 +1327,20 @@
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF 
-
-! fill hole in bathymetry
-         zmisfdep(:,:)=misfdep(:,:)
+ !
+ !! fill hole in bathymetry
          zmbathy (:,:)=mbathy (:,:)
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               ibtest = MAX(  zmbathy(ji-1,jj), zmbathy(ji+1,jj),   &
-                  &           zmbathy(ji,jj-1), zmbathy(ji,jj+1)  )
-               IF( ibtest < zmbathy(ji,jj) ) THEN
+               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) 
@@ -1221,42 +1358,94 @@
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF 
-        ! remove 1 cell pool of water stuck between ice shelf and bathymetry (need a 3D flood filling tools to do this properly)
-         DO jk = 1, jpk
-            WHERE (misfdep==0) misfdep=jpk
-            zmisfdep(:,:)=misfdep(:,:)
-            zmbathy (:,:)=mbathy (:,:)
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF ( jk .GE. zmisfdep(ji,jj) .AND. jk .LE. zmbathy(ji,jj) ) THEN
-                     IF ( (jk > zmbathy(ji,jj+1) .OR. jk < zmisfdep(ji,jj+1)) .AND. &
-                        & (jk > zmbathy(ji,jj-1) .OR. jk < zmisfdep(ji,jj-1)) .AND. &
-                        & (jk > zmbathy(ji+1,jj) .OR. jk < zmisfdep(ji+1,jj)) .AND. &
-                        & (jk > zmbathy(ji-1,jj) .OR. jk < zmisfdep(ji-1,jj))       ) THEN
-                        mbathy(ji,jj) = 0 ; misfdep(ji,jj) = jpk ; risfdep(ji,jj) = 0._wp ; bathy(ji,jj) = 0._wp
-                     ENDIF
-                  ENDIF
-               ENDDO
-            ENDDO
-            WHERE (misfdep==jpk) misfdep=0 
-            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
+ ! 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 
-
-      WHERE (mbathy(:,:) < misfdep(:,:))
-         misfdep(:,:) = 0
-         risfdep(:,:)  = 0._wp
-         mbathy(:,:)  = 0
-         bathy(:,:)   = 0._wp
+! 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 
@@ -1338,22 +1527,22 @@
             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) 
+               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) 
+               !       ... on ik 
+               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
+                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
+                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
+               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
+               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
+
+               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
+                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
                 ENDIF 
-                !       ... on ik / ik-1 
-                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
-                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
+               !       ... on ik / 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)
+               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
             ENDIF 
          END DO 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4747)
@@ -137,5 +137,5 @@
       READ  ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )
-      WRITE ( numond, namdyn_hpg )
+      IF(lwm) WRITE ( numond, namdyn_hpg )
       !
       IF(lwp) THEN                   ! Control print
@@ -397,7 +397,13 @@
       ! iniitialised to 0. zhpi zhpi 
       zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
-      
-      ! assume density of water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
+
+!==================================================================================     
+!=====Compute iceload and contribution of the half first wet layer ================= 
+!===================================================================================
+
+      ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
       ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
+
+      ! compute density of the water displaced by the ice shelf 
       zrhd = rhd ! save rhd
       DO jk = 1, jpk
@@ -434,5 +440,5 @@
             END DO
             IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &
-                               &                              * ( fsdepw(ji,jj,ikt) - gdept_1d(ikt-1) )
+                               &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
          END DO
       END DO
@@ -464,10 +470,12 @@
          END DO
       END DO
-      !
+!==================================================================================     
+!===== Compute partial cell contribution for the top cell ========================= 
+!==================================================================================
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
             iku = miku(ji,jj) ; 
             zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp
-            ze3wu  = fse3w(ji+1,jj,iku+1)-fse3w(ji,jj,iku+1) 
+            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
             ! u direction
             IF ( iku .GT. 1 ) THEN
@@ -478,12 +486,11 @@
                ! corrective term ( = 0 if z coordinate )
                zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)
-               ! zhpi will be added in interior loop and zuapint will be removed in the interior loop
+               ! zhpi will be added in interior loop
                ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap
                ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
                IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku)
 
-               ! zhpi will be added in interior loop and zuapint will be removed in the interior loop
-               ! case iku + 1 (remove the zphi term added in the interior loop)
-               zhpiint        =  zcoef0 / e1u(ji,jj)                                                             &    
+               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
+               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &    
                   &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          &
                   &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   &
@@ -492,7 +499,8 @@
                zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
             END IF
+               
             ! v direction
             ikv = mikv(ji,jj)
-            ze3wv  = fse3w(ji,jj+1,ikv+1)-fse3w(ji,jj,ikv+1) 
+            ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
             IF ( ikv .GT. 1 ) THEN
                ! case ikv
@@ -502,13 +510,13 @@
                ! corrective term ( = 0 if z coordinate )
                zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)
-               ! zhpi will be added in interior loop and zvapint will be removed in the interior loop
+               ! zhpi will be added in interior loop
                va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap
                ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
                IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv) 
-               ! zhpi will be added in interior loop and zvapint will be removed in the interior loop
-               ! case ikv + 1 (remove the zphj term added in the interior loop)
-               zhpjint        =  zcoef0 / e2v(ji,jj)                                                            &
+               
+               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
+               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              &
                   &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         &
-                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)  &
+                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    &
                   &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         &
                   &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  )
@@ -518,5 +526,8 @@
       END DO
 
-      ! interior value (2=<jk=<jpkm1)
+!==================================================================================     
+!===== Compute interior value ===================================================== 
+!==================================================================================
+
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -524,5 +535,5 @@
             DO jk = 2, jpkm1
                ! hydrostatic pressure gradient along s-surfaces
-               ! s-coordinate pressure gradient correction
+               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              &
                   &                            + zcoef0 / e1u(ji,jj)                                                           &
@@ -531,9 +542,12 @@
                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 &
                   &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   ) 
-             !corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
-               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+               ! s-coordinate pressure gradient correction
+               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
+               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    &
                   &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
                ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
 
+               ! hydrostatic pressure gradient along s-surfaces
+               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              &
                   &                            + zcoef0 / e2v(ji,jj)                                                           &
@@ -542,6 +556,7 @@
                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 &
                   &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   )
-             !corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
-               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
+               ! s-coordinate pressure gradient correction
+               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
+               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     &
                   &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
                ! add to the general momentum trend
@@ -550,6 +565,9 @@
          END DO
       END DO
-     
-      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
+
+!==================================================================================     
+!===== Compute bottom cell contribution (partial cell) ============================ 
+!==================================================================================
+
 # if defined key_vectopt_loop
          jj = 1
@@ -589,6 +607,6 @@
 # endif
       END DO
-
-
+     
+      ! set back to original density value into the ice shelf cell (maybe useless because it is masked)
       rhd = zrhd
       !
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4747)
@@ -220,10 +220,4 @@
       ENDIF
       !
-      IF( lk_vvl ) THEN     ! need it for (ISF)
-                     CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
-                     CALL iom_get( numror, jpdom_autoglo, 'fse3t  ', fse3t  (:,:,:) )
-                     CALL iom_get( numror, jpdom_autoglo, 'fse3w  ', fse3w  (:,:,:) )
-                     CALL iom_get( numror, jpdom_autoglo, 'fsdepw ', fsdepw (:,:,:) )
-      END IF
       CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields
       CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn      )
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 4747)
@@ -67,9 +67,9 @@
          imask                                ! temporary global workspace
       REAL(wp), DIMENSION(jpiglo,jpjglo) ::   &
-         zdta                   ! temporary data workspace
+         zdta, zdtaisf                     ! temporary data workspace
       REAL(wp) ::   zidom , zjdom          ! temporary scalars
 
       ! read namelist for ln_zco
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
 
       !!----------------------------------------------------------------------
@@ -109,9 +109,17 @@
       ENDIF
       CALL iom_close (inum)
+      
+      ! used to compute the land processor in case of not masked bathy file.
+      zdtaisf(:,:) = 0.0_wp
+      IF ( ln_isfcav ) THEN
+         CALL iom_open ( 'bathy_meter.nc', inum )   ! Meter bathy in case of partial steps
+         CALL iom_get ( inum, jpdom_unknown, 'isf_draft' , zdtaisf, kstart=(/jpizoom,jpjzoom/), kcount=(/jpiglo,jpjglo/) )
+      END IF
+      CALL iom_close (inum)
 
       ! land/sea mask over the global/zoom domain
 
       imask(:,:)=1
-      WHERE ( zdta(:,:) <= 0. ) imask = 0
+      WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0. ) imask = 0
 
       !  1. Dimension arrays for subdomains
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 4747)
@@ -130,6 +130,9 @@
                ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
                ! (ISF) case partial step top and bottom in adjacent cell in vertical
-               ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku))
-               ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv))
+               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
+               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
+               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
+               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
+               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
                !
                ! i- direction
@@ -189,12 +192,12 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku))
-               ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv))
-
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
+               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
+               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
+
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2
                ENDIF
 # if ! defined key_vectopt_loop
@@ -218,9 +221,9 @@
                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
                ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points
-               ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku))
-               ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv))
+               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
+               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
                IF( ze3wu >= 0._wp ) THEN 
                   pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
-                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 
+                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1 
                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
                                 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) &
@@ -267,7 +270,11 @@
                iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1
                ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1
-               ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1)
-               ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1)
                !
+               ! (ISF) case partial step top and bottom in adjacent cell in vertical
+               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
+               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
+               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
+               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
+               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
@@ -322,12 +329,12 @@
                iku = miku(ji,jj)
                ikv = mikv(ji,jj)
-               ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1)
-               ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1)
-
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
+               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
+               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
+
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
                ENDIF
 # if ! defined key_vectopt_loop
@@ -351,6 +358,6 @@
                iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
                ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1
-               ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1)
-               ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1)
+               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
+               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
                IF( ze3wu >= 0._wp ) THEN
                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 4726)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 4747)
@@ -167,4 +167,15 @@
                bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv
                !
+               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
+               IF ( miku(ji,jj) + 2 .LE. mbku(ji,jj) ) THEN
+                  bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   &
+                               &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) &
+                               &          * zecu * (1._wp - umask(ji,jj,1))
+               END IF
+               IF ( mikv(ji,jj) + 2 .LE. mbkv(ji,jj) ) THEN
+                  bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   &
+                               &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) &
+                               &          * zecv
+               END IF
                ! (ISF) ========================================================================
                ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points 
@@ -182,5 +193,15 @@
                tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1))
                ! (ISF) END ====================================================================
-
+               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
+               IF ( miku(ji,jj) + 2 .LE. mbku(ji,jj) ) THEN
+                  tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   &
+                               &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) &
+                               &          * zecu * (1._wp - umask(ji,jj,1))
+               END IF
+               IF ( mikv(ji,jj) + 2 .LE. mbkv(ji,jj) ) THEN
+                  tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   &
+                               &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) &
+                               &          * zecv * (1._wp - vmask(ji,jj,1))
+               END IF
             END DO
          END DO
