Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dommsk.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 6012)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 6012)
@@ -0,0 +1,132 @@
+MODULE dommsk
+   !!======================================================================
+   !!                       ***  MODULE dommsk   ***
+   !! Ocean initialization : domain land/sea masks, off-line case 
+   !!======================================================================
+   !! History :  3.3  ! 2010-10  (C. Ethe)  adapted from OPA_SRC/DOM/dommsk
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_msk        : compute land/ocean mask
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE lib_mpp         ! MPP library
+   USE in_out_manager  ! I/O manager
+   USE wrk_nemo  
+   USE lbclnk  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dom_msk    ! routine called by inidom.F90
+
+   REAL(wp)        :: rn_shlat   = 2.   ! type of lateral boundary condition on velocity
+   LOGICAL, PUBLIC :: ln_vorlat  = .false.   !  consistency of vorticity boundary condition 
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_msk
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk  ***
+      !!
+      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
+      !!
+      !! ** Method  :   The interior ocean/land mask is computed from tmask
+      !!              setting to zero the duplicated row and lines due to
+      !!              MPP exchange halos, est-west cyclic and north fold
+      !!              boundary conditions.
+      !!
+      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
+      !!               tpol     : ???
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
+      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
+      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
+      !
+      !!---------------------------------------------------------------------
+      
+      CALL wrk_alloc( jpi, jpj, imsk )
+      !
+      ! Interior domain mask (used for global sum)
+      ! --------------------
+      ssmask(:,:)  = tmask(:,:,1)
+      tmask_i(:,:) = tmask(:,:,1)
+      iif = jpreci                        ! thickness of exchange halos in i-axis
+      iil = nlci - jpreci + 1
+      ijf = jprecj                        ! thickness of exchange halos in j-axis
+      ijl = nlcj - jprecj + 1
+      !
+      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
+      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
+      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
+      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
+      !
+      !                                   ! north fold mask
+      tpol(1:jpiglo) = 1._wp
+      !                                
+      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
+      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
+      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
+         IF( mjg(ijl-1) == jpjglo-1 ) THEN
+            DO ji = iif+1, iil-1
+               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
+            END DO
+         ENDIF
+      ENDIF 
+      !
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
+      ! least 1 wet u point
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector loop
+            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+         END DO
+         DO ji = 1, jpim1      ! NO vector opt.
+            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
+               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
+         END DO
+      END DO
+      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( ssvmask, 'V', 1._wp )
+      CALL lbc_lnk( ssfmask, 'F', 1._wp )
+
+      ! 3. Ocean/land mask at wu-, wv- and w points 
+      !----------------------------------------------
+      wmask (:,:,1) = tmask(:,:,1) ! ????????
+      wumask(:,:,1) = umask(:,:,1) ! ????????
+      wvmask(:,:,1) = vmask(:,:,1) ! ????????
+      DO jk=2,jpk
+         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
+      END DO
+      !
+      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
+         imsk(:,:) = INT( tmask_i(:,:) )
+         WRITE(numout,*) ' tmask_i : '
+         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
+         WRITE (numout,*)
+         WRITE (numout,*) ' dommsk: tmask for each level'
+         WRITE (numout,*) ' ----------------------------'
+         DO jk = 1, jpk
+            imsk(:,:) = INT( tmask(:,:,jk) )
+            WRITE(numout,*)
+            WRITE(numout,*) ' level = ',jk
+            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
+         END DO
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, imsk )
+      !
+   END SUBROUTINE dom_msk
+   !!======================================================================
+END MODULE dommsk
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 6012)
@@ -489,12 +489,11 @@
          CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
 
-         ! Partial steps: before Horizontal DErivative
-         IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
-            &            CALL zps_hde    ( kt, jpts, pts, 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)                            &
-            &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-            &                                        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 first ocean level
+      ! Partial steps: before Horizontal DErivative
+      IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
+         &            CALL zps_hde    ( kt, jpts, pts, 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)                            &
+         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+         &                                        rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
 
          rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 6012)
@@ -636,15 +636,16 @@
 !!gm
 
-            IF( ln_zps .AND. .NOT. lk_c1d ) THEN      ! Partial steps: before horizontal gradient
-               IF(ln_isfcav) THEN                        ! ocean cavities: top and bottom cells (ISF)
-                  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,     &
-                     &                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-                     &                     grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )
-               ELSE                                      ! no ocean cavities: bottom cells
-                  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! 
-                     &                        rhd, gru , grv          )  ! of t, s, rd at the last ocean level
-               ENDIF
-            ENDIF
-            !
+            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
+               &  CALL zps_hde    ( kt, 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. .NOT. lk_c1d .AND.       ln_isfcav)      &
+               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF)
+               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level
+
+#if defined key_zdfkpp
+            CALL eos( tsn, rhd, fsdept_n(:,:,:) )                      ! Compute rhd
+!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd
+#endif
+
             DEALLOCATE( t_bkginc )
             DEALLOCATE( s_bkginc )
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 6012)
@@ -93,5 +93,5 @@
       ! -----------------------------------------------------------------------
 !!gm replace these lines :
-      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0
+      z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0
       IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
 !!gm   by :
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 6012)
@@ -194,5 +194,5 @@
                   DO ji = 1,jpi
                      ! Elevation
-                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)        
+                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)        
                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
@@ -358,6 +358,6 @@
                X1=ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj)
-               out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
+               out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj)
+               out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj)
             END DO
          END DO
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6012)
@@ -101,9 +101,5 @@
       IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
       ! Add ice shelf heat & salt input
-      IF( nn_isf .GE. 1 )  THEN
-          z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) )
-          z_frc_trd_s = z_frc_trd_s + glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) )
-      ENDIF
-
+      IF( ln_isf    ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) )
       ! Add penetrative solar radiation
       IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6012)
@@ -256,5 +256,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask              !: land/ocean mask of barotropic stream function
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level                (ISF)
@@ -390,5 +390,5 @@
 
       ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       &
-         &     tmask_i(jpi,jpj) , tmask_h(jpi, jpj),                                       & 
+         &     tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                       & 
          &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
          &     bmask(jpi,jpj)   ,                                                          &
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6012)
@@ -506,34 +506,4 @@
       CALL iom_close( inum )
       
-!!gm   THIS is TO BE REMOVED !!!!!!!
-
-! need to be define for the extended grid south of -80S
-! some point are undefined but you need to have e1 and e2 .NE. 0
-      WHERE (e1t==0.0_wp)
-         e1t=1.0e2
-      END WHERE
-      WHERE (e1v==0.0_wp)
-         e1v=1.0e2
-      END WHERE
-      WHERE (e1u==0.0_wp)
-         e1u=1.0e2
-      END WHERE
-      WHERE (e1f==0.0_wp)
-         e1f=1.0e2
-      END WHERE
-      WHERE (e2t==0.0_wp)
-         e2t=1.0e2
-      END WHERE
-      WHERE (e2v==0.0_wp)
-         e2v=1.0e2
-      END WHERE
-      WHERE (e2u==0.0_wp)
-         e2u=1.0e2
-      END WHERE
-      WHERE (e2f==0.0_wp)
-         e2f=1.0e2
-      END WHERE
-!!gm end
-       
     END SUBROUTINE hgr_read
     
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6012)
@@ -975,5 +975,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( 'Only vvl_zstar has been tested with ice shelf cavity' )
+      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
       IF(lwp) THEN                   ! Print the choice
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6012)
@@ -1258,4 +1258,14 @@
       END WHERE
 
+      ! remove very shallow ice shelf (less than ~ 10m if 75L)
+      WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1)
+         misfdep = 0; risfdep = 0.0_wp;
+         mbathy  = 0; bathy   = 0.0_wp;
+      END WHERE
+      WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp)
+         misfdep = 0; risfdep = 0.0_wp;
+         mbathy  = 0; bathy   = 0.0_wp;
+      END WHERE
+ 
 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
       icompt = 0 
@@ -1323,5 +1333,4 @@
                      zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
                          &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
- 
                      IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
                         IF (zbathydiff <= zrisfdepdiff) THEN
@@ -1369,5 +1378,5 @@
          END DO
  
- ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 
+ ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -1405,5 +1414,5 @@
             mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF
- ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 
+ ! point V misdep(ji,jj) == mbathy(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -1443,5 +1452,5 @@
          ENDIF 
  
- ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 
+ ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -1480,5 +1489,5 @@
          ENDIF 
  
- ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 
+ ! point U misfdep(ji,jj) == bathy(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -1749,4 +1758,8 @@
 ! end check compatibility ice shelf/bathy
       ! remove very shallow ice shelf (less than ~ 10m if 75L)
+      WHERE (risfdep(:,:) <= 10._wp)
+         misfdep = 1; risfdep = 0.0_wp;
+      END WHERE
+
       IF( icompt == 0 ) THEN 
          IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
@@ -1832,5 +1845,5 @@
                ENDIF 
             !       ... on ik / ik-1 
-               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) 
+               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90	(revision 6012)
@@ -20,5 +20,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean
+   USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean
    USE sbcrnf          ! river runoff 
    USE sbcisf          ! ice shelf
@@ -91,11 +91,11 @@
       END DO
       !
-      IF( ln_rnf                     )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field)
+      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field)
       !
-      IF( ln_divisf .AND. nn_isf > 0 )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
+      IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
       !
