Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 11931)
@@ -26,9 +26,7 @@
    !!-------------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
-   USE isf
    USE dom_oce        ! ocean space and time domain
    USE sbc_oce        ! Surface boundary condition: ocean fields
    USE sbcrnf         ! river runoffs
-   USE isfnxt
    USE phycst         ! physical constants
    USE dynadv         ! dynamics: vector invariant versus flux form
@@ -42,4 +40,6 @@
    USE trddyn         ! trend manager: dynamics
    USE trdken         ! trend manager: kinetic energy
+   USE isf       , ONLY: ln_isf     ! ice shelf
+   USE isfdynnxt , ONLY: isf_dynnxt ! ice shelf 
    !
    USE in_out_manager ! I/O manager
@@ -243,5 +243,7 @@
             END IF
             !
-            ! ice shelf melting
+            ! ice shelf melting (deal separatly as it can be in depth)
+            ! PM: we could probably define a generic subroutine to do the in depth correction
+            !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
             IF ( ln_isf ) CALL isf_dynnxt( kt, atfp * rdt )
             !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90	(revision 11931)
@@ -260,8 +260,9 @@
                &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   &
                &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
+
+            ! ice sheet coupling
+            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) sshb(:,:) = sshb(:,:) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
+
          ENDIF
-
-         ! ice sheet coupling
-         IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) sshb(:,:) = sshb(:,:) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
 
          sshn(:,:) = ssha(:,:)                              ! now <-- after
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11931)
@@ -14,5 +14,6 @@
    !!----------------------------------------------------------------------
 
-   USE in_out_manager, ONLY: wp, jpi,jpj, jpk, jpts ! I/O manager
+   USE par_oce       , ONLY: jpi, jpj, jpk
+   USE in_out_manager, ONLY: wp, jpts ! I/O manager
    USE lib_mpp       , ONLY: ctl_stop, mpp_sum      ! MPP library
    USE fldread        ! read input fields
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11931)
@@ -23,5 +23,5 @@
    !
    USE oce      , ONLY: tsn                    ! ocean tracers
-   USE dom_oce  , ONLY: jpi,jpj                ! ocean space and time domain
+   USE par_oce  , ONLY: jpi,jpj                ! ocean space and time domain
    USE phycst   , ONLY: grav,rau0,r1_rau0_rcp  ! physical constants
    USE eosbn2   , ONLY: l_useCT                ! l_useCT
@@ -90,5 +90,5 @@
       DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
          !
-         ! compute gammat every where (2d)
+         ! compute gammat everywhere (2d)
          ! useless if melt specified
          IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
@@ -101,5 +101,5 @@
             &                         zqhc   , zqoce, pqfwf  )
          !
-         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
+         ! define if we need to iterate
          SELECT CASE ( cn_gammablk )
          CASE ( 'spe','ad15' )
@@ -123,11 +123,11 @@
       END DO
       !
-      ! compute heat and water flux  (change signe directly in the melt subroutine)
+      ! compute heat and water flux ( > 0 out )
       pqfwf(:,:) = pqfwf(:,:) * mskisf_cav(:,:)
       zqoce(:,:) = zqoce(:,:) * mskisf_cav(:,:)
       zqhc (:,:) = zqhc(:,:)  * mskisf_cav(:,:)
       !
-      ! compute heat content flux
-      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2) ( > 0 out )
+      ! compute heat content flux ( > 0 out )
+      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
       !
       ! total heat flux ( >0 out )
@@ -163,10 +163,14 @@
       INTEGER :: ierr
       !!---------------------------------------------------------------------
-
+      !
+      !==============
       ! 0: allocation
+      !==============
       !
       CALL isf_alloc_cav()
       !
+      !==================
       ! 1: initialisation
+      !==================
       !
       ! top and bottom level of the 'top boundary layer'
@@ -179,10 +183,14 @@
       mskisf_cav(:,:)    = (1._wp - tmask(:,:,1)) * ssmask(:,:)
       !
+      !================
       ! 2: read restart
+      !================
       !
       ! read cav variable from restart
       IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
       !
+      !==========================================
       ! 3: specific allocation and initialisation (depending of scheme choice)
+      !==========================================
       !
       SELECT CASE ( TRIM(cn_isfcav_mlt) )
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 11931)
@@ -57,5 +57,8 @@
       !!---------------------------------------------------------------------
       !
-      ! compute velocity in the tbl if needed
+      !==========================================
+      ! 1.: compute velocity in the tbl if needed
+      !==========================================
+      !
       SELECT CASE ( cn_gammablk )
       CASE ( 'spe'  ) 
@@ -78,5 +81,8 @@
       END SELECT
       ! 
-      ! compute gamma
+      !==========================================
+      ! 2.: compute gamma
+      !==========================================
+      !
       SELECT CASE ( cn_gammablk )
       CASE ( 'spe'  ) ! gamma is constant (specified in namelist)
@@ -91,5 +97,8 @@
       END SELECT
       !
-      ! ouput exchange coeficient and tbl velocity
+      !==========================================
+      ! 3.: output and debug
+      !==========================================
+      !
       CALL iom_put('isfgammat', pgt(:,:))
       CALL iom_put('isfgammas', pgs(:,:))
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 11931)
@@ -3,5 +3,5 @@
    !!                       ***  MODULE  isfcavmlt  ***
    !! ice shelf module :  update surface ocean boundary condition under ice
-   !!                   shelf
+   !!                   shelves
    !!======================================================================
    !! History :  4.0  !  2019-09  (P. Mathiot) Original code
@@ -9,5 +9,5 @@
 
    !!----------------------------------------------------------------------
-   !!   isfcav_mlt    : update surface ocean boundary condition under ice shelf
+   !!   isfcav_mlt    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
    !!----------------------------------------------------------------------
 
@@ -23,5 +23,5 @@
    USE in_out_manager              ! I/O manager
    USE iom        , ONLY: iom_put  ! I/O library
-   USE fldread    , ONLY: fld_read !
+   USE fldread    , ONLY: fld_read, FLD, FLD_N !
    USE lib_fortran, ONLY: glob_sum !
    USE lib_mpp    , ONLY: ctl_stop !
