Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Biblio/Biblio.bib	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Biblio/Biblio.bib	(revision 5111)
@@ -271,4 +271,13 @@
   volume = {326},
   pages = {677--684}
+}
+
+@ARTICLE{Beckmann2003,
+  author = {A. Beckmann and H. Goosse},
+  title = {A parameterization of ice shelf-ocean interaction for climate models},
+  journal = OM
+  year = {2003}
+  volume = {5}
+  pages = {157--170}
 }
 
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 5111)
@@ -493,5 +493,6 @@
 $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step 
 bathymetry or $s$-coordinate (hybrid and partial step coordinates have not 
-yet been tested in NEMO v2.3). 
+yet been tested in NEMO v2.3). If using $z$-coordinate with partial step bathymetry
+(\np{ln\_zps}~=~true), ocean cavity beneath ice shelves can be open (\np{ln\_isfcav}~=~true).
 
 Contrary to the horizontal grid, the vertical grid is computed in the code and no 
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 5111)
@@ -627,5 +627,8 @@
 \eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of 
 the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 
-($e_{3w}$). 
+($e_{3w}$).
+ 
+$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}=true).
+This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}=true).
 
 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true)
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 5111)
@@ -1,6 +1,6 @@
 % ================================================================
-% Chapter � Surface Boundary Condition (SBC, ICB) 
-% ================================================================
-\chapter{Surface Boundary Condition (SBC, ICB) }
+% Chapter � Surface Boundary Condition (SBC, ISF, ICB) 
+% ================================================================
+\chapter{Surface Boundary Condition (SBC, ISF, ICB) }
 \label{SBC}
 \minitoc
@@ -48,6 +48,8 @@
 below ice-covered areas (using observed ice-cover or a sea-ice model) 
 (\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater 
-fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of a freshwater flux adjustment 
-in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the 
+fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of isf melting as lateral inflow (parametrisation) 
+(\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false) or as surface flux at the land-ice ocean interface
+(\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true); 
+the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the 
 transformation of the solar radiation (if provided as daily mean) into a diurnal 
 cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave 
@@ -60,4 +62,6 @@
 Finally, the different options that further modify the fluxes applied to the ocean are discussed.
 One of these is modification by icebergs (see \S\ref{ICB_icebergs}), which act as drifting sources of fresh water.
+An other exemple of modification is the ice shelf melting/freezing (see \S\ref{SBC_isf}), 
+which act as sources of fresh water due to ice shelf melting/freezing.
 
 
@@ -945,6 +949,52 @@
 
 %}
-
-
+% ================================================================
+%        Ice shelf melting
+% ================================================================
+\section   [Ice shelf melting (\textit{sbcisf})]
+                        {Ice shelf melting (\mdl{sbcisf})}
+\label{SBC_isf}
+%------------------------------------------namsbc_isf----------------------------------------------------
+\namdisplay{namsbc_isf}
+%--------------------------------------------------------------------------------------------------------
+Namelist variable in \ngn{namsbc}, \np{nn\_isf},  control the kind of ice shelf representation used. 
+\begin{description}
+\item[\np{nn\_isf}~=~1]
+The ice shelf cavity is represented. The fwf and heat flux are computed. 
+Full description, sensitivity and validation in preparation. 
+
+\item[\np{nn\_isf}~=~2]
+A parametrisation of isf is used. The ice shelf cavity is not represented. 
+The fwf is distributed along the ice shelf edge between the depth of the average grounding line (GL)
+(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}) as in (\np{nn\_isf}~=~3). 
+Furthermore the fwf is computed using the \citet{Beckmann2003} parametrisation of isf melting. 
+The effective melting length (\np{sn\_Leff\_isf}) is read in a file.
+
+\item[\np{nn\_isf}~=~3]
+A simple parametrisation of isf is used. The ice shelf cavity is not represented. 
+The fwf (\np{sn\_rnfisf}) is distributed along the ice shelf edge between the depth of the average grounding line (GL)
+(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}).
+Full description, sensitivity and validation in preparation.
+
+\item[\np{nn\_isf}~=~4]
+The ice shelf cavity is represented. However, the fwf (\np{sn\_fwfisf}) and heat flux (\np{sn\_qisf}) are 
+not computed but specified in file. 
+\end{description}
+
+\np{nn\_isf}~=~1 and \np{nn\_isf}~=~2 compute a melt rate based on the water masse properties, ocean velocities and depth.
+ This flux is thus highly dependent of the model resolution (horizontal and vertical), realism of the water masse onto the shelf ...
+
+\np{nn\_isf}~=~3 and \np{nn\_isf}~=~4 read the melt rate and heat flux in a file. You have total control of the fwf scenario.
+
+ This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too 
+coarse to have realistic melting or for sensitivity studies where you want to control your input. 
+Full description, sensitivity and validation in preparation. 
+
+There is 2 ways to apply the fwf to NEMO. The first possibility (\np{ln\_divisf}~=~false) applied the fwf
+ and heat flux directly on the salinity and temperature tendancy. The second possibility (\np{ln\_divisf}~=~true)
+ apply the fwf as for the runoff fwf (see \S\ref{SBC_rnf}). The mass/volume addition due to the ice shelf melting is,
+ at each relevant depth level, added to the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div} 
+(called from \mdl{divcur}). 
+%
 % ================================================================
 %        Handling of icebergs
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 5111)
@@ -830,5 +830,5 @@
 % Bottom Friction
 % ================================================================
