Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 4724)
@@ -104,5 +104,11 @@
          zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
       END DO
-      IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+      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
+      END IF
       !                                         
       zarho = SUM( area(:,:) * zbotpres(:,:) ) 
@@ -120,5 +126,11 @@
          zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
       END DO
-      IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+      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
+      END IF
       !    
       zarho = SUM( area(:,:) * zbotpres(:,:) ) 
@@ -145,6 +157,10 @@
       END DO
       IF( .NOT.lk_vvl ) THEN
-         ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
-         zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
+         DO ji=1,jpi
+            DO jj=1,jpj
+               ztemp = ztemp + SUM( zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) )
+               zsal  = zsal  + SUM( zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) )
+            END DO
+         END DO
       ENDIF
       IF( lk_mpp ) THEN  
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 4724)
@@ -196,8 +196,8 @@
               DO ji = 1,jpi
                 ! Elevation
-                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)        
+                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask_i(ji,jj)        
 #if defined key_dynspg_ts
-                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
-                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
+                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask_i(ji,jj)
+                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask_i(ji,jj)
 #endif
               END DO
@@ -294,6 +294,6 @@
                X1 = ana_amp(ji,jj,jh,1)
                X2 =-ana_amp(ji,jj,jh,2)
-               out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1)
-               out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1)
+               out_eta(ji,jj,jh       ) = X1 * tmask_i(ji,jj)
+               out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj)
             ENDDO
          ENDDO
@@ -328,6 +328,6 @@
                X1=ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_u(ji,jj,jh) = X1 * umask(ji,jj,1)
-               out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
+               out_u(ji,jj,jh) = X1 * umask_i(ji,jj)
+               out_u (ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
             ENDDO
          ENDDO
@@ -362,6 +362,6 @@
                X1=ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_v(ji,jj,jh)=X1 * vmask(ji,jj,1)
-               out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1)
+               out_v(ji,jj,jh)=X1 * vmask_i(ji,jj)
+               out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
             ENDDO
          ENDDO
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4724)
@@ -69,5 +69,5 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      INTEGER    ::   jk                          ! dummy loop indice
+      INTEGER    ::   jk, ji, jj                          ! dummy loop indice
       REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
       REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! -   -   -   -   -   -   -   - 
@@ -79,16 +79,24 @@
       REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     -
       REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     -
+      REAL(dp), DIMENSION(:,:), POINTER       ::   z2d0, z2d1
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')      
-
+      CALL wrk_alloc( jpi,jpj, z2d0, z2d1 )
+      tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ;
+      tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ;
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:)) * surf(:,:) ) ! volume fluxes
-      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
+      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes
+      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
       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
+      IF( nn_isf .GE. 1 )  THEN
+                        z_frc_trd_t = z_frc_trd_t +                      glob_sum( ( risf_tsc(:,:,jp_tem) - rdivisf * fwfisf(:,:) * -1.9 * r1_rau0 ) * surf(:,:) )
+                        z_frc_trd_s = z_frc_trd_s + (1.0_wp - rdivisf) * glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) )
+      ENDIF
 
       ! Add penetrative solar radiation
@@ -98,6 +106,13 @@
       !
       IF( .NOT. lk_vvl ) THEN
-         z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) )
-         z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) )
+         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)
+            ENDDO
+         ENDDO
+         z_wn_trd_t = - glob_sum( z2d0 ) 
+         z_wn_trd_s = - glob_sum( z2d1 )
       ENDIF
 
@@ -123,6 +138,13 @@
       ! heat & salt content variation (associated with ssh)
       IF( .NOT. lk_vvl ) THEN
-         z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) )
-         z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) )
+         z2d0=0.0_wp ; z2d1=0.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) ) 
+            ENDDO
+         ENDDO
+         z_ssh_hc = glob_sum( z2d0 ) 
+         z_ssh_sc = glob_sum( z2d1 ) 
       ENDIF
 