@@ -51,5 +51,4 @@
       !! ** Purpose    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
       !!
-      !!---------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
@@ -58,5 +57,4 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
-      !!---------------------------------------------------------------------
       !!---------------------------------------------------------------------
       !
@@ -110,5 +108,5 @@
       !!--------------------------------------------------------------------
       !
-      ! Calculate freezing temperature
+      ! Compute freezing temperature
       CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
       !
@@ -131,5 +129,5 @@
       !!----------------------------------------------------------------------
       !!
-      !!                          ***  ROUTINE isfcav_mlt_spe  ***
+      !!                          ***  ROUTINE isfcav_mlt_2eq  ***
       !!
       !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
@@ -144,5 +142,5 @@
       !!                 Tech.  Rep.  June,  Antarctic  Climate  &  Ecosystems  Cooperative  Research  Centre,  available  at:  
       !!                 http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
-      !!---------------------------------------------------------------------
+      !!
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
@@ -192,5 +190,5 @@
       !!                   MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), 
       !!                   Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016. 
-      !!---------------------------------------------------------------------
+      !!
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
@@ -308,4 +306,5 @@
       !
       CALL iom_put('isftfrz_cav', ztfrz * mskisf_cav(:,:) )
+      !
    END SUBROUTINE isfcav_mlt_oasis
 
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90	(revision 11931)
@@ -12,8 +12,8 @@
    !!   isfrst : read/write iceshelf variables in/from restart
    !!----------------------------------------------------------------------
-   USE isf                             ! ice shelf variable
+   USE isf                              ! ice shelf variable
    USE isfutils, ONLY : debug
-   USE lib_mpp, ONLY: mpp_sum, mpp_max ! mpp routine
-   USE domvvl , ONLY: dom_vvl_zgr      ! vertical scale factor interpolation
+   USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
+   USE domvvl  , ONLY: dom_vvl_zgr      ! vertical scale factor interpolation
    !
    USE oce            ! ocean dynamics and tracers
@@ -113,10 +113,10 @@
    SUBROUTINE isfcpl_rst_write(kt)
       !!---------------------------------------------------------------------
-      !!   isfrst_cpl_write : write icesheet coupling variables in restart
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT --------------------------------------
+      !!                   ***  ROUTINE iscpl_rst_write  ***
+      !! 
+      !! ** Purpose : write icesheet coupling variables in restart
+      !!
       !!-------------------------- IN  --------------------------------------
       INTEGER, INTENT(in) :: kt
-      !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
       !
@@ -136,11 +136,8 @@
       !!                   ***  ROUTINE iscpl_ssh  ***
       !! 
-      !! ** Purpose :   basic guess of ssh in new wet cell during coupling step
+      !! ** Purpose :   basic guess of ssh in new wet cell
       !! 
       !! ** Method  :   basic extrapolation from neigbourg cells
       !!
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT --------------------------------------
-      !!-------------------------- IN  --------------------------------------
       !!----------------------------------------------------------------------
       !!
@@ -344,6 +341,4 @@
       tsn(:,:,:,jp_sal) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
       !
-      tsb(:,:,:,:) = tsn(:,:,:,:)
-      !
       ! sanity check
       ! -----------------------------------------------------------------------------------------
@@ -366,8 +361,10 @@
       !!                   ***  ROUTINE iscpl_vol  ***
       !! 
-      !! ** Purpose :   
-      !!                
-      !! 
-      !! ** Method  :   
+      !! ** Purpose : compute the correction of the local divergence to apply  
+      !!              during the first time step after the coupling.
+      !!
+      !! ** Method  : - compute horizontal vol div. before/after coupling
+      !!              - compute vertical input
+      !!              - compute correction
       !!                
       !!----------------------------------------------------------------------
@@ -376,15 +373,17 @@
       INTEGER :: ikb, ikt
       !!
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b
-      !!----------------------------------------------------------------------
-      !
-      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln  ! vol flux div.         before/after coupling
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b  ! vertical scale factor before/after coupling
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b        ! mask                  before       coupling
+      !!----------------------------------------------------------------------
+      !
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b  , ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b  , ldxios = lrxios )
       !
-      ! get volume flux before coupling (>0 out)
+      ! 1.0: compute horizontal volume flux divergence difference before-after coupling
+      !
       DO jk = 1, jpk                                 ! Horizontal slab
+         ! 1.1: get volume flux before coupling (>0 out)
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -395,10 +394,10 @@
          ENDDO
          !
+         ! 1.2: get volume flux after coupling (>0 out)
          ! properly mask velocity 
          ! (velocity are still mask with old mask at this stage)
          un(:,:,jk) = un(:,:,jk) * umask(:,:,jk)
          vn(:,:,jk) = vn(:,:,jk) * vmask(:,:,jk)
-         !
-         ! get volume flux after coupling (>0 out)
+         ! compute volume flux divergence after coupling
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -409,10 +408,11 @@
          ENDDO
          !
-         ! get 3d volume flux difference (before - after cpl) (>0 out)
-         ! correction to add is _b - _n
+         ! 1.3: get 3d volume flux difference (before - after cpl) (>0 out)
+         !      correction to add is _b - _n
          risfcpl_vol(:,:,jk) = zqvolb(:,:,jk) - zqvoln(:,:,jk)
       END DO
       !
-      ! include the contribution of the vertical velocity in the volume flux correction
+      ! 2.0: include the contribution of the vertical velocity in the volume flux correction
+      !
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
@@ -426,10 +426,12 @@
       ENDDO
       !
-      ! mask volume flux divergence correction
+      CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
+      !
+      ! 3.0: set total correction (div, tra, ssh)
+      !
+      ! 3.1: mask volume flux divergence correction
       risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
       !
-      CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
-      !
-      ! get 3d tra increment to apply at the first time step
+      ! 3.2: get 3d tra increment to apply at the first time step
       ! temperature and salt content flux computed using local tsn 
       ! (very simple advection scheme)
