Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90	(revision 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90	(revision 11423)
@@ -906,14 +906,20 @@
       CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
       CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep          )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity
+      IF ( ln_isf ) THEN
+         IF (ln_isfcav_mlt) THEN
+            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity
+         END IF
+         IF (ln_isfpar_mlt) THEN
+            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity
+         END IF
+      END IF
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domvvl.F90	(revision 11423)
@@ -19,4 +19,5 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
+   USE isfcpl
    USE phycst          ! physical constant
    USE dom_oce         ! ocean space and time domain
@@ -813,4 +814,15 @@
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            !
+            ! coupling with an ice sheet model (grounding/calving line migration)
+            IF ( ln_iscpl ) THEN
+               !
+               ! extrapolate ssh over new water column if needed
+               CALL isfcpl_ssh()
+               !
+               ! everything need to compute over new water column
+               id1 = 0 ; id2 = 0 ; id3 = 0 ; id4 = 0 ; id5 = 0
+               !
+            END IF
             !                             ! --------- !
             !                             ! all cases !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/istate.F90	(revision 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/istate.F90	(revision 11423)
@@ -28,5 +28,4 @@
    USE dtauvd         ! data: U & V current             (dta_uvd routine)
    USE domvvl          ! varying vertical mesh
-   USE iscplrst        ! ice sheet coupling
    USE wet_dry         ! wetting and drying (needed for wad_istate)
    USE usrdef_istate   ! User defined initial state
@@ -86,5 +85,4 @@
          !                                    ! -------------------
          CALL rst_read                        ! Read the restart file
-         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90	(revision 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90	(revision 11423)
@@ -23,6 +23,4 @@
    USE sbcrnf          ! river runoff 
    USE isfhdiv         ! ice shelf
-   USE iscplhsb        ! ice sheet / ocean coupling
-   USE iscplini        ! ice sheet / ocean coupling
 #if defined key_asminc   
    USE asminc          ! Assimilation increment
@@ -100,5 +98,5 @@
       ! 
 #endif
-      IF( ln_isf )   CALL isf_hdiv( hdivn )      !==  ice shelf  ==!   (update hdivn field)
+      IF( ln_isf )   CALL isf_hdiv( kt, hdivn )      !==  ice shelf  ==!   (update hdivn field)
       !
       CALL lbc_lnk( 'divhor', hdivn, 'T', 1. )   !   (no sign change)
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplhsb.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplhsb.F90	(revision 11403)
+++ 	(revision )
@@ -1,324 +1,0 @@
-MODULE iscplhsb
-   !!======================================================================
-   !!                       ***  MODULE  iscplhsb  ***
-   !! Ocean forcing: ice sheet/ocean coupling (conservation)
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_alloc    : variable allocation
-   !!   iscpl_hsb      : compute and store the input of heat/salt/volume 
-   !!                    into the system due to the coupling process
-   !!   iscpl_div      : correction of divergence to keep volume conservation
-   !!----------------------------------------------------------------------
-   USE oce             ! global tra/dyn variable
-   USE dom_oce         ! ocean space and time domain
-   USE domwri          ! ocean space and time domain
-   USE domngb          ! 
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE iscplini        ! 
-   !
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE lbclnk          !
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_div   
-   PUBLIC   iscpl_cons       
-   !! * Substitutions  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_cons  ***
-      !! 
-      !! ** Purpose :   compute input into the system during the coupling step
-      !!                compute the correction term
-      !!                compute where the correction have to be applied
-      !! 
-      !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before
-      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx     !! corrective flux to have tracer conservation
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(out) :: pvol_flx    !! corrective flux to have volume conservation
-      REAL(wp),                     INTENT(in ) :: prdt_iscpl  !! coupling period 
-      !
-      INTEGER  ::   ji  , jj  , jk           ! loop index
-      INTEGER  ::   jip1, jim1, jjp1, jjm1
-      REAL(wp) ::   summsk, zsum , zsumn, zjip1_ratio  , zjim1_ratio, zdtem, zde3t, z1_rdtiscpl
-      REAL(wp) ::   zarea , zsum1, zsumb, zjjp1_ratio  , zjjm1_ratio, zdsal
-      REAL(wp), DIMENSION(jpi,jpj)        ::   zdssh   ! workspace
-      REAL(wp), DIMENSION(:), ALLOCATABLE ::   zlon, zlat
-      REAL(wp), DIMENSION(:), ALLOCATABLE ::   zcorr_vol, zcorr_tem, zcorr_sal
-      INTEGER , DIMENSION(:), ALLOCATABLE ::   ixpts, iypts, izpts, inpts
-      INTEGER :: jpts, npts
-      !!----------------------------------------------------------------------
-
-      ! get imbalance (volume heat and salt)
-      ! initialisation difference
-      zde3t = 0._wp   ;   zdsal = 0._wp   ;   zdtem = 0._wp
-
-      ! initialisation correction term
-      pvol_flx(:,:,:  ) = 0._wp
-      pts_flx (:,:,:,:) = 0._wp
-      
-      z1_rdtiscpl = 1._wp / prdt_iscpl 
-
-      ! mask tsn and tsb 
-      tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) *  tmask  (:,:,:)
-      tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) *  tmask  (:,:,:)
-
-      !==============================================================================
-      ! diagnose the heat, salt and volume input and compute the correction variable
-      !==============================================================================
-
-      ! 
-      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
-      
-      DO jk = 1,jpk-1
-         DO jj = 2,jpj-1
-            DO ji = fs_2,fs_jpim1
-               IF (tmask_h(ji,jj) == 1._wp) THEN
-
-                  ! volume differences
-                  zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
-
-                  ! heat diff
-                  zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-                  ! salt diff
-                  zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-               
-                  ! shh changes
-                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN 
-                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
-                     zdssh(ji,jj) = 0._wp
-                  END IF
-
-                  ! volume, heat and salt differences in each cell 
-                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t * z1_rdtiscpl
-                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * z1_rdtiscpl 
-                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * z1_rdtiscpl
-
-                  ! case where we close a cell: check if the neighbour cells are wet 
-                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
-
-                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
-
-                     zsum =   e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
-                       &    + e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk)
-
-                     IF ( zsum /= 0._wp ) THEN
-                        zjip1_ratio   = e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
-                        zjim1_ratio   = e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
-                        zjjp1_ratio   = e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
-                        zjjm1_ratio   = e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
-
-                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
-                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
-                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
-                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
-                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
-                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
-                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
-                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
-                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
-                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
-                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
-                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
-
-                        ! set to 0 the cell we distributed over neigbourg cells
-                        pvol_flx(ji,jj,jk       ) = 0._wp
-                        pts_flx (ji,jj,jk,jp_sal) = 0._wp
-                        pts_flx (ji,jj,jk,jp_tem) = 0._wp
-
-                     ELSE IF (zsum == 0._wp ) THEN
-                        ! case where we close a cell and no adjacent cell open
-                        ! check if the cell beneath is wet
-                        IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
-                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
-                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
-                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
-
-                           ! set to 0 the cell we distributed over neigbourg cells
-                           pvol_flx(ji,jj,jk       ) = 0._wp
-                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
-                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
-                        ELSE
-                        ! case no adjacent cell on the horizontal and on the vertical
-                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
-                                             ! open for narea== 1 (lwp=T)
-                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
-                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
-                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
-                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
-                           ENDIF
-                        ! We deal with these points later.
-                        END IF
-                     END IF
-                  END IF
-               END IF
-            END DO
-         END DO
-      END DO
-
-!!gm  ERROR !!!!
-!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
-!
-!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
-!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
-!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
-      STOP ' iscpl_cons:   please modify this module !'
-!!gm end
-      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point 
-      ! allocation and initialisation of the list of problematic point
-      ALLOCATE( inpts(jpnij) )
-      inpts(:) = 0
-
-      ! fill narea location with the number of problematic point
-      DO jk = 1,jpk-1
-         DO jj = 2,jpj-1
-            DO ji = fs_2,fs_jpim1
-               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
-                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
-                  inpts(narea) = inpts(narea) + 1 
-               END IF
-            END DO
-         END DO
-      END DO
-
-      ! build array of total problematic point on each cpu (share to each cpu)
-      CALL mpp_max('iscplhsb', inpts,jpnij) 
-
-      ! size of the new variable
-      npts  = SUM(inpts)    
-      
-      ! allocation of the coordinates, correction, index vector for the problematic points 
-      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
-      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
-      zcorr_vol(:) = -1.0e20_wp
-      zcorr_sal(:) = -1.0e20_wp
-      zcorr_tem(:) = -1.0e20_wp
-
-      ! fill new variable
-      jpts = SUM(inpts(1:narea-1))
-      DO jk = 1,jpk-1
-         DO jj = 2,jpj-1
-            DO ji = fs_2,fs_jpim1
-               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
-                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
-                  jpts = jpts + 1  ! positioning in the inpts vector for the area narea
-                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
-                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
-                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
-                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
-                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
-
-                  ! set flx to 0 (safer)
-                  pvol_flx(ji,jj,jk       ) = 0.0_wp
-                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
-                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
-               END IF
-            END DO
-         END DO
-      END DO
-
-      ! build array of total problematic point on each cpu (share to each cpu)
-      ! point coordinates
-      CALL mpp_max('iscplhsb', zlat ,npts)
-      CALL mpp_max('iscplhsb', zlon ,npts)
-      CALL mpp_max('iscplhsb', izpts,npts)
-
-      ! correction values 
-      CALL mpp_max('iscplhsb', zcorr_vol,npts)
-      CALL mpp_max('iscplhsb', zcorr_sal,npts)
-      CALL mpp_max('iscplhsb', zcorr_tem,npts)
-
-      ! put correction term in the closest cell          
-      DO jpts = 1,npts
-         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
-         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
-            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
-               jk = izpts(jpts)
-
-               IF (tmask_h(ji,jj) == 1._wp) THEN
-                  ! correct the vol_flx in the closest cell
-                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
-                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
-                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
-
-                  ! set correction to 0
-                  zcorr_vol(jpts) = 0.0_wp
-                  zcorr_sal(jpts) = 0.0_wp
-                  zcorr_tem(jpts) = 0.0_wp
-               END IF
-            END DO
-         END DO
-      END DO
-
-      ! deallocate variables 
-      DEALLOCATE(inpts)
-      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
-    
-      ! add contribution store on the hallo (lbclnk remove one of the contribution)
-      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
-      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
-      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
-
-!!gm  ERROR !!!!
-!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
-!
-!      ! compute sum over the halo and set it to 0.
-!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
-!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
-!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
-!!gm end
-
-      !
-   END SUBROUTINE iscpl_cons
-
-
-   SUBROUTINE iscpl_div( phdivn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE iscpl_div  ***
-      !!
-      !! ** Purpose :   update the horizontal divergenc
-      !!
-      !! ** Method  :
-      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
-      !!
-      !! ** Action  :   phdivn   increase by the iscpl correction term
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE iscpl_div
-
-END MODULE iscplhsb
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplini.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplini.F90	(revision 11403)
+++ 	(revision )
@@ -1,89 +1,0 @@
-MODULE iscplini
-   !!======================================================================
-   !!                       ***  MODULE  sbciscpl  ***
-   !! Ocean forcing:  ?????
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_init     : initialisation routine (namelist)
-   !!   iscpl_alloc    : allocation of correction variables
-   !!----------------------------------------------------------------------
-   USE oce             ! global tra/dyn variable
-   USE dom_oce         ! ocean space and time domain
-   !
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE in_out_manager  ! I/O manager
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_init      
-   PUBLIC   iscpl_alloc 
-   !                                 !!* namsbc_iscpl namelist *
-   LOGICAL , PUBLIC ::   ln_iscpl_hsb !:
-   INTEGER , PUBLIC ::   nn_fiscpl    !:
-   INTEGER , PUBLIC ::   nn_drown     !:
-   ! 
-   INTEGER , PUBLIC ::   nstp_iscpl   !:
-   REAL(wp), PUBLIC ::   rdt_iscpl    !: 
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdiv_iscpl   !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   htsc_iscpl   !:
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION iscpl_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE sbc_iscpl_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=iscpl_alloc )
-         !
-      CALL mpp_sum ( 'iscplini', iscpl_alloc )
-      IF( iscpl_alloc > 0 )   CALL ctl_warn('iscpl_alloc: allocation of arrays failed')
-   END FUNCTION iscpl_alloc
-
-
-   SUBROUTINE iscpl_init()
-      !!----------------------------------------------------------------------
-      !!----------------------------------------------------------------------
-      INTEGER ::   ios           ! Local integer output status for namelist read
-      NAMELIST/namsbc_iscpl/ nn_fiscpl, ln_iscpl_hsb, nn_drown
-      !!----------------------------------------------------------------------
-      !
-      nn_fiscpl = 0
-      ln_iscpl_hsb    = .FALSE.
-      REWIND( numnam_ref )              ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling
-      READ  ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_iscpl in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namsbc_iscpl in configuration namelist : Ice Sheet coupling
-      READ  ( numnam_cfg, namsbc_iscpl, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_iscpl in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namsbc_iscpl )
-      !
-      nstp_iscpl=MIN( nn_fiscpl, nitend-nit000+1 ) ! the coupling period have to be less or egal than the total number of time step
-      rdt_iscpl = nstp_iscpl * rn_rdt
-      !
-      IF (lwp) THEN
-         WRITE(numout,*) 'iscpl_rst:'
-         WRITE(numout,*) '~~~~~~~~~'
-         WRITE(numout,*) ' coupling     flag (ln_iscpl    )            = ', ln_iscpl
-         WRITE(numout,*) ' conservation flag (ln_iscpl_hsb)            = ', ln_iscpl_hsb
-         WRITE(numout,*) ' nb of stp for cons (rn_fiscpl  )            = ', nstp_iscpl
-         IF (nstp_iscpl .NE. nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified &
-            &                                           (larger than run length)'
-         WRITE(numout,*) ' coupling time step                       = ', rdt_iscpl
-         WRITE(numout,*) ' number of call of the extrapolation loop = ', nn_drown
-      ENDIF
-      !
-   END SUBROUTINE iscpl_init
-
-   !!======================================================================
-END MODULE iscplini
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplrst.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplrst.F90	(revision 11403)
+++ 	(revision )
@@ -1,407 +1,0 @@
-MODULE iscplrst
-   !!======================================================================
-   !!                       ***  MODULE  iscplrst  ***
-   !! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_stp          : step management
-   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
-   !!----------------------------------------------------------------------
-   USE oce             ! global tra/dyn variable
-   USE dom_oce         ! ocean space and time domain
-   USE domwri          ! ocean space and time domain
-   USE domvvl   , ONLY : dom_vvl_interpol
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE iscplini        ! ice sheet coupling: initialisation
-   USE iscplhsb        ! ice sheet coupling: conservation
-   !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE lbclnk          ! communication
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_stp          ! step management 
-   !!
-   !! * Substitutions  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE iscpl_stp
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_stp  ***
-      !! 
-      !! ** Purpose : compute initialisation
-      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
-      !!              compute correction term if needed
-      !! 
-      !!----------------------------------------------------------------------
-      INTEGER  ::   inum0
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zsmask_b
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztmask_b, zumask_b, zvmask_b
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t_b  , ze3u_b  , ze3v_b  
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepw_b
-      CHARACTER(20) :: cfile
-      !!----------------------------------------------------------------------
-      !
-      !                       ! 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, 'umask'  , zumask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:), ldxios = lrxios )  ! need to compute temperature correction
-      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
-      !
-      CALL iscpl_init()       ! read namelist
-      !                       ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
-      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
-      !
-      IF ( ln_iscpl_hsb ) THEN      ! compute correction if conservation needed
-         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
-         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
-      END IF
-         
-      !                       ! create  a domain file
-      IF( ln_meshmask .AND. ln_iscpl )   CALL dom_wri
-      !
-      IF ( ln_iscpl_hsb ) THEN
-         cfile='correction'
-         cfile = TRIM( cfile )
-         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE. )
-         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
-         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
-         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
-         CALL iom_close ( inum0 )
-      END IF
-      !
-      neuler = 0              ! next step is an euler time step
-      !
-      !                       ! set _b and _n variables equal
-      tsb (:,:,:,:) = tsn (:,:,:,:)
-      ub  (:,:,:)   = un  (:,:,:)
-      vb  (:,:,:)   = vn  (:,:,:)
-      sshb(:,:)     = sshn(:,:)
-      !
-      !                       ! set _b and _n vertical scale factor equal
-      e3t_b (:,:,:) = e3t_n (:,:,:)
-      e3u_b (:,:,:) = e3u_n (:,:,:)
-      e3v_b (:,:,:) = e3v_n (:,:,:)
-      !
-      e3uw_b (:,:,:) = e3uw_n (:,:,:)
-      e3vw_b (:,:,:) = e3vw_n (:,:,:)
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-      hu_b   (:,:)   = hu_n   (:,:)
-      hv_b   (:,:)   = hv_n   (:,:)
-      r1_hu_b(:,:)   = r1_hu_n(:,:)
-      r1_hv_b(:,:)   = r1_hv_n(:,:)
-      !
-   END SUBROUTINE iscpl_stp
-
-
-   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_rst_interpol  ***
-      !! 
-      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves 
-      !!                compute 2d fields of heat, salt and volume correction
-      !! 
-      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
-      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
-      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
-      !!
-      INTEGER :: ji, jj, jk, iz          !! loop index
-      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
-      !!
-      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
-      REAL(wp):: zdz, zdzm1, zdzp1
-      !!
-      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
-      REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn, ztrp
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
-      !!----------------------------------------------------------------------
-      !
-      !                 ! mask value to be sure
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
-      !
-      !                 ! compute wmask
-      zwmaskn(:,:,1) = tmask   (:,:,1)
-      zwmaskb(:,:,1) = ptmask_b(:,:,1)
-      DO jk = 2,jpk
-         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
-         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
-      END DO
-      !    
-      !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
-      sshb (:,:)=sshn(:,:)
-      zssh0(:,:)=sshn(:,:)
-      zsmask0(:,:) = psmask_b(:,:)
-      zsmask1(:,:) = psmask_b(:,:)
-      DO iz = 1, 10                 ! need to be tuned (configuration dependent) (OK for ISOMIP+)
-         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
-         DO jj = 2,jpj-1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               jip1=ji+1; jim1=ji-1;
-               jjp1=jj+1; jjm1=jj-1;
-               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
-               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
-                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
-                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
-                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
-                  zsmask1(ji,jj)=1._wp
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk_multi( 'iscplrst', sshn, 'T', 1., zsmask1, 'T', 1. )
-         zssh0   = sshn
-         zsmask0 = zsmask1
-      END DO
-      sshn(:,:) = sshn(:,:) * ssmask(:,:)
-
-!=============================================================================
-!PM: Is this needed since introduction of VVL by default?
-      IF ( .NOT.ln_linssh ) THEN
-      ! Reconstruction of all vertical scale factors at now time steps
-      ! =============================================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-         DO jk = 1,jpk
-            DO jj=1,jpj
-               DO ji=1,jpi
-                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
-                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   &
-                     &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-
-         ! Vertical scale factor interpolations
-         ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-         
-         ! t- and w- points depth
-         ! ----------------------
-         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-         gdepw_n(:,:,1) = 0.0_wp
-         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               DO jk = 2,mikt(ji,jj)-1
-                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
-               END DO
-               IF (mikt(ji,jj) > 1) THEN
-                  jk = mikt(ji,jj)
-                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
-               END IF
-               DO jk = mikt(ji,jj)+1, jpk
-                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
-               END DO
-            END DO
-         END DO
-
-      ! t-, u- and v- water column thickness
-      ! ------------------------------------
-         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         !                                        ! Inverse of the local depth
-         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
-         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
-
-      END IF
-
-!=============================================================================
-! compute velocity
-! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
-      DO jk = 1,jpk
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
-               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
-            END DO
-         END DO
-      END DO
-
-! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
-! compute barotropic velocity now and after 
-      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
-      zbub(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
-      zbvb(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
-      zbun(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
-      zbvn(:,:)   = SUM(ztrp,DIM=3)
-
-      ! new water column
-      zhu1=0.0_wp ;
-      zhv1=0.0_wp ;
-      DO jk  = 1,jpk
-        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      
-      ! compute correction      
-      zucorr = 0._wp
-      zvcorr = 0._wp
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
-               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
-            END IF
-            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
-               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
-            END IF
-         END DO
-      END DO 
-      
-      ! update velocity
-      DO jk  = 1,jpk
-         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
-         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
-      END DO
-
-!=============================================================================
-      ! compute temp and salt
-      ! compute new tn and sn if we open a new cell
-      tsb (:,:,:,:) = tsn(:,:,:,:)
-      zts0(:,:,:,:) = tsn(:,:,:,:)
-      ztmask1(:,:,:) = ptmask_b(:,:,:)
-      ztmask0(:,:,:) = ptmask_b(:,:,:)
-      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
-          DO jk = 1,jpk-1
-              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
-              DO jj = 2,jpj-1
-                 DO ji = fs_2,fs_jpim1
-                      jip1=ji+1; jim1=ji-1;
-                      jjp1=jj+1; jjm1=jj-1;
-                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
-                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 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) ) / summsk
-                         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) ) / summsk
-                         ztmask1(ji,jj,jk)=1
-                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
-                      !! vertical extrapolation if horizontal extrapolation failed
-                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
-                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
-                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 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))/summsk
-                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
-                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
-                            ztmask1(ji,jj,jk)=1._wp
-                         END IF
-                      END IF
-                  END DO
-              END DO
-          END DO
-          
-          CALL lbc_lnk_multi( 'iscplrst', tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
-
-          ! update
-          zts0(:,:,:,:) = tsn(:,:,:,:)
-          ztmask0 = ztmask1
-
-      END DO
-
-      ! mask new tsn field
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
-
-      ! 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
-      IF (.NOT.ln_linssh) THEN
-         DO jk = 2,jpk-1
-            DO jj = 1,jpj
-               DO ji = 1,jpi
-                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
-                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
-                     !compute weight
-                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
-                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
-                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - 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
-
-      ! closed pool
-      ! -----------------------------------------------------------------------------------------
-      ! case we open a cell but no neigbour cells available to get an estimate of T and S
-      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
-         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
-         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
-         umask(:,:,:) = 0._wp
-         vmask(:,:,:) = 0._wp
-      END WHERE
-      
-      ! set mbkt and mikt to 1 in thiese location
-      WHERE (SUM(tmask,dim=3) == 0)
-         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
-         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
-      END WHERE
-      ! -------------------------------------------------------------------------------------------
-      ! compute new tn and sn if we close cell 
-      ! nothing to do
-      ! 
-   END SUBROUTINE iscpl_rst_interpol
-
-   !!======================================================================
-END MODULE iscplrst
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11423)
@@ -22,45 +22,40 @@
    PRIVATE
 