-      IF( ln_iscpl  .AND. ln_hsb     )   CALL iscpl_div( hdivn )   !==  ice sheet  ==!   (update hdivn field)
+      IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn )  !==  ice sheet  ==!   (update hdivn field)
       !
-      CALL lbc_lnk( hdivn, 'T', 1. )                       !==  lateral boundary cond.  ==!   (no sign change)
+      CALL lbc_lnk( hdivn, 'T', 1. )                !==  lateral boundary cond.  ==!   (no sign change)
       !
       IF( nn_timing == 1 )  CALL timing_stop('div_hor')
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6012)
@@ -45,4 +45,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -130,4 +131,10 @@
       INTEGER ::   ioptio = 0      ! temporary integer
       INTEGER ::   ios             ! Local integer output status for namelist read
+      !!
+      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
+      REAL(wp) ::   znad
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop, zrhd ! hypothesys on isf density
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_isf  ! density at bottom of ISF
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  ziceload     ! density at bottom of ISF
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
@@ -190,6 +197,51 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
       ! 
-      ! initialisation of ice load
-      riceload(:,:)=0.0
+      ! initialisation of ice shelf load
+      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0
+      IF (       ln_isfcav ) THEN
+         CALL wrk_alloc( jpi,jpj, 2,  ztstop) 
+         CALL wrk_alloc( jpi,jpj,jpk, zrhd  )
+         CALL wrk_alloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'  
+
+         ! To use density and not density anomaly
+         znad=1._wp
+         
+         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
+         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
+
+         ! compute density of the water displaced by the ice shelf 
+         DO jk = 1, jpk
+            CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk))
+         END DO
+      
+         ! compute rhd at the ice/oce interface (ice shelf side)
+         CALL eos(ztstop,risfdep,zrhdtop_isf)
+
+         ! Surface value + ice shelf gradient
+         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
+         ! divided by 2 later
+         ziceload = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt=mikt(ji,jj)
+               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
+               DO jk=2,ikt-1
+                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
+                     &                              * (1._wp - tmask(ji,jj,jk))
+               END DO
+               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
+                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
+            END DO
+         END DO
+         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
+
+         CALL wrk_dealloc( jpi,jpj, 2,  ztstop) 
+         CALL wrk_dealloc( jpi,jpj,jpk, zrhd  )
+         CALL wrk_dealloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+      END IF
       !
    END SUBROUTINE dyn_hpg_init
@@ -447,5 +499,5 @@
    SUBROUTINE hpg_isf( kt )
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE hpg_sco  ***
+      !!                  ***  ROUTINE hpg_isf  ***
       !!
       !! ** Method  :   s-coordinate case. Jacobian scheme.
@@ -466,57 +518,26 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd
+      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj, 2, ztstop) 
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
-      CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
-      !
-     IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
-      ENDIF
-
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,  2, ztstop) 
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj)
+      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce )
+      !
       ! Local constant initialization
       zcoef0 = - grav * 0.5_wp
+   
       ! To use density and not density anomaly
-!      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
-!      ELSE                 ;     znad = 0._wp         ! Fixed volume
-!      ENDIF
       znad=1._wp
+
       ! iniitialised to 0. zhpi zhpi 
       zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
 
-      ! Partial steps: top & bottom before horizontal gradient
-!jc      CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi,   & 
-!jc               &       rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-!jc               &      grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )
-
-!==================================================================================     
-!=====Compute iceload and contribution of the half first wet layer ================= 
-!===================================================================================
-
-      ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
-      ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
-
-      ! compute density of the water displaced by the ice shelf 
-      zrhd = rhd ! save rhd
-      DO jk = 1, jpk
-           zdept(:,:)=gdept_1d(jk)
-           CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))
-      END DO
-      WHERE ( tmask(:,:,:) == 1._wp)
-        rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd
-      END WHERE
-      
-      ! compute rhd at the ice/oce interface (ice shelf side)
-      CALL eos(ztstop,risfdep,zrhdtop_isf)
-
       ! compute rhd at the ice/oce interface (ocean side)
+      ! usefull to reduce residual current in the test case ISOMIP with no melting
       DO ji=1,jpi
         DO jj=1,jpj
@@ -527,22 +548,8 @@
       END DO
       CALL eos(ztstop,risfdep,zrhdtop_oce)
-      !
-      ! Surface value + ice shelf gradient
-      ! 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
-         DO ji = 1, jpi   ! vector opt.
-            ikt=mikt(ji,jj)
-            ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
-            DO jk=2,ikt-1
-               ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
-                  &                              * (1._wp - tmask(ji,jj,jk))
-            END DO
-            IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &
-                               &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
-         END DO
-      END DO
-      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)
+
+!==================================================================================     
+!===== Compute surface value ===================================================== 
+!==================================================================================
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -550,14 +557,14 @@
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
             ! 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  ) )   &
-               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
-               &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              ) 
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    &
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   &
-               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
-               &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             ) 
+            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) )   &
+               &                                  - 0.5_wp * fse3w(ji,jj,ikt)                                         &
+               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
+               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            ) 
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j)                                    &
+               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
+               &                                  - 0.5_wp * fse3w(ji,jj,ikt)                                         & 
+               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
+               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
@@ -570,141 +577,34 @@
          END DO
       END DO
-!==================================================================================     
-!===== Compute partial cell contribution for the top cell ========================= 
-!==================================================================================
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            iku = miku(ji,jj) ; 
-            zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp
-            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-            ! u direction
-            IF ( iku .GT. 1 ) THEN
-               ! case iku
-               zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            &
-                      &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       &
-                      &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
-               ! 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
-               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)
-
-               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
-               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &    
-                  &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          &
-                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   &
-                  &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          &
-                  &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   )
-               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
-            END IF
-               
-            ! v direction
-            ikv = mikv(ji,jj)
-            ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-            IF ( ikv .GT. 1 ) THEN
-               ! case ikv
-               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
-               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) 
-               
-               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
-               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              &
-                  &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         &
-                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    &
-                  &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         &
-                  &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  )
-               zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
-            END IF
-         END DO
-      END DO
 
 !==================================================================================     
 !===== Compute interior value ===================================================== 
 !==================================================================================
-
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            iku=miku(ji,jj); ikv=mikv(ji,jj)
-            DO jk = 2, jpkm1
+      ! interior value (2=<jk=<jpkm1)
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                ! hydrostatic pressure gradient along s-surfaces
-               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              &
-                  &                            + zcoef0 / e1u(ji,jj)                                                           &
-                  &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   &
-                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   ) 
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
+                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
+                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
+                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
                ! s-coordinate pressure gradient correction
-               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
-               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    &
-                  &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
-               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
-
-               ! hydrostatic pressure gradient along s-surfaces
-               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
-               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              &
-                  &                            + zcoef0 / e2v(ji,jj)                                                           &
-                  &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   &
-                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 &
-                  &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   )
-               ! s-coordinate pressure gradient correction
-               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
-               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     &
-                  &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
+               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)
+               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)
                ! add to the general momentum trend
-               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
+               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
-!==================================================================================     
-!===== Compute bottom cell contribution (partial cell) ============================ 
-!==================================================================================
-
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-            iku = mbku(ji,jj)
-            ikv = mbkv(ji,jj)
-
-            IF (iku .GT. 1) THEN
-               ! remove old value (interior case)
-               zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   &
-                     &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj)
-               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) 
-               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap
-            END IF
-            ! v direction
-            IF (ikv .GT. 1) THEN
-               ! remove old value (interior case)
-               zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   &
-                     &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj)
-               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)   
-               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
-            END IF
-         END DO
-      END DO
-     
-      ! set back to original density value into the ice shelf cell (maybe useless because it is masked)
-      rhd = zrhd
-      !
-      CALL wrk_dealloc( jpi,jpj,2, ztstop)
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
-      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
+      !
+      CALL wrk_dealloc( jpi,jpj,2  , ztstop)
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj)
+      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce )
       !
    END SUBROUTINE hpg_isf
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6012)
@@ -95,4 +95,5 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikt          ! local integers
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
@@ -220,15 +221,15 @@
                ! Add volume filter correction: compatibility with tracer advection scheme
                ! => time filter + conservation correction (only at the first level)
-               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
+               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
                   fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
-                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
+                                 &                                         - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
                ELSE                     ! if ice shelf melting
                   DO jj = 1,jpj
                      DO ji = 1,jpi
-                        jk = mikt(ji,jj)
-                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
+                        ikt = mikt(ji,jj)
+                        fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0                     &
                                           &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
                                           &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
-                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
+                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,ikt)
                      END DO
                   END DO
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 6012)
@@ -231,6 +231,6 @@
       IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme' )
-      IF( ln_dynspg_ts .AND. ln_isfcav )   &
-           &   CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' )
+      IF( ln_dynspg_exp .AND. ln_isfcav )   &
+           &   CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' )
       !
       IF( ln_dynspg_exp)   nspg =  0
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6012)
@@ -145,4 +145,5 @@
       INTEGER  ::   ji, jj, jk, jn   		! dummy loop indices
       INTEGER  ::   ikbu, ikbv, noffset    	! local integers
+      INTEGER  ::   iktu, iktv             	! local integers
       REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf  	! local scalars
       REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      -
@@ -384,6 +385,6 @@
       DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
          DO ji = fs_2, fs_jpim1
-             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
-             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
+             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
+             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
           END DO
       END DO 
@@ -413,4 +414,29 @@
       zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
       zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