@@ -282,5 +304,5 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jk   ! 
+     INTEGER ::   jk,ji,jj   ! 
      INTEGER ::   id1   ! local integers
      !!----------------------------------------------------------------------
@@ -322,6 +344,10 @@
           frc_s = 0.d0                                           ! salt content   -    -   -    -        
           IF( .NOT. lk_vvl ) THEN
-             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
+             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
              frc_wn_t = 0.d0                                       ! initial heat content misfit due to free surface
              frc_wn_s = 0.d0                                       ! initial salt content misfit due to free surface
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 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4724)
@@ -312,5 +312,5 @@
       !                                                ! --------------------------------------------- !
 
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask_i(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask_i(:,:) )
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * lmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - lmask(:,:) )
       DO jk = 1, jpkm1
          ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
@@ -472,5 +472,5 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - lmask )
          DO jk = 1, jpkm1
             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
@@ -998,4 +998,5 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
+      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0)   CALL ctl_stop( 'vvl_ztilde, vvl_layer, vvl_ztilde_as_zstar, vvl_zstar_at_eqtor not tested with ice shelf cavity (only vvl_zstar was tested)' )
 
       IF(lwp) THEN                   ! Print the choice
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 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4724)
@@ -299,5 +299,5 @@
       ENDIF
 
-! need to be like this to compute the pressure gradient with ISF
+! 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
@@ -1023,6 +1023,7 @@
       END WHERE
 	 
-
+! basic check for the compatibility of bathy and icedep. I think it should be offline because it is not perfect and cannot solved all the situation
       icompt = 0 
+! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
       DO jl = 1, 10 
         IF( lk_mpp ) THEN
@@ -1058,5 +1059,5 @@
             bathy(:,:)   = 0._wp
          END WHERE
-
+! Case where bathy and icedep compatible but not the level variable mbathy/micedep because of partial cell condition
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -1075,20 +1076,25 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
+               ! T point
                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
+               ! V point
                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
+               ! V point -1
                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
+               ! U point
                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
+               ! U point -1
                IF( zmicedep(ji-1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN
                   mbathy(ji,jj) = zmbathy(ji,jj) + 1 
@@ -1217,5 +1223,5 @@
             mbathy(:,:) = INT( zbathy(:,:) )
         ENDIF 
-        ! remove pool of water stuck between ice shelf and bathymetry
+        ! remove 1 cell pool of water stuck between ice shelf and bathymetry (need a 3D flood filling tools to do this properly)
         DO jk = 1, jpk
         WHERE (micedep==0) micedep=jpk
@@ -1344,5 +1350,5 @@
                 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 
+                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
                    e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
                 ENDIF 
@@ -1398,5 +1404,4 @@
          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) )
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4724)
@@ -29,4 +29,5 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
+   USE sbc_oce         ! surface variable (only for the flag with ice shelf)
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
@@ -176,4 +177,5 @@
       IF( ln_hpg_prj )   ioptio = ioptio + 1
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
+      IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 )   CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity. Comparison in a GYRE simulation with bump in the middle show similar result than hpg_zps' )
       !
    END SUBROUTINE dyn_hpg_init
@@ -395,7 +397,8 @@
       ! iniitialised to 0. zhpi zhpi 
       zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
-
+      
+      ! assume density of water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
       ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
-      zrhd = rhd
+      zrhd = rhd ! save rhd
       DO jk = 1, jpk
            zdept(:,:)=gdept_1d(jk)
@@ -403,10 +406,11 @@
       END DO
       WHERE ( tmask(:,:,:) == 1._wp)
-        rhd(:,:,:) = zrhd(:,:,:)
+        rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd
       END WHERE
       
-
+      ! compute rhd at the ice/oce interface (ice shelf side)
       CALL eos(ztstop,icedep,zrhdtop_isf)
 
