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 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4666)
@@ -175,4 +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 
 
    !! All coordinates
@@ -250,7 +251,12 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask              !: land/ocean mask of barotropic stream function
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                              !: land/ocean mask of barotropic stream function
+
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   micedep                 !: top first ocean level                (ISF)
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icedep, iceh            !: Iceshelf draft, iceshef freeboard    (ISF)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lmask                   !: surface domain T-point mask 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
@@ -379,10 +385,16 @@
          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) )
 
-      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     &
-         &     tmask_i(jpi,jpj) , bmask(jpi,jpj) ,                     &
+      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      &
+         &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), &
+         &     bmask(jpi,jpj)   ,                                                       &
          &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )
 
+! (ISF) Allocation of basic array   
+      ALLOCATE( micedep(jpi,jpj) , icedep(jpi,jpj), iceh(jpi,jpj),      &
+         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           &
+         &     mikf(jpi,jpj), lmask(jpi,jpj), STAT=ierr(10) )
+
       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     & 
-         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(10) )
+         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) )
 
 #if defined key_noslip_accurate
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 4666)
@@ -111,6 +111,6 @@
       END DO
       !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
+      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
 
                              CALL dom_stp      ! time step
@@ -159,5 +159,5 @@
       READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namrun )
+      WRITE ( numond, namrun )
       !
       IF(lwp) THEN                  ! control print
@@ -241,5 +241,5 @@
       READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namdom )
+      WRITE ( numond, namdom )
 
       IF(lwp) THEN
@@ -303,5 +303,5 @@
       READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
 906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namcla )
+      WRITE( numond, namcla )
 
       IF(lwp) THEN
@@ -327,5 +327,5 @@
       READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
 908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namnc4 )
+      WRITE( numond, namnc4 )
 
       IF(lwp) THEN                        ! control print
@@ -365,24 +365,24 @@
       !
       IF(lk_mpp) THEN
-         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
-         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
-         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
-         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
+         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
+         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
+         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
       ELSE
-         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )    
-         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )    
-         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )    
-         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )    
-
-         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
+         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )    
+         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )    
+         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )    
+         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )    
+
+         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
          iimi1 = iloc(1) + nimpp - 1
          ijmi1 = iloc(2) + njmpp - 1
-         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
+         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
          iimi2 = iloc(1) + nimpp - 1
          ijmi2 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
+         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
          iima1 = iloc(1) + nimpp - 1
          ijma1 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
+         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
          iima2 = iloc(1) + nimpp - 1
          ijma2 = iloc(2) + njmpp - 1
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 4666)
@@ -638,4 +638,31 @@
       CALL iom_close( inum )
       
+! need to be define for the extended grid south of -80S
+! some point are undefined but you need to have e1 and e2 .NE. 0
+      WHERE (e1t==0.0_wp)
+         e1t=1.0e2
+      END WHERE
+      WHERE (e1v==0.0_wp)
+         e1v=1.0e2
+      END WHERE
+      WHERE (e1u==0.0_wp)
+         e1u=1.0e2
+      END WHERE
+      WHERE (e1f==0.0_wp)
+         e1f=1.0e2
+      END WHERE
+      WHERE (e2t==0.0_wp)
+         e2t=1.0e2
+      END WHERE
+      WHERE (e2v==0.0_wp)
+         e2v=1.0e2
+      END WHERE
+      WHERE (e2u==0.0_wp)
+         e2u=1.0e2
+      END WHERE
+      WHERE (e2f==0.0_wp)
+         e2f=1.0e2
+      END WHERE
+       
     END SUBROUTINE hgr_read
     
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 4666)
@@ -184,4 +184,17 @@
          END DO  
       END DO  