+      !       
+      !                                         ! Add top stress contribution from baroclinic velocities:      
+      IF (ln_bt_fw) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               iktu = miku(ji,jj)
+               iktv = mikv(ji,jj)
+               zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
+               zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
+            END DO
+         END DO
+      ELSE
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               iktu = miku(ji,jj)
+               iktv = mikv(ji,jj)
+               zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
+               zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      ! Note that the "unclipped" top friction parameter is used even with explicit drag
+      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * tfrua(:,:) * zwx(:,:)
+      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * tfrva(:,:) * zwy(:,:)
       !       
       IF (ln_bt_fw) THEN                        ! Add wind forcing
@@ -544,8 +570,8 @@
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
                DO ji = 2, fs_jpim1   ! Vector opt.
-                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)     &
+                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
                      &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
                      &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
-                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)     &
+                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
                      &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
                      &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
@@ -607,5 +633,5 @@
             END DO
          END DO
-         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
          CALL lbc_lnk( ssha_e, 'T',  1._wp )
 
@@ -622,10 +648,10 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1      ! NO Vector Opt.
-                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)  &
-                     &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
-                     &              +   e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
-                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)  &
-                     &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
-                     &              +   e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
+                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
+                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
+                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
+                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
+                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
+                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
                END DO
             END DO
@@ -661,8 +687,8 @@
             DO jj = 2, jpjm1                            
                DO ji = 2, jpim1
-                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e1e2u(ji  ,jj)    &
+                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
                      &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
                      &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
-                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e1e2v(ji  ,jj  )  &
+                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
                      &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
                      &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
@@ -736,4 +762,8 @@
          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
          !
+         ! Add top stresses:
+         zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:)
+         zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
+         !
          ! Surface pressure trend:
          DO jj = 2, jpjm1
@@ -755,5 +785,5 @@
                             &                                 + zu_trd(ji,jj)   &
                             &                                 + zu_frc(ji,jj) ) & 
-                            &   ) * umask(ji,jj,1)
+                            &   ) * ssumask(ji,jj)
 
                   va_e(ji,jj) = (                                 vn_e(ji,jj)   &
@@ -761,14 +791,14 @@
                             &                                 + zv_trd(ji,jj)   &
                             &                                 + zv_frc(ji,jj) ) &
-                            &   ) * vmask(ji,jj,1)
-               END DO
-            END DO
-
-         ELSE					 	! Flux form
+                            &   ) * ssvmask(ji,jj)
+               END DO
+            END DO
+
+         ELSE                                         ! Flux form
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
 
-                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
-                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
+                  zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj))
+                  zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj))
 
                   ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
@@ -791,6 +821,6 @@
             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
-            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
-            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
+            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
+            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
             !
          ENDIF
@@ -827,5 +857,5 @@
             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
-         ELSE                     			   ! Sum transports
+         ELSE                                              ! Sum transports
             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
@@ -881,6 +911,6 @@
          END DO
          ! Save barotropic velocities not transport:
-         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
-         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
+         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
+         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
       ENDIF
       !
@@ -897,5 +927,5 @@
       !
       IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
-         IF ( Agrif_NbStepint().EQ.0 ) THEN
+         IF ( Agrif_NbStepint() == 0 ) THEN
             ub2_i_b(:,:) = 0.e0
             vb2_i_b(:,:) = 0.e0
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 6012)
@@ -90,6 +90,6 @@
          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) )
-               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(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   
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 6012)
@@ -112,10 +112,12 @@
       !!
       INTEGER  ::   ji , jj , jk    ! dummy loop indices
-      INTEGER  ::   ii0, ii1, iku   ! temporary integer
-      INTEGER  ::   ij0, ij1, ikv   ! temporary integer
+      INTEGER  ::   ii0, ii1        ! temporary integer
+      INTEGER  ::   ij0, ij1        ! temporary integer
       REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars
       REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
       REAL(wp) ::   zck, zfk,      zbw             !   -      -
+      REAL(wp) ::   zdepu, zdepv                   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zslpml_hmlpu, zslpml_hmlpv
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr
@@ -126,4 +128,5 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
+      CALL wrk_alloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv )
 
       zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
@@ -149,4 +152,12 @@
                zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
                zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 
+               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj)
             END DO
          END DO
@@ -171,4 +182,22 @@
       ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
       !
+      IF ( ln_isfcav ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)       &
+                  &                                  - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji+1,jj  ) ) )
+               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)       &
+                  &                                  - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji  ,jj+1) ) )
+            END DO
+         END DO
+      ELSE
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
+               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
+            END DO
+         END DO
+      END IF
+
       DO jk = 2, jpkm1                            !* Slopes at u and v points
          DO jj = 2, jpjm1
@@ -186,12 +215,14 @@
                zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
                zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
-               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
-                  &                   + zfi  * uslpml(ji,jj)                                                     &
-                  &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
-                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk)
-               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
-                  &                   + zfj  * vslpml(ji,jj)                                                     &
-                  &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
-                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,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) )                            &
+                                - ( 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) )                              &
+                                - ( risfdep(ji,jj)    + risfdep(ji,jj+1)    ) - fse3v(ji,jj,mikv(ji,jj)) )
+               !
+               zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps )                                     &
+                  &                      + zfi  * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk)
+               zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps )                                     &
+                  &                      + zfj  * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk)
 !!gm  modif to suppress omlmask.... (as in Griffies case)
 !               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
@@ -265,7 +296,7 @@
                   &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) * e2t(ji,jj)
                zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
-                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk)
+                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * wmask (ji,jj,jk)
                zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
-                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk)
+                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * wmask (ji,jj,jk)
                !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
                !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
@@ -274,7 +305,7 @@
                !                                        ! 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
-               zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )
-               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
-               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
+               zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp )
+               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
+               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
 
 !!gm  modif to suppress omlmask....  (as in Griffies operator)
@@ -340,5 +371,4 @@
       CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
 
-
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
@@ -347,4 +377,5 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
+      CALL wrk_dealloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv )
       !
       IF( nn_timing == 1 )  CALL timing_stop('ldf_slp')
@@ -486,5 +517,5 @@
                   !
                   jk = nmln(ji,jj+jp) + 1
-                  IF( jk .GT. mbkt(ji,jj+jp) ) THEN  !ML reaches bottom
+                  IF( jk >  mbkt(ji,jj+jp) ) THEN  !ML reaches bottom
                      ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp
                   ELSE
@@ -699,5 +730,5 @@
             zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           &
                &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj)
-            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)           &
+            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)             &
                &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask(ji,jj,ik)
             zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           &
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 6012)
@@ -44,8 +44,8 @@
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
+   LOGICAL , PUBLIC ::   ln_isf         !: ice shelf melting
    LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS      
    LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice)
    INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3)
-   INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4) 
    INTEGER , PUBLIC ::   nn_ice_embd    !: flag for levitating/embedding sea-ice in the ocean
    !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect)
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6012)
@@ -35,36 +35,36 @@
    ! public in order to be able to output then 
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2]
    REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
-   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence 
-   INTEGER , PUBLIC ::   nn_isfblk                   !: 
-   INTEGER , PUBLIC ::   nn_gammablk                 !:
-   LOGICAL , PUBLIC ::   ln_conserve                 !:
-   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient
-   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient 
-   REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence
+   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified  
+   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
+   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
+   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
+   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
 
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
-   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m]
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl 
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-
-   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
+   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base
+
+   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K]
+   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s]
+   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3]
+   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C]
+   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg]
 
 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
-   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
-   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
-   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
-   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
+   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
+   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
+   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read
+   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read
    
    !! * Substitutions
@@ -77,22 +77,21 @@
 CONTAINS
  
-   SUBROUTINE sbc_isf(kt)
+  SUBROUTINE sbc_isf(kt)
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE sbc_isf  ***
-      !!---------------------------------------------------------------------
-      INTEGER, INTENT(in)          ::   kt         ! ocean time step
-      !
-      INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
-      INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
-      REAL(wp)                     ::   rmin
-      REAL(wp)                     ::   zhk
-      REAL(wp)                     ::   zt_frz, zpress
-      CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
-      CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
-      CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
-      INTEGER           ::   ios           ! Local integer output status for namelist read
-      !!
-      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
-         &                 sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !!                  ***  ROUTINE sbc_isf  ***
+      !!
+      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
+      !!              melting and freezing 
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt                   ! ocean time step
+      !
+      INTEGER               :: ji, jj               ! loop index
+      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
       !!---------------------------------------------------------------------
       !
@@ -100,133 +99,20 @@
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
          !                                      ! ====================== !