-\section  [Bottom Friction (\textit{zdfbfr})]   {Bottom Friction (\mdl{zdfbfr} module)}
+\section  [Bottom and top Friction (\textit{zdfbfr})]   {Bottom Friction (\mdl{zdfbfr} module)}
 \label{ZDF_bfr}
 
@@ -837,5 +837,8 @@
 %--------------------------------------------------------------------------------------------------------------
 
-Options are defined through the  \ngn{nambfr} namelist variables.
+Options to defined the top and bottom friction are defined through the  \ngn{nambfr} namelist variables.
+The top friction is activated only if the ice shelf cavities are opened (\np{ln\_isfcav}~=~true).
+As the friction processes at the top and bottom are the same, only the bottom friction is described in details.
+
 Both the surface momentum flux (wind stress) and the bottom momentum 
 flux (bottom friction) enter the equations as a condition on the vertical 
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/nambfr
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/nambfr	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/nambfr	(revision 5111)
@@ -5,9 +5,19 @@
                            !                              = 2 : nonlinear friction
    rn_bfri1    =    4.e-4  !  bottom drag coefficient (linear case)
-   rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case)
+   rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
+   rn_bfri2_max =   1.e-1  !  max. bottom drag coefficient (non linear case and ln_loglayer=T)
    rn_bfeb2    =    2.5e-3 !  bottom turbulent kinetic energy background  (m2/s2)
-   rn_bfrz0    =    3.e-3  ! bottom roughness for loglayer bfr coeff 
+   rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T
    ln_bfr2d    = .false.   !  horizontal variation of the bottom friction coef (read a 2D mask file )
    rn_bfrien   =    50.    !  local multiplying factor of bfr (ln_bfr2d=T)
+   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case)
+   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
+   rn_tfri2_max =   1.e-1  !  max. top drag coefficient (non linear case and ln_loglayer=T)
+   rn_tfeb2    =    0.0    !  top turbulent kinetic energy background  (m2/s2)
+   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T
+   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file )
+   rn_tfrien   =    50.    !  local multiplying factor of tfr (ln_tfr2d=T)
+
    ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true)
+   ln_loglayer = .false.   !  logarithmic formulation (non linear case)
 /
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namdyn_hpg
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namdyn_hpg	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namdyn_hpg	(revision 5111)
@@ -5,4 +5,5 @@
    ln_hpg_zps  = .true.    !  z-coordinate - partial steps (interpolation)
    ln_hpg_sco  = .false.   !  s-coordinate (standard jacobian formulation)
+   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to ice shelf cavity
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc	(revision 5111)
@@ -19,4 +19,9 @@
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave
    ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf)
