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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90	(revision 11541)
@@ -64,4 +64,5 @@
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
       REAL(wp) ::   zraur, zdep   ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) :: ztmp
       !!----------------------------------------------------------------------
       !
@@ -85,4 +86,11 @@
          END DO  
       END DO
+
+      ztmp = 0.0
+      DO jk = jpkm1,1,-1
+         ztmp(:,:) = ztmp(:,:) + hdivn(:,:,jk) * e3t_n(:,:,jk)
+      END DO
+      CALL debug('divhor 1',ztmp)
+
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN
@@ -105,4 +113,10 @@
       CALL lbc_lnk( 'divhor', hdivn, 'T', 1. )   !   (no sign change)
       !
+      ztmp = 0.0
+      DO jk = jpkm1,1,-1
+         ztmp(:,:) = ztmp(:,:) + hdivn(:,:,jk) * e3t_n(:,:,jk)
+      END DO
+      CALL debug('divhor 2',ztmp)
+      !
       IF( ln_timing )   CALL timing_stop('div_hor')
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 11541)
@@ -244,5 +244,5 @@
             !
             ! ice shelf melting
-            IF ( ln_isf ) CALL isf_dynnxt( zcoef )
+            IF ( ln_isf ) CALL isf_dynnxt( kt, atfp * rdt )
             !
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90	(revision 11541)
@@ -382,5 +382,5 @@
 !!gm  Is it correct to do so ?   I think so...
       
-      
+      !
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
@@ -631,10 +631,11 @@
       !                                   !* Right-Hand-Side of the barotropic ssh equation
       !                                   ! -----------------------------------------------
-      !                                         ! Surface net water flux and rivers
+      !                                         ! Surface net water flux, rivers and ice shelves
       IF (ln_bt_fw) THEN
          zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
       ELSE
          zztmp = r1_rau0 * r1_2
-         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
+         zssh_frc(:,:) = zztmp * (  emp(:,:)        + emp_b(:,:)                    &
+                &                 - rnf(:,:)        - rnf_b(:,:)                    &
                 &                 + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)             &
                 &                 + fwfisf_par(:,:) + fwfisf_par_b(:,:)             )
@@ -645,10 +646,17 @@
       ENDIF
       !
-      IF( ll_isfcpl_cons ) THEN
-         zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:)
-      END IF
-      !
-      IF( ln_isfcpl .AND. ln_rstart .AND. kt == nit000 ) THEN
-         zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:)
+      !                                         ! ice sheet coupling
+      IF ( ln_isf .AND. ln_isfcpl ) THEN
+         !
+         ! ice sheet coupling
+         IF( ln_rstart .AND. kt == nit000 ) THEN
+            zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:)
+         END IF
+         !
+         ! conservation option
+         IF( ln_isfcpl_cons ) THEN
+            zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:)
+         END IF
+         !
       END IF
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90	(revision 11541)
@@ -261,4 +261,8 @@
                &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
          ENDIF
+
+         ! ice sheet coupling
+         IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) sshb(:,:) = sshb(:,:) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
+
          sshn(:,:) = ssha(:,:)                              ! now <-- after
       ENDIF
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11541)
@@ -71,5 +71,4 @@
    ! 0.1 -------- ice shelf cavity parameter --------------
    LOGICAL , PUBLIC            :: l_isfoasis
-   REAL(wp), PUBLIC            :: r1_Lfusisf               !: 1/rLfusisf
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   risfload                    !: ice shelf load
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   fwfisf_oasis
@@ -97,6 +96,4 @@
    !
    ! 0.4 -------- coupling namelist parameter -------------
-   LOGICAL , PUBLIC                                        :: ll_isfcpl      !:
-   LOGICAL , PUBLIC                                        :: ll_isfcpl_cons !:
    INTEGER , PUBLIC                                        ::   nstp_iscpl   !:
    REAL(wp), PUBLIC                                        ::   rdt_iscpl    !: 
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11541)
@@ -1,22 +1,22 @@
 MODULE isfcav
    !!======================================================================
-   !!                       ***  MODULE  sbcisf  ***
-   !! Surface module :  update surface ocean boundary condition under ice
-   !!                   shelf
+   !!                       ***  MODULE  isfcav  ***
+   !! Ice shelf cavity module :  update ice shelf melting under ice
+   !!                            shelf
    !!======================================================================
    !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