+      ! compute rhd at the ice/oce interface (ocean side)
       DO ji=1,jpi
         DO jj=1,jpj
@@ -419,5 +423,5 @@
       !
       ! Surface value + ice shelf gradient
-      ! compute pressure (used to compute hpgi/j for all the level from 1 to miku/v)
+      ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
       ziceload = 0._wp
       DO jj = 1, jpj
@@ -433,10 +437,11 @@
          END DO
       END DO
-      ! compute zp from first level to first wet cell (blue and purple area)
+      riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
+      ! compute zp from z=0 to first T wet point (correction due to zps not yet applied)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1)
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
-            ! we assume ISF is in isostatic equilibrium with a density equal to the reference density
+            ! we assume ISF is in isostatic equilibrium
             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    &
                &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   &
@@ -449,5 +454,5 @@
                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
                &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             ) 
-            ! s-coordinate pressure gradient correction
+            ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
@@ -471,9 +476,10 @@
                       &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       &
                       &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
- 
-               IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku)
+               ! corrective term ( = 0 if z coordinate )
                zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)
                ! zhpi will be added in interior loop and zuapint will be removed in the interior loop
                ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap
+               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
+               IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku)
 
                ! zhpi will be added in interior loop and zuapint will be removed in the interior loop
@@ -488,16 +494,15 @@
             ! v direction
             ikv = mikv(ji,jj)
-            !ze3wv  = 0.5 * (e3w_0(ikv+1) / e3t_0(ikv)) * ( fse3t(ji,jj+1,ikv)-fse3t(ji,jj,ikv) )
             ze3wv  = fse3w(ji,jj+1,ikv+1)-fse3w(ji,jj,ikv+1) 
             IF ( ikv .GT. 1 ) THEN
                ! case ikv
-         !      ze3wv             =  - (fsde3w(ji,jj+1,ikv) - fsde3w(ji,jj,ikv) )
                zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            &
                      &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           &
                      &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
-
+               ! corrective term ( = 0 if z coordinate )
                zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)
                ! zhpi will be added in interior loop and zvapint will be removed in the interior loop
                va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap
+               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
                IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv) 
                ! zhpi will be added in interior loop and zvapint will be removed in the interior loop
@@ -526,5 +531,5 @@
                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 &
                   &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   ) 
-               ! corrective term               
+             !corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
                   &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
@@ -537,5 +542,5 @@
                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 &
                   &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   )
-               !
+             !corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
                   &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
@@ -563,4 +568,5 @@
                ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap
                ! put new value
+               ! -zpshpi to avoid double contribution of the partial step in the top layer 
                zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj)
                zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
@@ -574,4 +580,5 @@
                va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap
                ! put new value
+               ! -zpshpj to avoid double contribution of the partial step in the top layer 
                zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj)
                zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 4724)
@@ -246,4 +246,6 @@
       IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' )
+      IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. nn_isf .NE. 0 )   &
+           &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' )
       !
       IF( lk_esopa     )   nspg = -1
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4724)
@@ -111,5 +111,5 @@
       ! 
       z1_rau0 = 0.5_wp * r1_rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * lmask(:,:)
 
 #if ! defined key_dynspg_ts
@@ -291,5 +291,5 @@
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
-         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * lmask(:,:)
          sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 4724)
@@ -70,9 +70,9 @@
       ! and ssh which is used to calculate gradients
 
-      uo_e(:,:) = 0._wp ;   uo_e(1:jpi, 1:jpj) = ssu_m(:,:)
-      vo_e(:,:) = 0._wp ;   vo_e(1:jpi, 1:jpj) = ssv_m(:,:)
+      uo_e(:,:) = 0._wp ;   uo_e(1:jpi, 1:jpj) = ssu_m(:,:) * umask(:,:,1)
+      vo_e(:,:) = 0._wp ;   vo_e(1:jpi, 1:jpj) = ssv_m(:,:) * vmask(:,:,1)
       ff_e(:,:) = 0._wp ;   ff_e(1:jpi, 1:jpj) = ff   (:,:)