-   PUBLIC   isf_alloc, isf_alloc_par, isf_alloc_cav
-
-   ! public in order to be able to output then 
-
-   LOGICAL, PUBLIC :: ln_isfpar_mlt                  !: logical for the computation of melt inside the cavity
-   LOGICAL, PUBLIC :: ln_isfcav_mlt                  !: logical for the use of ice shelf parametrisation
-   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
-   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient    []
-   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient    []
-   REAL(wp), PUBLIC ::   rn_htbl                     !: Losch top boundary layer thickness [m]
-   CHARACTER(LEN=256), PUBLIC :: cn_isfload          !: ice shelf load computation method
-   CHARACTER(LEN=256), PUBLIC :: cn_gammablk         !: gamma formulation
-   CHARACTER(LEN=256), PUBLIC :: cn_isfcav_mlt, cn_isfpar_mlt !: melt formulation (cavity/param)
-   TYPE(FLD_N), PUBLIC                                  :: sn_isfcav_fwf   !: information about the isf melting file to be read
-   TYPE(FLD_N), PUBLIC                                  :: sn_isfpar_fwf   !: information about the isf melting file to be read
-   TYPE(FLD_N), PUBLIC                                  :: sn_isfpar_zmax  !: information about the grounding line depth file to be read
-   TYPE(FLD_N), PUBLIC                                  :: sn_isfpar_zmin  !: information about the calving   line depth file to be read
-   TYPE(FLD_N), PUBLIC                                  :: sn_isfpar_Leff  !: information about the effective length     file to be read
-
-   LOGICAL, PUBLIC :: l_isfcpl
-   !
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt_par, misfkt_cav   !: Level of ice shelf base
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkb_par, misfkb_cav   !: Level of ice shelf base
-   !
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mskisf_par, mskisf_cav   !: Level of ice shelf base
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfload                     !: ice shelf load
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl_0                  !: thickness of tbl (initial value)  [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl_cav, rhisf_tbl_par !: thickness of tbl                  [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rfrac_tbl_cav, rfrac_tbl_par !: fraction of the deepest cell affect by isf tbl  []
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf0_tbl_par
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fwfisf_cpl
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fwfisf_par, fwfisf_par_b !: net fwf from the ice shelf        [kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fwfisf_cav, fwfisf_cav_b !: net fwf from the ice shelf        [kg/m2/s]
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_cav_tsc_b, risf_cav_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s]  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_par_tsc_b, risf_par_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s]  
-
+   PUBLIC   isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl
+   !
+   !-------------------------------------------------------
+   ! 0 :              namelist parameter
+   !-------------------------------------------------------
+   !
+   ! 0.1 -------- ice shelf cavity parameter --------------
+   CHARACTER(LEN=256), PUBLIC :: cn_isfdir
+   CHARACTER(LEN=256), PUBLIC :: cn_isfload      !: ice shelf load computation method
+   !
+   ! 0.2 -------- ice shelf cavity melt namelist parameter -------------
+   LOGICAL           , PUBLIC :: ln_isfcav_mlt   !: logical for the use of ice shelf parametrisation
+   REAL(wp)          , PUBLIC :: rn_gammat0      !: temperature exchange coeficient    []
+   REAL(wp)          , PUBLIC :: rn_gammas0      !: salinity    exchange coeficient    []
+   REAL(wp)          , PUBLIC :: rn_htbl         !: Losch top boundary layer thickness [m]
+   CHARACTER(LEN=256), PUBLIC :: cn_gammablk     !: gamma formulation
+   CHARACTER(LEN=256), PUBLIC :: cn_isfcav_mlt   !: melt formulation (cavity/param)
+   TYPE(FLD_N)       , PUBLIC :: sn_isfcav_fwf   !: information about the isf melting file to be read
+   !
+   ! 0.3 -------- ice shelf param. melt namelist parameter -------------
+   LOGICAL           , PUBLIC :: ln_isfpar_mlt   !: logical for the computation of melt inside the cavity
+   CHARACTER(LEN=256), PUBLIC :: cn_isfpar_mlt   !: melt formulation (cavity/param)
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_fwf   !: information about the isf melting file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmax  !: information about the grounding line depth file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmin  !: information about the calving   line depth file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_Leff  !: information about the effective length     file to be read
+   !
+   ! 0.4 -------- coupling namelist parameter -------------
+   LOGICAL, PUBLIC :: ln_isfcpl      !:
+   LOGICAL, PUBLIC :: ln_isfcpl_cons !:
+   INTEGER, PUBLIC :: nn_drown       !:
+   !
+   !-------------------------------------------------------
+   ! 1 :              ice shelf parameter
+   !-------------------------------------------------------
+   !
    REAL(wp), PARAMETER, PUBLIC :: rLfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
    REAL(wp), PARAMETER, PUBLIC :: rcpisf = 2000.0_wp       !: specific heat of ice shelf             [J/kg/K]
@@ -68,20 +63,44 @@
    REAL(wp), PARAMETER, PUBLIC :: rhoisf = 920.0_wp        !: volumic mass of ice shelf              [kg/m3]
    REAL(wp), PARAMETER, PUBLIC :: rtsurf = -20.0           !: surface temperature                    [C]
-   REAL(wp), PARAMETER, PUBLIC :: risf_eps = 1.e-20_wp        
-
-   REAL(wp), PUBLIC            :: risf_lamb1, risf_lamb2, risf_lamb3  ! freezing point linearization coeficient
-
+   !
+   !-------------------------------------------------------
+   ! 2 :              ice shelf global variables
+   !-------------------------------------------------------
+   !
+   ! 0.1 -------- ice shelf cavity parameter --------------
+   LOGICAL , PUBLIC            :: l_isfoasis
    REAL(wp), PUBLIC            :: r1_Lfusisf               !: 1/rLfusisf
-
-   TYPE(FLD)  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf
-   TYPE(FLD)  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfpar_fwf
-
-!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
-   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
-   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
-   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
-   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   risfload                    !: ice shelf load
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   fwfisf_oasis
+   !
+   ! 0.2 -------- ice shelf cavity melt namelist parameter -------------
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: mskisf_cav                    !:
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: misfkt_cav   , misfkb_cav     !: 
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rhisf_tbl_cav, rfrac_tbl_cav  !: 
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: fwfisf_cav   , fwfisf_cav_b   !: before and now net fwf from the ice shelf        [kg/m2/s]
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     :: sf_isfcav_fwf                 !:
+   !
+   REAL(wp) , PUBLIC                                      :: risf_lamb1, risf_lamb2, risf_lamb3  ! freezing point linearization coeficient
+   !
+   ! 0.3 -------- ice shelf param. melt namelist parameter -------------
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: mskisf_par                    !:
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: misfkt_par   , misfkb_par     !:
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rhisf_tbl_par, rfrac_tbl_par  !: 
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: fwfisf_par   , fwfisf_par_b   !: before and now net fwf from the ice shelf        [kg/m2/s]
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_par_tsc , risf_par_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     :: sf_isfpar_fwf                 !:
+   !
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rhisf0_tbl_par                !: thickness of tbl (initial value)  [m]
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: risfLeff                      !:
+   !
+   ! 0.4 -------- coupling namelist parameter -------------
+   LOGICAL , PUBLIC                                        :: ll_isfcpl      !:
+   LOGICAL , PUBLIC                                        :: ll_isfcpl_cons !:
+   INTEGER , PUBLIC                                        ::   nstp_iscpl   !:
+   REAL(wp), PUBLIC                                        ::   rdt_iscpl    !: 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   risfcpl_vol, risfcpl_cons_vol  !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   risfcpl_tsc, risfcpl_cons_tsc  !:
+   !
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -121,4 +140,5 @@
       CALL mpp_sum ( 'isf', ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
+      !
    END SUBROUTINE isf_alloc_par
 
@@ -145,4 +165,59 @@
       ierr = ierr + ialloc
       !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
+      !
+   END SUBROUTINE isf_alloc_cav
+
+   SUBROUTINE isf_alloc_cpl()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_alloc_cpl  ***
+      !!
+      !! ** Purpose : 
+      !!
+      !! ** Method  : 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      ierr = 0
+      !
+      ALLOCATE( risfcpl_tsc(jpi,jpj,jpk,jpts)      , risfcpl_vol(jpi,jpj,jpk)      ,             &
+         &      risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_warn('STOP','isfcpl: failed to allocate arrays.')
+      !
+   END SUBROUTINE isf_alloc_cpl
+
+   SUBROUTINE isf_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_alloc  ***
+      !!
+      !! ** Purpose : 
+      !!
+      !! ** Method  : 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      !
+      ierr = 0       ! set to zero if no array to be allocated
+      !
+      ALLOCATE(fwfisf_par(jpi,jpj)  , fwfisf_par_b(jpi,jpj), &
+         &     fwfisf_cav(jpi,jpj)  , fwfisf_cav_b(jpi,jpj), &
+         &     fwfisf_oasis(jpi,jpj),            STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(risfload(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
       ALLOCATE( mskisf_cav(jpi,jpj), STAT=ialloc)
       ierr = ierr + ialloc
@@ -150,38 +225,6 @@
       CALL mpp_sum ( 'isf', ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
-   END SUBROUTINE isf_alloc_cav
-
-   SUBROUTINE isf_alloc()
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_alloc  ***
-      !!
-      !! ** Purpose : 
-      !!
-      !! ** Method  : 
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER :: ierr, ialloc
-      !!----------------------------------------------------------------------
-      !
-      ierr = 0       ! set to zero if no array to be allocated
-      !
-      ALLOCATE(fwfisf_par(jpi,jpj), fwfisf_par_b(jpi,jpj), &
-         &     fwfisf_cav(jpi,jpj), fwfisf_cav_b(jpi,jpj), &
-         &     fwfisf_cpl(jpi,jpj),            STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risfload(jpi,jpj), STAT=ialloc)
-      ierr = ierr + ialloc
-      !
-      CALL mpp_sum ( 'isf', ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
 
    END SUBROUTINE isf_alloc
-   !
+
 END MODULE isf
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11423)
@@ -17,10 +17,9 @@
    USE isftbl         !
    USE isfcavmlt
-   USE isfgammats
+   USE isfcavgam
    USE isfdiags
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
-   USE eosbn2         ! equation of state
-   USE zdfdrg         ! vertical physics: top/bottom drag coef.
+   USE eosbn2         ! l_useCT
    !
    USE in_out_manager ! I/O manager
@@ -166,5 +165,5 @@
          ALLOCATE( sf_isfcav_fwf(1), STAT=ierr )
          ALLOCATE( sf_isfcav_fwf(1)%fnow(jpi,jpj,1), sf_isfcav_fwf(1)%fdta(jpi,jpj,1,2) )
-         CALL fld_fill( sf_isfcav_fwf, (/ sn_isfcav_fwf /), cn_dirisf, 'isf_cav_init', 'read fresh water flux isf data', 'namisf' )
+         CALL fld_fill( sf_isfcav_fwf, (/ sn_isfcav_fwf /), cn_isfdir, 'isf_cav_init', 'read fresh water flux isf data', 'namisf' )
 
          IF(lwp) WRITE(numout,*)
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 11423)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 11423)
@@ -0,0 +1,231 @@
+MODULE isfcavgam
+   !!======================================================================
+   !!                       ***  MODULE  isfgammats  ***
+   !! Ice shelf module :  compute exchange coeficient at the ice/ocean interface
+   !!======================================================================
+   !! History :  4.1  !  (P. Mathiot) original
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isf_gammats       : compute exchange coeficient gamma 
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE isf
+   USE isftbl
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state
+   USE zdfdrg         ! vertical physics: top/bottom drag coef.
+   USE iom
+   USE in_out_manager ! I/O manager
+   !
+   IMPLICIT NONE
+   !
+   PRIVATE
+   !
+   PUBLIC   isfcav_gammats
+
+   !!----------------------------------------------------------------------
+   !! 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
+   !
+   !!-----------------------------------------------------------------------------------------------------
+   !!                              PUBLIC SUBROUTINES 
+   !!-----------------------------------------------------------------------------------------------------
+   !
+   SUBROUTINE isfcav_gammats( pttbl, pstbl, pqoce, pqfwf, pgt, pgs )
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange for heat and fwf flux
+      !!
+      !! ** Method     : select the gamma formulation
+      !!                 3 method available (cst, AD15 and HJ99)
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pgt  , pgs      ! gamma t and gamma s 
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf    ! isf heat and fwf
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl    ! top boundary layer tracer
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)                :: zutbl, zvtbl    ! top boundary layer velocity
+      !!---------------------------------------------------------------------
+      !
+      ! compute velocity in the tbl if needed
+      SELECT CASE ( cn_gammablk )
+      CASE ( 'spe'  ) 
+         ! gamma is constant (specified in namelist)
+         ! nothing to do
+      CASE ('ad15', 'hj99')
+         CALL isf_tbl(un(:,:,:) ,zutbl(:,:),'U')
+         CALL isf_tbl(vn(:,:,:) ,zvtbl(:,:),'V')
+         CALL iom_put('utbl',zutbl(:,:))
+         CALL iom_put('vtbl',zvtbl(:,:))
+      CASE DEFAULT
+         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
+      END SELECT
+      ! 
+      ! compute gamma
+      SELECT CASE ( cn_gammablk )
+      CASE ( 'spe'  ) ! gamma is constant (specified in namelist)
+         pgt(:,:) = rn_gammat0
+         pgs(:,:) = rn_gammas0
+      CASE ( 'ad15' ) ! gamma is proportional to u*
+         CALL gammats_AD15 (              zutbl, zvtbl, rCdU_top, r_ke0_top,               pgt, pgs )
+      CASE ( 'hj99' ) ! gamma depends of stability of boundary layer and u*
+         CALL gammats_HJ99 (pttbl, pstbl, zutbl, zvtbl, rCdU_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
+      CASE DEFAULT
+         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
+      END SELECT
+      !
+      ! ouput exchange coeficient and tbl velocity
+      CALL iom_put('isfgammat', pgt(:,:))
+      CALL iom_put('isfgammas', pgs(:,:))
+      !
+   END SUBROUTINE isfcav_gammats
+   !
+   !!-----------------------------------------------------------------------------------------------------
+   !!                              PRIVATE SUBROUTINES 
+   !!-----------------------------------------------------------------------------------------------------
+   !
+   SUBROUTINE gammats_AD15(putbl, pvtbl, pCd, pke2, &   ! <<== in
+      &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange coefficient 
+      !!
+      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
+      !!
+      !! ** Reference  : Jenkins et al., 2010, JPO, p2298-2312
+      !!                 Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pgt, pgs     ! gammat and gammas [m/s]
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl ! velocity in the losch top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd          ! drag coefficient
+      REAL(wp),                     INTENT(in   ) :: pke2         ! background velocity
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar
+      !!---------------------------------------------------------------------
+      !
+      ! compute ustar (AD15 eq. 27)
+      zustar(:,:) = SQRT( pCd(:,:) * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) )
+      !
+      ! Compute gammats
+      pgt(:,:) = zustar(:,:) * rn_gammat0
+      pgs(:,:) = zustar(:,:) * rn_gammas0
+      !
+      ! output ustar
+      CALL iom_put('isfustar',zustar(:,:))
+      !
+   END SUBROUTINE gammats_AD15
+
+   SUBROUTINE gammats_HJ99( pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
+      &                                                            pgt  , pgs    )  ! ==>> out gammats [m/s]
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange coefficient 
+      !!
+      !! ** Method     : gamma is velocity dependent and stability dependent
+      !!
+      !! ** Reference  : Holland and Jenkins, 1999, JPO, p1787-1800
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pgt, pgs     ! gammat and gammas
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp),                     INTENT(in   ) :: pke2           ! background velocity squared
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf   ! surface heat flux and fwf flux
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd            ! drag coeficient
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl   ! velocity in the losch top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl   ! tracer   in the losch top boundary layer
+      !!---------------------------------------------------------------------
+      INTEGER  :: ji, jj                     ! loop index
+      INTEGER  :: ikt                        ! local integer
+      REAL(wp) :: zdku, zdkv                 ! U, V shear 
+      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
+      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
+      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
+      REAL(wp) :: zhmax                      ! limitation of mol
+      REAL(wp) :: zetastar                   ! stability parameter
+      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence 
+      REAL(wp) :: zcoef                      ! temporary coef
+      REAL(wp) :: zdep
+      REAL(wp) :: zeps = 1.0e-20_wp    
+      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
+      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
+      REAL(wp), DIMENSION(2) :: zts, zab
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
+      !!---------------------------------------------------------------------
+      !
+      ! compute ustar
+      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + r_ke0_top ) )
+      !
+      ! output ustar
+      CALL iom_put('isfustar',zustar(:,:))
+      !
+      ! compute Pr and Sc number (eq ??)
+      zPr =   13.8_wp
+      zSc = 2432.0_wp
+      !
+      ! compute gamma mole (eq ??)
+      zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
+      zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
+      !
+      ! compute gamma
+      DO ji = 2, jpi
+         DO jj = 2, jpj
+            ikt = mikt(ji,jj)
+
+            IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
+               pgt = rn_gammat0
+               pgs = rn_gammas0
+            ELSE
+               ! compute Rc number (as done in zdfric.F90)
+!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
+               zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
+               !                                            ! shear of horizontal velocity
+               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
+                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
+               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
+                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+               !                                            ! richardson number (minimum value set to zero)
+               zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )
+
+               ! compute bouyancy 
+               zts(jp_tem) = pttbl(ji,jj)
+               zts(jp_sal) = pstbl(ji,jj)
+               zdep        = gdepw_n(ji,jj,ikt)
+               !
+               CALL eos_rab( zts, zdep, zab )
+               !
+               ! compute length scale (Eq ??)
+               zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
+               !
+               ! compute Monin Obukov Length
+               ! Maximum boundary layer depth (Eq ??)
+               zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
+               !
+               ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
+               zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+               !
+               ! compute eta* (stability parameter) (Eq ??)
+               zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
+               !
+               ! compute the sublayer thickness (Eq ??)
+               zhnu = 5 * znu / zustar(ji,jj)
+               !
+               ! compute gamma turb (Eq ??)
+               zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
+               &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
+               !
+               ! compute gammats
+               pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+               pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
+            END IF
+         END DO
+      END DO
+
+   END SUBROUTINE gammats_HJ99
+
+END MODULE isfcavgam
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 11423)
@@ -198,5 +198,5 @@
             ! solve the 2nd order equation to find zsfrz
             zaqe  = risf_lamb1 * (zeps1 + zeps3)