@@ -438,5 +440,5 @@
       risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_sal)
       !
-      ! ssh correction (for dynspg_ts)
+      ! 3.3: ssh correction (for dynspg_ts)
       risfcpl_ssh(:,:) = 0.0
       DO jk = 1,jpk
@@ -459,23 +461,32 @@
       !!----------------------------------------------------------------------
       !
-      TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts 
-      !
-      INTEGER  ::   ji   , jj  , jk           ! loop index
-      INTEGER  ::   jip1 , jim1, jjp1, jjm1
-      INTEGER  ::   iig  , ijg, ik
-      INTEGER  ::   istart, iend, jisf
-      INTEGER  ::   nisfg , ingb, ifind
-      INTEGER, DIMENSION(jpnij) :: nisfl
+      TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
+      !
+      INTEGER  ::   ji   , jj  , jk  , jproc          ! loop index
+      INTEGER  ::   jip1 , jim1, jjp1, jjm1           ! dummy indices
+      INTEGER  ::   iig  , ijg, ik                    ! dummy indices
+      INTEGER  ::   jisf                              ! start, end and current position in the increment array
+      INTEGER  ::   ingb, ifind                       ! 0/1 target found or need to be found 
+      INTEGER  ::   ingb, ifind, nisfl_area           ! global number of cell concerned by the wet->dry case 
+      INTEGER, DIMENSION(jpnij) :: nisfl              ! local  number of cell concerned by the wet->dry case
       !
       REAL(wp) ::   z1_sum, z1_rdtiscpl
-      REAL(wp) ::   zdtem, zdsal, zdvol, zratio
-      REAL(wp) ::   zlon , zlat
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b    !! mask before
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b      !! scale factor before
-      !!----------------------------------------------------------------------
+      REAL(wp) ::   zdtem, zdsal, zdvol, zratio       ! tem, sal, vol increment
+      REAL(wp) ::   zlon , zlat                       ! target location  
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b    ! mask before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b      ! scale factor before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b      ! scale factor before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b      ! scale factor before
+      !!----------------------------------------------------------------------
+
+      !==============================================================================
+      ! 1.0: initialisation
+      !==============================================================================
 
       ! get restart variable
-      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
-      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:), ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b(:,:,:), ldxios = lrxios   ) ! need to extrapolate T/S
+      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:)  , ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'tn'     , zt_b(:,:,:)    , ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'sn'     , zs_b(:,:,:)    , ldxios = lrxios )
 
       ! compute run length
@@ -483,21 +494,18 @@
       rdt_iscpl   = nstp_iscpl * rn_rdt
       z1_rdtiscpl = 1._wp / rdt_iscpl 
+
       IF (lwp) WRITE(numout,*) '            nb of stp for cons  = ', nstp_iscpl
       IF (lwp) WRITE(numout,*) '            coupling time step  = ', rdt_iscpl
 
-      ! mask tsn and tsb 
-      tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) * ztmask_b(:,:,:)
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) *  tmask  (:,:,:)
-      tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) * ztmask_b(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) *  tmask  (:,:,:)
-
-      !==============================================================================
-      ! diagnose the heat, salt and volume input and compute the correction variable
-      !==============================================================================
-
+      ! initialisation correction
       risfcpl_cons_vol = 0.0
       risfcpl_cons_ssh = 0.0
       risfcpl_cons_tsc = 0.0
 
+      !==============================================================================
+      ! 2.0: diagnose the heat, salt and volume input and compute the correction variable
+      !      for case where we wet a cell or cell still wet (no change in cell status)
+      !==============================================================================
+
       DO jk = 1,jpk-1
          DO jj = nldj,nlej
@@ -508,15 +516,16 @@
 
                ! heat diff
-               zdtem = tsn(ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk)   &
-                     - tsb(ji,jj,jk,jp_tem) * ze3t_b(ji,jj,jk)
+               zdtem = tsn (ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                     - zt_b(ji,jj,jk)        * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
 
                ! salt diff
-               zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk)   &
-                     - tsb(ji,jj,jk,jp_sal) * ze3t_b(ji,jj,jk)
+               zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                     - zs_b(ji,jj,jk)       * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
             
                ! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
-               risfcpl_cons_vol(ji,jj,jk)        =   zdvol * e1e2t(ji,jj) * z1_rdtiscpl
-               risfcpl_cons_tsc(ji,jj,jk,jp_sal) = - zdsal * e1e2t(ji,jj) * z1_rdtiscpl
-               risfcpl_cons_tsc(ji,jj,jk,jp_tem) = - zdtem * e1e2t(ji,jj) * z1_rdtiscpl
+               ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
+               risfcpl_cons_vol(ji,jj,jk)        = (   zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk)        ) * z1_rdtiscpl
+               risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
+               risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
 
             END DO
@@ -524,9 +533,12 @@
       END DO
       !
-      ! redistribute on valid point the vol/heat/salt removed during the coupling (ie when we dry a cell)
-      ! where we dry a cell get the number of point
+      !==============================================================================
+      ! 3.0: diagnose the heat, salt and volume input and compute the correction variable
+      !      for case where we close a cell
+      !==============================================================================
+      !
       ! compute the total number of point receiving a correction increment for each processor
       ! local
-      nisfl=0
+      nisfl(:)=0
       DO jk = 1,jpk-1
          DO jj = nldj,nlej
@@ -539,16 +551,14 @@
       !
       ! global 
-      CALL mpp_sum('isfcpl',nisfl)
-      nisfg  = SUM(nisfl           )
-      istart = SUM(nisfl(1:narea-1))
-      iend   = SUM(nisfl(1:narea  ))
+      CALL mpp_sum('isfcpl',nisfl  )
       !
       ! allocate list of point receiving correction
-      ALLOCATE(zisfpts(nisfg))
+      ALLOCATE(zisfpts(nisfl(narea)))
+      !
       zisfpts(:) = isfcons(0,0,0,-HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), 0)
       !
       ! start computing the correction and fill zisfpts
       ! local
-      jisf = istart
+      jisf = 0
       DO jk = 1,jpk-1
          DO jj = nldj,nlej