-         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
-         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
-
-         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
-         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
-         IF(lwm) WRITE ( numond, namsbc_isf )
-
-         IF ( lwp ) WRITE(numout,*)
-         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
-         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
-         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
-         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
-         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
-         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
-         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
-         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
-         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
-         IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
-         IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
-         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
-         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
-            rdivisf = 1._wp
-         ELSE
-            rdivisf = 0._wp
-         END IF
-         !
-         ! Allocate public variable
-         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
-         !
-         ! initialisation
-         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
-         risf_tsc(:,:,:)  = 0._wp
-         !
-         ! define isf tbl tickness, top and bottom indice
-         IF      (nn_isf == 1) THEN
-            rhisf_tbl(:,:) = rn_hisf_tbl
-            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
-         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
-            ALLOCATE( sf_rnfisf(1), STAT=ierror )
-            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-
-            !: read effective lenght (BG03)
-            IF (nn_isf == 2) THEN
-               ! Read Data and save some integral values
-               CALL iom_open( sn_Leff_isf%clname, inum )
-               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
-               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
-               CALL iom_close(inum)
-               !
-               risfLeff = risfLeff*1000           !: convertion in m
-            END IF
-
-           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
-            CALL iom_open( sn_depmax_isf%clname, inum )
-            cvarhisf = TRIM(sn_depmax_isf%clvar)
-            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
-            CALL iom_close(inum)
-            !
-            CALL iom_open( sn_depmin_isf%clname, inum )
-            cvarzisf = TRIM(sn_depmin_isf%clvar)
-            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
-            CALL iom_close(inum)
-            !
-            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
-
-           !! compute first level of the top boundary layer
-           DO ji = 1, jpi
-              DO jj = 1, jpj
-                  jk = 2
-                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
-                  misfkt(ji,jj) = jk-1
-               END DO
-            END DO
-
-         ELSE IF ( nn_isf == 4 ) THEN
-            ! as in nn_isf == 1
-            rhisf_tbl(:,:) = rn_hisf_tbl
-            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
-            
-            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
-            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
-            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
-            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
-         END IF
-         
-         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
-         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               ikt = misfkt(ji,jj)
-               ikb = misfkt(ji,jj)
-               ! thickness of boundary layer at least the top level thickness
-               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
-
-               ! determine the deepest level influenced by the boundary layer
-               ! test on tmask useless ?????
-               DO jk = ikt, mbkt(ji,jj)
-                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
-               END DO
-               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
-               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
-               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
-
-               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
-               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
-            END DO
-         END DO
+         CALL sbc_isf_init
+      !                                         ! ---------------------------------------- !
+      ELSE                                      !          Swap of forcing fields          !
+         !                                      ! ---------------------------------------- !
+         fwfisf_b  (:,:  ) = fwfisf  (:,:  )    ! Swap the ocean forcing fields except at nit000
+         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)    ! where before fields are set at the end of the routine
          !
       END IF
 
-      !                                            ! ---------------------------------------- !
-      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
-         !                                         ! ---------------------------------------- !
-         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
-         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
-         !
-      ENDIF
-
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
-
-
-         ! compute salf and heat flux
-         IF (nn_isf == 1) THEN
-            ! realistic ice shelf formulation
+         ! allocation
+         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
+
+         ! compute salt and heat flux
+         SELECT CASE ( nn_isf )
+         CASE ( 1 )    ! realistic ice shelf formulation
             ! compute T/S/U/V for the top boundary layer
             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
@@ -242,49 +128,43 @@
             CALL sbc_isf_cav (kt)
 
-         ELSE IF (nn_isf == 2) THEN
-            ! Beckmann and Goosse parametrisation 
+         CASE ( 2 )    ! Beckmann and Goosse parametrisation 
             stbl(:,:)   = soce
             CALL sbc_isf_bg03(kt)
 
-         ELSE IF (nn_isf == 3) THEN
-            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
+         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
-            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
+            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting) 
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux
             stbl(:,:)   = soce
 
-         ELSE IF (nn_isf == 4) THEN
-            ! specified fwf and heat flux forcing beneath the ice shelf
+         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
             CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
-            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
-            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
-            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
+            fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting)
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux
             stbl(:,:)   = soce
 
-         END IF
+         END SELECT
+
          ! compute tsc due to isf
-         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
-!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
-         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
-         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
+         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
+         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
+         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj))
+            END DO
+         END DO
+         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
          
-         ! salt effect already take into account in vertical advection
-         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
-
-         ! output
-         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
-         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
-
-         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
-         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
- 
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
+         risf_tsc(:,:,jp_sal) = 0.0_wp
+
          ! lbclnk
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
          CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
-         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
-         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
-
-         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.)
+         CALL lbc_lnk(qisf(:,:)           ,'T',1.)
+
+         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
             IF( ln_rstart .AND.    &                     ! Restart: read in restart file
                  & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
@@ -297,6 +177,12 @@
                risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
             END IF
-         ENDIF
+         END IF
          ! 
+         ! output
+         CALL iom_put('qisf'  , qisf)
+         CALL iom_put('fwfisf', fwfisf)
+
+         ! deallocation
+         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
       END IF
       !  
@@ -317,79 +203,205 @@
                &    STAT= sbc_isf_alloc )
          !
-         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
+         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc )
          IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
          !
-      ENDIF
+      END IF
   END FUNCTION
 
-
-   SUBROUTINE sbc_isf_bg03(kt)
-      !!==========================================================================
-      !!                 *** SUBROUTINE sbcisf_bg03  ***
-      !! add net heat and fresh water flux from ice shelf melting
-      !! into the adjacent ocean using the parameterisation by
-      !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
-      !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
-      !!  (hereafter BG)
-      !!==========================================================================
-      !!----------------------------------------------------------------------
-      !!   sbc_isf_bg03      : routine called from sbcmod
-      !!----------------------------------------------------------------------
-      !!
-      !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
-      !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
-      !!
+  SUBROUTINE sbc_isf_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_init  ***
+      !!
+      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER               :: ji, jj, jk           ! loop index
+      INTEGER               :: ik                   ! current level index
+      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer
+      INTEGER               :: inum, ierror
+      INTEGER               :: ios                  ! Local integer output status for namelist read
+      REAL(wp)              :: zhk
+      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
+      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale
+      !!----------------------------------------------------------------------
+      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
+                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_isf )
+
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
+      IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
+      IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
+      IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
+      IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
+      IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
+      IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
+      !
+      ! Allocate public variable
+      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
+      !
+      ! initialisation
+      qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp
+      risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp
+      !
+      ! define isf tbl tickness, top and bottom indice
+      SELECT CASE ( nn_isf )
+      CASE ( 1 ) 
+         rhisf_tbl(:,:) = rn_hisf_tbl
+         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+
+      CASE ( 2 , 3 )
+         ALLOCATE( sf_rnfisf(1), STAT=ierror )
+         ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+
+         !  read effective lenght (BG03)
+         IF (nn_isf == 2) THEN
+            CALL iom_open( sn_Leff_isf%clname, inum )
+            cvarLeff = TRIM(sn_Leff_isf%clvar)
+            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
+            CALL iom_close(inum)
+            !
+            risfLeff = risfLeff*1000.0_wp           !: convertion in m
+         END IF
+
+         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
+         CALL iom_open( sn_depmax_isf%clname, inum )
+         cvarhisf = TRIM(sn_depmax_isf%clvar)
+         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
+         CALL iom_close(inum)
+         !
+         CALL iom_open( sn_depmin_isf%clname, inum )
+         cvarzisf = TRIM(sn_depmin_isf%clvar)
+         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
+         CALL iom_close(inum)
+         !
+         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
+
+         !! compute first level of the top boundary layer
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+                ik = 2
+                DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
+                misfkt(ji,jj) = ik-1
+            END DO
+         END DO
+
+      CASE ( 4 ) 
+         ! as in nn_isf == 1
+         rhisf_tbl(:,:) = rn_hisf_tbl
+         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+         
+         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
+         ALLOCATE( sf_fwfisf(1), STAT=ierror )
+         ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+
+      END SELECT
+         
+      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
+
+      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ikt = misfkt(ji,jj)
+            ikb = misfkt(ji,jj)
+            ! thickness of boundary layer at least the top level thickness
+            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
+
+            ! determine the deepest level influenced by the boundary layer
+            DO jk = ikt+1, mbkt(ji,jj)
+               IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+            END DO
+            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
+            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
+            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
+
+            zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
+            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+         END DO
+      END DO
+
+  END SUBROUTINE sbc_isf_init
+
+  SUBROUTINE sbc_isf_bg03(kt)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_bg03  ***
+      !!
+      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
+      !!          into the adjacent ocean
+      !!
+      !! ** Method  :   See reference
+      !!
+      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!         (hereafter BG)
       !! History :
-      !!      !  06-02  (C. Wang) Original code
+      !!         06-02  (C. Wang) Original code
       !!----------------------------------------------------------------------
       INTEGER, INTENT ( in ) :: kt
       !
-    INTEGER :: ji, jj, jk, jish  !temporary integer
-    INTEGER :: ijkmin
-    INTEGER :: ii, ij, ik 
-    INTEGER :: inum
-
-    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
-    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
-    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
-    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
-    
-    !!----------------------------------------------------------------------
-    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
-     !
-    DO ji = 1, jpi
-       DO jj = 1, jpj
-          ik = misfkt(ji,jj)
-          !! Initialize arrays to 0 (each step)
-          zt_sum = 0.e0_wp
-          IF ( ik .GT. 1 ) THEN
-    ! 3. -----------the average temperature between 200m and 600m ---------------------
-             DO jk = misfkt(ji,jj),misfkb(ji,jj)
-             ! Calculate freezing temperature
-                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, fsdept(ji,jj,ik)) 
-                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
-             ENDDO
-             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
-    
-    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
-          ! For those corresponding to zonal boundary    
-             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
-                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
+      INTEGER  :: ji, jj, jk ! dummy loop index
+      INTEGER  :: ik         ! current level
+      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
+      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
+      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
+      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
+      !!----------------------------------------------------------------------
+
+      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
+      !
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+            ik = misfkt(ji,jj)
+            !! Initialize arrays to 0 (each step)
+            zt_sum = 0.e0_wp
+            IF ( ik > 1 ) THEN
+               ! 1. -----------the average temperature between 200m and 600m ---------------------
+               DO jk = misfkt(ji,jj),misfkb(ji,jj)
+                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
+                  ! after verif with UNESCO, wrong sign in BG eq. 2
+                  ! Calculate freezing temperature
+                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
+                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
+               END DO
+               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
+               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
+               ! For those corresponding to zonal boundary    
+               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
+                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
              
