Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90	(revision 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90	(revision 11553)
@@ -26,4 +26,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean variables
+   USE isfstp
    USE dom_oce        ! domain: ocean
    USE sbc_oce        ! surface boundary condition: ocean
@@ -448,4 +449,8 @@
 #endif
       !
+      ! need to know ln_isfcpl for istate (restart extrapolation)
+      ! and for domvvl_init (ssh extrapolation when recompute e3 variable
+      CALL isf_nam
+      !
    END SUBROUTINE dom_nam
 
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domvvl.F90	(revision 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domvvl.F90	(revision 11553)
@@ -822,5 +822,5 @@
                CALL isfcpl_ssh()
                !
-               ! everything need to compute over new water column
+               ! e3t need to be computed where ice shelf draft changed (new ht_0)
                id1 = 0 ; id2 = 0 ; id3 = 0 ; id4 = 0 ; id5 = 0
                !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/IOM/restart.F90	(revision 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/IOM/restart.F90	(revision 11553)
@@ -19,4 +19,6 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers 
+   USE isf             ! ice shelf
+   USE isfcpl          ! ice shelf coupling
    USE dom_oce         ! ocean space and time domain
    USE sbc_ice         ! only lk_si3 
@@ -282,4 +284,7 @@
       ENDIF
       !
+      ! ice sheet coupling: extrapolation of restart to fill new wet cell and compute divergence correction
+      IF ( ln_isfcpl ) CALL isfcpl_init()
+      !
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
          tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
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 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11553)
@@ -182,9 +182,16 @@
       ierr = 0
       !
-      ALLOCATE( risfcpl_ssh(jpi,jpj), risfcpl_tsc(jpi,jpj,jpk,jpts)      , risfcpl_vol(jpi,jpj,jpk)      ,             &
-         &      risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) ,             &
-         &      risfcpl_cons_tsc_b(jpi,jpj,jpk,jpts) , risfcpl_cons_vol_b(jpi,jpj,jpk),          &
-                risfcpl_cons_ssh_b(jpi,jpj), risfcpl_cons_ssh(jpi,jpj), STAT=ialloc )
-      ierr = ierr + ialloc
+      ALLOCATE( risfcpl_ssh(jpi,jpj), risfcpl_tsc(jpi,jpj,jpk,jpts), risfcpl_vol(jpi,jpj,jpk), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      risfcpl_tsc(:,:,:,:) = 0.0 ; risfcpl_vol(:,:,:) = 0.0 ; risfcpl_ssh(:,:) = 0.0
+
+      IF ( ln_isfcpl_cons) THEN
+         ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) ,risfcpl_cons_ssh(jpi,jpj), STAT=ialloc )
+         ierr = ierr + ialloc
+         !
+         risfcpl_cons_tsc(:,:,:,:) = 0.0 ; risfcpl_cons_vol(:,:,:) = 0.0 ; risfcpl_cons_ssh(:,:) = 0.0
+         !
+      END IF
       !
       CALL mpp_sum ( 'isf', ierr )
@@ -226,4 +233,13 @@
       CALL mpp_sum ( 'isf', ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
+      !
+      ! initalisation of fwf and tsc array to 0
+      risfload(:,:)       = 0.0_wp
+      fwfisf_oasis(:,:)   = 0.0_wp
+      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
+      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
+      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
+      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
+      !
 
    END SUBROUTINE isf_alloc
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 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11553)
@@ -192,4 +192,7 @@
       END SELECT
       !
+      ! compute mask
+      mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
+      !
    END SUBROUTINE isf_cav_init
 
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 11553)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90	(revision 11553)
@@ -0,0 +1,660 @@
+MODULE isfcpl
+   !!======================================================================
+   !!                       ***  MODULE  isfcpl  ***
+   !!
+   !! iceshelf coupling module : module managing the coupling between NEMO and an ice sheet model
+   !!
+   !!======================================================================
+   !! History :  4.1  !  2019-07  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfrst : read/write iceshelf variables in/from restart
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE isf            ! ice shelf variable
+   USE isfutils       ! debuging
+   USE lib_mpp        ! mpp routine
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   !
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isfcpl_rst_write, isfcpl_init  ! iceshelf restart read and write 
+   PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons
+
+   TYPE isfcons
+      INTEGER :: ii     ! i global
+      INTEGER :: jj     ! j global
+      INTEGER :: kk     ! k level
+      REAL(wp):: dvol   ! volume increment
+      REAL(wp):: dsal   ! salt increment
+      REAL(wp):: dtem   ! heat increment
+      REAL(wp):: lon    ! lon
+      REAL(wp):: lat    ! lat
+      INTEGER :: ngb    ! 0/1 (valid location or not (ie on halo or no neigbourg))
+   END TYPE
+   !
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE isfcpl_init()
+      ! 
+      CALL isf_alloc_cpl()
+      !
+      ! extrapolation ssh
+      CALL isfcpl_ssh()
+      !
+      ! extrapolation tracer properties
+      CALL isfcpl_tra()
+      !
+      ! correction of the horizontal divergence and associated temp. and salt content flux
+      CALL isfcpl_vol()
+      !
+      ! apply the 'conservation' method
+      IF ( ln_isfcpl_cons ) CALL isfcpl_cons()
+      !
+      ! mask velocity properly (mask used in restart not compatible with new mask)
+      un(:,:,:) = un(:,:,:,) * umask(:,:,:)
+      vn(:,:,:) = vn(:,:,:,) * vmask(:,:,:)
+      !
+      ! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
+      !           - decide how to manage thickness level change in conservation
+      !
+   END SUBROUTINE isfcpl_init
+   ! 
+   SUBROUTINE isfcpl_rst_write(kt)
+      !!---------------------------------------------------------------------
+      !!   isfrst_cpl_write : write icesheet coupling variables in restart
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT --------------------------------------
+      !!-------------------------- IN  --------------------------------------
+      INTEGER, INTENT(in) :: kt
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      IF( lwxios ) CALL iom_swap( cwxios_context )
+      CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n , ldxios = lwxios )
+      IF( lwxios ) CALL iom_swap( cxios_context )
+      !
+   END SUBROUTINE isfcpl_rst_write
+
+   SUBROUTINE isfcpl_ssh ()
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_ssh  ***
+      !! 
+      !! ** Purpose :   basic guess of ssh in new wet cell during coupling step
+      !! 
+      !! ** Method  :   basic extrapolation from neigbourg cells
+      !!
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT --------------------------------------
+      !!-------------------------- IN  --------------------------------------
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER :: ji, jj, jd          !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1
+      !!
+      REAL(wp):: zsummsk
+      REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
+      !!----------------------------------------------------------------------
+      !
+      CALL iom_get( numror, jpdom_autoglo, 'ssmask'  , zssmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+
+      ! compute new ssh if we open a full water column 
+      ! rude average of the closest neigbourgs (e1e2t not taking into account)
+      !
+      zssh(:,:)     = sshn(:,:)
+      zssmask0(:,:) = zssmask_b(:,:)
+      !
+      DO jd = 1, nn_drown
+         !
+         zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
+         DO jj = 2,jpj-1
+            DO ji = 2,jpi-1
+               jip1=ji+1; jim1=ji-1;
+               jjp1=jj+1; jjm1=jj-1;
+               !
+               zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
+               !
+               IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
+                  sshn(ji,jj)=( zssh(jip1,jj)*zssmask0(jip1,jj)     &
+                  &           + zssh(jim1,jj)*zssmask0(jim1,jj)     &
+                  &           + zssh(ji,jjp1)*zssmask0(ji,jjp1)     &
+                  &           + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
+                  zssmask_b(ji,jj) = 1._wp
+               ENDIF
+            END DO
+         END DO
+         !
+         zssh(:,:) = sshn(:,:)
+         zssmask0(:,:) = zssmask_b(:,:)
+         !
+         CALL lbc_lnk_multi( 'iscplrst', zssh, 'T', 1., zssmask0, 'T', 1. )
+         !
+      END DO
+      !
+      ! update sshn
+      sshn(:,:) = zssh(:,:) * ssmask(:,:)
+      !
+      ! force to start on an euler time step
+      neuler = 0
+      !
+      sshb(:,:) = sshn(:,:)
+      !
+   END SUBROUTINE isfcpl_ssh
+
+   SUBROUTINE isfcpl_tra ()
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_tra  ***
+      !! 
+      !! ** Purpose :   compute new tn, sn in case of evolving geometry of ice shelves 
+      !! 
+      !! ** Method  :   tn, sn : basic extrapolation from neigbourg cells
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
+      !REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
+      !!
+      INTEGER :: ji, jj, jk, jd          !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
+      !!
+      REAL(wp):: zsummsk
+      REAL(wp):: zdz, zdzm1, zdzp1
+      !!
+      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
+      !!----------------------------------------------------------------------
+      ! 
+      IF(lwp) WRITE(numout,*) ' isfrst_cpl_read: read ice sheet coupling restart variable ' 
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) ''
+      !
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+      !CALL iom_get( numror, jpdom_autoglo, 'wmask'  , zwmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+      !CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
+      !
+      ! 
+      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
+      !PM: Is this IF needed since change to VVL by default
+      !bugged : to be corrected (PM)
+      ! back up original t/s/mask
+      !tsb (:,:,:,:) = tsn(:,:,:,:)
+      !
+     ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
+
+!      IF (.NOT.ln_linssh) THEN
+!         DO jk = 2,jpk-1
+!            DO jj = 1,jpj
+!               DO ji = 1,jpi
+!                  IF (wmask(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ztmask_b(ji,jj,1)==0._wp) ) THEN
+!
+!                     !compute weight
+!                     zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw_n(ji,jj,jk+1))
+!                     zdzm1 = MAX(0._wp,gdepw_n(ji,jj,jk  ) - pdepw_b(ji,jj,jk  ))
+!                     zdz   = e3t_n(ji,jj,jk) - zdzp1 - zdzm1 ! if isf : e3t = gdepw_n(ji,jj,jk+1)- gdepw_n(ji,jj,jk)
+!
+!                     IF (zdz .LT. 0._wp) THEN
+!                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
+!                     END IF
+!
+!                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
+!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem) &
+!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/e3t_n(ji,jj,jk)
+!
+!                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
+!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal) &
+!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/e3t_n(ji,jj,jk)
+!
+!                  END IF
+!               END DO
+!            END DO
+!         END DO
+!      END IF
+
+      zts0(:,:,:,:)  = tsn(:,:,:,:)
+      ztmask0(:,:,:) = ztmask_b(:,:,:)
+      ztmask1(:,:,:) = ztmask_b(:,:,:)
+      !
+      ! iterate the extrapolation processes nn_drown times
+      DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
+         DO jk = 1,jpk-1
+            !
+            ! define new wet cell
+            zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
+            !
+            DO jj = 2,jpj-1
+               DO ji = 2,jpi-1
+                  jip1=ji+1; jim1=ji-1;
+                  jjp1=jj+1; jjm1=jj-1;
+                  !
+                  ! check if a wet neigbourg cell is present
+                  zsummsk = ztmask0(jip1,jj  ,jk) + ztmask0(jim1,jj  ,jk) &
+                          + ztmask0(ji  ,jjp1,jk) + ztmask0(ji  ,jjm1,jk)
+                  !
+                  ! if neigbourg wet cell available at the same level
+                  IF ( zdmask(ji,jj) == 1._wp  .AND. zsummsk /= 0._wp ) THEN
+                     !
+                     ! horizontal basic extrapolation
+                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1) * ztmask0(jip1,jj  ,jk) &
+                     &               + zts0(jim1,jj  ,jk,1) * ztmask0(jim1,jj  ,jk) &
+                     &               + zts0(ji  ,jjp1,jk,1) * ztmask0(ji  ,jjp1,jk) &
+                     &               + zts0(ji  ,jjm1,jk,1) * ztmask0(ji  ,jjm1,jk) ) / zsummsk
+                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2) * ztmask0(jip1,jj  ,jk) &
+                     &               + zts0(jim1,jj  ,jk,2) * ztmask0(jim1,jj  ,jk) &
+                     &               + zts0(ji  ,jjp1,jk,2) * ztmask0(ji  ,jjp1,jk) &
+                     &               + zts0(ji  ,jjm1,jk,2) * ztmask0(ji  ,jjm1,jk) ) / zsummsk
+                     !
+                     ! update mask for next pass
+                     ztmask1(ji,jj,jk)=1
+                     !
+                  ! in case no neigbourg wet cell available at the same level
+                  ! check if a wet cell is available below
+                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
+                     !
+                     ! vertical extrapolation if horizontal extrapolation failed
+                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
+                     !
+                     ! check if a wet neigbourg cell is present
+                     zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
+                     IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
+                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
+                        &               + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
+                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
+                        &               + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
+                        !
+                        ! update mask for next pass
+                        ztmask1(ji,jj,jk)=1._wp
+                     END IF
+                  END IF
+               END DO
+            END DO
+         END DO
+         !
+         ! update temperature and salinity and mask
+         zts0(:,:,:,:)  = tsn(:,:,:,:)
+         ztmask0(:,:,:) = ztmask1(:,:,:)
+         !
+         CALL lbc_lnk_multi( 'iscplrst', zts0(:,:,:,jp_tem), 'T', 1., zts0(:,:,:,jp_sal), 'T', 1., ztmask0, 'T', 1.)
+         !
+      END DO  ! nn_drown
+      !
+      ! mask new tsn field
+      tsn(:,:,:,jp_tem) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
+      tsn(:,:,:,jp_sal) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
+      !
+      ! sanity check
+      ! -----------------------------------------------------------------------------------------
+      ! case we open a cell but no neigbour cells available to get an estimate of T and S
+      DO jk = 1,jpk-1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               IF (tmask(ji,jj,jk) == 1._wp .AND. tsn(ji,jj,jk,2) == 0._wp)              &
+                  &   CALL ctl_stop('STOP', 'failing to fill all new weet cell,     &
+                  &                          try increase nn_drown or activate XXXX &
+                  &                         in your domain cfg computation'         )
+            END DO
+         END DO
+      END DO
+      ! 
+   END SUBROUTINE isfcpl_tra
+
+   SUBROUTINE isfcpl_vol
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_vol  ***
+      !! 
+      !! ** Purpose :   
+      !!                
+      !! 
+      !! ** Method  :   
+      !!                
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER :: ji, jj, jk 
+      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
+      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)
+      DO jk = 1, jpk                                 ! Horizontal slab
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                  &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
+                  &                * ztmask_b(ji,jj,jk)
+            END DO
+         ENDDO
+         !
+         ! 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)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * e3u_n(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                  &                 + e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * e3v_n(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
+                  &               * tmask(ji,jj,jk)
+            END DO
+         ENDDO
+         !
+         ! 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
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            !
+            ikt = mikt(ji,jj)
+            IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN
+               risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1))  ! test sign
+            ENDIF
+            !
+         END DO
+      ENDDO
+      !
+      ! mask volume flux divergence correction
+      risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
+      !
+      CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
+      CALL lbc_lnk( 'iscpl', zqvolb, 'T', 1. )
+      CALL lbc_lnk( 'iscpl', zqvoln, 'T', 1. )
+      DO jk = 1,jpk
+      CALL debug('zqvolb',zqvolb(:,:,jk))
+      END DO
+      DO jk = 1,jpk
+      CALL debug('zqvoln',zqvoln(:,:,jk))
+      END DO
+      !
+      ! get 3d tra increment to apply at the first time step
+      ! temperature and salt content flux computed using local tsn 
+      ! (very simple advection scheme)
+      ! (>0 out)
+      risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_tem)
+      risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_sal)
+      !
+      ! ssh correction (for dynspg_ts)
+      risfcpl_ssh(:,:) = 0.0
+      DO jk = 1,jpk
+         risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
+      END DO
+
+   END SUBROUTINE isfcpl_vol
+
+   SUBROUTINE isfcpl_cons
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_cons  ***
+      !! 
+      !! ** Purpose :   compute the corrective increment in volume/salt/heat to put back the vol/heat/salt
+      !!                removed or added during the coupling processes (wet or dry new cell)
+      !! 
+      !! ** Method  :   - compare volume/heat/salt before and after
+      !!                - look for the closest wet cells (share amoung neigbourgs if there are)
+      !!                - build the correction increment to applied at each time step
+      !!                
+      !!----------------------------------------------------------------------
+      !
+      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
+      !
+      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
+      !!----------------------------------------------------------------------
+
+      ! 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 )
+
+      ! compute run length
+      nstp_iscpl  = nitend - nit000 + 1
+      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
+      !==============================================================================
+
+      risfcpl_cons_vol = 0.0
+      risfcpl_cons_ssh = 0.0
+      risfcpl_cons_tsc = 0.0
+
+      DO jk = 1,jpk-1
+         DO jj = nldj,nlej
+            DO ji = nldi,nlei
+
+               ! volume diff
+               zdvol = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+
+               ! 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)
+
+               ! 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)
+            
+               ! 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
+
+            END DO
+         END DO
+      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
+      ! compute the total number of point receiving a correction increment for each processor
+      ! local
+      nisfl=0
+      DO jk = 1,jpk-1
+         DO jj = nldj,nlej
+            DO ji = nldi,nlei
+               jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
+               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      ! global 
+      CALL mpp_sum('isfcpl',nisfl)
+      nisfg  = SUM(nisfl           )
+      istart = SUM(nisfl(1:narea-1))
+      iend   = SUM(nisfl(1:narea  ))
+      !
+      ! allocate list of point receiving correction
+      ALLOCATE(zisfpts(nisfg))
+      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
+      DO jk = 1,jpk-1
+         DO jj = nldj,nlej
+            DO ji = nldi,nlei
+               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
+
+                  jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
+
+                  zdvol = risfcpl_cons_vol(ji,jj,jk       )
+                  zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
+                  zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
+
+                  IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
+                     ! spread correction amoung neigbourg wet cells (horizontal direction first)
+                     ! as it is a rude correction corner and lateral cell have the same weight
+                     !
+                     z1_sum =  1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
+                     !
+                     ! lateral cells
+                     IF (tmask(jip1,jj  ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj  , jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jim1,jj  ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj  , jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(ji  ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji  , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(ji  ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji  , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     !
+                     ! corner  cells
+                     IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     !
+                  ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
+                     ! spread correction amoung neigbourg wet cells (vertical direction)
+                     CALL update_isfpts(zisfpts, jisf, ji  , jj  , jk+1, zdvol, zdsal, zdtem, 1., 0)
+                  ELSE
+                     ! need to find where to put correction in later on
+                     CALL update_isfpts(zisfpts, jisf, ji  , jj  , jk  , zdvol, zdsal, zdtem, 1., 1)
+                  END IF
+               END IF
+            END DO
+         END DO
+      END DO
+      !
+      ! 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
+      !
+      ! mask (>0 out)
+      risfcpl_cons_vol(:,:,:       ) = risfcpl_cons_vol(:,:,:       ) * tmask(:,:,:)
+      risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
+      risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
+      !
+      ! 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.)
+      !
+      ! ssh correction (for dynspg_ts)
+      DO jk = 1,jpk
+         risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
+      END DO
+      risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
+      !
+   END SUBROUTINE isfcpl_cons
+   !
+   SUBROUTINE update_isfpts(sisfpts, kpts, ki, kj, kk, pdvol, pdsal, pdtem, pratio, kfind)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE update_isfpts  ***
+      !!
+      !! ** Purpose : if a cell become dry, we need to put the corrective increment elsewhere
+      !!
+      !! ** Action  : update the list of point
+      !!
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      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 :: ifind
+      !!----------------------------------------------------------------------
+      ! 
+      ! increment position
+      kpts = kpts + 1
+      !
+      ! define if we need to look for closest valid wet cell (no neighbours or neigbourg on halo)
+      IF ( PRESENT(kfind) ) THEN
+         ifind = kfind
+      ELSE
+         ifind = ( 1 - tmask_h(ki,kj) ) * tmask(ki,kj,kk)
+      END IF
+      !
+      ! update isfpts structure
+      sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
+      !
+   END SUBROUTINE update_isfpts
+   !
+   SUBROUTINE get_correction( ki, kj, kk, zlon, zlat, pvolinc, psalinc, pteminc, kfind)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE get_correction  ***
+      !!
+      !! ** Action : - Find the closest valid cell if needed (wet and not on the halo)
+      !!             - Scale the correction depending of pratio (case where multiple wet neigbourgs)
+      !!             - Fill the correction array
+      !!
+      !!----------------------------------------------------------------------
+      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 :: jj, ji, iig, ijg
+      !!----------------------------------------------------------------------
+      !
+      ! define global indice of correction location
+      iig = ki ; ijg = kj
+      IF ( kfind == 1 ) CALL dom_ngb( zlon, zlat, 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
+            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
+            risfcpl_cons_tsc(ji,jj,kk,jp_tem) =  risfcpl_cons_tsc(ji,jj,kk,jp_tem) + pteminc
+         END DO
+      END DO
+
+   END SUBROUTINE get_correction
+
+END MODULE isfcpl
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90	(revision 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90	(revision 11553)
@@ -73,8 +73,4 @@
       !!----------------------------------------------------------------------
       !
-      ! print
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '   ice shelf case: set the ice-shelf load following isomip case'
-      !
       znad = 1._wp                     !- To use density and not density anomaly
       !
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 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 11553)
@@ -39,5 +39,5 @@
    PRIVATE
 
-   PUBLIC   isf_stp, isf_stp_init  ! routine called in sbcmod and divhor
+   PUBLIC   isf_stp, isf_init, isf_nam  ! routine called in sbcmod and divhor
 
    !!----------------------------------------------------------------------
@@ -114,5 +114,5 @@
    END SUBROUTINE isf_stp
 
-   SUBROUTINE isf_stp_init
+   SUBROUTINE isf_init
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isfstp_init  ***
@@ -127,36 +127,64 @@
       !!              - call cav/param init routine
       !!----------------------------------------------------------------------
-      INTEGER               :: inum, ierror
-      INTEGER               :: ios                  ! Local integer output status for namelist read
-      INTEGER               :: ikt, ikb
-      INTEGER               :: ji, jj
-      !!----------------------------------------------------------------------
-      NAMELIST/namisf/ ln_isf       ,                                                                               & 
-         &             ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf,   &
-         &             ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff, &
-         &             ln_isfcpl    , nn_drown, ln_isfcpl_cons,                                                     &
-         &             cn_isfload   , cn_isfdir
-      !!----------------------------------------------------------------------
+      !
+      ! constrain: l_isfoasis need to be known
+      !
+      ! Read namelist
+      CALL isf_nam()
       !
       ! Allocate public array
       CALL isf_alloc()
       !
-      ! initalisation of fwf and tsc array to 0
-      riceload(:,:)       = 0.0_wp
-      fwfisf_oasis(:,:)   = 0.0_wp
-      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
-      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
-      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
-      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
-      !
-      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
-      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist', lwp )
-      !
-      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
-      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namisf )
-      !
+      ! check option compatibility
+      CALL isf_ctl()
+      !
+      ! compute ice shelf load
+      IF ( ln_isfcav ) CALL isf_load( risfload )
+      !
+      ! terminate routine now if no ice shelf melt formulation specify
+      IF ( ln_isf ) THEN
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation melt in the cavity
+         IF ( ln_isfcav_mlt ) THEN
+            !
+            ! initialisation  of cav variable
+            CALL isf_cav_init()
+            !
+            ! read cav variable from restart
+            IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
+            !
+         END IF
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation parametrised melt
+         IF ( ln_isfpar_mlt ) THEN
+            !
+            ! initialisation  of par variable
+            CALL isf_par_init()
+            !
+            ! read par variable from restart
+            IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
+            !
+         END IF
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation ice sheet coupling
+         IF( ln_isfcpl ) THEN 
+            !
+            ! prepare writing restart
+            IF( lwxios ) CALL iom_set_rstw_var_active('ssmask')
+            IF( lwxios ) CALL iom_set_rstw_var_active('tmask')
+            IF( lwxios ) CALL iom_set_rstw_var_active('e3t_n')
+            IF( lwxios ) CALL iom_set_rstw_var_active('e3u_n')
+            IF( lwxios ) CALL iom_set_rstw_var_active('e3v_n')
+            !
+         END IF
+         !
+      END IF
+         
+  END SUBROUTINE isf_init
+
+  SUBROUTINE isf_ctl()!
       IF (lwp) THEN
          WRITE(numout,*)
@@ -166,4 +194,5 @@
          !
          WRITE(numout,*) '   ice shelf cavity (open or parametrised)  ln_isf = ', ln_isf
+         WRITE(numout,*)
          !
          IF ( ln_isf ) THEN
@@ -181,5 +210,5 @@
             WRITE(numout,*) ''
             !
-            WRITE(numout,*) '      ice shelf melt parametrisation          ln_isfpar_mlt    = ', ln_isfpar_mlt
+            WRITE(numout,*) '      ice shelf melt parametrisation          ln_isfpar_mlt   = ', ln_isfpar_mlt
             IF ( ln_isfpar_mlt ) THEN
                WRITE(numout,*) '         isf parametrisation formulation         cn_isfpar_mlt   = ', TRIM(cn_isfpar_mlt)
@@ -187,8 +216,8 @@
             WRITE(numout,*) ''
             !
-            WRITE(numout,*) '      Coupling to an ice sheet model          ln_isfcpl        = ', ln_isfcpl
+            WRITE(numout,*) '      Coupling to an ice sheet model          ln_isfcpl       = ', ln_isfcpl
             IF ( ln_isfcpl ) THEN
-               WRITE(numout,*) '         conservation activated ln_isfcpl_cons           = ', ln_isfcpl_cons
-               WRITE(numout,*) '            number of call of the extrapolation loop = ', nn_drown
+               WRITE(numout,*) '         conservation activated ln_isfcpl_cons     = ', ln_isfcpl_cons
+               WRITE(numout,*) '         number of call of the extrapolation loop  = ', nn_drown
             ENDIF
             WRITE(numout,*) ''
@@ -208,21 +237,10 @@
       END IF
       !
+
       !---------------------------------------------------------------------------------------------------------------------
-      ! initialisation ice shelf load
-      IF ( ln_isfcav ) THEN
-         !
-         ! compute ice shelf mask
-         mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
-         !
-         ! compute ice shelf load
-         CALL isf_load( risfload )
-         !
-      END IF
-      !
-      !---------------------------------------------------------------------------------------------------------------------
-      ! sanity check
+      ! sanity check  ! issue ln_isfcav not yet known as well as l_isfoasis  => move this call in isf_stp ?
       ! melt in the cavity without cavity
-      IF ( ln_isfcav_mlt .AND. (.NOT. ln_isfcav) ) &
-         &   CALL ctl_stop('ice shelf melt in the cavity activated (ln_isfcav_mlt) but no cavity detected in domcfg (ln_isfcav), STOP' )
+       IF ( ln_isfcav_mlt .AND. (.NOT. ln_isfcav) ) &
+          &   CALL ctl_stop('ice shelf melt in the cavity activated (ln_isfcav_mlt) but no cavity detected in domcfg (ln_isfcav), STOP' )
       !
       ! ice sheet coupling without cavity
@@ -233,5 +251,7 @@
          &   CALL ctl_stop( 'The coupling between NEMO and an ice sheet model with the conservation option does not work with the linssh option' )
       !
-      IF ( l_isfoasis ) THEN
+      IF ( l_isfoasis .AND. .NOT. ln_isf ) CALL ctl_stop( ' OASIS send ice shelf fluxes to NEMO but NEMO does not have the isf module activated' )
+      !
+      IF ( l_isfoasis .AND. ln_isf ) THEN
          !
          CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
@@ -251,72 +271,35 @@
             IF ( TRIM(cn_isfpar_mlt) == 'oasis' .AND. TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis and cn_isfcav_mlt = oasis not coded' )
          END IF
-      END IF
-      !
-      ! terminate routine now if no ice shelf melt formulation specify
-      IF ( ln_isf ) THEN
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation melt in the cavity
-         IF ( ln_isfcav_mlt ) THEN
-            !
-            ! initialisation  of cav variable
-            CALL isf_cav_init()
-            !
-            ! read cav variable from restart
-            IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
-            !
-         END IF
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation parametrised melt
-         IF ( ln_isfpar_mlt ) THEN
-            !
-            ! initialisation  of par variable
-            CALL isf_par_init()
-            !
-            ! read par variable from restart
-            IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
-            !
-         END IF
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation ice sheet coupling
-         IF( ln_isfcpl ) THEN
-            ! 
-            ! prepare writing restart
-            IF( lwxios ) CALL iom_set_rstw_var_active('ssmask')
-            IF( lwxios ) CALL iom_set_rstw_var_active('tmask')
-            IF( lwxios ) CALL iom_set_rstw_var_active('e3t_n')
-            IF( lwxios ) CALL iom_set_rstw_var_active('e3u_n')
-            IF( lwxios ) CALL iom_set_rstw_var_active('e3v_n')
-            !
-            IF( ln_rstart ) THEN
-               !
-               CALL isf_alloc_cpl()
-               !
-               ! extrapolation tracer properties
-               CALL isfcpl_tra()
-               !
-               ! correction of the horizontal divergence and associated temp. and salt content flux
-               CALL isfcpl_vol()
-               !
-               ! apply the 'conservation' method
-               IF ( ln_isfcpl_cons ) CALL isfcpl_cons()
-               !
-               ! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
-               !           - decide how to manage thickness level change in conservation
-               !
-               tsb    (:,:,:,:) = tsn (:,:,:,:)
-               sshb   (:,:)     = sshn(:,:)
-               ub     (:,:,:)   = un(:,:,:)
-               vb     (:,:,:)   = vn(:,:,:)
-               !
-            END IF
-            !
-         END IF
-         !
-      END IF
-         
-  END SUBROUTINE isf_stp_init
+         !
+      END IF
+   END SUBROUTINE isf_ctl
+   !
+   SUBROUTINE isf_nam
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_nam  ***
+      !!
+      !! ** Purpose :   Read ice shelf namelist cfg and ref
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER               :: ios                  ! Local integer output status for namelist read
+      !!----------------------------------------------------------------------
+      NAMELIST/namisf/ ln_isf       ,                                                                               & 
+         &             ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf,   &
+         &             ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff, &
+         &             ln_isfcpl    , nn_drown, ln_isfcpl_cons,                                                     &
+         &             cn_isfload   , cn_isfdir
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namisf )
+
+   END SUBROUTINE isf_nam
+   !!
    !!======================================================================
 END MODULE isfstp
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 11553)
@@ -474,5 +474,5 @@
       srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE.
 
-      IF( srcv(jpr_isf)%laction .AND. ln_isf ) THEN
+      IF( srcv(jpr_isf)%laction ) THEN
          l_isfoasis = .TRUE.  ! -> isf fwf comes from oasis
          IF(lwp) WRITE(numout,*)
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90	(revision 11549)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90	(revision 11553)
@@ -435,4 +435,5 @@
          RETURN                                       ! end of initialization
       ENDIF
+      !
       
                            CALL  istate_init    ! ocean initial state (Dynamics and tracers)
@@ -453,5 +454,4 @@
       !                                      ! Active tracers
       IF( ln_traqsr    )   CALL tra_qsr_init      ! penetrative solar radiation qsr
-                           CALL isf_stp_init      ! ice shelf
                            CALL tra_bbc_init      ! bottom heat flux
                            CALL tra_bbl_init      ! advective (and/or diffusive) bottom boundary layer scheme
@@ -476,4 +476,7 @@
       !                                      ! Icebergs
                            CALL icb_init( rdt, nit000)   ! initialise icebergs instance
+
+                                                ! ice shelf
+                           CALL isf_init
 
       !                                      ! Misc. options