@@ -593,26 +603,49 @@
       !
       ! share data among all processes because for some point we need to find the closest wet point (could be on other process)
-      DO jisf = 1,nisfg
-         !
-         ! indices (conversion to global indices and sharing)
-         iig = zisfpts(jisf)%ii       ; ijg = zisfpts(jisf)%jj       ; ik = zisfpts(jisf)%kk
-         CALL mpp_max('isfcpl',iig)   ; CALL mpp_max('isfcpl',ijg)   ; CALL mpp_max('isfcpl',ik)
-         !
-         ! data
-         zdvol = zisfpts(jisf)%dvol   ; zdsal = zisfpts(jisf)%dsal   ; zdtem = zisfpts(jisf)%dtem
-         CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
-         !
-         ! location
-         zlat = zisfpts(jisf)%lat     ; zlon = zisfpts(jisf)%lon
-         CALL mpp_max('isfcpl',zlat)  ; CALL mpp_max('isfcpl',zlon)
-         !
-         ! find flag
-         ingb = zisfpts(jisf)%ngb
-         CALL mpp_max('isfcpl',ingb)
-         !
-         ! fill the 3d correction array
-         CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
-         !
-      END DO
+      DO jproc=1,jpnij
+         ! 
+         ! share total number of isf point treated for proc jproc
+         IF (jproc==narea) THEN
+            nisfl_area=nisfl(jproc)
+         ELSE
+            nisfl_area=0
+         END IF
+         CALL mpp_max('isfcpl',nisfl_area)
+         !
+         DO jisf = 1,nisfl_area
+            !
+            IF (jproc==narea) THEN
+               ! indices (conversion to global indices and sharing)
+               iig = zisfpts(jisf)%ii       ; ijg = zisfpts(jisf)%jj       ; ik = zisfpts(jisf)%kk
+               !
+               ! data
+               zdvol = zisfpts(jisf)%dvol   ; zdsal = zisfpts(jisf)%dsal   ; zdtem = zisfpts(jisf)%dtem
+               !
+               ! location
+               zlat = zisfpts(jisf)%lat     ; zlon = zisfpts(jisf)%lon
+               !
+               ! find flag
+               ingb = zisfpts(jisf)%ngb
+            ELSE
+               iig  =0   ; ijg  =0   ; ik   =0  
+               zdvol=-HUGE(1.0) ; zdsal=-HUGE(1.0) ; zdtem=-HUGE(1.0)
+               zlat =-HUGE(1.0) ; zlon =-HUGE(1.0)   
+               ingb = 0
+            END IF
+            !
+            ! share data (need synchronisation of data as get_correction call a global com)
+            CALL mpp_max('isfcpl',iig)   ; CALL mpp_max('isfcpl',ijg)   ; CALL mpp_max('isfcpl',ik)
+            CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
+            CALL mpp_max('isfcpl',zlat)  ; CALL mpp_max('isfcpl',zlon)
+            CALL mpp_max('isfcpl',ingb)
+            !
+            ! fill the 3d correction array
+            CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
+         END DO
+      END DO
+      !
+      !==============================================================================
+      ! 4.0: finalisation and compute ssh equivalent of the volume correction
+      !==============================================================================
       !
       ! mask (>0 out)
@@ -622,5 +655,6 @@
       !
       ! add lbclnk
-      CALL lbc_lnk_multi( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1., risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1., risfcpl_cons_vol, 'T', 1.)
+      CALL lbc_lnk_multi( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1., risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1., &
+         &                            risfcpl_cons_vol(:,:,:)       , 'T', 1.)
       !
       ! ssh correction (for dynspg_ts)
@@ -643,9 +677,12 @@
       !!----------------------------------------------------------------------
       TYPE(isfcons), DIMENSION(:), INTENT(inout) :: sisfpts
-      INTEGER,      INTENT(inout) :: kpts
-      !!----------------------------------------------------------------------
-      INTEGER,      INTENT(in   ) :: ki, kj, kk
-      INTEGER,      INTENT(in   ), OPTIONAL :: kfind
-      REAL(wp),     INTENT(in   ) :: pdvol, pdsal, pdtem, pratio
+      INTEGER,                     INTENT(inout) :: kpts
+      !!----------------------------------------------------------------------
+      INTEGER,      INTENT(in   )           :: ki, kj, kk                  !    target location (kfind=0) 
+      !                                                                    ! or source location (kfind=1)
+      INTEGER,      INTENT(in   ), OPTIONAL :: kfind                       ! 0  target cell already found
+      !                                                                    ! 1  target to be determined
+      REAL(wp),     INTENT(in   )           :: pdvol, pdsal, pdtem, pratio ! vol/sal/tem increment 
+      !                                                                    ! and ratio in case increment span over multiple cells.
       !!----------------------------------------------------------------------
       INTEGER :: ifind
@@ -667,5 +704,5 @@
    END SUBROUTINE update_isfpts
    !
-   SUBROUTINE get_correction( ki, kj, kk, zlon, zlat, pvolinc, psalinc, pteminc, kfind)
+   SUBROUTINE get_correction( ki, kj, kk, plon, plat, pvolinc, psalinc, pteminc, kfind)
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE get_correction  ***
@@ -676,7 +713,7 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in)           :: ki, kj, kk, kfind        ! target point
-      REAL(wp), INTENT(in)           :: zlon, zlat
-      REAL(wp), INTENT(in)           :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
+      INTEGER , INTENT(in) :: ki, kj, kk, kfind        ! target point indices
+      REAL(wp), INTENT(in) :: plon, plat               ! target point lon/lat
+      REAL(wp), INTENT(in) :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
       !!----------------------------------------------------------------------
       INTEGER :: jj, ji, iig, ijg
@@ -685,10 +722,10 @@
       ! define global indice of correction location
       iig = ki ; ijg = kj
-      IF ( kfind == 1 ) CALL dom_ngb( zlon, zlat, iig, ijg,'T', kk)
+      IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
       !
       ! fill the correction array
       DO jj = mj0(ijg),mj1(ijg)
          DO ji = mi0(iig),mi1(iig)