-   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
    !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
+   !!            4.1  !  2019-09  (P. Mathiot) Split ice shelf cavity and ice shelf parametrisation
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   sbc_isf       : update sbc under ice shelf
+   !!   isf_cav       : update ice shelf melting under ice shelf
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
-   USE isf            !
-   USE isftbl         !
-   USE isfcavmlt
-   USE isfcavgam
-   USE isfdiags
+   USE isf            ! ice shelf public variables
+   USE isftbl         ! ice shelf top boundary layer properties
+   USE isfcavmlt      ! ice shelf melt formulation
+   USE isfcavgam      ! ice shelf melt exchange coeficient
+   USE isfdiags       ! ice shelf diags
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -26,5 +26,5 @@
    USE iom            ! I/O library
    USE fldread        ! read input field at current time step
-   USE lbclnk         !
+   USE lbclnk         ! lbclnk
 
    IMPLICIT NONE
@@ -43,14 +43,16 @@
    SUBROUTINE isf_cav( kt, ptsc, pqfwf )
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE sbc_isf_cav  ***
+      !!                     ***  ROUTINE isf_cav  ***
       !!
       !! ** Purpose :   handle surface boundary condition under ice shelf
       !!
-      !! ** Method  : -
+      !! ** Method  :   based on Mathiot et al. (2017)
       !!
-      !! ** Action  :   utau, vtau : remain unchanged
-      !!                taum, wndm : remain unchanged
-      !!                qns        : update heat flux below ice shelf
-      !!                emp, emps  : update freshwater flux below ice shelf
+      !! ** Action  :   - compute geometry of the Losch top bournary layer (see Losch et al. 2008)
+      !!                - depending on the chooses option
+      !!                   - compute temperature/salt in the tbl
+      !!                   - compute exchange coeficient
+      !!                   - compute heat and fwf fluxes
+      !!                   - output
       !!---------------------------------------------------------------------
       !!-------------------------- OUT --------------------------------------
@@ -68,8 +70,4 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl
       !!---------------------------------------------------------------------
-      !
-      ! compute misfkb_par, rhisf_tbl
-      rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
-      CALL isf_tbl_lvl( ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
       !
       ! compute T/S/U/V for the top boundary layer
@@ -118,5 +116,5 @@
             END IF
          END SELECT
-
+         !
       END DO
       !
@@ -145,7 +143,7 @@
    SUBROUTINE isf_cav_init
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_diags_2dto3d ***
+      !!                  ***  ROUTINE isf_cav_init ***
       !!
-      !! ** Purpose : 
+      !! ** Purpose : initialisation of variable needed to compute melt under an ice shelf
       !!
       !!----------------------------------------------------------------------
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 11541)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE  isfgammats  ***
-   !! Ice shelf module :  compute exchange coeficient at the ice/ocean interface
+   !! Ice shelf gamma module :  compute exchange coeficient at the ice/ocean interface
    !!======================================================================
    !! History :  4.1  !  (P. Mathiot) original
@@ -8,5 +8,5 @@
 
    !!----------------------------------------------------------------------
-   !!   isf_gammats       : compute exchange coeficient gamma 
+   !!   isfcav_gammats       : compute exchange coeficient gamma 
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 11541)
@@ -1,5 +1,5 @@
 MODULE isfcavmlt
    !!======================================================================
-   !!                       ***  MODULE  isfcav_mlt  ***
+   !!                       ***  MODULE  isfcavmlt  ***
    !! ice shelf module :  update surface ocean boundary condition under ice
    !!                   shelf
@@ -9,9 +9,8 @@
 
    !!----------------------------------------------------------------------
-   !!   isfcav_mlt    : 
+   !!   isfcav_mlt    : update surface ocean boundary condition under ice shelf
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
    USE isf            ! ice shelf public variables
-   USE isfutils
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -42,4 +41,6 @@
       &                           pqhc, pqoce, pqfwf  )
       !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt  ***
       !!
       !! ** Purpose    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
@@ -79,7 +80,10 @@
 ! -------------------------------------------------------------------------------------------------------
 