-             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                  
-             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
-             !add to salinity trend
-          ELSE
-             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
-          END IF
-       END DO
-    END DO
-    !
-    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                  
+               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
+               !add to salinity trend
+            ELSE
+               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
+            END IF
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
       !
   END SUBROUTINE sbc_isf_bg03
 
-
-   SUBROUTINE sbc_isf_cav( kt )
+  SUBROUTINE sbc_isf_cav( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_isf_cav  ***
@@ -406,155 +418,119 @@
       INTEGER, INTENT(in)          ::   kt         ! ocean time step
       !
-      LOGICAL :: ln_isomip = .true.
-      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
-      REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d 
-      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      INTEGER  ::   nit
       REAL(wp) ::   zlamb1, zlamb2, zlamb3
       REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
       REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
-      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
-      REAL(wp) ::   zgammat, zgammas
-      REAL(wp) ::   zeps = 1.e-20_wp        !==   Local constant initialization   ==!
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
-      INTEGER  ::   ierror     ! return error code
-      LOGICAL  ::   lit=.TRUE.
-      INTEGER  ::   nit
+      REAL(wp) ::   zeps = 1.e-20_wp        
+      REAL(wp) ::   zerr
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz
+      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
+      LOGICAL  ::   lit
       !!---------------------------------------------------------------------
-      !
-      ! coeficient for linearisation of tfreez
-      zlamb1=-0.0575
-      zlamb2=0.0901
-      zlamb3=-7.61e-04
+      ! coeficient for linearisation of potential tfreez
+      ! Crude approximation for pressure (but commonly used)
+      zlamb1 =-0.0573_wp
+      zlamb2 = 0.0832_wp
+      zlamb3 =-7.53e-08_wp * grav * rau0
       IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
       !
-      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
-
-      zcfac=0.0_wp 
-      IF (ln_conserve)  zcfac=1.0_wp
-      zpress(:,:)=0.0_wp
-      zgammat2d(:,:)=0.0_wp
-      zgammas2d(:,:)=0.0_wp
-      !
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! Crude approximation for pressure (but commonly used)
-            ! 1e-04 to convert from Pa to dBar
-            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
-            !
-         END DO
+      CALL wrk_alloc( jpi,jpj, zfrz  , zgammat, zgammas  )
+      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
+
+      ! initialisation
+      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
+      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp    
+      zfwflx (:,:) = 0.0_wp
+
+      ! compute ice shelf melting
+      nit = 1 ; lit = .TRUE.
+      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
+         SELECT CASE ( nn_isfblk )
+         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
+            ! Calculate freezing temperature
+            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
+
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+            
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
+                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
+               END DO
+            END DO
+
+            ! Compute heat flux and upward fresh water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! compute coeficient to solve the 2nd order equation
+                  zeps1 = rcp*rau0*zgammat(ji,jj)
+                  zeps2 = rlfusisf*rau0*zgammas(ji,jj)
+                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
+                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
+                  zeps6 = zeps4-ttbl(ji,jj)
+                  zeps7 = zeps4-tsurf
+                  zaqe  = zlamb1 * (zeps1 + zeps3)
+                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
+                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
+                  zcqe  = zeps2*stbl(ji,jj)
+                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
+
+                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
+                  ! compute s freeze
+                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
+                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
+
+                  ! compute t freeze (eq. 22)
+                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
+  
+                  ! zfwflx is upward water flux
+                  ! zhtflx is upward heat flux (out of ocean)
+                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
+                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
+                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
+               END DO
+            END DO
+
+            ! compute heat and water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         END SELECT
+
+         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
+         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
+         ELSE                           
+            ! check total number of iteration
+            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+            ELSE                 ; nit = nit + 1
+            END IF
+
+            ! compute error between 2 iterations
+            ! if needed save gammat and compute zhtflx_b for next iteration
+            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
+            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
+            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
+            END IF
+         END IF
       END DO
-
-! Calculate in-situ temperature (ref to surface)
-      zti(:,:)=tinsitu( ttbl, stbl, zpress )
-! Calculate freezing temperature
-      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
-
-      
-      zhtflx=0._wp ; zfwflx=0._wp
-      IF (nn_isfblk == 1) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (mikt(ji,jj) > 1 ) THEN
-                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
-                  DO WHILE ( lit )
-! compute gamma
-                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
-! zhtflx is upward heat flux (out of ocean)
-                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
-! zwflx is upward water flux
-                     zfwflx = - zhtflx/lfusisf
-! test convergence and compute gammat
-                     IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
-
-                     nit = nit + 1
-                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-
-! save gammat and compute zhtflx_b
-                     zgammat2d(ji,jj)=zgammat
-                     zhtflx_b = zhtflx
-                  END DO
-
-                  qisf(ji,jj) = - zhtflx
-! For genuine ISOMIP protocol this should probably be something like
-                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
-               ELSE
-                  fwfisf(ji,jj) = 0._wp
-                  qisf(ji,jj)   = 0._wp
-               END IF
-            !
-            END DO
-         END DO
-
-      ELSE IF (nn_isfblk == 2 ) THEN
-
-! More complicated 3 equation thermodynamics as in MITgcm
-         DO jj = 2, jpj
-            DO ji = 2, jpi
-               IF (mikt(ji,jj) > 1 ) THEN
-                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
-                  DO WHILE ( lit )
-                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
-
-                     zeps1=rcp*rau0*zgammat
-                     zeps2=lfusisf*rau0*zgammas
-                     zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)
-                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
-                     zeps6=zeps4-zti(ji,jj)
-                     zeps7=zeps4-tsurf
-                     zaqe=zlamb1 * (zeps1 + zeps3)
-                     zaqer=0.5/MIN(zaqe,-zeps)
-                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
-                     zcqe=zeps2*stbl(ji,jj)
-                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
-! Presumably zdis can never be negative because gammas is very small compared to gammat
-                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
-                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
-                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
-  
-! zfwflx is upward water flux (MAX usefull if kappa = 0
-                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) )
-! zhtflx is upward heat flux (out of ocean)
-! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
-                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
-! zwflx is upward water flux
-! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
-                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
-! test convergence and compute gammat
-                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
-
-                     nit = nit + 1
-                     IF (nit .GE. 51) THEN
-                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
-                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
-! save gammat and compute zhtflx_b
-                     zgammat2d(ji,jj)=zgammat
-                     zgammas2d(ji,jj)=zgammas
-                     zhtflx_b = zhtflx
-
-                  END DO
-! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
-                  qisf(ji,jj) = - zhtflx 
-! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
-                  fwfisf(ji,jj) = zfwflx 
-               ELSE
-                  fwfisf(ji,jj) = 0._wp
-                  qisf(ji,jj)   = 0._wp
-               ENDIF
-               !
-            END DO
-         END DO
-      ENDIF
-      ! lbclnk
-      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
-      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
-      ! output
-      CALL iom_put('isfgammat', zgammat2d)
-      CALL iom_put('isfgammas', zgammas2d)
-      !
-      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
+      !
+      CALL iom_put('isfgammat', zgammat)
+      CALL iom_put('isfgammas', zgammas)
+      ! 
+      CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  )
+      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
@@ -562,6 +538,5 @@
    END SUBROUTINE sbc_isf_cav
 
-
-   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
+   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat flux
@@ -572,10 +547,10 @@
       !!                Jenkins et al., 2010, JPO, p2298-2312
       !!---------------------------------------------------------------------
-      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
-      INTEGER , INTENT(in)    :: ji,jj
-      LOGICAL , INTENT(inout) :: lit
-
-      INTEGER  :: ikt                 ! loop index
-      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
+      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs
+      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf
+      !
+      INTEGER  :: ikt                        
+      INTEGER  :: ji, jj                     ! loop index
+      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity
       REAL(wp) :: zdku, zdkv                 ! U, V shear 
       REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
@@ -587,169 +562,211 @@
       REAL(wp) :: zcoef                      ! temporary coef
       REAL(wp) :: zdep
-      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
-      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
-      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
-      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
+      REAL(wp) :: zeps = 1.0e-20_wp    
+      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
+      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
       REAL(wp), DIMENSION(2) :: zts, zab
       !!---------------------------------------------------------------------