-            ! correct the vol_flx in the closest cell
+            ! correct the vol_flx and corresponding heat/salt flx in the closest cell
             risfcpl_cons_vol(ji,jj,kk)        =  risfcpl_cons_vol(ji,jj,kk       ) + pvolinc
             risfcpl_cons_tsc(ji,jj,kk,jp_sal) =  risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdiags.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdiags.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdiags.F90	(revision 11931)
@@ -46,5 +46,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac              ! thickness of the tbl and fraction of last cell affected by the tbl
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqfwf, pqoce, pqlat, pqhc ! 2d var to map in 3d
-      CHARACTER(LEN=256), INTENT(in) :: cdisf                               ! parametrisation or interactive melt
+      CHARACTER(LEN=3), INTENT(in) :: cdisf                                 ! parametrisation or interactive melt
       !!---------------------------------------------------------------------
       CHARACTER(LEN=256) :: cvarqfwf  , cvarqoce  , cvarqlat  , cvarqhc
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdynnxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdynnxt.F90	(revision 11931)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdynnxt.F90	(revision 11931)
@@ -0,0 +1,90 @@
+MODULE isfdynnxt
+   !!=========================================================================
+   !!                       ***  MODULE  isfnxt  ***
+   !! Ice shelf update: compute the dynnxt ice shelf contribution
+   !!=========================================================================
+   !! History :  OPA  !  2019-09  (P. Mathiot)  Original code
+   !!-------------------------------------------------------------------------
+  
+   !!-------------------------------------------------------------------------
+   !!   isfnxt       : aplly correction needed for the ice shelf to ensure conservation
+   !!-------------------------------------------------------------------------
+
+   USE isf
+
+   USE phycst , ONLY: r1_rau0                ! physical constant
+   USE dom_oce, ONLY: e3t_b, e3t_n, r1_e1e2t ! time and space domain
+
+   USE in_out_manager
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isf_dynnxt
+
+CONTAINS
+
+   SUBROUTINE isf_dynnxt ( kt, pcoef )
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_dynnxt  ***
+      !!
+      !! ** Purpose : compute the ice shelf volume filter correction for cavity, param, ice sheet coupling case
+      !!
+      !!-------------------------- OUT -------------------------------------
+      INTEGER ,                     INTENT(in   ) :: kt
+      !
+      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
+      !!--------------------------------------------------------------------
+      INTEGER :: jk  ! loop index
+      !!--------------------------------------------------------------------
+      !
+      ! ice shelf cavity
+      IF ( ln_isfcav_mlt ) CALL isf_dynnxt_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, pcoef)
+      !
+      ! ice shelf parametrised
+      IF ( ln_isfpar_mlt ) CALL isf_dynnxt_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, pcoef)
+      !
+      IF ( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
+         DO jk = 1, jpkm1
+            e3t_b(:,:,jk) =   e3t_b(:,:,jk) - pcoef * risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
+         END DO
+      END IF
+      !
+   END SUBROUTINE isf_dynnxt
+
+   SUBROUTINE isf_dynnxt_mlt ( ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, pcoef )
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_dynnxt_mlt  ***
+      !!
+      !! ** Purpose : compute the ice shelf volume filter correction for cavity or param
+      !!
+      !!-------------------------- IN  -------------------------------------
+      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot     ! top and bottom level of tbl
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl    ! fraction of bottom cell included in tbl, tbl thickness
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfwf , pfwf_b   ! now/before fwf
+      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
+      !!----------------------------------------------------------------------
+      INTEGER :: ji,jj,jk
+      REAL(wp), DIMENSION(jpi,jpj) :: zfwfinc
+      !!----------------------------------------------------------------------
+      !
+      ! compute fwf conservation correction
+      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:) * r1_rau0
+      !
+      ! add the increment in the tbl
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( ktop(ji,jj) <= jk .AND. jk < kbot(ji,jj)  ) THEN
+                  e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zfwfinc(ji,jj) * e3t_n(ji,jj,jk)
+               ELSEIF ( jk == kbot(ji,jj) ) THEN
+                  e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zfwfinc(ji,jj) * e3t_n(ji,jj,jk) * pfrac(ji,jj)
+               ENDIF
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE isf_dynnxt_mlt
+
+END MODULE isfdynnxt
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 11931)
@@ -46,10 +46,4 @@
             ! conservation option
             IF ( ln_isfcpl_cons ) CALL isf_hdiv_cpl(risfcpl_cons_vol, phdiv)