-   SUBROUTINE isfcav_mlt_spe(kt, pstbl,              &  ! <<== in
+   SUBROUTINE isfcav_mlt_spe(kt, pstbl,          &  ! <<== in
       &                      pqhc , pqoce, pqfwf )  ! ==>> out
       !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt_spe  ***
+      !!
       !! ** Purpose    : - read ice shelf melt from forcing file
       !!                 - compute ocea-ice heat flux (assuming it is equal to latent heat)
@@ -109,7 +113,10 @@
    END SUBROUTINE isfcav_mlt_spe
 
-   SUBROUTINE isfcav_mlt_2eq(pgt, pttbl, pstbl,       &  ! <<== in
-      &                       pqhc , pqoce, pqfwf )  ! ==>> out
-      !!----------------------------------------------------------------------
+   SUBROUTINE isfcav_mlt_2eq(pgt , pttbl, pstbl, &  ! <<== in
+      &                      pqhc, pqoce, pqfwf  )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt_spe  ***
+      !!
       !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
       !!
@@ -141,5 +148,5 @@
       !
       ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
-      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * zthd(:,:) * r1_Lfusisf  ! fresh water flux  ( > 0 out )
+      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * zthd(:,:) / rLfusisf    ! fresh water flux  ( > 0 out )
       pqoce(:,:) = - pqfwf(:,:) * rLfusisf                         ! ocea-ice flux     ( > 0 out )
       pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                 ! heat content flux ( > 0 out )
@@ -150,7 +157,10 @@
    END SUBROUTINE isfcav_mlt_2eq
 
-   SUBROUTINE isfcav_mlt_3eq(pgt, pgs , pttbl, pstbl, &
-      &                       pqhc, pqoce, pqfwf  )
-      !!----------------------------------------------------------------------
+   SUBROUTINE isfcav_mlt_3eq(pgt, pgs , pttbl, pstbl, &  ! <<== in
+      &                           pqhc, pqoce, pqfwf  )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt_3eq  ***
+      !!
       !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using the 3 equation formulation 
       !!
@@ -232,7 +242,8 @@
    END SUBROUTINE isfcav_mlt_3eq
 
-   SUBROUTINE isfcav_mlt_oasis(kt, pstbl,              &  ! <<== in
+   SUBROUTINE isfcav_mlt_oasis(kt, pstbl,          &  ! <<== in
       &                        pqhc , pqoce, pqfwf )  ! ==>> out
       !!----------------------------------------------------------------------
+      !!                          ***  ROUTINE isfcav_mlt_oasis  ***
       !!
       !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
@@ -282,11 +293,3 @@
    END SUBROUTINE isfcav_mlt_oasis
 
-   !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
-   ! compute sfrz (eq 26)
-   !END SUBROUTINE isfmlt_3eq_frz_ktm1
-
 END MODULE isfcavmlt
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 11541)
@@ -20,9 +20,8 @@
       !!                  ***  SUBROUTINE isf_hdiv  ***
       !!       
-      !! ** Purpose :   
+      !! ** Purpose :   update the horizontal divergence with the ice shelf contribution
+      !!                (parametrisation, explicit, ice sheet coupling conservation
+      !!                 increment)
       !!
-      !! ** Method  :   
-      !!
-      !! ** Action  :   phdiv   decreased by the fwf inflow (isf melt in this case)
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdiv   ! horizontal divergence
@@ -38,14 +37,17 @@
          IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, phdiv)
          !
-      END IF
-      !
-      ! ice sheet coupling contribution (if conservation needed)
-      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)
+         ! ice sheet coupling contribution 
+         IF ( ln_isfcpl ) THEN
+            !
+            ! correct divergence only for the first time step
+            IF ( kt == nit000   ) CALL isf_hdiv_cpl(risfcpl_vol       , phdiv)
+            IF ( kt == nit000+1 ) CALL isf_hdiv_cpl(risfcpl_vol*0.5_wp, phdiv)
+            !
+            ! correct divergence every time step to remove any trend due to coupling
+            ! conservation option
+            IF ( ln_isfcpl_cons ) CALL isf_hdiv_cpl(risfcpl_cons_vol, phdiv)
+            !
+         END IF
          !
       END IF
@@ -57,11 +59,10 @@
       !!                  ***  SUBROUTINE sbc_isf_div  ***
       !!       
-      !! ** Purpose :   update the horizontal divergence with the runoff inflow
+      !! ** Purpose :   update the horizontal divergence with the ice shelf inflow
       !!