-            zaqer = 0.5_wp/MIN(zaqe,-risf_eps)
+            zaqer = 0.5_wp/MIN(zaqe,-zeps)
             zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
             zcqe  = zeps2*pstbl(ji,jj)
@@ -248,5 +248,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
       !!--------------------------------------------------------------------
-      REAL(wp)                     :: zfwf_fld, zfwf_cpl                 ! total fwf in the forcing fields (pattern) and from the cpl interface (amount)
+      REAL(wp)                     :: zfwf_fld, zfwf_oasis               ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
       REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
       REAL(wp), DIMENSION(jpi,jpj) :: zfwf                               ! 2d fwf map after scaling
@@ -266,8 +266,8 @@
       !
       ! compute glob sum from atm->oce ice shelf fwf
-      zfwf_cpl = glob_sum('isfcav_mlt', fwfisf_cpl(:,:))
+      zfwf_oasis = glob_sum('isfcav_mlt', fwfisf_oasis(:,:))
       !
       ! scale fwf
-      zfwf(:,:) = zfwf(:,:) * zfwf_cpl / zfwf_fld
+      zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
       ! 
       ! define fwf and qoce
@@ -281,4 +281,5 @@
    !SUBROUTINE isfmlt_3eq_frz_ktm1
    ! compute tfrz based on sfrz value at kt-1 (need to be SAVED local array)