+      
+      ! (ISF) define barotropic mask and mask the ice shelf point
+      lmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
+      
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( REAL( micedep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
+                  tmask(ji,jj,jk) = 0._wp
+               END IF
+            END DO  
+         END DO  
+      END DO  
 
 !!gm  ????
@@ -207,5 +220,5 @@
       ! Interior domain mask (used for global sum)
       ! --------------------
-      tmask_i(:,:) = tmask(:,:,1)
+      tmask_i(:,:) = lmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
       iif = jpreci                         ! ???
       iil = nlci - jpreci + 1
@@ -250,12 +263,26 @@
          END DO
       END DO
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector loop
+            umask_i(ji,jj)  = lmask(ji,jj) * lmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            vmask_i(ji,jj)  = lmask(ji,jj) * lmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+         END DO
+         DO ji = 1, jpim1      ! NO vector opt.
+            fmask_i(ji,jj) =  lmask(ji,jj  ) * lmask(ji+1,jj  )   &
+               &            * lmask(ji,jj+1) * lmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
+         END DO
+      END DO
       CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
       CALL lbc_lnk( vmask, 'V', 1._wp )
       CALL lbc_lnk( fmask, 'F', 1._wp )
+      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( vmask_i, 'V', 1._wp )
+      CALL lbc_lnk( fmask_i, 'F', 1._wp )
 
 
       ! 4. ocean/land mask for the elliptic equation
       ! --------------------------------------------
-      bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point
+      bmask(:,:) = lmask(:,:)       ! elliptic equation is written at t-point
       !
       !                               ! Boundary conditions
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 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4666)
@@ -169,10 +169,21 @@
       fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
       fsdepw_b(:,:,1) = 0.0_wp
-      DO jk = 2, jpk
-         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
-         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
-         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
-         fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk)
-         fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1)
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            DO jk = 1,mikt(ji,jj)-1
+               fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+               fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
+               fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk)
+               fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk)
+            END DO
+            DO jk = mikt(ji,jj), jpk
+               fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+               fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)
+               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
+            END DO
+         END DO
       END DO
 
@@ -185,6 +196,6 @@
          hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
       END DO
-      hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
-      hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
+      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
+      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
 
       ! Restoring frequencies for z_tilde coordinate
@@ -293,5 +304,5 @@
       !                                                ! --------------------------------------------- !
 
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask_i(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask_i(:,:) )
       DO jk = 1, jpkm1
          ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
@@ -313,5 +324,5 @@
             zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
+         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
 
          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
@@ -338,5 +349,5 @@
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          END IF
@@ -463,5 +474,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk)
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -531,6 +542,6 @@
       END DO
       !                                        ! Inverse of the local depth
-      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
-      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
+      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
 
       CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
@@ -830,5 +841,5 @@
                DO jk=1,jpk
                   fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
-                      &                            / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) &
+                      &                            / ( ht_0(:,:) + 1._wp - tmask_i(:,:) ) * tmask(:,:,jk) &
                       &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
                END DO
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 4666)
@@ -132,14 +132,34 @@
       
       CALL dom_uniq( zprw, 'T' )
-      zprt = tmask(:,:,1) * zprw                               !    ! unique point mask
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            jk=mikt(ji,jj) 
+            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+         END DO
+      END DO                             !    ! unique point mask
       CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'U' )
-      zprt = umask(:,:,1) * zprw
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            jk=miku(ji,jj) 
+            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+         END DO
+      END DO
       CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'V' )
-      zprt = vmask(:,:,1) * zprw
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            jk=mikv(ji,jj) 
+            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+         END DO
+      END DO
       CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'F' )
-      zprt = fmask(:,:,1) * zprw
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            jk=mikf(ji,jj) 
+            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+         END DO
+      END DO
       CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 )  
 
@@ -168,6 +188,10 @@
       
       ! note that mbkt is set to 1 over land ==> use surface tmask
-      zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp )
+      zprt(:,:) = lmask(:,:) * REAL( mbkt(:,:) , wp )
       CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
+      zprt(:,:) = lmask(:,:) * REAL( mikt(:,:) , wp )
+      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points
+      zprt(:,:) = lmask(:,:) * REAL( icedep(:,:) , wp )
+      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points
             
       IF( ln_sco ) THEN                                         ! s-coordinate
@@ -203,6 +227,6 @@
             DO jj = 1,jpj   
                DO ji = 1,jpi
-                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
-                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
+                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * lmask(ji,jj)
+                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * lmask(ji,jj)
                END DO
             END DO
@@ -228,6 +252,6 @@
             DO jj = 1,jpj   
                DO ji = 1,jpi
-                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
-                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
+                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * lmask(ji,jj)
+                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * lmask(ji,jj)
                END DO
             END DO
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 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4666)
@@ -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
@@ -102,4 +103,6 @@
       !
       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
       !!----------------------------------------------------------------------
       !
@@ -145,4 +148,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 +298,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
+! 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!
@@ -451,4 +467,26 @@
          END DO
          !