-      !! ** Method  :   
-      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the 
-      !!                          divergence and expressed in m/s
+      !! ** Method  :   pfwf is positive (outflow) and expressed as kg/m2/s
+      !!                increase the divergence
       !!
-      !! ** Action  :   phdivn   decreased by the runoff inflow
+      !! ** Action  :   phdivn   increased by the ice shelf outflow
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
@@ -73,5 +74,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt, ikb 
-      REAL(wp), DIMENSION(jpi,jpj) :: zqvol,ztmp
+      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
       !!----------------------------------------------------------------------
       !
@@ -79,5 +80,5 @@
       !
       ! compute integrated divergence correction
-      zqvol(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rau0 / phtbl(:,:)
+      zhdiv(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rau0 / phtbl(:,:)
       !
       ! update divergence at each level affected by ice shelf top boundary layer
@@ -88,8 +89,8 @@
             ! level fully include in the ice shelf boundary layer
             DO jk = ikt, ikb - 1
-               phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zqvol(ji,jj)
+               phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zhdiv(ji,jj)
             END DO
             ! level partially include in ice shelf boundary layer 
-            phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zqvol(ji,jj) * pfrac(ji,jj)
+            phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zhdiv(ji,jj) * pfrac(ji,jj)
          END DO
       END DO
@@ -101,9 +102,12 @@
       !!                  ***  SUBROUTINE isf_hdiv_cpl  ***
       !!       
-      !! ** Purpose :   
+      !! ** Purpose :   update the horizontal divergence with the ice shelf 
+      !!                coupling conservation increment
       !!
-      !! ** Method  :   
+      !! ** Method  :   pqvol is positive (outflow) and expressed as m3/s
+      !!                increase the divergence
       !!
-      !! ** Action  :  
+      !! ** Action  :   phdivn   increased by the ice shelf outflow
+      !!
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90	(revision 11541)
@@ -80,5 +80,5 @@
       !
       !                                !- assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
-      zts_top(:,:,jp_tem) = -1.9_wp   ;   zts_top(:,:,jp_sal) = 34.4_wp
+      zts_top(:,:,jp_tem) = -1.9_wp   ;   zts_top(:,:,jp_sal) = 35.4_wp
       !
       DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfnxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfnxt.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfnxt.F90	(revision 11541)
@@ -8,8 +8,9 @@
   
    !!-------------------------------------------------------------------------
-   !!   dyn_nxt       : obtain the next (after) horizontal velocity
+   !!   isfnxt       : aplly correction need for the ice shelf
    !!-------------------------------------------------------------------------
 
    USE isf
+   USE phycst
    USE dom_oce
    USE in_out_manager
@@ -23,5 +24,5 @@
 CONTAINS
 
-   SUBROUTINE isf_dynnxt ( pcoef )
+   SUBROUTINE isf_dynnxt ( kt, pcoef )
       !!--------------------------------------------------------------------
       !!                  ***  ROUTINE isf_dynnxt  ***
@@ -31,7 +32,10 @@
       !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
-      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
+      INTEGER ,                     INTENT(in   ) :: kt
+      !
+      REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt
       !!-------------------------- IN  -------------------------------------
       !!--------------------------------------------------------------------
+      INTEGER :: jk  ! loop index
       !!--------------------------------------------------------------------
       !
@@ -41,4 +45,10 @@
       ! ice shelf parametrised
       IF ( ln_isfpar_mlt ) CALL isf_dynnxt_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, pcoef)
+      !
+      IF ( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
+         DO jk = 1, jpkm1
+            e3t_b(:,:,jk) =   e3t_b(:,:,jk) - pcoef * risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
+         END DO
+      END IF
       !
    END SUBROUTINE isf_dynnxt
@@ -63,5 +73,5 @@
       !
       ! compute fwf conservation correction
-      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:)
+      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:) * r1_rau0
       !
       ! add the increment in the tbl
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 11541)
@@ -8,16 +8,16 @@
    !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
    !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