-      !
-      IF( nn_gammablk == 0 ) THEN
-      !! gamma is constant (specified in namelist)
-         gt = rn_gammat0
-         gs = rn_gammas0
-         lit = .FALSE.
-      ELSE IF ( nn_gammablk == 1 ) THEN
-      !! gamma is assume to be proportional to u* 
-      !! WARNING in case of Losh 2008 tbl parametrization, 
-      !! you have to used the mean value of u in the boundary layer) 
-      !! not yet coded
-      !! Jenkins et al., 2010, JPO, p2298-2312
-         ikt = mikt(ji,jj)
-      !! Compute U and V at T points
-   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
-   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
-          zut = utbl(ji,jj)
-          zvt = vtbl(ji,jj)
-
-      !! compute ustar
-         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
-      !! Compute mean value over the TBL
-
-      !! Compute gammats
-         gt = zustar * rn_gammat0
-         gs = zustar * rn_gammas0
-         lit = .FALSE.
-      ELSE IF ( nn_gammablk == 2 ) THEN
-      !! gamma depends of stability of boundary layer
-      !! WARNING in case of Losh 2008 tbl parametrization, 
-      !! you have to used the mean value of u in the boundary layer) 
-      !! not yet coded
-      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
-      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+      CALL wrk_alloc( jpi,jpj, zustar )
+      !
+      SELECT CASE ( nn_gammablk )
+      CASE ( 0 ) ! gamma is constant (specified in namelist)
+         !! ISOMIP formulation (Hunter et al, 2006)
+         pgt(:,:) = rn_gammat0
+         pgs(:,:) = rn_gammas0
+
+      CASE ( 1 ) ! gamma is assume to be proportional to u*
+         !! Jenkins et al., 2010, JPO, p2298-2312
+         !! Adopted by Asay-Davis et al. (2015)
+
+         !! compute ustar (eq. 24)
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+
+         !! Compute gammats
+         pgt(:,:) = zustar(:,:) * rn_gammat0
+         pgs(:,:) = zustar(:,:) * rn_gammas0
+      
+      CASE ( 2 ) ! gamma depends of stability of boundary layer
+         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
+         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+         !! compute ustar
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+
+         !! compute Pr and Sc number (can be improved)
+         zPr =   13.8_wp
+         zSc = 2432.0_wp
+
+         !! compute gamma mole
+         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
+         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
+
+         !! compute gamma
+         DO ji=2,jpi
+            DO jj=2,jpj
                ikt = mikt(ji,jj)
 
-      !! Compute U and V at T points
-               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
-               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
-
-      !! compute ustar
-               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
-               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
-                 gt = rn_gammat0
-                 gs = rn_gammas0
+               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think
+                  pgt = rn_gammat0
+                  pgs = rn_gammas0
                ELSE
-      !! compute Rc number (as done in zdfric.F90)
-               zcoef = 0.5 / fse3w(ji,jj,ikt)
-               !                                            ! shear of horizontal velocity
-               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
-                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
-               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
-                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
-               !                                            ! richardson number (minimum value set to zero)
-               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
-
-      !! compute Pr and Sc number (can be improved)
-               zPr =   13.8
-               zSc = 2432.0
-
-      !! compute gamma mole
-               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
-               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
-
-      !! compute bouyancy 
-               zts(jp_tem) = ttbl(ji,jj)
-               zts(jp_sal) = stbl(ji,jj)
-               zdep        = fsdepw(ji,jj,ikt)
-               !
-               CALL eos_rab( zts, zdep, zab )
+                  !! compute Rc number (as done in zdfric.F90)
+                  zcoef = 0.5_wp / fse3w(ji,jj,ikt)
+                  !                                            ! shear of horizontal velocity
+                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
+                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
+                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
+                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+                  !                                            ! richardson number (minimum value set to zero)
+                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
+
+                  !! compute bouyancy 
+                  zts(jp_tem) = ttbl(ji,jj)
+                  zts(jp_sal) = stbl(ji,jj)
+                  zdep        = fsdepw(ji,jj,ikt)
                   !
-      !! compute length scale 
-               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      !! compute Monin Obukov Length
-               ! Maximum boundary layer depth
-               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
-               ! Compute Monin obukhov length scale at the surface and Ekman depth:
-               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
-               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
-
-      !! compute eta* (stability parameter)
-               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
-
-      !! compute the sublayer thickness
-               zhnu = 5 * znu / zustar
-      !! compute gamma turb
-               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
-               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
-
-      !! compute gammats
-               gt = zustar / (zgturb + zgmolet)
-               gs = zustar / (zgturb + zgmoles)
+                  CALL eos_rab( zts, zdep, zab )
+                  !
+                  !! compute length scale 
+                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+                  !! compute Monin Obukov Length
+                  ! Maximum boundary layer depth
+                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp
+                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
+                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+
+                  !! compute eta* (stability parameter)
+                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp)))
+
+                  !! compute the sublayer thickness
+                  zhnu = 5 * znu / zustar(ji,jj)
+
+                  !! compute gamma turb
+                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
+                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
+
+                  !! compute gammats
+                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
                END IF
-      END IF
-      !
-   END SUBROUTINE
-
-
-   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
+            END DO
+         END DO
+         CALL lbc_lnk(pgt(:,:),'T',1.)
+         CALL lbc_lnk(pgs(:,:),'T',1.)
+      END SELECT
+      CALL wrk_dealloc( jpi,jpj, zustar )
+      !
+   END SUBROUTINE sbc_isf_gammats
+
+   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE sbc_isf_tbl  ***
       !!
-      !! ** Purpose : compute mean T/S/U/V in the boundary layer 
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
-      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
+      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin
+      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout
+      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out
+      !
+      REAL(wp) :: ze3, zhk
+      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl
+
+      INTEGER :: ji, jj, jk                  ! loop index
+      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl
+      !!----------------------------------------------------------------------
+      ! allocation
+      CALL wrk_alloc( jpi,jpj, zhisf_tbl)
       
-      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
-
-      REAL(wp) :: ze3, zhk
-      REAL(wp), DIMENSION(:,:), POINTER :: zikt
-
-      INTEGER :: ji,jj,jk
-      INTEGER :: ikt,ikb
-      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
-
-      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
-      CALL wrk_alloc( jpi,jpj, zikt )
-
-      ! get first and last level of tbl
-      mkt(:,:) = misfkt(:,:)
-      mkb(:,:) = misfkb(:,:)
-
-      varout(:,:)=0._wp
-      DO jj = 2,jpj
-         DO ji = 2,jpi
-            IF (ssmask(ji,jj) == 1) THEN
-               ikt = mkt(ji,jj)
-               ikb = mkb(ji,jj)
+      ! initialisation
+      pvarout(:,:)=0._wp
+   
+      SELECT CASE ( cd_ptin )
+      CASE ( 'U' ) ! compute U in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt+1, mbku(ji,jj)
+                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3u_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+      
+      CASE ( 'V' ) ! compute V in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt+1, mbkv(ji,jj)
+                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3v_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+
+      CASE ( 'T' ) ! compute T in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = misfkt(ji,jj)
+               ikb = misfkb(ji,jj)
 
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
                   ze3 = fse3t_n(ji,jj,jk)
-                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
-                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
-                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
-                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
-                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
                zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
-               IF (cptin == 'T') &
-                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
-               IF (cptin == 'U') &
-                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
-               IF (cptin == 'V') &
-                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
-            END IF
-         END DO
-      END DO
-
-      CALL wrk_dealloc( jpi,jpj, mkt, mkb )      
-      CALL wrk_dealloc( jpi,jpj, zikt ) 
-
-      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
-      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+      END SELECT
+
+      ! mask mean tbl value
+      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
+
+      ! deallocation
+      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )      
       !
    END SUBROUTINE sbc_isf_tbl
@@ -768,16 +785,15 @@
       !! ** Action  :   phdivn   decreased by the runoff inflow
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
+      ! 
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt, ikb 
-      INTEGER  ::   nk_isf
-      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
-      REAL(wp)     ::   zfact     ! local scalar
+      REAL(wp) ::   zhk
+      REAL(wp) ::   zfact     ! local scalar
       !!----------------------------------------------------------------------
       !
       zfact   = 0.5_wp
       !
-      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
+      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -788,7 +804,6 @@
 
                ! determine the deepest level influenced by the boundary layer
-               ! test on tmask useless ?????
-               DO jk = ikt, mbkt(ji,jj)
-                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+               DO jk = ikt+1, mbkt(ji,jj)
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) <  rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -796,10 +811,11 @@
                r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
 
-               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
+               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
                ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
             END DO
          END DO
-      END IF ! vvl case
-      !
+      END IF 
+      !
+      !==   ice shelf melting distributed over several levels   ==!
       DO jj = 1,jpj
          DO ji = 1,jpi
@@ -809,91 +825,13 @@
                DO jk = ikt, ikb - 1
                   phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
-                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
+                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
                END DO
                ! level partially include in ice shelf boundary layer 
                phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
-                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
-            !==   ice shelf melting mass distributed over several levels   ==!
+                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
          END DO
       END DO
       !
    END SUBROUTINE sbc_isf_div