-      ua_e(:,:) = 0._wp ;   ua_e(1:jpi, 1:jpj) = utau (:,:)
-      va_e(:,:) = 0._wp ;   va_e(1:jpi, 1:jpj) = vtau (:,:)
+      ua_e(:,:) = 0._wp ;   ua_e(1:jpi, 1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
+      va_e(:,:) = 0._wp ;   va_e(1:jpi, 1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
 
       CALL lbc_lnk_e( uo_e, 'U', -1._wp, 1, 1 )
@@ -93,5 +93,5 @@
       !! so fudge some numbers all the way around the boundary
 
-      ssh_e(:,:) = 0._wp ;   ssh_e(1:jpi, 1:jpj) = ssh_m(:,:)
+      ssh_e(:,:) = 0._wp ;   ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)
       ssh_e(0    ,    :) = ssh_e(1  ,  :)
       ssh_e(jpi+1,    :) = ssh_e(jpi,  :)
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 4724)
@@ -257,4 +257,7 @@
          END DO
       END DO
+      utau(:,:) = utau(:,:) * umask(:,:,1)
+      vtau(:,:) = vtau(:,:) * vmask(:,:,1)
+      taum(:,:) = taum(:,:) * tmask(:,:,1)
       CALL lbc_lnk( taum, 'T', 1. )
 
@@ -264,4 +267,5 @@
 !CDIR COLLAPSE
       wndm(:,:) = sf(jp_wndm)%fnow(:,:,1)
+      wndm(:,:) = wndm(:,:) * tmask(:,:,1)
 
       !------------------------------------------------!
@@ -270,5 +274,5 @@
       
       CALL blk_clio_qsr_oce( qsr )
-
+      qsr(:,:) = qsr(:,:) * tmask(:,:,1) ! no shortwave radiation into the ocean beneath ice shelf
       !------------------------!
       !   Other ocean fluxes   !
@@ -376,4 +380,5 @@
          &     - zqla(:,:)             * pst(:,:) * zcevap                &   ! remove evap.   heat content at SST in Celcius
          &     + sf(jp_prec)%fnow(:,:,1) * sf(jp_tair)%fnow(:,:,1) * zcprec   ! add    precip. heat content at Tair in Celcius
+      qns(:,:) = qns(:,:) * tmask(:,:,1)
       ! NB: if sea-ice model, the snow precip are computed and the associated heat is added to qns (see blk_ice_clio)
 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 4724)
@@ -233,11 +233,14 @@
          ! Interpolate utau, vtau into the grid_V and grid_V
          !-------------------------------------------------
-
+      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1
                utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
-               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )
+               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )        &
+               &                 * MAX(tmask(ji,jj,1),tmask(ji+1,jj  ,1))
                vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
-               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )
+               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )        &
+               &                 * MAX(tmask(ji,jj,1),tmask(ji  ,jj+1,1))
             END DO
          END DO
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 4724)
@@ -156,4 +156,5 @@
             END DO
          END DO
+         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
          CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4724)
@@ -179,6 +179,6 @@
          !                                           !----------------!
          !
-         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
-         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
+         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point
+         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)                    ! (C-grid dynamics :  U- & V-points as the ocean)
          !
          t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 4724)
@@ -126,6 +126,6 @@
             DO jj = 2, jpj
                DO ji = 2, jpi   ! NO vector opt. possible
-                  u_oce(ji,jj) = 0.5_wp * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) ) * tmu(ji,jj)
-                  v_oce(ji,jj) = 0.5_wp * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) ) * tmu(ji,jj)
+                  u_oce(ji,jj) = 0.5_wp * ( ssu_m(ji-1,jj  ) * umask(ji-1,jj  ,1) + ssu_m(ji-1,jj-1) * umask(ji-1,jj-1,1) ) * tmu(ji,jj)
+                  v_oce(ji,jj) = 0.5_wp * ( ssv_m(ji  ,jj-1) * vmask(ji  ,jj-1,1) + ssv_m(ji-1,jj-1) * vmask(ji-1,jj-1,1) ) * tmu(ji,jj)
                END DO
             END DO