-   !!            4.0  !  2019-09  (P. Mathiot) Restructuration
+   !!            4.1  !  2019-09  (P. Mathiot) Restructuration
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   isf_par       : compute ice shelf melt using a prametrisation of ice shelf cavities
+   !!   isfpar       : compute ice shelf melt using a prametrisation of ice shelf cavities
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
-   USE isf
-   USE isfparmlt
-   USE isftbl
-   USE isfdiags
-   USE isfutils
+   USE isf            ! ice shelf
+   USE isfutils       !
+   USE isfparmlt      ! ice shelf parametrisation
+   USE isftbl         ! ice shelf depth average
+   USE isfdiags       ! ice shelf diagnostics
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -45,10 +45,6 @@
       !!                     ***  ROUTINE isf_par ***      
       !!
-      !! ** Purpose :   
+      !! ** Purpose : compute the heat and fresh water due to ice shelf melting/freezing using a parametrisation 
       !!
-      !! ** Method  : 
-      !!
-      !! ** Action  :   
-      !!                
       !!---------------------------------------------------------------------
       !!-------------------------- OUT --------------------------------------
@@ -60,8 +56,4 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh
       !!---------------------------------------------------------------------
-      !
-      ! compute misfkb_par, rhisf_tbl
-      rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
-      CALL isf_tbl_lvl( ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
       !
       ! compute heat content, latent heat and melt fluxes (2d)
@@ -94,7 +86,6 @@
       !!                  ***  ROUTINE isf_par_init  ***
       !!
-      !! ** Purpose : 
+      !! ** Purpose : initialisation of the variable needed for the parametrisation of ice shelf melt
       !!
-      !! ** Method  :  
       !!----------------------------------------------------------------------
       INTEGER               :: ierr
@@ -125,5 +116,5 @@
       !
       ! compute ktop
-      CALL isftbl_ktop(ztblmin, misfkt_par)
+      CALL isf_tbl_ktop(ztblmin, misfkt_par)
       !
       ! initial tbl thickness
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90	(revision 11541)
@@ -2,6 +2,6 @@
    !!======================================================================
    !!                       ***  MODULE  isfparmlt  ***
-   !! Surface module :  update surface ocean boundary condition under ice
-   !!                   shelf
+   !! Ice shelf parametrisation module :  update surface ocean boundary condition under ice
+   !!                   shelf using an ice shelf melt parametrisation
    !!======================================================================
    !! History :  4.0  ! original code
@@ -9,6 +9,6 @@
 
    USE oce            ! ocean dynamics and tracers
-   USE isf
-   USE isftbl
+   USE isf            ! ice shelf
+   USE isftbl         ! ice shelf depth average
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -17,6 +17,6 @@
    USE in_out_manager ! I/O manager
    USE iom            ! I/O library
-   USE fldread
-   USE lib_fortran
+   USE fldread        !
+   USE lib_fortran    !
 
    IMPLICIT NONE
@@ -143,5 +143,5 @@
       ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
       pqoce(:,:) =   rau0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:)
-      pqfwf(:,:) = - pqoce(:,:) * r1_Lfusisf           ! derived from the latent heat flux
+      pqfwf(:,:) = - pqoce(:,:) / rLfusisf             ! derived from the latent heat flux
       pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
       !
@@ -158,5 +158,5 @@
    SUBROUTINE isfpar_mlt_oasis(kt, pqhc , pqoce, pqfwf )
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE isfpar_oasis  ***
+      !!                  ***  ROUTINE isfpar_mlt_oasis  ***
       !!
       !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
@@ -173,5 +173,5 @@
       !!--------------------------------------------------------------------
       INTEGER                           :: jk                            ! loop index
-      REAL(wp)                          :: zfwf_fld, zfwf_oasis            ! 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
@@ -182,6 +182,5 @@
       CALL fld_read ( kt, nn_fsbc, sf_isfpar_fwf   )
       !
-      ! compute ptfrz
-      ! 1. ------------Mean freezing point
+      ! 1. ------------Mean freezing point (needed for heat content flux)
       DO jk = 1,jpk
          CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
@@ -189,4 +188,5 @@
       CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
       !
+      ! 2. ------------Scale isf melt pattern with total amount from oasis
       ! ice shelf 2d map
       zfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)
@@ -203,5 +203,5 @@
       zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
       ! 
-      ! define fwf and qoce
+      ! i3. -----------Define fwf and qoce
       ! ocean heat flux is assume to be equal to the latent heat
       pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out )
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 11541)
@@ -75,4 +75,5 @@
          !
          ! compute tbl lvl/h
+         rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
          CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
          !