+   ! => should reduce error due to linarisation
    ! compute qfwf (eq 24)
    ! compute zqoce, zqlat, zqcon, zqhc
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfgammats.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfgammats.F90	(revision 11403)
+++ 	(revision )
@@ -1,231 +1,0 @@
-MODULE isfgammats
-   !!======================================================================
-   !!                       ***  MODULE  isfgammats  ***
-   !! Ice shelf module :  compute exchange coeficient at the ice/ocean interface
-   !!======================================================================
-   !! History :  4.1  !  (P. Mathiot) original
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   isf_gammats       : compute exchange coeficient gamma 
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE isf
-   USE isftbl
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constants
-   USE eosbn2         ! equation of state
-   USE zdfdrg         ! vertical physics: top/bottom drag coef.
-   USE iom
-   USE in_out_manager ! I/O manager
-   !
-   IMPLICIT NONE
-   !
-   PRIVATE
-   !
-   PUBLIC   isfcav_gammats
-
-   !!----------------------------------------------------------------------
-   !! 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
-   !
-   !!-----------------------------------------------------------------------------------------------------
-   !!                              PUBLIC SUBROUTINES 
-   !!-----------------------------------------------------------------------------------------------------
-   !
-   SUBROUTINE isfcav_gammats( pttbl, pstbl, pqoce, pqfwf, pgt, pgs )
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : compute the coefficient echange for heat and fwf flux
-      !!
-      !! ** Method     : select the gamma formulation
-      !!                 3 method available (cst, AD15 and HJ99)
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pgt  , pgs      ! gamma t and gamma s 
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf    ! isf heat and fwf
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl    ! top boundary layer tracer
-      !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)                :: zutbl, zvtbl    ! top boundary layer velocity
-      !!---------------------------------------------------------------------
-      !
-      ! compute velocity in the tbl if needed
-      SELECT CASE ( cn_gammablk )
-      CASE ( 'spe'  ) 
-         ! gamma is constant (specified in namelist)
-         ! nothing to do
-      CASE ('ad15', 'hj99')
-         CALL isf_tbl(un(:,:,:) ,zutbl(:,:),'U')
-         CALL isf_tbl(vn(:,:,:) ,zvtbl(:,:),'V')
-         CALL iom_put('utbl',zutbl(:,:))
-         CALL iom_put('vtbl',zvtbl(:,:))
-      CASE DEFAULT
-         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
-      END SELECT
-      ! 
-      ! compute gamma
-      SELECT CASE ( cn_gammablk )
-      CASE ( 'spe'  ) ! gamma is constant (specified in namelist)
-         pgt(:,:) = rn_gammat0
-         pgs(:,:) = rn_gammas0
-      CASE ( 'ad15' ) ! gamma is proportional to u*
-         CALL gammats_AD15 (              zutbl, zvtbl, rCdU_top, r_ke0_top,               pgt, pgs )
-      CASE ( 'hj99' ) ! gamma depends of stability of boundary layer and u*
-         CALL gammats_HJ99 (pttbl, pstbl, zutbl, zvtbl, rCdU_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
-      CASE DEFAULT
-         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
-      END SELECT
-      !
-      ! ouput exchange coeficient and tbl velocity
-      CALL iom_put('isfgammat', pgt(:,:))
-      CALL iom_put('isfgammas', pgs(:,:))
-      !
-   END SUBROUTINE isfcav_gammats
-   !
-   !!-----------------------------------------------------------------------------------------------------
-   !!                              PRIVATE SUBROUTINES 
-   !!-----------------------------------------------------------------------------------------------------
-   !
-   SUBROUTINE gammats_AD15(putbl, pvtbl, pCd, pke2, &   ! <<== in
-      &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : compute the coefficient echange coefficient 
-      !!
-      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
-      !!
-      !! ** Reference  : Jenkins et al., 2010, JPO, p2298-2312
-      !!                 Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pgt, pgs     ! gammat and gammas [m/s]
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl ! velocity in the losch top boundary layer
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd          ! drag coefficient
-      REAL(wp),                     INTENT(in   ) :: pke2         ! background velocity
-      !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) :: zustar
-      !!---------------------------------------------------------------------
-      !
-      ! compute ustar (AD15 eq. 27)
-      zustar(:,:) = SQRT( pCd(:,:) * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) )
-      !
-      ! Compute gammats
-      pgt(:,:) = zustar(:,:) * rn_gammat0
-      pgs(:,:) = zustar(:,:) * rn_gammas0
-      !
-      ! output ustar
-      CALL iom_put('isfustar',zustar(:,:))
-      !
-   END SUBROUTINE gammats_AD15
-
-   SUBROUTINE gammats_HJ99( pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
-      &                                                            pgt  , pgs    )  ! ==>> out gammats [m/s]
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : compute the coefficient echange coefficient 
-      !!
-      !! ** Method     : gamma is velocity dependent and stability dependent
-      !!
-      !! ** Reference  : Holland and Jenkins, 1999, JPO, p1787-1800
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pgt, pgs     ! gammat and gammas
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp),                     INTENT(in   ) :: pke2           ! background velocity squared
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf   ! surface heat flux and fwf flux
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd            ! drag coeficient
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl   ! velocity in the losch top boundary layer
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl   ! tracer   in the losch top boundary layer
-      !!---------------------------------------------------------------------
-      INTEGER  :: ji, jj                     ! loop index
-      INTEGER  :: ikt                        ! local integer
-      REAL(wp) :: zdku, zdkv                 ! U, V shear 
-      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
-      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
-      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
-      REAL(wp) :: zhmax                      ! limitation of mol
-      REAL(wp) :: zetastar                   ! stability parameter
-      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence 
-      REAL(wp) :: zcoef                      ! temporary coef
-      REAL(wp) :: zdep
-      REAL(wp) :: zeps = 1.0e-20_wp    
-      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
-      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
-      REAL(wp), DIMENSION(2) :: zts, zab
-      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
-      !!---------------------------------------------------------------------
-      !
-      ! compute ustar
-      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + r_ke0_top ) )
-      !
-      ! output ustar
-      CALL iom_put('isfustar',zustar(:,:))
-      !
-      ! compute Pr and Sc number (eq ??)
-      zPr =   13.8_wp
-      zSc = 2432.0_wp
-      !
-      ! compute gamma mole (eq ??)
-      zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
-      zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
-      !
-      ! compute gamma
-      DO ji = 2, jpi
-         DO jj = 2, jpj
-            ikt = mikt(ji,jj)
-
-            IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
-               pgt = rn_gammat0
-               pgs = rn_gammas0
-            ELSE
-               ! compute Rc number (as done in zdfric.F90)
-!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
-               zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
-               !                                            ! shear of horizontal velocity
-               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
-                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
-               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
-                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
-               !                                            ! richardson number (minimum value set to zero)
-               zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )
-
-               ! compute bouyancy 
-               zts(jp_tem) = pttbl(ji,jj)
-               zts(jp_sal) = pstbl(ji,jj)
-               zdep        = gdepw_n(ji,jj,ikt)
-               !
-               CALL eos_rab( zts, zdep, zab )
-               !
-               ! compute length scale (Eq ??)
-               zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
-               !
-               ! compute Monin Obukov Length
-               ! Maximum boundary layer depth (Eq ??)
-               zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
-               !
-               ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
-               zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
-               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
-               !
-               ! compute eta* (stability parameter) (Eq ??)
-               zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
-               !
-               ! compute the sublayer thickness (Eq ??)
-               zhnu = 5 * znu / zustar(ji,jj)
-               !
-               ! compute gamma turb (Eq ??)
-               zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
-               &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
-               !
-               ! compute gammats
-               pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
-               pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
-            END IF
-         END DO
-      END DO
-
-   END SUBROUTINE gammats_HJ99
-
-END MODULE isfgammats
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 11423)
@@ -4,4 +4,5 @@
    USE isf
    USE phycst