-            !
-            IF ( ln_isfdebug ) THEN
-               CALL debug('isfdiv: phdiv'      ,phdiv(:,:,:))
-               CALL debug('isfdiv: risfcpl_vol',risfcpl_vol(:,:,:))
-               CALL debug('isfdiv: fwfisf     ',fwfisf_cav(:,:)+fwfisf_cav_b(:,:))
-            END IF
             !
          END IF
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfnxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfnxt.F90	(revision 11930)
+++ 	(revision )
@@ -1,94 +1,0 @@
-MODULE isfnxt
-   !!=========================================================================
-   !!                       ***  MODULE  isfnxt  ***
-   !! Ice shelf update: compute the dynnxt ice shelf contribution
-   !!=========================================================================
-   !! History :  OPA  !  2019-09  (P. Mathiot)  Original code
-   !!-------------------------------------------------------------------------
-  
-   !!-------------------------------------------------------------------------
-   !!   isfnxt       : aplly correction needed for the ice shelf to ensure conservation
-   !!-------------------------------------------------------------------------
-
-   USE isf
-
-   USE phycst , ONLY: r1_rau0                ! physical constant
-   USE dom_oce, ONLY: e3t_b, e3t_n, r1_e1e2t ! time and space domain
-
-   USE in_out_manager
-
-   IMPLICIT NONE
-
-   PRIVATE
-
-   PUBLIC isf_dynnxt
-
-CONTAINS
-
-   SUBROUTINE isf_dynnxt ( kt, pcoef )
-      !!--------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_dynnxt  ***
-      !!
-      !! ** Purpose : compute the ice shelf volume filter correction for cavity, param, ice sheet coupling case
-      !!
-      !!--------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      INTEGER ,                     INTENT(in   ) :: kt
-      !
-      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
-      !!-------------------------- IN  -------------------------------------
-      !!--------------------------------------------------------------------
-      INTEGER :: jk  ! loop index
-      !!--------------------------------------------------------------------
-      !
-      ! ice shelf cavity
-      IF ( ln_isfcav_mlt ) CALL isf_dynnxt_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, pcoef)
-      !
-      ! ice shelf parametrised
-      IF ( ln_isfpar_mlt ) CALL isf_dynnxt_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, pcoef)
-      !
-      IF ( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
-         DO jk = 1, jpkm1
-            e3t_b(:,:,jk) =   e3t_b(:,:,jk) - pcoef * risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
-         END DO
-      END IF
-      !
-   END SUBROUTINE isf_dynnxt
-
-   SUBROUTINE isf_dynnxt_mlt ( ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, pcoef )
-      !!--------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_dynnxt_mlt  ***
-      !!
-      !! ** Purpose : compute the ice shelf volume filter correction for cavity or param
-      !!
-      !!--------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      !!-------------------------- IN  -------------------------------------
-      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot     ! top and bottom level of tbl
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl    ! fraction of bottom cell included in tbl, tbl thickness
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfwf , pfwf_b   ! now/before fwf
-      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
-      !!----------------------------------------------------------------------
-      INTEGER :: ji,jj,jk
-      REAL(wp), DIMENSION(jpi,jpj) :: zfwfinc
-      !!----------------------------------------------------------------------
-      !
-      ! compute fwf conservation correction
-      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:) * r1_rau0
-      !
-      ! add the increment in the tbl
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( ktop(ji,jj) <= jk .AND. jk < kbot(ji,jj)  ) THEN
-                  e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zfwfinc(ji,jj) * e3t_n(ji,jj,jk)
-               ELSEIF ( jk == kbot(ji,jj) ) THEN
-                  e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zfwfinc(ji,jj) * e3t_n(ji,jj,jk) * pfrac(ji,jj)
-               ENDIF
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE isf_dynnxt_mlt
-
-END MODULE isfnxt
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 11931)
@@ -14,5 +14,4 @@
    !!   isfpar       : compute ice shelf melt using a prametrisation of ice shelf cavities
    !!----------------------------------------------------------------------
-   !USE oce            ! ocean dynamics and tracers
    USE isf            ! ice shelf
    !
@@ -23,5 +22,6 @@
    USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine
    !
-   USE dom_oce  , ONLY: jpi,jpj, bathy ! ocean space and time domain
+   USE dom_oce  , ONLY: bathy          ! ocean space and time domain
+   USE par_oce  , ONLY: jpi,jpj        ! ocean space and time domain
    USE phycst   , ONLY: r1_rau0_rcp    ! physical constants
    !
@@ -62,13 +62,13 @@
       CALL isfpar_mlt( kt, zqhc, zqoce, pqfwf  )
       !
-      ! compute heat and water flux  (change signe directly in the melt subroutine)
+      ! compute heat and water flux ( > 0 out )
       pqfwf(:,:) = pqfwf(:,:) * mskisf_par(:,:)
       zqoce(:,:) = zqoce(:,:) * mskisf_par(:,:)
       zqhc (:,:) = zqhc(:,:)  * mskisf_par(:,:)
       !
-      ! compute heat content flux
+      ! compute heat content flux ( > 0 out )
       zqlat(:,:) = pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
       !
-      ! total heat flux
+      ! total heat flux ( > 0 out )
       zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
       !
@@ -107,5 +107,5 @@
       CALL isf_alloc_par()
       !
-      ! par
+      ! initialisation
       misfkt_par(:,:)     = 1         ; misfkb_par(:,:)       = 1         
       rhisf_tbl_par(:,:)  = 1e-20     ; rfrac_tbl_par(:,:)    = 0.0_wp
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90	(revision 11931)
@@ -18,5 +18,5 @@
    USE in_out_manager              ! I/O manager
    USE iom        , ONLY: iom_put  ! I/O library
-   USE fldread    , ONLY: fld_read !
+   USE fldread    , ONLY: fld_read, FLD, FLD_N !
    USE lib_fortran, ONLY: glob_sum !
    USE lib_mpp    , ONLY: ctl_stop !
@@ -49,5 +49,5 @@
       !!                        1 : Specified melt flux
       !!                        2 : Beckmann & Goose parameterization
-      !!----------------------------------------------------------------------
+      !!
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqfwf, pqoce, pqhc  ! fresh water, ice-ocean heat and heat content fluxes
@@ -81,5 +81,4 @@
       !!              data read into a forcing files.
       !!
-      !!----------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf, pqoce  ! fresh water and ice-ocean heat fluxes
@@ -105,4 +104,5 @@
       pqoce(:,:) =   pqfwf(:,:) * rLfusisf             ! ocean/ice shelf flux assume to be equal to latent heat flux
       pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
+      !
       CALL iom_put('isftfrz_par', ztfrz )
       !
@@ -155,4 +155,6 @@
       ! output thermal driving
       CALL iom_put('isfthermald_par',( ztfrz(:,:) - ztavg(:,:) ) * mskisf_par(:,:))
+      !
+      ! output freezing point used to define the thermal driving and heat content fluxes
       CALL iom_put('isftfrz_par', ztfrz )
       !
@@ -206,5 +208,5 @@
       zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
       ! 
-      ! i3. -----------Define fwf and qoce
+      ! 3. -----------Define fwf and qoce
       ! ocean heat flux is assume to be equal to the latent heat
       pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out )
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90	(revision 11931)
@@ -11,5 +11,5 @@
    !!----------------------------------------------------------------------
    !
-   USE dom_oce, ONLY: jpi,jpj,jpk,jpts ! time and space domain
+   USE par_oce, ONLY: jpi,jpj,jpk,jpts ! time and space domain
    !
    USE in_out_manager ! I/O manager