-
-
-   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_init  ***
-      !!
-      !! ** Purpose :   Compute the in-situ temperature [Celcius]
-      !!
-      !! ** Method  :   
-      !!
-      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
-      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
-!      REAL(wp) :: fsatg
-!      REAL(wp) :: pfps, pfpt, pfphp 
-      REAL(wp) :: zt, zs, zp, zh, zq, zxk
-      INTEGER  :: ji, jj            ! dummy loop indices
-      !
-      CALL wrk_alloc( jpi,jpj, pti  )
-      ! 
-      DO jj=1,jpj
-         DO ji=1,jpi
-            zh = ppress(ji,jj)
-! Theta1
-            zt = ptem(ji,jj)
-            zs = psal(ji,jj)
-            zp = 0.0
-            zxk= zh * fsatg( zs, zt, zp )
-            zt = zt + 0.5 * zxk
-            zq = zxk
-! Theta2
-            zp = zp + 0.5 * zh
-            zxk= zh*fsatg( zs, zt, zp )
-            zt = zt + 0.29289322 * ( zxk - zq )
-            zq = 0.58578644 * zxk + 0.121320344 * zq
-! Theta3
-            zxk= zh * fsatg( zs, zt, zp )
-            zt = zt + 1.707106781 * ( zxk - zq )
-            zq = 3.414213562 * zxk - 4.121320344 * zq
-! Theta4
-            zp = zp + 0.5 * zh
-            zxk= zh * fsatg( zs, zt, zp )
-            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
-         END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj, pti  )
-      !
-   END FUNCTION tinsitu
-
-
-   FUNCTION fsatg( pfps, pfpt, pfphp )
-      !!----------------------------------------------------------------------
-      !!                 ***  FUNCTION fsatg  ***
-      !!
-      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
-      !!
-      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
-      !! 
-      !! ** units      :   pressure        pfphp    decibars
-      !!                   temperature     pfpt     deg celsius (ipts-68)
-      !!                   salinity        pfps     (ipss-78)
-      !!                   adiabatic       fsatg    deg. c/decibar
-      !!----------------------------------------------------------------------
-      REAL(wp) :: pfps, pfpt, pfphp 
-      REAL(wp) :: fsatg
-      !
-      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
-        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
-        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
-        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
-        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
-      !
-    END FUNCTION fsatg
-    !!======================================================================
+   !!======================================================================
 END MODULE sbcisf
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6012)
@@ -90,5 +90,5 @@
       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
-         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
+         &             ln_ssr    , ln_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
          &             nn_lsm    , nn_limflx , nn_components, ln_cpl
       INTEGER  ::   ios
@@ -143,5 +143,5 @@
          WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
          WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
-         WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf
+         WRITE(numout,*) '              iceshelf melting                           ln_isf      = ', ln_isf
          WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr
          WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb
@@ -181,10 +181,10 @@
 
       !                          ! Checks:
-      IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf 
-         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_isf arrays' )
+      IF( .NOT. ln_isf ) THEN                      ! variable initialisation if no ice shelf 
+         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
          fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
          risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
-         rdivisf       = 0.0_wp
       END IF
+
       IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero
 
@@ -378,5 +378,5 @@
       IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs
 
-      IF( nn_isf   /= 0  )   CALL sbc_isf( kt )                    ! compute iceshelves
+      IF( ln_isf         )   CALL sbc_isf( kt )                   ! compute iceshelves
 
       IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 6012)
@@ -142,10 +142,10 @@
          IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
             rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+            CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
                rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
             END WHERE
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-               ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study)
-               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp
+               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp
             END WHERE
          ELSE                                                        ! use SST as runoffs temperature
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 6012)
@@ -243,6 +243,6 @@
       ENDIF
       IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities
-         IF(  ln_traadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF
-            & ln_traadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' )
+         IF(  ln_traadv_cen .AND. nn_cen_v == 4    .OR.   &                            ! NO 4th order with ISF
+            & ln_traadv_fct .AND. nn_fct_v == 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' )
       ENDIF
       !
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 6012)
@@ -103,4 +103,5 @@
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  ::  ikt
       INTEGER  ::  ierr             ! local integer
       REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
@@ -228,5 +229,4 @@
             DO jj = 1, jpjm1              ! bottom correction (partial bottom cell)
                DO ji = 1, fs_jpim1   ! vector opt.
-!!gm  the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere
                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)          
                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
@@ -255,15 +255,4 @@
             ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
             ENDIF
-!!gm I don't understand why we should need this.... since wmask is used instead of tmask
-!         IF ( ln_isfcav ) THEN
-!            DO jj = 1, jpj
-!               DO ji = 1, jpi   ! vector opt.
-!                  ikt = mikt(ji,jj) ! surface level
-!                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)
-!                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)
-!               END DO
-!            END DO
-!         END IF
-!!gm 
 
             DO jj = 1 , jpjm1            !==  Horizontal fluxes
@@ -272,9 +261,9 @@
                   zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   !
-                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
-                     &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
-                  !
-                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
-                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
+                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
+                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. )
+                  !
+                  zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
+                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. )
                   !
                   zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
@@ -317,7 +306,7 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   !
-                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
+                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
-                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
+                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
                      !
@@ -343,5 +332,5 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
@@ -358,5 +347,5 @@
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
-                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)
+                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk) / fse3w(ji,jj,jk)
                      END DO
                   END DO
@@ -366,5 +355,5 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk)                      &
                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 6012)
@@ -28,5 +28,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbcrnf          ! river runoffs
-   USE sbcisf          ! ice shelf melting/freezing
+   USE sbcisf          ! ice shelf melting
    USE zdf_oce         ! ocean vertical mixing
    USE domvvl          ! variable volume
@@ -262,13 +262,9 @@
          ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
          ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
-         IF (nn_isf .GE. 1) THEN 
-            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
-         ELSE
-            ll_isf = .FALSE.
-         END IF
+         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
       ELSE                          
          ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
          ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
-         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
+         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting
       ENDIF
       !
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6012)
@@ -119,5 +119,5 @@
       INTEGER  ::   ikt, ikb 
       REAL(wp) ::   zfact, z1_e3t, zdep
-      REAL(wp) ::   zalpha, zhk
+      REAL(wp) ::   zt_frz, zpress
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -218,5 +218,5 @@
       !----------------------------------------
       !
-      IF( nn_isf > 0 ) THEN
+      IF( ln_isf ) THEN
          zfact = 0.5_wp
          DO jj = 2, jpj
@@ -227,20 +227,18 @@
    
                ! level fully include in the ice shelf boundary layer
-               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
                ! sign - because fwf sign of evapo (rnf sign of precip)
                DO jk = ikt, ikb - 1
                ! 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)) * r1_hisf_tbl(ji,jj)
-                  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)
+                  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) )             &
+                     &           * r1_hisf_tbl(ji,jj)
                END DO
    
                ! level partially include in ice shelf boundary layer 
                ! 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)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
-               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) * ralpha(ji,jj) 
+               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) )             &
+                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
+
             END DO
          END DO
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 6012)
@@ -191,20 +191,19 @@
       !
    END SUBROUTINE zps_hde
-
-
-   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi,                                   &
-      &                              prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv ,   &
-      &                                   pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE zps_hde  ***
+   !
+   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
+      &                          prd, pgru, pgrv, pgrui, pgrvi )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE zps_hde_isf  ***
       !!                    
       !! ** Purpose :   Compute the horizontal derivative of T, S and rho
       !!      at u- and v-points with a linear interpolation for z-coordinate
-      !!      with partial steps.
+      !!      with partial steps for top (ice shelf) and bottom.
       !!
       !! ** Method  :   In z-coord with partial steps, scale factors on last 
       !!      levels are different for each grid point, so that T, S and rd 
       !!      points are not at the same depth as in z-coord. To have horizontal
-      !!      gradients again, we interpolate T and S at the good depth : 
+      !!      gradients again, we interpolate T and S at the good depth :
+      !!      For the bottom case:
       !!      Linear interpolation of T, S   
       !!         Computation of di(tb) and dj(tb) by vertical interpolation:
@@ -235,31 +234,22 @@
       !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
       !!
+      !!      For the top case (ice shelf): As for the bottom case but upside down
+      !!
       !! ** Action  : compute for top and bottom interfaces
       !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points
       !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points
-      !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl)
-      !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl)
-      !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   )           ::  kt                ! ocean time-step index
-      INTEGER                              , INTENT(in   )           ::  kjpt              ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta               ! 4D tracers fields
-      !                                                              !!  u-point ! v-point !
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu    , pgtv    ! bottom GRADh( ptra )  
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui   , pgtvi   ! top    GRADh( ptra )
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd               ! 3D density anomaly fields
-      !                                                              !!  u-point ! v-point !
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru    , pgrv    ! bottom GRADh( prd  )
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru    , pmrv    ! bottom SUM  ( prd  )
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu    , pgzv    ! bottom GRADh( z    ) 
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru  , pge3rv  ! bottom GRADh( prd  ) weighted by e3w
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui   , pgrvi   ! top    GRADh( prd  ) 
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui   , pmrvi   ! top    SUM  ( prd  ) 
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui   , pgzvi   ! top    GRADh( z    ) 
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui , pge3rvi ! top    GRADh( prd  ) weighted by e3w
+      !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts 
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF)
+      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 ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top)
       !
       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
+      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             ! 
@@ -277,45 +267,143 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
-               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
+
+               iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
+               ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
+               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku)
+               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv)
+               !
+               ! i- direction
+               IF( ze3wu >= 0._wp ) THEN      ! case 1
+                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
+                  ! interpolated values of tracers
+                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of  tracers
+                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
+               ELSE                           ! case 2
+                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
+                  ! interpolated values of tracers
+                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
+               ENDIF
+               !
+               ! j- direction
+               IF( ze3wv >= 0._wp ) THEN      ! case 1
+                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
+                  ! interpolated values of tracers
+                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
+               ELSE                           ! case 2
+                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
+                  ! interpolated values of tracers
+                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
+               ENDIF
+
+            END DO
+         END DO
+         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
+         !
+      END DO
+
+      ! horizontal derivative of density anomalies (rd)
+      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
+         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
+         !
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               iku = mbku(ji,jj)
+               ikv = mbkv(ji,jj)
+               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku)
+               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv)
+               !
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)    ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)    ! -     -      case 2
+               ENDIF
+
+            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
+         CALL eos( zti, zhi, zri )
+         CALL eos( ztj, zhj, zrj )
+
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
+            DO ji = 1, jpim1
+
+               iku = mbku(ji,jj)
+               ikv = mbkv(ji,jj)
+               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku)
+               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv)
+
+               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
+               ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
+               ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
+               ENDIF
+
+            END DO
+         END DO
+
+         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
+         !
+      END IF
+      !
+      !     !==  (ISH)  compute grui and gruvi  ==!
+      !
+      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1
+               ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1
+               !
                ! (ISF) case partial step top and bottom in adjacent cell in vertical
                ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