+   USE in_out_manager
 
    IMPLICIT NONE
@@ -13,5 +14,5 @@
 CONTAINS
 
-   SUBROUTINE isf_hdiv( phdiv )
+   SUBROUTINE isf_hdiv( kt, phdiv )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE isf_hdiv  ***
@@ -25,4 +26,5 @@
       REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdiv   ! horizontal divergence
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
       !
       ! ice shelf cavity contribution
@@ -33,5 +35,13 @@
       !
       ! ice sheet coupling contribution (if conservation needed)
-      !IF ( ln_iscpl_hsb  ) CALL isf_hdiv_cpl(hdiv_iscpl, phdivn)
+      IF ( ll_isfcpl ) THEN
+         !
+         ! correct divergence only for the first time step
+         IF ( kt == nit000 ) CALL isf_hdiv_cpl(risfcpl_vol, phdiv)
+         !
+         ! correct divergence every time step to remove any trend due to coupling
+         IF ( ll_isfcpl_cons ) CALL isf_hdiv_cpl(risfcpl_cons_vol, phdiv)
+         !
+      END IF
       !
    END SUBROUTINE isf_hdiv
@@ -81,3 +91,15 @@
    END SUBROUTINE isf_hdiv_mlt
 
+   SUBROUTINE isf_hdiv_cpl(pqvol, phdiv)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pqvol
+   
+   INTEGER :: jk
+      
+      DO jk=1,jpk 
+         phdiv(:,:,jk) =  phdiv(:,:,jk) + pqvol(:,:,jk) * r1_e1e2t(:,:) / e3t_n(:,:,jk)
+      END DO
+
+   END SUBROUTINE
+
 END MODULE isfhdiv
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90	(revision 11423)
@@ -90,26 +90,26 @@
       !
       !                                !- Surface value + ice shelf gradient