+         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
+         IF( cp_cfg == "isomip" ) THEN 
+	        ! 
+	        icedep(:,:)=200.e0 
+	        micedep(:,:)=1 
+	        ij0 = 1 ; ij1 = 40 
+	        DO jj = mj0(ij0), mj1(ij1) 
+	           icedep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
+                END DO 
+            WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp 
+	        ! 
+         ELSEIF ( cp_cfg == "isomip2" ) THEN
+         ! 
+            icedep(:,:)=0.e0
+            micedep(:,:)=1
+            ij0 = 1 ; ij1 = 40
+            DO jj = mj0(ij0), mj1(ij1)
+               icedep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
+            END DO
+            WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp
+         END IF
+         !
          !                                            ! ================ !
       ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
@@ -492,5 +530,14 @@
             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close( inum )
-            !                                                
+            !  
+            icedep(:,:)=0._wp         
+            micedep(:,:)=1             
+            IF ( nn_isf == 1 .OR. nn_isf == 4 ) THEN
+               CALL iom_open ( 'isf_draft_meter.nc', inum ) 
+               CALL iom_get  ( inum, jpdom_data, 'isf_draft', icedep )
+               CALL iom_close( inum )
+               WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp
+            END IF
+            !       
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
                !
@@ -783,4 +830,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 micedep 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( micedep(:,:) , 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 +953,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   ! (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) ::   eps=0.99            ! small offset to avoid large pool in case bathy slightly greater than icedep
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zbathy, zmask   ! 3D workspace (ISH)
+      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmicedep   ! 3D workspace (ISH)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
       !!---------------------------------------------------------------------
@@ -875,4 +972,6 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zprt )
+      CALL wrk_alloc( jpi, jpj, zbathy, zmask)
+      CALL wrk_alloc( jpi, jpj, zmbathy, zmicedep)
       !
       IF(lwp) WRITE(numout,*)
@@ -906,4 +1005,260 @@
          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
       END DO