+   nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf)
+                           !  0 =no isf                  1 = presence of ISF
+                           !  2 = bg03 parametrisation   3 = rnf file for isf
+                           !  4 = ISF fwf specified
+                           !  option 1 and 4 need ln_isfcav = .true. (domzgr)
    ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
    nn_fwb      = 3         !  FreshWater Budget: =0 unchecked
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc_isf
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc_isf	(revision 5111)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc_isf	(revision 5111)
@@ -0,0 +1,30 @@
+!-----------------------------------------------------------------------
+&namsbc_isf    !  Top boundary layer (ISF)
+!-----------------------------------------------------------------------
+!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation !
+!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  !
+! nn_isf == 4
+   sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+! nn_isf == 3
+   sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+! nn_isf == 2 and 3
+   sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+! nn_isf == 2
+   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+! for all case
+   ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?)
+! only for nn_isf = 1 or 2
+   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula
+   rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula
+! only for nn_isf = 1
+   nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??)
+   rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008)
+                          ! 0 => thickness of the tbl = thickness of the first wet cell
+   ln_conserve = .true.   ! conservative case (take into account meltwater advection)
+   nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s)
+                          ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
+                          !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate
+                          ! 2 = velocity and stability dependent Gamma    Holland et al. 1999
+/
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namzgr
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namzgr	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namzgr	(revision 5111)
@@ -5,3 +5,4 @@
    ln_zps      = .true.    !  z-coordinate - partial steps   (T/F)
    ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F)
+   ln_isfcav   = .false.   !  ice shelf cavity               (T/F)
 /
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5111)
@@ -83,5 +83,5 @@
    ln_zps      = .true.    !  z-coordinate - partial steps   (T/F)
    ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F)
-   ln_isfcav   = .false.   !  ice shelf cavity
+   ln_isfcav   = .false.   !  ice shelf cavity               (T/F)
 /
 !-----------------------------------------------------------------------
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 5111)
@@ -105,9 +105,13 @@
       END DO
       IF( .NOT.lk_vvl ) THEN
-         DO ji=1,jpi
-            DO jj=1,jpj
-               zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
-            END DO
-         END DO
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+               END DO
+            END DO
+         ELSE
+            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END IF
       END IF
       !                                         
@@ -127,9 +131,13 @@
       END DO
       IF( .NOT.lk_vvl ) THEN
-         DO ji=1,jpi
-            DO jj=1,jpj
-               zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
-            END DO
-         END DO
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+               END DO
+            END DO
+         ELSE
+            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END IF
       END IF
       !    
@@ -157,10 +165,15 @@
       END DO
       IF( .NOT.lk_vvl ) THEN
-         DO ji=1,jpi
-            DO jj=1,jpj
-               ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
-               zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
-            END DO
-         END DO
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
+                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
+               END DO
+            END DO
+         ELSE
+            ztemp = ztemp + zarea_ssh(:,:) * tsn(:,:,1,jp_tem) 
+            zsal  = zsal  + zarea_ssh(:,:) * tsn(:,:,1,jp_sal) 
+         END IF
       ENDIF
       IF( lk_mpp ) THEN  
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5111)
@@ -96,8 +96,8 @@
       z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                               ! heat fluxes
       z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                               ! salt fluxes
-      ! Add runoff heat & salt input
+      ! Add runoff    heat & salt input
       IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) )
       IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
-      ! Add geothermal ice shelf
+      ! Add ice shelf heat & salt input
       IF( nn_isf .GE. 1 )  THEN
           z_frc_trd_t = z_frc_trd_t &
@@ -112,11 +112,15 @@
       !
       IF( .NOT. lk_vvl ) THEN
-         z2d0=0.0_wp ; z2d1=0.0_wp
-         DO ji=1,jpi
-            DO jj=1,jpj
-              z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem)
-              z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal)
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem)
+                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal)
+               ENDDO
             ENDDO
-         ENDDO
+         ELSE
+            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem)
+            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal)
+         END IF
          z_wn_trd_t = - glob_sum( z2d0 ) 
          z_wn_trd_s = - glob_sum( z2d1 )
@@ -144,11 +148,15 @@
       ! heat & salt content variation (associated with ssh)
       IF( .NOT. lk_vvl ) THEN
-         z2d0 = 0._wp   ;   z2d1 = 0._wp
-         DO ji = 1, jpi
-            DO jj = 1, jpj
-              z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
-              z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+         IF ( ln_isfcav ) THEN
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
+                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+               END DO
             END DO
-         END DO
+         ELSE
+            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 
+            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
+         END IF
          z_ssh_hc = glob_sum( z2d0 ) 
          z_ssh_sc = glob_sum( z2d1 ) 
@@ -277,10 +285,15 @@
           frc_s = 0._wp                                           ! salt content   -    -   -    -        
           IF( .NOT. lk_vvl ) THEN