-      risfload = 0._wp                       ! compute pressure due to ice shelf load 
+      risfload(:,:) = 0._wp                       ! compute pressure due to ice shelf load 
       DO jj = 1, jpj                         ! (used to compute hpgi/j for all the level from 1 to miku/v)
          DO ji = 1, jpi                      ! divided by 2 later
             ikt = mikt(ji,jj)
             !
-            ! top layer of the ice shelf
-            risfload(ji,jj) = risfload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1)
-            !
-            ! core layers of the ice shelf
-            DO jk = 2, ikt-1
-               risfload(ji,jj) = risfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk)
-            END DO
-            !
-            ! deepest part of the ice shelf (between deepest T point and ice/ocean interface
-            risfload(ji,jj) = risfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
-               &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) )
-            !
+            IF ( ikt > 1 ) THEN
+               !
+               ! top layer of the ice shelf
+               risfload(ji,jj) = risfload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1)
+               !
+               ! core layers of the ice shelf
+               DO jk = 2, ikt-1
+                  risfload(ji,jj) = risfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk)
+               END DO
+               !
+               ! deepest part of the ice shelf (between deepest T point and ice/ocean interface
+               risfload(ji,jj) = risfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
+                  &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) )
+               !
+            END IF
          END DO
       END DO
-      !
-      ! mask it 
-      risfload(:,:) = risfload(:,:) * mskisf_cav(:,:)  ! need to be saved for diaar5
       !
    END SUBROUTINE isf_load_isomip
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 11423)
@@ -146,5 +146,5 @@
          ALLOCATE( sf_isfpar_fwf(1), STAT=ierr )
          ALLOCATE( sf_isfpar_fwf(1)%fnow(jpi,jpj,1), sf_isfpar_fwf(1)%fdta(jpi,jpj,1,2) )