@@ -94,4 +95,5 @@
          !
          ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
          CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
          !
@@ -140,4 +142,5 @@
       CALL isf_alloc()
       !
+      ! initalisation of fwf and tsc array to 0
       riceload(:,:)       = 0.0_wp
       fwfisf_oasis(:,:)   = 0.0_wp
@@ -166,5 +169,5 @@
          IF ( ln_isf ) THEN
             WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt
-            IF ( ln_isfcav ) THEN
+            IF ( ln_isfcav_mlt) THEN
                WRITE(numout,*) '         melt formulation                        cn_isfcav_mlt   = ', TRIM(cn_isfcav_mlt)
                WRITE(numout,*) '         thickness of the top boundary layer     rn_htbl     = ', rn_htbl
@@ -184,5 +187,13 @@
             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,*) ''
+            !
          ELSE
+
             IF ( ln_isfcav ) THEN
                WRITE(numout,*) ''
@@ -192,10 +203,4 @@
          END IF
 
-         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
-         !
          IF (ln_isfcav) WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload)
          WRITE(numout,*) ''
@@ -251,7 +256,5 @@
       IF ( ln_isf ) THEN
          !
-         ! initialisation useful variable
-         r1_Lfusisf =  1._wp / rLfusisf
-         !
+         !---------------------------------------------------------------------------------------------------------------------
          ! initialisation melt in the cavity
          IF ( ln_isfcav_mlt ) THEN
@@ -276,48 +279,41 @@
             !
          END IF
-      END IF
-      !
-      !---------------------------------------------------------------------------------------------------------------------
-      ! initialisation ice sheet coupling
-      !
-      ll_isfcpl     = .FALSE.
-      ll_isfcpl_cons= .FALSE.
-      !
-      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('wmask')
-         !IF( lwxios ) CALL iom_set_rstw_var_active('gdepw_n')
-         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()
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation ice sheet coupling
+         IF( ln_isfcpl ) THEN
+            ! 
+            ! prepare writing restart
+            IF( lwxios ) CALL iom_set_rstw_var_active('ssmask')
+            IF( lwxios ) CALL iom_set_rstw_var_active('tmask')
+            IF( lwxios ) CALL iom_set_rstw_var_active('e3t_n')
+            IF( lwxios ) CALL iom_set_rstw_var_active('e3u_n')
+            IF( lwxios ) CALL iom_set_rstw_var_active('e3v_n')
+            !
+            IF( ln_rstart ) THEN
+               !
+               CALL isf_alloc_cpl()
+               !
+               ! extrapolation tracer properties
+               CALL isfcpl_tra()
+               !
+               ! correction of the horizontal divergence and associated temp. and salt content flux
+               CALL isfcpl_vol()
+               !
+               ! apply the 'conservation' method
+               IF ( ln_isfcpl_cons ) CALL isfcpl_cons()
+               !
+               ! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
+               !           - decide how to manage thickness level change in conservation
+               !
+               tsb    (:,:,:,:) = tsn (:,:,:,:)
+               sshb   (:,:)     = sshn(:,:)
+               ub     (:,:,:)   = un(:,:,:)
+               vb     (:,:,:)   = vn(:,:,:)
+               !
             END IF
             !
-            ! 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 IF
          
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 11541)
@@ -21,5 +21,5 @@
    PRIVATE
 
-   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isftbl_ktop, isftbl_kbot
+   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
 
 CONTAINS
@@ -30,4 +30,9 @@
       !!
       !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
+      !!
+      !! ** Method : Average properties over a specific thickness
+      !!
+      !! ** Reference : inspired from : Losch, Modeling ice shelf cavities in a z coordinate ocean general circulation model
+      !!                https://doi.org/10.1029/2007JC004368 , 2008
       !!
       !!--------------------------------------------------------------------
@@ -103,5 +108,5 @@
    SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
       !!--------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_tbl_lvl  ***
+      !!                  ***  ROUTINE isf_tbl_avg  ***
       !!
       !! ** Purpose : compute mean property in the boundary layer
@@ -145,6 +150,7 @@
       !!                  ***  ROUTINE isf_tbl_lvl  ***
       !!
-      !! ** Purpose : - compute bottom level fully included in the top boundary layer
+      !! ** Purpose : - compute bottom level off the top boundary layer
       !!              - thickness of the top boundary layer