-             DO ji=1,jpi
-                DO jj=1,jpj
-                   ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
-                   ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+             IF ( ln_isfcav ) THEN
+                DO ji=1,jpi
+                   DO jj=1,jpj
+                      ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
+                      ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                   ENDDO
                 ENDDO
-             ENDDO
+             ELSE
+                ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
+                ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+             END IF
              frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
              frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5111)
@@ -298,15 +298,17 @@
       ENDIF
 
+!      IF ( ln_isfcav ) THEN
 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
-      DO jk = 1, jpkm1
-         e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
-      END DO
-      e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
-
-      DO jk = 2, jpk
-         e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
-      END DO
-      e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
+         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)) 
+!      END IF
 
 !!gm BUG in s-coordinate this does not work!
@@ -473,13 +475,12 @@
          ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
          IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 
-	        ! 
-	        risfdep(:,:)=200.e0 
-	        misfdep(:,:)=1 
-	        ij0 = 1 ; ij1 = 40 
-	        DO jj = mj0(ij0), mj1(ij1) 
-	           risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
-                END DO 
+            risfdep(:,:)=200.e0 
+            misfdep(:,:)=1 
+            ij0 = 1 ; ij1 = 40 
+            DO jj = mj0(ij0), mj1(ij1) 
+               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
+            END DO 
             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
-	        ! 
+         ! 
          ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
          ! 
@@ -536,5 +537,5 @@
             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close( inum )
-            !  
+            !                                                
             risfdep(:,:)=0._wp         
             misfdep(:,:)=1             
@@ -963,11 +964,10 @@
       !!----------------------------------------------------------------------
       !!
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
       INTEGER  ::   ik, it           ! temporary integers
-      INTEGER  ::   id, jd, nprocd
       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, zmin       ! Maximum and minimum depth
+      REAL(wp) ::   zmax             ! Maximum depth
       REAL(wp) ::   zdiff            ! temporary scalar
       REAL(wp) ::   zrefdep          ! temporary scalar
@@ -1046,6 +1046,6 @@
                      &                             * ((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) ) 
+                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
@@ -1202,5 +1202,5 @@
      
       ! Compute gdep3w_0 (vertical sum of e3w)
-      IF ( ln_isfcav ) THEN
+      IF ( ln_isfcav ) THEN ! if cavity
          WHERE (misfdep == 0) misfdep = 1
          DO jj = 1,jpj
@@ -1335,10 +1335,10 @@
             misfdep(jpi,:) = misfdep(  2  ,:) 
          ENDIF
- 
+
          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
             mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
             mbathy(jpi,:) = mbathy(  2  ,:)
          ENDIF
- 
+
          ! split last cell if possible (only where water column is 2 cell or less)
          DO jk = jpkm1, 1, -1
@@ -1358,5 +1358,5 @@
             END WHERE
          END DO
- 
+
  
  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
@@ -1639,5 +1639,5 @@
                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == 0 ) THEN
+               IF( ibtest == 0 .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
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5111)
@@ -140,5 +140,5 @@
             &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
             &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level
-         IF( ln_zps .AND. ln_isfcav)                                       &
+         IF( ln_zps .AND.       ln_isfcav)                                 &
             &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
             &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5111)
@@ -127,4 +127,6 @@
       IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
           CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
+      IF( ln_dynzad_zts .AND. ln_isfcav )   &
+          CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' )
 
       !                               ! Set nadv
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5111)
@@ -167,4 +167,7 @@
                            & the standard jacobian formulation hpg_sco or &
                            & the pressure jacobian formulation hpg_prj')
+
+      IF( ln_hpg_isf .AND. .NOT. ln_isfcav )   &
+           &   CALL ctl_stop( ' hpg_isf not available is ln_isfcav = false ' )
       !
       !                               ! Set nhpg from ln_hpg_... flags
@@ -187,4 +190,7 @@
       IF( ( .NOT. ln_hpg_isf ) .AND. ln_isfcav )   &
           &  CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
+      ! 
+      ! initialisation of ice load
+      riceload(:,:)=0.0
       !
    END SUBROUTINE dyn_hpg_init
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5111)
@@ -90,17 +90,30 @@
          DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
             DO ji = fs_2, fs_jpim1        ! vector opt.
-               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) ) !* wumask(ji,jj,jk)
-               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) ) !* wvmask(ji,jj,jk)
+               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
+               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
             END DO  
          END DO   
       END DO