-         CALL fld_fill( sf_isfpar_fwf, (/ sn_isfpar_fwf /), cn_dirisf, 'isf_par_init', 'read fresh water flux isf data', 'namisf' )
+         CALL fld_fill( sf_isfpar_fwf, (/ sn_isfpar_fwf /), cn_isfdir, 'isf_par_init', 'read fresh water flux isf data', 'namisf' )
 
          IF(lwp) WRITE(numout,*)
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90	(revision 11423)
@@ -162,5 +162,5 @@
       !!
       !! ** Purpose    : - read ice shelf melt from forcing file => pattern
-      !!                 - total amount of fwf is given by sbccpl (fwfisf_cpl)
+      !!                 - total amount of fwf is given by sbccpl (fwfisf_oasis)
       !!                 - scale fwf and compute heat fluxes
       !!
@@ -172,5 +172,5 @@
       !!--------------------------------------------------------------------
       INTEGER                           :: jk                            ! loop index
-      REAL(wp)                          :: zfwf_fld, zfwf_cpl            ! total fwf in the forcing fields (pattern) and from the cpl interface (amount)
+      REAL(wp)                          :: zfwf_fld, zfwf_oasis            ! total fwf in the forcing fields (pattern) and from the cpl interface (amount)
       REAL(wp), DIMENSION(jpi,jpj)      :: ztfrz                         ! tbl freezing temperature
       REAL(wp), DIMENSION(jpi,jpj)      :: zfwf                          ! 2d fwf map after scaling
@@ -195,8 +195,8 @@
       !
       ! compute glob sum from atm->oce ice shelf fwf
-      zfwf_cpl = glob_sum('isfcav_mlt', fwfisf_cpl(:,:))
+      zfwf_oasis = glob_sum('isfcav_mlt', fwfisf_oasis(:,:))
       !
       ! scale fwf
-      zfwf(:,:) = zfwf(:,:) * zfwf_cpl / zfwf_fld
+      zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
       ! 
       ! define fwf and qoce
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90	(revision 11423)
@@ -20,5 +20,5 @@
    PRIVATE
 
-   PUBLIC isfrst_read, isfrst_write  ! iceshelf restart read and write 
+   PUBLIC isfrst_read, isfrst_write ! iceshelf restart read and write 
 
    !!----------------------------------------------------------------------
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 11423)
@@ -31,4 +31,5 @@
    USE isfcav         ! ice shelf cavity
    USE isfload        ! ice shelf load
+   USE isfcpl         ! isf variables
    USE isf            ! isf variables
 
@@ -87,6 +88,6 @@
          ! before time step 
          IF ( kt /= nit000 ) THEN 
-            risf_par_tsc_b (:,:) = risf_par_tsc (:,:)
-            fwfisf_par_b(:,:)    = fwfisf_par(:,:)
+            risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
+            fwfisf_par_b  (:,:)   = fwfisf_par  (:,:)
          END IF
          !
@@ -100,4 +101,8 @@
          IF (lrst_oce) CALL isfrst_write(kt, 'par', risf_par_tsc, fwfisf_par)
          !
+      END IF
+
+      IF ( ln_isfcpl ) THEN
+         IF (lrst_oce) CALL isfcpl_rst_write(kt)
       END IF
       !
@@ -122,6 +127,8 @@
       INTEGER               :: ji, jj
       !!----------------------------------------------------------------------
-      NAMELIST/namisf/ 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
+      NAMELIST/namisf/ 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
       !!----------------------------------------------------------------------
       !
@@ -130,5 +137,5 @@
       !
       riceload(:,:)       = 0.0_wp
-      fwfisf_cpl(:,:)     = 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
@@ -148,29 +155,40 @@
       IF(lwm) WRITE ( numond, namisf )
       !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'isf_init : ice shelf initialisation'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-      IF(lwp) WRITE(numout,*) '   Namelist namisf :'