-               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-               !
+               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku)
+               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 
+
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
-                  ! interpolated values of tracers
-                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
+                  zmaxu = ze3wu / fse3w(ji+1,jj,ikup1)
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
+               ELSE                           ! case 2
+                  zmaxu = - ze3wu / fse3w(ji,jj,ikup1)
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
                   ! gradient of  tracers
-                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
-               ELSE                           ! case 2
-                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
-                  ! interpolated values of tracers
-                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
-                  ! gradient of tracers
-                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
+                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
                ENDIF
                !
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
-                  ! interpolated values of tracers
-                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
-                  ! gradient of tracers
-                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
-               ELSE                           ! case 2
-                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
-                  ! interpolated values of tracers
-                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
-                  ! gradient of tracers
-                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
+                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikvp1)
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
+               ELSE                           ! case 2
+                  zmaxv =  - ze3wv / fse3w(ji,jj,ikvp1)
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
+               ENDIF
+
+            END DO
+         END DO
+         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
          !
       END DO
@@ -323,22 +411,21 @@
       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
          !
-         pgru  (:,:)=0._wp   ;   pgrv  (:,:) = 0._wp
-         pgzu  (:,:)=0._wp   ;   pgzv  (:,:) = 0._wp 
-         pmru  (:,:)=0._wp   ;   pmru  (:,:) = 0._wp 
-         pge3ru(:,:)=0._wp   ;   pge3rv(:,:) = 0._wp 
-         !
-         DO jj = 1, jpjm1                 ! depth of the partial step level
-            DO ji = 1, jpim1
-               iku = mbku(ji,jj)
-               ikv = mbkv(ji,jj)
-               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
-               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-               !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2
-               ENDIF
+         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp;
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               iku = miku(ji,jj)
+               ikv = mikv(ji,jj)
+               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku)
+               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 
+               !
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)    ! -     -      case 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)    ! -     -      case 2
+               ENDIF
+
             END DO
          END DO
@@ -346,168 +433,24 @@
          CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
-
+         !
          DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
-               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
-               ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points
-               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
-               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-               IF( ze3wu >= 0._wp ) THEN 
-                  pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku)
-                  pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
-                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1 
-                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
-                                * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) &
-                                   - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
-               ELSE  
-                  pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu)
-                  pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2
-                  pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2
-                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
-                                * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) &
-                                   -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN
-                  pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 
-                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1
-                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1
-                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
-                                * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) &
-                                   - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
-               ELSE 
-                  pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv)
-                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2
-                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2
-                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
-                                * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) &
-                                   -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions
-         !
-      END IF
-      !
-      !     !==  (ISH)  compute grui and gruvi  ==!
-      !
-      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1
-               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1
-               !
-               ! (ISF) case partial step top and bottom in adjacent cell in vertical
-               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
-               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
-               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
-               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-               ! i- direction
-               IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1)
-                  ! interpolated values of tracers
-                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) )
-                  ! gradient of tracers
-                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
-               ELSE                           ! case 2
-                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1)
-                  ! interpolated values of tracers
-                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) )
-                  ! gradient of  tracers
-                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
-               ENDIF
-               !
-               ! j- direction
-               IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1)
-                  ! interpolated values of tracers
-                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) )
-                  ! gradient of tracers
-                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
-               ELSE                           ! case 2
-                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1)
-                  ! interpolated values of tracers
-                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) )
-                  ! gradient of tracers
-                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
-               ENDIF
-            END DO!!
-         END DO!!
-         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
-         !
-      END DO
-
-      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
-         !
-         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
-         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ;
-         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ;
-         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ;
-         !
-         DO jj = 1, jpjm1        ! depth of the partial step level
-            DO ji = 1, jpim1
-               iku = miku(ji,jj)
-               ikv = mikv(ji,jj)
-               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-               !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
-               ENDIF
-            END DO
-         END DO
-         !
-         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
-         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
-         !
-         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
-            DO ji = 1, jpim1
-               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
-               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1
-               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
-               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-               IF( ze3wu >= 0._wp ) THEN
-                 pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)
-                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1
-                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1 
-                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  &
-                    &           * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
-                    &              - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
-               ELSE
-                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
-                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2
-                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
-                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   &
-                    &           * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
-                    &              -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN
-                 pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 
-                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1
-                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
-                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  & 
-                     &           * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
-                     &                - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
-                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco
-               ELSE
-                 pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)
-                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2
-                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
-                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   &
-                    &           * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
-                    &              -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions
+
+               iku = miku(ji,jj) 
+               ikv = mikv(ji,jj) 
+               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku)
+               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv) 
+
+               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
+               ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1
+               ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2
+               ENDIF
+
+            END DO
+         END DO
+         CALL lbc_lnk( pgrui   , 'U', -1. ); CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions
          !
       END IF  
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6012)
@@ -31,5 +31,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer
 
    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
@@ -111,11 +111,10 @@
       END DO
       !
-      ! w-level of the turbocline
+      ! w-level of the turbocline and mixing layer (iom_use)
       imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
       DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
          DO jj = 1, jpj
             DO ji = 1, jpi
-               imkt = mikt(ji,jj)
-               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline 
+               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline 
             END DO
          END DO
@@ -127,12 +126,18 @@
             iikn = nmln(ji,jj)
             imkt = mikt(ji,jj)
-            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth 
-            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth
-            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
+            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO
       IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
-         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
-         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
+         IF ( iom_use("mldr10_1") ) THEN
+            IF( .NOT. ln_isfcav ) CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
+            IF(       ln_isfcav ) CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
+         END IF
+         IF ( iom_use("mldkz5") ) THEN
+            IF( .NOT. ln_isfcav ) CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
+            IF(       ln_isfcav ) CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
+         END IF
       ENDIF
       
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 6012)
@@ -60,7 +60,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtsu, gtsv   !: horizontal gradient of T, S bottom u-point 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aru , arv    
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzu , gzv    
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3ru, ge3rv   !: horizontal gradient of T, S and rd at top v-point  
 
    !! (ISF) interpolated gradient (only used for ice shelf case) 
@@ -68,7 +65,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtui, gtvi   !: horizontal gradient of T, S and rd at top u-point 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   grui, grvi   !: horizontal gradient of T, S and rd at top v-point  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   arui, arvi   !: horizontal average  of rd          at top v-point  
-   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
@@ -115,12 +109,6 @@
          &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
          &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
-         &     aru(jpi,jpj)      , arv(jpi,jpj)      ,                     &
-         &     gzu(jpi,jpj)      , gzv(jpi,jpj)      ,                     &
          &     gru(jpi,jpj)      , grv(jpi,jpj)      ,                     &
-         &     ge3ru(jpi,jpj)    , ge3rv(jpi,jpj)    ,                     &
          &     gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts),                     &
-         &     arui(jpi,jpj)     , arvi(jpi,jpj)     ,                     &
-         &     gzui(jpi,jpj)     , gzvi(jpi,jpj)     ,                     &
-         &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     &
          &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     &
          &     riceload(jpi,jpj),                             STAT=ierr(2) )
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6012)
@@ -166,8 +166,6 @@
 
          IF( ln_zps .AND.       ln_isfcav)                               &
-            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,   &    ! Partial steps for top cell (ISF)
-            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-            &                                   grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
-
+            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
          IF( ln_traldf_triad ) THEN 
                          CALL ldf_slp_triad( kstp )                       ! before slope for triad operator
@@ -186,5 +184,4 @@
       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
                          CALL wzv           ( kstp )  ! now cross-level velocity 
-
 !!gm : why also here ????
       IF(ln_sto_eos  )   CALL sto_pts( tsn )                             ! Random T/S fluctuations
@@ -196,11 +193,10 @@
 !!                                         but ensures reproductible results
 !!                                         with previous versions using split-explicit free surface          
-            IF( ln_zps .AND. .NOT. ln_isfcav)   &                           ! Partial steps: bottom before horizontal gradient
-               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &  ! of t, s, rd at the last ocean level
-               &                                          rhd, gru , grv    )
-            IF( ln_zps .AND.       ln_isfcav)   &                           ! Partial steps: top & bottom before horizontal gradient
-               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,   & 
-               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                               grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )
+            IF( ln_zps .AND. .NOT. ln_isfcav )                               &
+               &            CALL zps_hde    ( kstp, jpts, tsn, 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 )                                          &
+               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
 !!jc: fs simplification
                             
Index: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 6011)
+++ branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 6012)
@@ -246,5 +246,4 @@
    END SUBROUTINE trc_ini_state
 
-
    SUBROUTINE top_alloc
       !!----------------------------------------------------------------------