-      DO jj = 2, jpjm1              ! Surface and bottom values set to zero
-         DO ji = fs_2, fs_jpim1           ! vector opt.
-            zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
-            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
-            zwuw(ji,jj,jpk) = 0._wp
-            zwvw(ji,jj,jpk) = 0._wp
-         END DO  
-      END DO
+      !
+      ! Surface and bottom advective fluxes set to zero
+      IF ( ln_isfcav ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1           ! vector opt.
+               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
+               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
+               zwuw(ji,jj,jpk) = 0._wp
+               zwvw(ji,jj,jpk) = 0._wp
+            END DO
+         END DO
+      ELSE
+         DO jj = 2, jpjm1        
+            DO ji = fs_2, fs_jpim1           ! vector opt.
+               zwuw(ji,jj, 1 ) = 0._wp
+               zwvw(ji,jj, 1 ) = 0._wp
+               zwuw(ji,jj,jpk) = 0._wp
+               zwvw(ji,jj,jpk) = 0._wp
+            END DO  
+         END DO
+      END IF
 
       DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
@@ -196,9 +209,10 @@
          END DO
       END DO
-
-      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
+      !
+      ! Surface and bottom advective fluxes set to zero
+      DO jj = 2, jpjm1        
          DO ji = fs_2, fs_jpim1           ! vector opt.
-            zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
-            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
+            zwuw(ji,jj, 1 ) = 0._wp
+            zwvw(ji,jj, 1 ) = 0._wp
             zwuw(ji,jj,jpk) = 0._wp
             zwvw(ji,jj,jpk) = 0._wp
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5111)
@@ -157,6 +157,4 @@
          !
          !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
-         ! surface initialisation 
-         zdzr(:,:,1) = 0._wp 
          ! interior value
          DO jk = 2, jpkm1
@@ -169,6 +167,8 @@
                &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
          END DO
+         ! surface initialisation 
+         zdzr(:,:,1) = 0._wp 
          IF ( ln_isfcav ) THEN
-         ! if isf need to overwrite the interior value at at the first ocean point
+            ! if isf need to overwrite the interior value at at the first ocean point
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -185,14 +185,23 @@
          ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
          !
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji  ,jj)
-               IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji+1,jj)
-               IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj), hmlpt(ji+1,jj))
-               IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji  ,jj)
-               IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji,jj+1)
-               IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji,jj+1))
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
+                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp)
+                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp)
+                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
+                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp)
+                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp)
+               ENDDO
             ENDDO
-         ENDDO
+         ELSE
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
+                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
+               ENDDO
+            ENDDO
+         END IF
          DO jk = 2, jpkm1                            !* Slopes at u and v points
             DO jj = 2, jpjm1
@@ -208,20 +217,17 @@
                   zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  )
                   !                                      ! uslp and vslp output in zwz and zww, resp.
-                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
-                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
+                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj  ,jk) )
+                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) )
                   ! thickness of water column between surface and level k at u/v point
-                  zdepu = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                   &
-                             - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) )  &
-                             - fse3u(ji,jj,miku(ji,jj))                                         )
-                  zdepv = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                   &
-                             - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) &
-                             - fse3v(ji,jj,mikv(ji,jj))                                         )
-                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                              &
-                     &                 + zfi  * uslpml(ji,jj)                                                     &
-                     &                        * zdepu / MAX( zhmlpu(ji,jj), 5._wp )
+                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              &
+                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) )
+                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              &
+                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) )
+                  !
+                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          &
+                     &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj)
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk)
-                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                              &
-                     &                 + zfj  * vslpml(ji,jj)                                                     &
-                     &                        * zdepv / MAX( zhmlpv(ji,jj), 5._wp ) 
+                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          &
+                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj) 
                   zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk)
                   
@@ -276,8 +282,8 @@
                   uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
                      &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   &
-                     &                            *   umask(ji,jj,jk-1) !* umask(ji,jj,jk) * umask(ji,jj,jk+1)
+                     &                            *   umask(ji,jj,jk-1)
                   vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   &
                      &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   &
-                     &                            *   vmask(ji,jj,jk-1) !* vmask(ji,jj,jk) * vmask(ji,jj,jk+1)
+                     &                            *   vmask(ji,jj,jk-1)
                END DO
             END DO