+      ! (ISF) compute micedep
+      WHERE( icedep(:,:) == 0._wp ) ;   micedep(:,:) = 1   ! no ice shelf : set micedep to 1  
+      ELSEWHERE                     ;   micedep(:,:) = 2   ! iceshelf : initialize micedep to second level 
+      END WHERE  
+
+      ! Compute micedep 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 < icedep(:,:) .AND. icedep(:,:) >= zdepth )   micedep(:,:) = jk+1 
+      END DO 
+      WHERE (icedep <= e3t_1d(1) .AND. icedep .GT. 0._wp)
+         icedep = 0.
+         micedep= 1
+      END WHERE
+	 
+
+      icompt = 0 
+      DO jl = 1, 10 
+        IF( lk_mpp ) THEN
+           zbathy(:,:) = FLOAT( micedep(:,:) )
+           CALL lbc_lnk( zbathy, 'T', 1. )
+           micedep(:,:) = INT( zbathy(:,:) )
+           CALL lbc_lnk( icedep, '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 
+            micedep( 1 ,:) = micedep(jpim1,:)           ! local domain is cyclic east-west 
+            micedep(jpi,:) = micedep(  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) 
+            icedep = 0._wp
+            micedep= 0
+            bathy  = 0._wp
+         ENDWHERE
+
+         WHERE (bathy(:,:) < icedep(:,:)+eps)
+            micedep(:,:) = 0 
+            icedep(:,:)  = 0._wp
+            mbathy(:,:)  = 0
+            bathy(:,:)   = 0._wp
+         END WHERE
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF (bathy(ji,jj) .GT. icedep(ji,jj) .AND. mbathy(ji,jj) .LT. micedep(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
+               END IF
+            END DO
+         END DO
+
+
+         ! At least 2 levels for water thickness at T, U, and V point.
+        zmicedep(:,:)=micedep(:,:)
+        zmbathy (:,:)=mbathy (:,:)
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF( zmicedep(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 )
+               ENDIF
+               IF( zmicedep(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 )
+               ENDIF
+               IF( zmicedep(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 )
+               ENDIF
+               IF( zmicedep(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
+               IF( zmicedep(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
+            ENDDO
+         ENDDO
+  	IF( lk_mpp ) THEN 
+	   zbathy(:,:) = FLOAT( micedep(:,:) ) 
+ 	   CALL lbc_lnk( zbathy, 'T', 1. ) 
+ 	   micedep(:,:) = INT( zbathy(:,:) ) 
+	   CALL lbc_lnk( icedep, 'T', 1. ) 
+           CALL lbc_lnk( bathy, 'T', 1. )
+           zbathy(:,:) = FLOAT( mbathy(:,:) )
+           CALL lbc_lnk( zbathy, 'T', 1. )
+           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
+                  icedep(ji,jj)=0._wp
+                  micedep(ji,jj)=0
+               END IF
+            END DO
+         END DO
+  	IF( lk_mpp ) THEN 
+	   zbathy(:,:) = FLOAT( micedep(:,:) ) 
+ 	   CALL lbc_lnk( zbathy, 'T', 1. ) 
+ 	   micedep(:,:) = INT( zbathy(:,:) ) 
+	   CALL lbc_lnk( icedep, 'T', 1. ) 
+           CALL lbc_lnk( bathy, 'T', 1. )
+           zbathy(:,:) = FLOAT( mbathy(:,:) )
+           CALL lbc_lnk( zbathy, 'T', 1. )
+           mbathy(:,:) = INT( zbathy(:,:) )
+	ENDIF 
+
+         ! if single point on isf coast line
+         DO jk = 1, jpk
+         WHERE (micedep==0) micedep=jpk
+         zmask=0
+         WHERE (micedep .LE. jk) zmask=1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF (micedep(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
+                     icedep(ji,jj)=gdepw_1d(jk+1) ; micedep(ji,jj)=jk+1
+                     IF (micedep(ji,jj) .GT. mbathy(ji,jj)) micedep(ji,jj) = jpk
+                !     bathy(ji,jj)=0.  ; mbathy(ji,jj)=0
+                  !END IF
+               END IF
+               END IF
+            END DO
+         END DO
+         WHERE (micedep==jpk) 
+             micedep=0 ; icedep=0._wp ; mbathy=0 ; bathy=0._wp
+         END WHERE
+
+  	IF( lk_mpp ) THEN 
+	   zbathy(:,:) = FLOAT( micedep(:,:) ) 
+ 	   CALL lbc_lnk( zbathy, 'T', 1. ) 
+ 	   micedep(:,:) = INT( zbathy(:,:) ) 
+	   CALL lbc_lnk( icedep, '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 (micedep==0) micedep=jpk
+        zmicedep(:,:)=micedep(:,:)
+        zmbathy (:,:)=mbathy (:,:)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+                  ibtest=MIN(zmicedep(ji-1,jj), zmicedep(ji+1,jj), zmicedep(ji,jj-1), zmicedep(ji,jj+1))
+                  IF( ibtest > zmicedep(ji,jj)) THEN
+                     micedep(ji,jj) = ibtest
+                     icedep(ji,jj)  = gdepw_1d(ibtest)
+                  ENDIF
+            ENDDO
+         ENDDO
+         WHERE (micedep==jpk) 
+             micedep=0 ; icedep=0. ; mbathy=0 ; bathy=0
+         END WHERE
+         IF( lk_mpp ) THEN 
+            zbathy(:,:) = FLOAT( micedep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            micedep(:,:) = INT( zbathy(:,:) ) 
+            CALL lbc_lnk( icedep, '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
+        zmicedep(:,:)=micedep(:,:)
+        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
+                  mbathy(ji,jj) = ibtest
+                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
+               ENDIF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:) = FLOAT( micedep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            micedep(:,:) = INT( zbathy(:,:) ) 
+            CALL lbc_lnk( icedep, 'T', 1. ) 
+            CALL lbc_lnk( bathy, 'T', 1. )
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+        ENDIF 
+        ! remove pool of water stuck between ice shelf and bathymetry
+        DO jk = 1, jpk
+        WHERE (micedep==0) micedep=jpk
+        zmicedep(:,:)=micedep(:,:)
+        zmbathy (:,:)=mbathy (:,:)
+        DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF( jk .GE. zmicedep(ji,jj) .AND. jk .LE. zmbathy(ji,jj) ) THEN
+               IF( (jk > zmbathy(ji,jj+1) .OR. jk < zmicedep(ji,jj+1)) .AND. &
+                 & (jk > zmbathy(ji,jj-1) .OR. jk < zmicedep(ji,jj-1)) .AND. &
+                 & (jk > zmbathy(ji+1,jj) .OR. jk < zmicedep(ji+1,jj)) .AND. &
+                 & (jk > zmbathy(ji-1,jj) .OR. jk < zmicedep(ji-1,jj))       ) THEN
+                  mbathy(ji,jj) = 0 ; micedep(ji,jj) = jpk ; icedep(ji,jj) = 0._wp ; bathy(ji,jj) = 0._wp
+               ENDIF
+               ENDIF
+            ENDDO
+        ENDDO
+        WHERE (micedep==jpk) micedep=0 
+        IF( lk_mpp ) THEN 
+            zbathy(:,:) = FLOAT( micedep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            micedep(:,:) = INT( zbathy(:,:) ) 
+            CALL lbc_lnk( icedep, 'T', 1. ) 
+            CALL lbc_lnk( bathy, 'T', 1. )
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+        ENDIF 
+        ENDDO
+      END DO 
+
+      WHERE (mbathy(:,:) < micedep(:,:))
+         micedep(:,:) = 0
+         icedep(:,:)  = 0._wp
+         mbathy(:,:)  = 0
+         bathy(:,:)   = 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 +1293,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 +1328,52 @@
          END DO
       END DO
+      !
+      ! (ISF) Definition of e3t, u, v, w for ISF case
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+            ik = micedep(ji,jj) 
+            IF( ik > 1 ) THEN               ! ice shelf point only 
+                IF( icedep(ji,jj) < gdepw_1d(ik) )  icedep(ji,jj)= gdepw_1d(ik) 
+                gdepw_0(ji,jj,ik) = icedep(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 
+                   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 = micedep(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,*) ' icedep = ', icedep(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 +1394,16 @@
          END DO
       END DO
+      ! (ISF) define e3uw
+      DO jk = 2,jpk                          
+         DO jj = 1, jpjm1 
+            DO ji = 1, fs_jpim1   ! vector opt. 
+! (ISF)** NEEDS CHANGING TO SECOND OPTION FOR ICE SHELF BUT WILL CHANGE RESULTS WITHOUT ICE (DIFFER AT THE 1e-13 LEVEL) 
+               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 +1447,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 (micedep == 0) micedep = 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, micedep(ji,jj)
+               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+            END DO
+            IF (micedep(ji,jj) .GE. 2) gdep3w_0(ji,jj,micedep(ji,jj)) = icedep(ji,jj) + 0.5_wp * e3w_0(ji,jj,micedep(ji,jj))
+            DO jk = micedep(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 +1489,6 @@
       !
       CALL wrk_dealloc( jpi, jpj, jpk, zprt )
+      CALL wrk_dealloc( jpi, jpj, zmask, zbathy )
+      CALL wrk_dealloc( jpi, jpj, zmicedep, zmbathy )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 4666)
@@ -263,4 +263,10 @@
                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
                   ENDIF
+                  ik = mikt(ji,jj)
+                  IF( ik > 1 ) THEN
+                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
+                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
+                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
+                  END IF
                END DO
             END DO
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 4666)
@@ -76,5 +76,5 @@
       !
 
-      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' '
       IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
@@ -110,4 +110,9 @@
          ELSEIF( cp_cfg == 'gyre' ) THEN         
             CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
+        ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN
+            IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain' 
+            tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields 
+            tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:)
+            tsb(:,:,:,:)=tsn(:,:,:,:)  
          ELSE                                    ! Initial T-S, U-V fields read in files
             IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
@@ -129,6 +134,7 @@
          CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities
 #if ! defined key_c1d
-         IF( ln_zps )   CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  & ! zps: before hor. gradient
-            &                                       rhd, gru , grv  )   ! of t,s,rd at ocean bottom
+         IF( ln_zps )    CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient
+            &                                      rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             !
+            &                                      gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level
 #endif
          !   
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 4666)
@@ -50,5 +50,5 @@
    REAL(wp), PUBLIC ::   rau0     = 1026._wp         !: volumic mass of reference     [kg/m3]
 #else
-   REAL(wp), PUBLIC ::   rau0     = 1035._wp         !: volumic mass of reference     [kg/m3]
+   REAL(wp), PUBLIC ::   rau0     = 1028.4_wp       !: volumic mass of reference     [kg/m3]
 #endif
    REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg]