@@ -134,6 +134,6 @@
             !
          CASE( 'C' )                  !== C-grid ice dynamics :   U & V-points (same as ocean)
-            u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
-            v_oce(:,:) = ssv_m(:,:)
+            u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point
+            v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
             !
          END SELECT
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 4724)
@@ -32,4 +32,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE sbc_oce
 
 
@@ -201,4 +202,7 @@
       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 )  &
+      &   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity')
+
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' )
 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 4724)
@@ -233,9 +233,9 @@
                   zt_frz = -1.9 !tfreez1D( tsn(ji,jj,jk,jp_sal), zpress )
                ! compute trend
-                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                        &
-                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)        &
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
+                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
                      &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &
                      &           * r1_hisf_tbl(ji,jj)
-                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                        &
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
                      &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
                END DO
@@ -243,14 +243,14 @@
                ! level partially include in ice shelf boundary layer 
                zhk   = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
-               zalpha = rhisf_tbl(ji,jj) * ( 1 - zhk ) / fse3t(ji,jj,ikb)     ! proportion of bottom cell influenced by boundary layer
+               zalpha = rhisf_tbl(ji,jj) * ( 1._wp - zhk ) / fse3t(ji,jj,ikb)     ! proportion of bottom cell influenced by boundary layer
                ! compute tfreez for the temperature correction (we add water at freezing temperature)
                zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04
                zt_frz = -1.9 !tfreez1D( tsn(ji,jj,ikb,jp_sal), zpress )
                ! compute trend
-               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                         &
-                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)        &
+               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
+                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
                   &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
                   &              * r1_hisf_tbl(ji,jj) * zalpha
-               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                         &
+               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
                   &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * zalpha 
             END DO
@@ -261,7 +261,7 @@
                &                    'at it= ', kt,' date= ', ndastp
             IF(lwp) WRITE(numout,*) '~~~~'
-            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) )
-            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b', risf_tsc(:,:,jp_tem) )
-            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b', risf_tsc(:,:,jp_sal) )
+            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
+            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
+            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
          ENDIF
       END IF
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 4724)
@@ -81,6 +81,10 @@
       !!          di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~
       !!
-      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points
-      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 
+      !! ** Action  : compute for top and bottom interfaces
+      !!              - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points
+      !!              - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points
+      !!              - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl)
+      !!              - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl)
+      !!              - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points 
       !!----------------------------------------------------------------------
       !
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 4724)
@@ -188,4 +188,5 @@
          !
          CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+         CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition
          !
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 4724)
@@ -14,4 +14,5 @@
    !!----------------------------------------------------------------------
    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
@@ -117,4 +118,6 @@
       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 )   &
+         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' )
       !
       !                               ! ... Convection
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 4724)
@@ -57,4 +57,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzui, gzvi   !: horizontal gradient of z           at top v-point  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3rui, ge3rvi   !: horizontal gradient of T, S and rd at top v-point  
+   !! (ISF) ice load
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   riceload
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rke          !: kinetic energy
@@ -109,5 +111,6 @@
          &     gzui(jpi,jpj)     , gzvi(jpi,jpj)     ,                     &
          &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     &
-         &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                           STAT=ierr(2) )
+         &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     &
+         &     riceload(jpi,jpj),                             STAT=ierr(2) )
          !
       ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) )
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4704)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4724)
@@ -252,12 +252,12 @@
             &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level
       ELSE                                                  ! centered hpg  (eos then time stepping)
-         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
-                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
-            IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, 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
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
+                             CALL eos    ( tsb, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, 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
       ENDIF
 