@@ -292,5 +298,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   !                                  !* Local vertical density gradient evaluated from N^2
-                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk)
                   !                                  !* Slopes at w point
                   !                                        ! i- & j-gradient of density at w-points
@@ -308,5 +314,5 @@
                   zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  )
                   !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
-                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )    ! zfk=1 in the ML otherwise zfk=0
+                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
                   zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp )
                   zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) &
@@ -746,6 +752,8 @@
             DO ji = 1, jpi
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
-               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
+               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN
+                  omlmask(ji,jj,jk) = 1._wp
+               ELSE
+                  omlmask(ji,jj,jk) = 0._wp
                ENDIF
             END DO
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5111)
@@ -89,6 +89,9 @@
          !
          IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
-         !
-         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
+         IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
+         !
+         area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface
+         ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
+         ! and in case of no melt, it can generate HSSW.
          !
 #if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice
@@ -107,6 +110,6 @@
             z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
-            emp(:,:) = emp(:,:) - z_fwf 
-            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1)
+            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
          ENDIF
          !
@@ -139,6 +142,6 @@
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
             zcoef = fwfold * rcp
-            emp(:,:) = emp(:,:) + fwfold
-            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:)  ! account for change to the heat budget due to fw correction
+            emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1)
+            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
          ENDIF
          !
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 5111)
@@ -61,20 +61,16 @@
       !!---------------------------------------------------------------------
       
-      !                                        !* first wet T-, U-, V- ocean level (ISF) variables (T, S, depth, velocity)
+      !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zub(ji,jj)        = ub (ji,jj,miku(ji,jj))
-            zvb(ji,jj)        = vb (ji,jj,mikv(ji,jj))
             zts(ji,jj,jp_tem) = tsn(ji,jj,mikt(ji,jj),jp_tem)
             zts(ji,jj,jp_sal) = tsn(ji,jj,mikt(ji,jj),jp_sal)
          END DO
       END DO
+      zub(:,:)        = ub (:,:,1       )
+      zvb(:,:)        = vb (:,:,1       )
       !
       IF( lk_vvl ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zdep(ji,jj) = fse3t_n(ji,jj,mikt(ji,jj))
-            END DO
-         END DO
+         zdep(:,:) = fse3t_n(:,:,1)
       ENDIF
       !                                                   ! ---------------------------------------- !
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5111)
@@ -206,5 +206,5 @@
       IF( lk_esopa         )   ioptio =          1
 
-      IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck ) .AND. nn_isf .NE. 0 )  &
+      IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts ) .AND. ln_isfcav )  &
       &   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity')
 
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 5111)
@@ -147,15 +147,23 @@
          ! Surface value
          IF( lk_vvl ) THEN   
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable
-               END DO
-            END DO
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable
+                  END DO
+               END DO
+            ELSE
+               zwz(:,:,1) = 0.e0          ! volume variable
+            END IF
          ELSE                
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
-               END DO
-            END DO   
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                  END DO
+               END DO   
+            ELSE
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
+            END IF
          ENDIF
 
@@ -212,9 +220,13 @@
          END DO
          ! surface value
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zwz(ji,jj,mikt(ji,jj)) = 0.e0
-            END DO
-         END DO
+         IF ( ln_isfcav ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwz(ji,jj,mikt(ji,jj)) = 0.e0
+               END DO
+            END DO
+         ELSE
+            zwz(:,:,1) = 0.e0
+         END IF
          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
          CALL lbc_lnk( zwz, 'W',  1. )
@@ -360,8 +372,4 @@
 
          ! upstream tracer flux in the k direction
-         ! Surface value
-         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable
-         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
-         ENDIF
          ! Interior value
          DO jk = 2, jpkm1
@@ -374,4 +382,26 @@
             END DO
          END DO
+         ! Surface value
+         IF( lk_vvl ) THEN
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf
+                  END DO
+               END DO
+            ELSE
+               zwz(:,:,1) = 0.e0                              ! volume variable + no isf
+            END IF
+         ELSE
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf
+                  END DO
+               END DO
+            ELSE
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf
+            END IF
+         ENDIF
 
          ! total advective trend
@@ -583,9 +613,9 @@
 
       DO jk = 1, jpkm1
+         ikm1 = MAX(jk-1,1)
+         z2dtt = p2dt(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ikm1 = MAX(jk-1,1)
-               z2dtt = p2dt(jk)
-               
+
                ! search maximum in neighbourhood
                zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 5111)
@@ -290,4 +290,9 @@
       IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0
 
+      ! Initialisation of gtui/gtvi in case of no cavity
+      IF ( .NOT. ln_isfcav ) THEN
+         gtui(:,:,:) = 0.0_wp
+         gtvi(:,:,:) = 0.0_wp
+      END IF
       !                                        ! T & S profile (to be coded +namelist parameter
 
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5111)
@@ -187,5 +187,5 @@
             END DO
          END DO
-            ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
+         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
          zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2)
          zdkt (:,:,1) = zdk1t(:,:,1)