@@ -31,11 +31,12 @@
    SUBROUTINE isfrst_read(cdisf, ptsc, pfwf, ptsc_b, pfwf_b )
       !!---------------------------------------------------------------------
+      !!
       !!   isfrst_read : read iceshelf variables from restart
-      !!---------------------------------------------------------------------
+      !!
       !!-------------------------- OUT --------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) :: pfwf_b
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(  out) :: ptsc_b
       !!-------------------------- IN  --------------------------------------
-      CHARACTER(LEN=256)               , INTENT(in   ) :: cdisf
+      CHARACTER(LEN=3)                 , INTENT(in   ) :: cdisf
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: pfwf
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc
@@ -66,16 +67,14 @@
       ENDIF
 
-      CALL FLUSH(numout)
-
    END SUBROUTINE isfrst_read
    ! 
    SUBROUTINE isfrst_write(kt, cdisf, ptsc, pfwf )
       !!---------------------------------------------------------------------
+      !!
       !!   isfrst_write : write iceshelf variables in restart
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT --------------------------------------
+      !!
       !!-------------------------- IN  --------------------------------------
       INTEGER                          , INTENT(in   ) :: kt
-      CHARACTER(LEN=256)               , INTENT(in   ) :: cdisf
+      CHARACTER(LEN=3)                 , INTENT(in   ) :: cdisf
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: pfwf
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 11931)
@@ -26,4 +26,5 @@
    !
    USE lib_mpp, ONLY: ctl_stop, ctl_nam
+   USE fldread, ONLY: FLD, FLD_N
    USE in_out_manager ! I/O manager
    USE timing
@@ -62,7 +63,11 @@
       IF( ln_timing )   CALL timing_start('isf')
       !
+      !=======================================================================
+      ! 1.: compute melt and associated heat fluxes in the ice shelf cavities
+      !=======================================================================
+      !
       IF ( ln_isfcav_mlt ) THEN
          !
-         ! before time step 
+         ! 1.1: before time step 
          IF ( kt /= nit000 ) THEN 
             risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
@@ -70,16 +75,20 @@
          END IF
          !
-         ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
          rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
          CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
          !
-         ! compute ice shelf melt
+         ! 1.3: compute ice shelf melt
          CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav)
          !
       END IF
       ! 
+      !=================================================================================
+      ! 2.: compute melt and associated heat fluxes for not resolved ice shelf cavities
+      !=================================================================================
+      !
       IF ( ln_isfpar_mlt ) THEN
          !
-         ! before time step 
+         ! 2.1: before time step 
          IF ( kt /= nit000 ) THEN 
             risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
@@ -87,13 +96,17 @@
          END IF
          !
-         ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
-         ! by simplicity, we assume the top level where param applied do not change with time
+         ! 2.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         ! by simplicity, we assume the top level where param applied do not change with time (done in init part)
          rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
          CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
          !
-         ! compute ice shelf melt
+         ! 2.3: compute ice shelf melt
          CALL isf_par( kt, risf_par_tsc, fwfisf_par)
          !
       END IF
+      !
+      !==================================================================================
+      ! 3.: output specific restart variable in case of coupling with an ice sheet model
+      !==================================================================================
       !
       IF ( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write(kt)
@@ -230,5 +243,5 @@
       IF ( l_isfoasis .AND. ln_isf ) THEN
          !
-         CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
+         CALL ctl_stop( ' OASIS and ice shelf not tested' )
          !
          ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation 
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 11931)
@@ -37,5 +37,4 @@
       !!                https://doi.org/10.1029/2007JC004368 , 2008
       !!
-      !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) :: pvarout ! 2d average of pvarin
@@ -115,5 +114,4 @@
       !!              over a thickness phtbl. The bottom level is partially counted (pfrac).
       !!
-      !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pvarout      ! tbl property averaged over phtbl between level ktop and kbot
@@ -154,5 +152,4 @@
       !!              - fraction of the bottom level affected by the tbl
       !!
-      !!---------------------------------------------------------------------
       !!-------------------------- OUT --------------------------------------
       INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
@@ -209,5 +206,4 @@
       !! ** Purpose : compute bottom level of the isf top boundary layer
       !!
-      !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
@@ -246,5 +242,4 @@
       !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
       !!
-      !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       INTEGER,  DIMENSION(jpi,jpj), INTENT(  out) :: ktop        ! top level affected by the ice shelf parametrisation
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 11931)
@@ -14,5 +14,5 @@
    USE iom           , ONLY: iom_open, iom_get, iom_close, jpdom_data ! read input file
    USE lib_fortran   , ONLY: glob_sum, glob_min, glob_max             ! compute global value
-   USE dom_oce       , ONLY: jpi,jpj,jpk                              ! domain size
+   USE par_oce       , ONLY: jpi,jpj,jpk                              ! domain size
    USE in_out_manager, ONLY: wp, lwp, numout                          ! miscelenious
 
@@ -35,5 +35,4 @@
       !! ** Purpose : read input file
       !!
-      !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pvar          ! output variable
@@ -57,6 +56,4 @@
       !! ** Purpose : add debug print
       !!
-      !!--------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
       !!-------------------------- IN  -------------------------------------
       CHARACTER(LEN=256)          , INTENT(in   ) :: cdtxt
@@ -82,6 +79,4 @@
       !! ** Purpose : add debug print
       !!
-      !!--------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
       !!-------------------------- IN  -------------------------------------
       CHARACTER(LEN=256)              , INTENT(in   ) :: cdtxt
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traisf.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traisf.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traisf.F90	(revision 11931)
@@ -68,15 +68,10 @@
          ! ensure 0 trend due to unconservation of the ice shelf coupling
          IF ( ln_isfcpl_cons ) CALL tra_isf_cpl(risfcpl_cons_tsc, tsa)
-
-         IF ( ln_isfdebug ) THEN
-            CALL debug('tra_isf: risfcpl_tsc T',risfcpl_tsc(:,:,1))
-            CALL debug('tra_isf: risfcpl_tsc S',risfcpl_tsc(:,:,2))
-         END IF
          !
       END IF
       !
       IF ( ln_isfdebug ) THEN