-      !
-      IF(lwp) WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt
-      IF ( ln_isfcav ) THEN
-         IF(lwp) WRITE(numout,*) '         melt formulation                        cn_isfcav_mlt   = ', cn_isfcav_mlt
-         IF(lwp) WRITE(numout,*) '         thickness of the top boundary layer     rn_htbl     = ', rn_htbl
-         IF(lwp) WRITE(numout,*) '         gamma formulation                       cn_gammablk = ', cn_gammablk 
-         IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN 
-            IF(lwp) WRITE(numout,*) '            gammat coefficient                      rn_gammat0  = ', rn_gammat0  
-            IF(lwp) WRITE(numout,*) '            gammas coefficient                      rn_gammas0  = ', rn_gammas0  
-            IF(lwp) WRITE(numout,*) '            top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top 
-         END IF
-      END IF
-      !
-      IF(lwp) WRITE(numout,*) ''
-      !
-      IF(lwp) WRITE(numout,*) '      ice shelf melt parametrisation          ln_isfpar_mlt    = ', ln_isfpar_mlt
-      IF ( ln_isfpar_mlt ) THEN
-         IF(lwp) WRITE(numout,*) '         isf parametrisation formulation         cn_isfpar_mlt   = ', cn_isfpar_mlt
-      END IF
-      IF(lwp) WRITE(numout,*) ''
-      !
+      IF (lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'isf_init : ice shelf initialisation'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namisf :'
+         !
+         WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt
+         IF ( ln_isfcav ) THEN
+            WRITE(numout,*) '         melt formulation                        cn_isfcav_mlt   = ', TRIM(cn_isfcav_mlt)
+            WRITE(numout,*) '         thickness of the top boundary layer     rn_htbl     = ', rn_htbl
+            WRITE(numout,*) '         gamma formulation                       cn_gammablk = ', TRIM(cn_gammablk) 
+            IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN 
+               WRITE(numout,*) '            gammat coefficient                      rn_gammat0  = ', rn_gammat0  
+               WRITE(numout,*) '            gammas coefficient                      rn_gammas0  = ', rn_gammas0  
+               WRITE(numout,*) '            top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top 
+            END IF
+         END IF
+         WRITE(numout,*) ''
+         !
+         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)
+         END IF
+         WRITE(numout,*) ''
+         !
+         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
+         ENDIF
+         !
+         WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload)
+
+      END IF
+      !
+      !---------------------------------------------------------------------------------------------------------------------
       ! sanity check
       ! melt in the cavity without cavity
@@ -178,11 +196,14 @@
          &   CALL ctl_stop('ice shelf melt in the cavity activated (ln_isfcav_mlt) but no cavity detected in domcfg (ln_isfcav), STOP' )
       !
-      IF ( ln_cpl ) THEN
+      IF ( ln_isfcpl .AND. ln_isfcpl_cons .AND. ln_linssh ) &
+         &   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
          !
          CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
          !
          ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation 
-         IF ( ln_isfcav_mlt .AND. TRIM(cn_isfcav_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis is the only option availble with ln_cpl' )
-         IF ( ln_isfpar_mlt .AND. TRIM(cn_isfpar_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis is the only option availble with ln_cpl' )
+         IF ( ln_isfcav_mlt .AND. TRIM(cn_isfcav_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis is the only option availble if fwf send by oasis' )
+         IF ( ln_isfpar_mlt .AND. TRIM(cn_isfpar_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis is the only option availble if fwf send by oasis' )
          !
          ! oasis melt computation not tested (coded but not tested)
@@ -198,4 +219,5 @@
       END IF
       !
+      !---------------------------------------------------------------------------------------------------------------------
       ! initialisation ice shelf load
       IF ( ln_isfcav ) THEN
@@ -212,4 +234,7 @@
       r1_Lfusisf =  1._wp / rLfusisf
       !
+      ll_isfcpl     = .FALSE.
+      ll_isfcpl_cons= .FALSE.
+      !
       ! initialisation melt in the cavity
       IF ( ln_isfcav_mlt ) THEN
@@ -223,4 +248,5 @@
       END IF
       !
+      !---------------------------------------------------------------------------------------------------------------------
       ! initialisation parametrised melt
       IF ( ln_isfpar_mlt ) THEN
@@ -234,4 +260,42 @@
       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
+            !
+            ll_isfcpl = .TRUE.
+            !
+            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 )  THEN
+               ll_isfcpl_cons = .TRUE.
+               CALL isfcpl_cons()
+            END IF
+            !
+            ! Need to include in the cpl cons the isfrst_cpl_div contribution
+            ! decide how to manage thickness level change in conservation
+            !
+            tsb    (:,:,:,:) = tsn (:,:,:,:)
+            sshb   (:,:)     = sshn(:,:)
+            !
+         END IF
+      END IF
+         
   END SUBROUTINE isf_stp_init
    !!======================================================================
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 11423)
@@ -36,5 +36,5 @@
    USE eosbn2         ! 
    USE sbcrnf  , ONLY : l_rnfcpl
-   USE isf     , ONLY : l_isfcpl, fwfisf_cpl
+   USE isf     , ONLY : l_isfoasis, fwfisf_oasis
 #if defined key_cice
    USE ice_domain_size, only: ncat
@@ -475,5 +475,5 @@
 
       IF( srcv(jpr_isf)%laction .AND. ln_isf ) THEN
-         l_isfcpl             = .TRUE.                      ! -> no need to read isf in sbcisf
+         l_isfoasis = .TRUE.  ! -> isf fwf comes from oasis
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   iceshelf received from oasis '
@@ -1409,5 +1409,5 @@
          ! ice shelf fwf
          IF( srcv(jpr_isf)%laction )  THEN
-            fwfisf_cpl(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)  
+            fwfisf_oasis(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)  
          END IF
         
@@ -1713,5 +1713,5 @@
       ENDIF
       IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
-        fwfisf_cpl(:,:) = - frcv(jpr_isf)%z3(:,:,1)  
+        fwfisf_oasis(:,:) = - frcv(jpr_isf)%z3(:,:,1)  
       ENDIF
 
@@ -1752,5 +1752,5 @@
       ENDIF
       IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
-        fwfisf_cpl(:,:) = - frcv(jpr_isf)%z3(:,:,1)
+        fwfisf_oasis(:,:) = - frcv(jpr_isf)%z3(:,:,1)
       ENDIF
       !
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 11403)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traisf.F90	(revision 11423)
@@ -16,5 +16,4 @@
    USE isf            ! Ice shelf variable
    USE isfutils       !
-   USE iscplini       ! Ice sheet coupling
    !
    USE in_out_manager ! I/O manager
@@ -56,5 +55,8 @@
       !
       ! ice sheet coupling case
-      !CALL tra_isf_cpl()
+      IF ( ll_isfcpl .AND. kt == nit000 ) CALL tra_isf_cpl(risfcpl_tsc, tsa)
+      !
+      ! ice sheet coupling case (
+      IF ( ll_isfcpl_cons ) CALL tra_isf_cpl(risfcpl_cons_tsc, tsa)
       !
       IF( ln_timing )   CALL timing_stop('tra_isf')
@@ -76,5 +78,4 @@
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: phtbl, pfrac
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc , ptsc_b
-      REAL(wp), DIMENSION(jpi,jpj) :: ztc
       !!----------------------------------------------------------------------
       INTEGER                      :: ji,jj,jk  ! loop index   
@@ -107,23 +108,24 @@
    END SUBROUTINE tra_isf_mlt
    !
-   !SUBROUTINE tra_isf_cpl
-       !
-!      !----------------------------------------
-!      !        Ice Sheet coupling imbalance correction to have conservation
-!      !----------------------------------------
-!      !
-!      IF( ln_iscpl .AND. ln_iscpl_hsb) THEN         ! input of heat and salt due to river runoff 
-!         DO jk = 1,jpk
-!            DO jj = 2, jpj 
-!               DO ji = fs_2, fs_jpim1
-!                  zdep = 1._wp / e3t_n(ji,jj,jk) 
-!                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
-!                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
-!               END DO  
-!            END DO  
-!         END DO
-!      ENDIF
-!      !
-!   END SUBROUTINE tra_isf_cpl
+   SUBROUTINE tra_isf_cpl( ptsc, ptsa )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_isf_cpl  ***
+      !!
+      !! *** Action :: Update tsa with the ice shelf coupling trend 
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: ptsa
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) :: ptsc  !( >0 out )
+      !!----------------------------------------------------------------------
+      INTEGER :: jk
+      !!----------------------------------------------------------------------
+      !
+      DO jk = 1,jpk
+         ptsa(:,:,jk,jp_tem) = ptsa(:,:,jk,jp_tem) - ptsc(:,:,jk,jp_tem) * r1_e1e2t(:,:) / e3t_n(:,:,jk)
+         ptsa(:,:,jk,jp_sal) = ptsa(:,:,jk,jp_sal) - ptsc(:,:,jk,jp_sal) * r1_e1e2t(:,:) / e3t_n(:,:,jk)
+      END DO
+      !
+   END SUBROUTINE tra_isf_cpl
    !
 END MODULE traisf
