Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5215)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5216)
@@ -969,5 +969,5 @@
       !!
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
-      INTEGER  ::   ik, it           ! temporary integers
+      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
       LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
@@ -1093,4 +1093,5 @@
          e3vw_0(:,:,jk) = e3w_1d(jk)
       END DO
+
       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
          DO jj = 1, jpjm1
@@ -1105,20 +1106,22 @@
       IF ( ln_isfcav ) THEN
       ! (ISF) define e3uw (adapted for 2 cells in the water column)
-      ! Need to test if the modification of only mikt and mbkt levels is enough
-         DO jk = 2,jpk                          
-            DO jj = 1, jpjm1 
-               DO ji = 1, fs_jpim1   ! vector opt. 
-                  e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) &
-                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) )
-                  e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) &
-                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) )
-               END DO 
-            END DO 
+         DO jj = 2, jpjm1 
+            DO ji = 2, fs_jpim1   ! vector opt. 
+               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
+               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
+               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
+               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
+            END DO
          END DO
       END IF
-      
+
       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
       !
+
       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
@@ -1263,4 +1266,16 @@
          risfdep(:,:) = 0. ; misfdep(:,:) = 1
       END WHERE
+
+      ! remove very shallow ice shelf (less than ~ 10m if 75L)
+      IF ( cp_cfg .NE. "isomip" ) THEN
+         WHERE (misfdep(:,:) <= 10 .AND. misfdep(:,:) .GT. 1)
+            misfdep = 0; risfdep = 0.0_wp;
+            mbathy  = 0; bathy   = 0.0_wp;
+         END WHERE
+         WHERE (bathy(:,:) <= 30.0 .AND. gphit < -60)
+            misfdep = 0; risfdep = 0.0_wp;
+            mbathy  = 0; bathy   = 0.0_wp;
+         END WHERE
+      END IF
  
 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
@@ -1787,5 +1802,5 @@
                ENDIF 
             !       ... on ik / ik-1 
-               e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
+               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5215)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5216)
@@ -585,5 +585,5 @@
 
       !! compute ustar
-         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) )
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
 
       !! Compute gammats
@@ -596,5 +596,5 @@
       !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
       !! compute ustar
-         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) )
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
 
       !! compute Pr and Sc number (can be improved)
@@ -678,50 +678,101 @@
 
       REAL(wp) :: ze3, zhk
-      REAL(wp), DIMENSION(:,:), POINTER :: zikt
+      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl
 
       INTEGER :: ji,jj,jk
       INTEGER :: ikt,ikb
-      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
-
-      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
-      CALL wrk_alloc( jpi,jpj, zikt )
+
+      CALL wrk_alloc( jpi,jpj, zhisf_tbl)
 
       ! get first and last level of tbl
-      mkt(:,:) = misfkt(:,:)
-      mkb(:,:) = misfkb(:,:)
 
       varout(:,:)=0._wp
-      DO jj = 2,jpj
-         DO ji = 2,jpi
-            ikt = mkt(ji,jj)
-            ikb = mkb(ji,jj)
+      IF (cptin == 'U') THEN
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
+            ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt))
+
+            ! determine the deepest level influenced by the boundary layer
+            ! test on tmask useless ?????
+               DO jk = ikt, mbku(ji,jj)
+                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
 
             ! level fully include in the ice shelf boundary layer
-            DO jk = ikt, ikb - 1
-               ze3 = fse3t_n(ji,jj,jk)
-               IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
-               IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
-                  &                                                       * r1_hisf_tbl(ji,jj) * ze3
-               IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
-                  &                                                       * r1_hisf_tbl(ji,jj) * ze3
-            END DO
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3u_n(ji,jj,jk)
+                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
 
             ! level partially include in ice shelf boundary layer 
-            zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
-            IF (cptin == 'T') &
-                &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
-            IF (cptin == 'U') &
-                &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
-            IF (cptin == 'V') &
-                &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
-         END DO
-      END DO
+               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj))
+            END DO
+         END DO
+         CALL lbc_lnk(varout,'T',-1.)
+      END IF
+
+      IF (cptin == 'V') THEN
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
+           ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt))
+
+            ! determine the deepest level influenced by the boundary layer
+            ! test on tmask useless ?????
+               DO jk = ikt, mbkv(ji,jj)
+                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+
+            ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3v_n(ji,jj,jk)
+                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+            ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1))
+            END DO
+         END DO
+         CALL lbc_lnk(varout,'T',-1.)
+      END IF
+
+      IF (cptin == 'T') THEN
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = misfkt(ji,jj)
+               ikb = misfkb(ji,jj)
+
+            ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3t_n(ji,jj,jk)
+                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
+               END DO
+
+            ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
+               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+      END IF
       varout(:,:) = varout(:,:) * ssmask(:,:)
 
-      CALL wrk_dealloc( jpi,jpj, mkt, mkb )      
-      CALL wrk_dealloc( jpi,jpj, zikt ) 
-
-      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
-      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
+      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )      
 
    END SUBROUTINE sbc_isf_tbl
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5215)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5216)
@@ -271,6 +271,6 @@
                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 )    ! if level first is a p-step, ik.m1=1
-               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
-               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
+               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku)
+               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv)
                !
                ! i- direction
@@ -319,12 +319,12 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
-               ze3wv = gdept_0(ji,jj+1,ikv) - 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
+               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku)
+               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(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
                ENDIF
 
@@ -343,6 +343,6 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
-               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
+               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku)
+               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv)
 
                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
@@ -370,6 +370,6 @@
                ! 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,jj,iku) - gdept_0(ji+1,jj,iku)
-               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 
+               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku)
+               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 
 
                ! i- direction
@@ -417,13 +417,13 @@
                iku = miku(ji,jj)
                ikv = mikv(ji,jj)
-               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
-               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,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
+               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku)
+               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,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
                ENDIF
 
@@ -442,6 +442,6 @@
                iku = miku(ji,jj) 
                ikv = mikv(ji,jj) 
-               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
-               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 
+               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku)
+               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 
 
                IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