@@ -200,6 +200,6 @@
          END IF
 
-            ! 2. Horizontal fluxes
-            ! --------------------   
+         ! 2. Horizontal fluxes
+         ! --------------------   
          DO jk = 1, jpkm1
             DO jj = 1 , jpjm1
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5111)
@@ -82,5 +82,5 @@
       !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
       !!
-      !! ** Action  : compute for top and bottom interfaces
+      !! ** Action  : compute for top interfaces
       !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points
       !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points
@@ -91,9 +91,9 @@
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts 
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom)
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom)
       !
       INTEGER  ::   ji, jj, jn      ! Dummy loop indices
-      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points
-      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars
+      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points
+      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
       REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             ! 
@@ -172,5 +172,5 @@
             END DO
          END DO
-         
+
          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
          ! step and store it in  zri, zrj for each  case
@@ -193,5 +193,5 @@
             END DO
          END DO
-         CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions
+         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
          !
       END IF
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5111)
@@ -120,13 +120,18 @@
                   zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
                   zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max)
-! (ISF)
-                  ikbt = mikt(ji,jj)
-! JC: possible WAD implementation should modify line below if layers vanish
-                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
-                  ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
-                  ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max)
-
                END DO
             END DO
+! (ISF)
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ikbt = mikt(ji,jj)
+! JC: possible WAD implementation should modify line below if layers vanish
+                     ztmp = (1-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
+                     ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                     ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max)
+                  END DO
+               END DO
+            END IF
          !   
          ELSE
@@ -152,42 +157,50 @@
                !
                ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
-               IF ( miku(ji,jj) + 2 .GE. 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 .GE. 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 * (1._wp - vmask(ji,jj,1))
-               END IF
-               ! (ISF) ========================================================================
-               ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points 
-               ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points)
-               !
-               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
-                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
-               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
-                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
-               !
-               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 )
-               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 )
-               !
-               tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1))
-               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 .GE. 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 .GE. 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))
+               IF ( ln_isfcav ) THEN
+                  IF ( miku(ji,jj) + 1 .GE. 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) + 1 .GE. 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 * (1._wp - vmask(ji,jj,1))
+                  END IF
                END IF
             END DO
          END DO
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  ! (ISF) ========================================================================
+                  ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points 
+                  ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points)
+                  !
+                  zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
+                     &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
+                  zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
+                     &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
+              !
+                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 )
+                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 )
+              !
+                  tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1))
+                  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) + 1 .GE. 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) + 1 .GE. 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
+         END IF
          !
          CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 5111)
@@ -14,5 +14,4 @@
    !!----------------------------------------------------------------------
    USE par_oce         ! mesh and scale factors
-   USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test)
    USE ldftra_oce      ! ocean active tracers: lateral physics
    USE ldfdyn_oce      ! ocean dynamics lateral physics
@@ -118,5 +117,5 @@
       IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   &
          &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' )
-      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0 )   &
+      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav )   &
          &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' )
       !
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 5111)
@@ -336,6 +336,6 @@
                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
                   &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
-                  &           / (  fse3uw_n(ji,jj,jk)         &
-                  &              * fse3uw_b(ji,jj,jk) )
+                  &                            / (  fse3uw_n(ji,jj,jk)               &
+                  &                              *  fse3uw_b(ji,jj,jk)  )
                avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
                   &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
Index: /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 5110)
+++ /branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 5111)
@@ -143,7 +143,7 @@
  
       tra(:,:,:,:) = 0._wp
-      IF( ln_zps .AND. .NOT. lk_c1d )   &              ! Partial steps: before horizontal gradient of passive
+      IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav )   &              ! Partial steps: before horizontal gradient of passive
         &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient
-      IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) &
+      IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   &
         &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level
 