-         CALL debug('tra_isf: tsa T'        ,tsa(:,:,:,1))
-         CALL debug('tra_isf: tsa S'        ,tsa(:,:,:,2))
+         CALL debug('tra_isf: tsa T', tsa(:,:,:,1))
+         CALL debug('tra_isf: tsa S', tsa(:,:,:,2))
       END IF
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 11931)
@@ -341,4 +341,5 @@
                   IF( ll_isf ) THEN
                      !
+                     ! melt in the cavity
                      IF ( ln_isfcav_mlt ) THEN
                         ! level fully include in the Losch_2008 ice shelf boundary layer
@@ -357,4 +358,6 @@
                         END IF
                      END IF
+                     !
+                     ! parametrised melt (cavity closed)
                      IF ( ln_isfpar_mlt ) THEN
                         ! level fully include in the Losch_2008 ice shelf boundary layer
@@ -374,7 +377,13 @@
                      END IF
                      !
-                     IF (ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN    ! risfcpl_vol_n = 0 and risfcpl_vol_b = risfcpl_vol
-                        ztc_f  = ztc_f  + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj)
-                        ze3t_f = ze3t_f - zfact1 * risfcpl_vol(ji,jj,jk   ) * r1_e1e2t(ji,jj)
+                     ! ice sheet coupling correction
+                     IF ( ln_isfcpl ) THEN
+                        !
+                        ! at kt = nit000,  risfcpl_vol_n = 0 and risfcpl_vol_b = risfcpl_vol so contribution nul
+                        IF ( ln_rstart .AND. kt == nit000+1 ) THEN
+                           ztc_f  = ztc_f  + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj)
+                           ze3t_f = ze3t_f - zfact1 * risfcpl_vol(ji,jj,jk   ) * r1_e1e2t(ji,jj)
+                        END IF
+                        !
                      END IF
                      !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 11931)
@@ -241,5 +241,5 @@
                          CALL tra_sbc       ( kstp )  ! surface boundary condition
       IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
-      IF( ln_isf     ) CALL tra_isf ( kstp )  ! ice shelf heat flux
+      IF( ln_isf     )   CALL tra_isf       ( kstp )  ! ice shelf heat flux
       IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
       IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isf.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isf.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isf.F90	(revision 11931)
@@ -22,5 +22,5 @@
    PRIVATE
 
-   PUBLIC   isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl
+   PUBLIC   isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl, isf_dealloc_cpl
    !
    !-------------------------------------------------------
@@ -175,7 +175,5 @@
       !!                  ***  ROUTINE isf_alloc_cpl  ***
       !!
-      !! ** Purpose : 
-      !!
-      !! ** Method  : 
+      !! ** Purpose : allocate array use for the ice sheet coupling
       !!
       !!----------------------------------------------------------------------
@@ -202,11 +200,28 @@
    END SUBROUTINE isf_alloc_cpl
 
+   SUBROUTINE isf_dealloc_cpl()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_dealloc_cpl  ***
+      !!
+      !! ** Purpose : de-allocate useless public 3d array used for ice sheet coupling
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      ierr = 0
+      !
+      DEALLOCATE( risfcpl_ssh, risfcpl_tsc, risfcpl_vol, STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop('STOP','isfcpl: failed to deallocate arrays.')
+      !
+   END SUBROUTINE isf_dealloc_cpl
+
    SUBROUTINE isf_alloc()
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isf_alloc  ***
       !!
-      !! ** Purpose : 
-      !!
-      !! ** Method  : 
+      !! ** Purpose : allocate array used for the ice shelf cavity (cav and par)
       !!
       !!----------------------------------------------------------------------
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isfstp.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 11931)
@@ -97,5 +97,10 @@
       END IF
       !
-      IF ( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write(kt)
+      IF ( ln_isfcpl ) THEN
+         ! after step nit000 + 2 we do not need anymore the risfcpl_ arrays
+         IF ( kt == nit000 + 2 ) CALL isf_dealloc_cpl()
+
+         IF ( lrst_oce ) CALL isfcpl_rst_write(kt)
+      END IF
       !
       IF( ln_timing )   CALL timing_stop('isf')
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 11930)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 11931)
@@ -118,11 +118,22 @@
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
             z_fwf = glob_sum( 'sbcfwb',  e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) - snwice_fmass(:,:) ) )
+            !
+            ! correction for ice sheet coupling testing (ie remove the excess through the surface)
+            ! test impact on the melt as conservation correction made in depth
+            ! test conservation level as sbcfwb is conserving
+            ! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
+            IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
+               z_fwf = z_fwf + glob_sum( 'sbcfwb',  e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rau0 )
+            END IF
+            !
             z_fwf = z_fwf / area
             zcoef = z_fwf * rcp
-            emp(:,:) = emp(:,:) - z_fwf            * tmask(:,:,1) ! (Eq. 34 AD2015)
-            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! could be sst_m if we don't want any bouyancy fluxes
-            sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! could be sss_m if we don't want any bouyancy fluxes
+            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1) ! (Eq. 34 AD2015)
+            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
+            sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
             !qns(:,:) = qns(:,:) + zcoef * ( -1.9 ) * tmask(:,:,1) ! (Eq. 35 AD2015) ! could be sst_m if we don't want any bouyancy fluxes
             !sfx(:,:) = sfx(:,:) + z_fwf * ( 33.8 ) * tmask(:,:,1) ! (Eq. 36 AD2015) ! could be sss_m if we don't want any bouyancy fluxes
+            !qns(:,:) = qns(:,:) + zcoef * ( -1.0 ) * tmask(:,:,1) ! use for ISOMIP+ coupling sanity check (keep ssh cst while playing with cpl conservation option)
+            !sfx(:,:) = sfx(:,:) + z_fwf * ( 34.2 ) * tmask(:,:,1) ! use for ISOMIP+ coupling sanity check (keep ssh cst while playing with cpl conservation option)
          ENDIF
          !