+      !!              - fraction of the bottom level affected by the tbl
       !!
       !!---------------------------------------------------------------------
@@ -180,5 +186,5 @@
       !
       ! get ktbl
-      CALL isftbl_kbot(ktop, phtbl, pe3, kbot)
+      CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
       !
       ! get pfrac
@@ -197,7 +203,7 @@
    END SUBROUTINE isf_tbl_lvl
    !
-   SUBROUTINE isftbl_kbot(ktop, phtbl, pe3, kbot)
-      !!--------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_tbl_lvl  ***
+   SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_bot  ***
       !!
       !! ** Purpose : compute bottom level of the isf top boundary layer
@@ -232,9 +238,9 @@
       END DO
       !
-   END SUBROUTINE isftbl_kbot
-      !
-   SUBROUTINE isftbl_ktop(pdep, ktop)
-      !!--------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_tbl_lvl  ***
+   END SUBROUTINE isf_tbl_kbot
+      !
+   SUBROUTINE isf_tbl_ktop(pdep, ktop)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_top  ***
       !!
       !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
@@ -263,5 +269,5 @@
       END DO
       !
-   END SUBROUTINE isftbl_ktop
+   END SUBROUTINE isf_tbl_ktop
 
 END MODULE isftbl
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 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traisf.F90	(revision 11541)
@@ -48,21 +48,18 @@
       IF( ln_timing )   CALL timing_start('tra_isf')
       !
-      IF (ln_isf) THEN
-         !
-         ! cavity case
-         IF ( ln_isfcav_mlt ) CALL tra_isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, risf_cav_tsc_b, tsa)
-         !
-         ! parametrisation case
-         IF ( ln_isfpar_mlt ) CALL tra_isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, risf_par_tsc_b, tsa)
-         !
-      END IF
+      ! cavity case
+      IF ( ln_isfcav_mlt ) CALL tra_isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, risf_cav_tsc_b, tsa)
+      !
+      ! parametrisation case
+      IF ( ln_isfpar_mlt ) CALL tra_isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, risf_par_tsc_b, tsa)
       !
       ! ice sheet coupling case
-      IF ( ll_isfcpl ) THEN
+      IF ( ln_isfcpl ) THEN
          !
-         IF ( kt == nit000 ) CALL tra_isf_cpl(risfcpl_tsc, tsa)
+         IF ( kt == nit000  ) CALL tra_isf_cpl(risfcpl_tsc       , tsa)
+         IF ( kt == nit000+1) CALL tra_isf_cpl(risfcpl_tsc*0.5_wp, tsa)
          !
          ! ensure 0 trend due to unconservation of the ice shelf coupling
-         IF ( ll_isfcpl_cons ) CALL tra_isf_cpl(risfcpl_cons_tsc, tsa)
+         IF ( ln_isfcpl_cons ) CALL tra_isf_cpl(risfcpl_cons_tsc, tsa)
          !
       END IF
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 11541)
@@ -340,4 +340,5 @@
                   ! ice shelf
                   IF( ll_isf ) THEN
+                     !
                      IF ( ln_isfcav_mlt ) THEN
                         ! level fully include in the Losch_2008 ice shelf boundary layer
@@ -372,4 +373,10 @@
                         END IF
                      END IF
+                     !
+                     IF (ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN    ! risfcpl_vol_n = 0 and risfcpl_vol_b = risfcpl_vol
+                        ztc_f  = ztc_f  + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj)
+                        ze3t_f = ze3t_f - zfact1 * risfcpl_vol(ji,jj,jk   ) * r1_e1e2t(ji,jj)
+                     END IF
+                     !
                   END IF
                   !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 11529)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 11541)
@@ -114,5 +114,5 @@
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
       IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-      IF( ln_isf     .OR. ln_isfcpl )   CALL isf_stp ( kstp )
+      IF( ln_isf     )   CALL isf_stp ( kstp )
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
 
@@ -242,5 +242,5 @@
                          CALL tra_sbc       ( kstp )  ! surface boundary condition
       IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
-      IF( ln_isf .OR. ln_isfcpl ) CALL tra_isf ( kstp )  ! ice shelf heat flux
+      IF( ln_isf     ) CALL tra_isf ( kstp )  ! ice shelf heat flux
       IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
       IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme
