Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/dom_oce.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/dom_oce.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/dom_oce.F90	(revision 11403)
@@ -157,4 +157,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m)
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep, bathy
 
    !!----------------------------------------------------------------------
@@ -271,6 +272,7 @@
          &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  )
          !
-         !
-      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) )
+      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  ) 
+         !
+      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) )
          !
       ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domzgr.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domzgr.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domzgr.F90	(revision 11403)
@@ -71,5 +71,6 @@
       INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices
       !
-      INTEGER  ::   jk                  ! dummy loop index
+      INTEGER  ::   ji,jj,jk            ! dummy loop index
+      INTEGER  ::   ikt, ikb            ! top/bot index
       INTEGER  ::   ioptio, ibat, ios   ! local integer
       REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
@@ -138,6 +139,15 @@
       !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
       CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1
-      
-
+      !
+      !                                ! ice shelf draft and bathymetry
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ikt = mikt(ji,jj)
+            ikb = mbkt(ji,jj)
+            bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1)
+            risfdep(ji,jj) = gdepw_0(ji,jj,ikt  )
+         END DO
+      END DO
+      !
       !                                ! deepest/shallowest W level Above/Below ~10m
 !!gm BUG in s-coordinate this does not work!
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynhpg.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynhpg.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynhpg.F90	(revision 11403)
@@ -31,5 +31,6 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
-   USE isf             ! ice shelf
+   USE isf             ! ice shelf  (risfload variable)
+   USE isfload         ! ice shelf  (isf_load routine )
    USE sbc_oce         ! surface variable (only for the flag with ice shelf)
    USE dom_oce         ! ocean space and time domain
@@ -214,46 +215,5 @@
       ENDIF
       !                          
-      IF ( .NOT. ln_isfcav ) THEN     !--- no ice shelf load
-         riceload(:,:) = 0._wp
-         !
-      ELSE                            !--- set an ice shelf load
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ice shelf case: set the ice-shelf load'
-         ALLOCATE( zts_top(jpi,jpj,jpts) , zrhd(jpi,jpj,jpk) , zrhdtop_isf(jpi,jpj) , ziceload(jpi,jpj) ) 
-         !
-         znad = 1._wp                     !- To use density and not density anomaly
-         !
-         !                                !- 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
-         !
-         DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
-            CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) )
-         END DO
-         !
-         !                                !- compute rhd at the ice/oce interface (ice shelf side)
-         CALL eos( zts_top , risfdep, zrhdtop_isf )
-         !
-         !                                !- Surface value + ice shelf gradient
-         ziceload = 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)
-               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))
-               DO jk = 2, ikt-1
-                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &
-                     &                              * (1._wp - tmask(ji,jj,jk))
-               END DO
-               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
-                  &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) )
-            END DO
-         END DO
-         riceload(:,:) = ziceload(:,:)  ! need to be saved for diaar5
-         !
-         DEALLOCATE( zts_top , zrhd , zrhdtop_isf , ziceload ) 
-      ENDIF
-      !
    END SUBROUTINE dyn_hpg_init
-
 
    SUBROUTINE hpg_zco( kt )
@@ -572,5 +532,5 @@
       !!         ua = ua - 1/e1u * zhpi
       !!         va = va - 1/e2v * zhpj
-      !!      iceload is added and partial cell case are added to the top and bottom
+      !!      iceload is added
       !!      
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
@@ -617,10 +577,10 @@
                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &
                &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            ) 
+               &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            ) 
             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    &
                &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
                &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            ) 
+               &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11403)
@@ -11,5 +11,5 @@
 
    !!----------------------------------------------------------------------
-   !!   isfmlt       : compute iceshelf melt and heat flux
+   !!   isf          : define and allocate ice shelf variables
    !!----------------------------------------------------------------------
 
@@ -32,4 +32,5 @@
    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)
@@ -47,14 +48,15 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mskisf_par, mskisf_cav   !: Level of ice shelf base
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl_cav, rhisf_tbl_par !: thickness of tbl  [m]
+   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(:,:)   ::   risfdep                  !: Iceshelf draft                              (ISF)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   bathy                    !: Bathymetry (needed for isf tbl definition)  (ISF)
    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]  
@@ -91,5 +93,5 @@
    SUBROUTINE isf_alloc_par()
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isfmlt_alloc  ***
+      !!                  ***  ROUTINE isf_alloc_par  ***
       !!
       !! ** Purpose : 
@@ -105,4 +107,16 @@
       ierr = ierr + ialloc
       !
+      ALLOCATE(misfkt_par(jpi,jpj), misfkb_par(jpi,jpj), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rfrac_tbl_par(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rhisf_tbl_par(jpi,jpj), rhisf0_tbl_par(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( mskisf_par(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
       CALL mpp_sum ( 'isf', ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
@@ -111,5 +125,5 @@
    SUBROUTINE isf_alloc_cav()
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isfmlt_alloc  ***
+      !!                  ***  ROUTINE isf_alloc_cav  ***
       !!
       !! ** Purpose : 
@@ -122,4 +136,16 @@
       ierr = 0       ! set to zero if no array to be allocated
       !
+      ALLOCATE(misfkt_cav(jpi,jpj), misfkb_cav(jpi,jpj), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rfrac_tbl_cav(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rhisf_tbl_cav(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( mskisf_cav(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
       CALL mpp_sum ( 'isf', ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
@@ -128,5 +154,5 @@
    SUBROUTINE isf_alloc()
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isfmlt_alloc  ***
+      !!                  ***  ROUTINE isf_alloc  ***
       !!
       !! ** Purpose : 
@@ -140,26 +166,16 @@
       ierr = 0       ! set to zero if no array to be allocated
       !
-      ALLOCATE(misfkt_par(jpi,jpj), misfkb_par(jpi,jpj),             &
-         &     misfkt_cav(jpi,jpj), misfkb_cav(jpi,jpj), STAT=ialloc )
+      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(fwfisf_par(jpi,jpj), fwfisf_par_b(jpi,jpj),             &
-         &     fwfisf_cav(jpi,jpj), fwfisf_cav_b(jpi,jpj), STAT=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),             &
-         &     risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )
+      ALLOCATE(risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )
       ierr = ierr + ialloc
       !
-      ALLOCATE( rfrac_tbl_cav(jpi,jpj), rfrac_tbl_par(jpi,jpj), STAT=ialloc)
-      ierr = ierr + ialloc
-      !
-      ALLOCATE( rhisf_tbl_par(jpi,jpj), rhisf_tbl_cav(jpi,jpj), STAT=ialloc)
-      ierr = ierr + ialloc
-      !
-      ALLOCATE( mskisf_cav(jpi,jpj), mskisf_par(jpi,jpj), STAT=ialloc)
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risfdep(jpi,jpj), bathy(jpi,jpj), STAT=ialloc)
+      ALLOCATE(risfload(jpi,jpj), STAT=ialloc)
       ierr = ierr + ialloc
       !
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11403)
@@ -156,9 +156,8 @@
       ! allocation isfcav gamtisf, gamsisf, 
       CALL isf_alloc_cav()
-
-      ! initialisation
-      mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
       !
-      misfkt_cav(:,:) = mikt(:,:)
+      ! cav
+      misfkt_cav(:,:)     = mikt(:,:) ; misfkb_cav(:,:)       = 1
+      rhisf_tbl_cav(:,:)  = 0.0_wp    ; rfrac_tbl_cav(:,:)    = 0.0_wp
       !
       SELECT CASE ( TRIM(cn_isfcav_mlt) )
Index: MO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav_gam.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav_gam.F90	(revision 11402)
+++ 	(revision )
@@ -1,230 +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
-      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. (2015)
-      !!---------------------------------------------------------------------
-      !!-------------------------- 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: MO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav_mlt.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav_mlt.F90	(revision 11402)
+++ 	(revision )
@@ -1,215 +1,0 @@
-MODULE isfcavmlt
-   !!======================================================================
-   !!                       ***  MODULE  isfcav_mlt  ***
-   !! ice shelf module :  update surface ocean boundary condition under ice
-   !!                   shelf
-   !!======================================================================
-   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   isfcav_mlt    : 
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE isf
-   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 in_out_manager ! I/O manager
-   USE iom            ! I/O library
-   USE fldread        ! read input field at current time step
-   USE lbclnk         !
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   isfcav_mlt
-
-   !!----------------------------------------------------------------------
-   !! 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 SUBROUTINE ----------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
-
-   SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl,   &
-      &                           pqhc, pqoce, pqfwf    )
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : 
-      !!
-      !! ** Method     : 
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
-      !!-------------------------- IN  -------------------------------------
-      INTEGER, INTENT(in) :: kt
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
-      !!---------------------------------------------------------------------
-      !!---------------------------------------------------------------------
-      !
-      ! compute latent heat and melt (2d)
-      SELECT CASE ( cn_isfcav_mlt )
-      CASE ( 'spe' )   ! ice shelf melt specified (read input file, and heat fluxes derived from
-         CALL isfmlt_spe( kt, pstbl,                  &
-            &                  pqhc, pqoce, pqfwf    )
-      CASE ( '2eq' )   !  ISOMIP  formulation (2 equations) for volume flux (Hunter et al., 2006)
-         CALL isfmlt_2eq( pgt, pttbl, pstbl,          &
-            &                  pqhc , pqoce, pqfwf    )
-      CASE ( '3eq' )   ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
-         CALL isfmlt_3eq( pgt, pgs , pttbl, pstbl,   &
-            &                  pqhc, pqoce, pqfwf    )
-      CASE ( 'oasis' ) ! fwf pass trough oasis
-         !CALL isfmlt_oasis( kt, pstbl,               &
-         !   &                    zqhc, zqoce, pqfwf  )
-      CASE DEFAULT
-         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
-      END SELECT
-      !
-   END SUBROUTINE isfcav_mlt
-
-! -------------------------------------------------------------------------------------------------------
-! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
-
-   SUBROUTINE isfmlt_spe(kt, pstbl,                  &  ! <<== in
-      &                       pqhc, pqoce, pqfwf    )  ! ==>> out
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : 
-      !!
-      !! ** Method     : 
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
-      !!-------------------------- IN  -------------------------------------
-      INTEGER                     , INTENT(in   ) :: kt
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl         ! salinity in top boundary layer
-      !!--------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
-      !!--------------------------------------------------------------------
-      !
-      ! Calculate freezing temperature
-      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
-      !
-      ! read input file
-      CALL fld_read ( kt, nn_fsbc, sf_isfcav_fwf )
-      !
-      ! define fwf and qoce
-      pqfwf(:,:) = -sf_isfcav_fwf(1)%fnow(:,:,1)            ! fwf
-      ! ocean heat flux is assume to be equal to the latent heat
-      pqoce(:,:) =  pqfwf(:,:) * rLfusisf             ! ocean heat flux
-      pqhc (:,:) =  pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux
-      !
-   END SUBROUTINE isfmlt_spe
-
-   SUBROUTINE isfmlt_2eq(pgt, pttbl, pstbl,          &  ! <<== in
-      &                       pqhc , pqoce, pqfwf    )  ! ==>> out
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : 
-      !!
-      !! ** Method     : 
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt           ! temperature exchange coeficient
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
-      !!--------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
-      !!--------------------------------------------------------------------
-      !
-      ! Calculate freezing temperature
-      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
-      !
-      ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
-      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * ( pttbl(:,:)-ztfrz(:,:) ) * r1_Lfusisf  ! fresh water flux ( > 0 out )
-      pqoce(:,:) = - pqfwf(:,:) * rLfusisf                                         ! ocea-ice flux (assume to be equal to latent heat flux) ( > 0 out )
-      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                                 ! heat content flux ( > 0 out )
-      !
-      ! output thermal driving
-      CALL iom_put('isfthermald_cav',( pttbl(:,:)-ztfrz(:,:) ))
-      !
-   END SUBROUTINE isfmlt_2eq
-
-   SUBROUTINE isfmlt_3eq(pgt, pgs , pttbl, pstbl,   &
-      &                       pqhc, pqoce, pqfwf    )
-      !!----------------------------------------------------------------------
-      !! ** Purpose    : 
-      !!
-      !! ** Method     : 
-      !!---------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! temperature exchange coeficient
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
-      !!--------------------------------------------------------------------
-      REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
-      REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
-      REAL(wp) :: zeps = 1.e-20
-      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz             ! freezing point
-      REAL(wp), DIMENSION(jpi,jpj) :: zqcon             ! conductive flux through the ice shelf
-      !
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      !!--------------------------------------------------------------------
-      !
-      ! compute upward heat flux zhtflx and upward water flux zwflx
-      ! Resolution of a 2d equation from equation 24, 25 and 26
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            ! compute coeficient to solve the 2nd order equation
-            zeps1 = rau0_rcp*pgt(ji,jj)
-            zeps2 = rLfusisf*rau0*pgs(ji,jj)
-            zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps)
-            zeps4 = risf_lamb2+risf_lamb3*risfdep(ji,jj)
-            zeps6 = zeps4-pttbl(ji,jj)
-            zeps7 = zeps4-rtsurf
-            !
-            ! solve the 2nd order equation to find zsfrz
-            zaqe  = risf_lamb1 * (zeps1 + zeps3)
-            zaqer = 0.5_wp/MIN(zaqe,-risf_eps)
-            zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
-            zcqe  = zeps2*pstbl(ji,jj)
-            zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
-            !
-            ! Presumably zdis can never be negative because gammas is very small compared to gammat
-            zsfrz=(-zbqe-SQRT(zdis))*zaqer
-            IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer  ! check this if this if is needed
-            !
-            ! compute t freeze (eq. 25)
-            ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
-            !
-            ! compute the upward water and heat flux (eq. 24 and eq. 26)
-            ! ocean heat content flux added later on.
-            pqfwf(ji,jj) = rau0 * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux (> 0 out)
-            pqoce(ji,jj) = pgt(ji,jj) * rau0_rcp * ( pttbl(ji,jj) - ztfrz(ji,jj) )        ! ocean-ice heat flux (> 0 out)
-            pqhc (ji,jj) = pqfwf(ji,jj) * ztfrz(ji,jj) * rcp                              ! heat content   flux (> 0 out)
-            zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf ) ! to be check
-            !
-         END DO
-      END DO
-      !
-      ! output conductive heat flux through the ice
-      CALL iom_put('qconisf', zqcon)
-      !
-      ! output thermal driving
-      CALL iom_put('isfthermald_cav',( pttbl(:,:) - ztfrz(:,:) ))
-      !
-   END SUBROUTINE isfmlt_3eq
-
-   !SUBROUTINE isfmlt_3eq_frz_ktm1
-   ! compute tfrz based on sfrz value at kt-1 (need to be SAVED local array)
-   ! 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/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 11403)
@@ -0,0 +1,288 @@
+MODULE isfcavmlt
+   !!======================================================================
+   !!                       ***  MODULE  isfcav_mlt  ***
+   !! ice shelf module :  update surface ocean boundary condition under ice
+   !!                   shelf
+   !!======================================================================
+   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfcav_mlt    : 
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE isf            ! ice shelf public variables
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE fldread        ! read input field at current time step
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   isfcav_mlt
+
+   !!----------------------------------------------------------------------
+   !! 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 SUBROUTINE ----------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl, &
+      &                           pqhc, pqoce, pqfwf  )
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose    : compute ice shelf fwf/heqt fluxes 
+      !!
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER, INTENT(in) :: kt
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
+      !!---------------------------------------------------------------------
+      !!---------------------------------------------------------------------
+      !
+      ! compute latent heat and melt (2d)
+      SELECT CASE ( cn_isfcav_mlt )
+      CASE ( 'spe' )   ! ice shelf melt specified (read input file, and heat fluxes derived from
+         CALL isfmlt_spe( kt, pstbl,               &
+            &                  pqhc, pqoce, pqfwf  )
+      CASE ( '2eq' )   !  ISOMIP  formulation (2 equations) for volume flux (Hunter et al., 2006)
+         CALL isfmlt_2eq( pgt, pttbl, pstbl,       &
+            &                  pqhc , pqoce, pqfwf )
+      CASE ( '3eq' )   ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
+         CALL isfmlt_3eq( pgt, pgs , pttbl, pstbl, &
+            &                  pqhc, pqoce, pqfwf  )
+      CASE ( 'oasis' ) ! fwf pass trough oasis
+         CALL isfmlt_oasis( kt, pstbl,              &
+            &                   pqhc, pqoce, pqfwf  )
+      CASE DEFAULT
+         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
+      END SELECT
+      !
+   END SUBROUTINE isfcav_mlt
+
+! -------------------------------------------------------------------------------------------------------
+! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE isfmlt_spe(kt, pstbl,              &  ! <<== in
+      &                      pqhc , pqoce, pqfwf )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : - read ice shelf melt from forcing file
+      !!                 - compute ocea-ice heat flux (assuming it is equal to latent heat)
+      !!                 - compute heat content flux
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
+      !!--------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
+      !!--------------------------------------------------------------------
+      !
+      ! Calculate freezing temperature
+      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+      !
+      ! read input file
+      CALL fld_read ( kt, nn_fsbc, sf_isfcav_fwf )
+      !
+      ! define fwf and qoce
+      ! ocean heat flux is assume to be equal to the latent heat
+      pqfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)      ! fwf                ( >0 out)
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
+      !
+   END SUBROUTINE isfmlt_spe
+
+   SUBROUTINE isfmlt_2eq(pgt, pttbl, pstbl,       &  ! <<== in
+      &                       pqhc , pqoce, pqfwf )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
+      !!
+      !! ** Method     : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux.
+      !!                 From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being :
+      !!                   qfwf  = Gammat * Rau0 * Cp * ( Tw - Tfrz ) / Lf 
+      !!                   qhoce = qlat
+      !!                   qhc   = qfwf * Cp * Tfrz
+      !!
+      !! ** Reference  : Hunter,  J.  R.:  Specification  for  test  models  of  ice  shelf  cavities,  
+      !!                 Tech.  Rep.  June,  Antarctic  Climate  &  Ecosystems  Cooperative  Research  Centre,  available  at:  
+      !!                 http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt           ! temperature exchange coeficient
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
+      !!--------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
+      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
+      !!--------------------------------------------------------------------
+      !
+      ! Calculate freezing temperature
+      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+      !
+      ! thermal driving
+      zthd (:,:) = pttbl(:,:) - ztfrz(:,:)
+      !
+      ! 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 )
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf                         ! ocea-ice flux     ( > 0 out )
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                 ! heat content flux ( > 0 out )
+      !
+      ! output thermal driving
+      CALL iom_put('isfthermald_cav', zthd )
+      !
+   END SUBROUTINE isfmlt_2eq
+
+   SUBROUTINE isfmlt_3eq(pgt, pgs , pttbl, pstbl, &
+      &                       pqhc, pqoce, pqfwf  )
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using the 3 equation formulation 
+      !!
+      !! ** Method     : The melt rate is determined considering the heat balance, the salt balance
+      !!                 at the phase change interface and a linearisation of the equation of state.
+      !!
+      !! ** Reference  : - Holland, D. M. and Jenkins, A.,
+      !!                   Modeling Thermodynamic Ice-Ocean Interactions at the Base of an Ice Shelf,
+      !!                   J. Phys. Oceanogr., 29, 1999.
+      !!                 - Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, 
+      !!                   R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, 
+      !!                   P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.:
+      !!                   Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: 
+      !!                   MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), 
+      !!                   Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016. 
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs          ! heat/salt exchange coeficient
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl        ! mean temperature and salinity in top boundary layer
+      !!--------------------------------------------------------------------
+      REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7       ! dummy local scalar for quadratic equation resolution
+      REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac     ! dummy local scalar for quadratic equation resolution
+      REAL(wp) :: zeps = 1.e-20
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing point
+      REAL(wp), DIMENSION(jpi,jpj) :: zqcon         ! conductive flux through the ice shelf
+      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      !!--------------------------------------------------------------------
+      !
+      ! compute upward heat flux zhtflx and upward water flux zwflx
+      ! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            ! compute coeficient to solve the 2nd order equation
+            zeps1 = rau0_rcp*pgt(ji,jj)
+            zeps2 = rLfusisf*rau0*pgs(ji,jj)
+            zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps)
+            zeps4 = risf_lamb2+risf_lamb3*risfdep(ji,jj)
+            zeps6 = zeps4-pttbl(ji,jj)
+            zeps7 = zeps4-rtsurf
+            !
+            ! solve the 2nd order equation to find zsfrz
+            zaqe  = risf_lamb1 * (zeps1 + zeps3)
+            zaqer = 0.5_wp/MIN(zaqe,-risf_eps)
+            zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
+            zcqe  = zeps2*pstbl(ji,jj)
+            zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
+            !
+            ! Presumably zdis can never be negative because gammas is very small compared to gammat
+            zsfrz=(-zbqe-SQRT(zdis))*zaqer
+            IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer  ! check this if this if is needed
+            !
+            ! compute t freeze (eq. 25)
+            ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
+            !
+            ! thermal driving
+            zthd(ji,jj) = pttbl(ji,jj) - ztfrz(ji,jj)
+            !
+            ! compute the upward water and heat flux (eq. 24 and eq. 26)
+            pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out)
+            pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out)
+            pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out)
+            !
+            zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )                                  ! conductive flux through the ice (> 0 out)
+            !
+         END DO
+      END DO
+      !
+      ! output conductive heat flux through the ice
+      CALL iom_put('qconisf', zqcon)
+      !
+      ! output thermal driving
+      CALL iom_put('isfthermald_cav', zthd)
+      !
+   END SUBROUTINE isfmlt_3eq
+
+   SUBROUTINE isfmlt_oasis(kt, pstbl,              &  ! <<== in
+      &                        pqhc , pqoce, pqfwf )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
+      !!
+      !! ** Purpose    : - read ice shelf melt from forcing file => pattern
+      !!                 - total amount of fwf is given by sbccpl (fwfisf_cpl)
+      !!                 - scale fwf and compute heat fluxes
+      !!
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
+      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), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
+      REAL(wp), DIMENSION(jpi,jpj) :: zfwf                               ! 2d fwf map after scaling
+      !!--------------------------------------------------------------------
+      !
+      ! Calculate freezing temperature
+      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+      !
+      ! read input file
+      CALL fld_read ( kt, nn_fsbc, sf_isfcav_fwf )
+      !
+      ! ice shelf 2d map
+      zfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)
+      !
+      ! compute glob sum from input file
+      zfwf_fld = glob_sum('isfcav_mlt', zfwf(:,:))
+      !
+      ! compute glob sum from atm->oce ice shelf fwf
+      zfwf_cpl = glob_sum('isfcav_mlt', fwfisf_cpl(:,:))
+      !
+      ! scale fwf
+      zfwf(:,:) = zfwf(:,:) * zfwf_cpl / zfwf_fld
+      ! 
+      ! define fwf and qoce
+      ! ocean heat flux is assume to be equal to the latent heat
+      pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out)
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
+      !
+   END SUBROUTINE isfmlt_oasis
+
+   !SUBROUTINE isfmlt_3eq_frz_ktm1
+   ! compute tfrz based on sfrz value at kt-1 (need to be SAVED local array)
+   ! 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/isfgammats.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfgammats.F90	(revision 11403)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfgammats.F90	(revision 11403)
@@ -0,0 +1,231 @@
+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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 11403)
@@ -1,15 +1,13 @@
 MODULE isfhdiv
 
-USE oce
-USE dom_oce
-USE isf
-USE isfutils
-USE phycst
+   USE dom_oce
+   USE isf
+   USE phycst
 
-IMPLICIT NONE
+   IMPLICIT NONE
 
-PRIVATE
+   PRIVATE
 
-PUBLIC isf_hdiv
+   PUBLIC isf_hdiv
 
 CONTAINS
@@ -59,5 +57,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt, ikb 
-      REAL(wp), DIMENSION(jpi,jpj) :: zqvol
+      REAL(wp), DIMENSION(jpi,jpj) :: zqvol,ztmp
       !!----------------------------------------------------------------------
       !
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 11403)
@@ -0,0 +1,117 @@
+MODULE isfload
+   !!======================================================================
+   !!                       ***  MODULE  isfload  ***
+   !! isfload module :  compute ice shelf load (needed for the hpg)
+   !!======================================================================
+   !! History :  4.1  !  2019-09  (P. Mathiot) original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfload      : compute ice shelf load
+   !!----------------------------------------------------------------------
+
+   USE isf             ! ice shelf variables
+   USE dom_oce         ! vertical scale factor
+   USE in_out_manager  ! 
+   USE eosbn2          ! eos routine
+   USE lib_mpp         ! ctl_stop routine
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isf_load
+
+CONTAINS
+
+   SUBROUTINE isf_load ( pisfload )
+      !!--------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_load  ***
+      !!
+      !! ** Purpose : compute the ice shelf load
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pisfload
+      !!-------------------------- IN  -------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      ! quality test: ice shelf in a stratify/uniform ocean should not drive any flow.
+      !               the smaller the residual flow is, the better it is.
+      !
+      ! ice shelf cavity
+      SELECT CASE ( cn_isfload )
+      CASE ( 'isomip' )
+         CALL isf_load_isomip ( pisfload )
+      CASE DEFAULT
+         CALL ctl_stop('STOP','method cn_isfload to compute ice shelf load does not exist (isomip), check your namelist')
+      END SELECT
+      !
+   END SUBROUTINE isf_load
+
+   SUBROUTINE isf_load_isomip( pisfload )
+      !!--------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_load  ***
+      !!
+      !! ** Purpose : compute the ice shelf load
+      !!
+      !! ** Method  : The ice shelf is assumed to be in hydro static equilibrium
+      !!              in water at -1.9 C and 34.4 PSU. Weight of the ice shelf is
+      !!              integrated from top to bottom.
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pisfload
+      !!-------------------------- IN  -------------------------------------
+      !!--------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk
+      INTEGER  :: ikt
+      REAL(wp)                          :: znad        ! 
+      REAL(wp), DIMENSION(jpi,jpj)      :: zrhdtop_isf ! water density    displaced by the ice shelf (at the interface)
+      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top     ! water properties displaced by the ice shelf   
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: zrhd        ! water density    displaced by the ice shelf
+      !!----------------------------------------------------------------------
+      !
+      ! print
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '   ice shelf case: set the ice-shelf load following isomip case'
+      !
+      znad = 1._wp                     !- To use density and not density anomaly
+      !
+      !                                !- 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
+      !
+      DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
+         CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) )
+      END DO
+      !
+      !                                !- compute rhd at the ice/oce interface (ice shelf side)
+      CALL eos( zts_top , risfdep, zrhdtop_isf )
+      !
+      !                                !- Surface value + ice shelf gradient
+      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) )
+            !
+         END DO
+      END DO
+      !
+      ! mask it 
+      risfload(:,:) = risfload(:,:) * mskisf_cav(:,:)  ! need to be saved for diaar5
+      !
+   END SUBROUTINE isf_load_isomip
+
+END MODULE isfload
Index: MO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfmlt.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfmlt.F90	(revision 11402)
+++ 	(revision )
@@ -1,202 +1,0 @@
-MODULE isfmlt
-   !!======================================================================
-   !!                       ***  MODULE  sbcisf  ***
-   !! Surface module :  compute iceshelf melt and heat flux
-   !!======================================================================
-   !! 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 param/explicit ice shelf and re-organisation
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   isfmlt       : compute iceshelf melt and heat flux
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constants
-   USE eosbn2         ! equation of state
-   USE sbc_oce        ! surface boundary condition: ocean fields
-   USE zdfdrg         ! vertical physics: top/bottom drag coef.
-   !
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O library
-   USE fldread        ! read input field at current time step
-   USE lbclnk         !
-   USE lib_fortran    ! glob_sum
-   !
-   USE isfrst         ! iceshelf restart
-   USE isftbl         ! ice shelf boundary layer
-   USE isfpar         ! ice shelf parametrisation
-   USE isfcav         ! ice shelf cavity
-   USE isf            ! isf variables
-
-   IMPLICIT NONE
-
-   PRIVATE
-
-   PUBLIC   isf_mlt, isf_mlt_init  ! routine called in sbcmod and divhor
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
- 
-  SUBROUTINE isf_mlt( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_mlt  ***
-      !!
-      !! ** Purpose : 
-      !!
-      !! ** Method  : 
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !!---------------------------------------------------------------------
-      !
-      IF ( ln_isfcav_mlt ) THEN
-         !
-         ! before time step 
-         IF ( kt /= nit000 ) THEN 
-            risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
-            fwfisf_cav_b(:,:)      = fwfisf_cav(:,:)
-         END IF
-         !
-         ! compute tbl lvl/h
-         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
-         !
-         ! compute ice shelf melt
-         CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav)
-         !
-         ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
-         IF (lrst_oce) CALL isfrst_write(kt, 'cav', risf_cav_tsc, fwfisf_cav)
-         !
-      END IF
-      ! 
-      IF ( ln_isfpar_mlt ) THEN
-         !
-         ! before time step 
-         IF ( kt /= nit000 ) THEN 
-            risf_par_tsc_b (:,:) = risf_par_tsc (:,:)
-            fwfisf_par_b(:,:)    = fwfisf_par(:,:)
-         END IF
-         !
-         ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
-         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
-         !
-         ! compute ice shelf melt
-         CALL isf_par( kt, risf_par_tsc, fwfisf_par)
-         !
-         ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
-         IF (lrst_oce) CALL isfrst_write(kt, 'par', risf_par_tsc, fwfisf_par)
-         !
-      END IF
-      !
-   END SUBROUTINE isf_mlt
-
-   SUBROUTINE isf_mlt_init
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isfmlt_init  ***
-      !!
-      !! ** Purpose :
-      !!
-      !! ** Method  :
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER               :: inum, ierror
-      INTEGER               :: ios                  ! Local integer output status for namelist read
-      INTEGER               :: ikt, ikb
-      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
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
-      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
-      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namisf )
-      !
-      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,*) ''
-      !
-      ! Allocate public array
-      CALL isf_alloc()
-      !
-      ! initialisation
-      ! cav
-      misfkt_cav(:,:)     = mikt(:,:) ; misfkb_cav(:,:)       = mikt(:,:)
-      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
-      rhisf_tbl_cav(:,:)  = 1e-20     ; rfrac_tbl_cav(:,:)    = 0.0_wp
-      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
-      !
-      mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
-      !
-      ! par
-      misfkt_par(:,:)     = 1         ; misfkb_par(:,:)       = 1         
-      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
-      rhisf_tbl_par(:,:)  = 1e-20     ; rfrac_tbl_par(:,:)    = 0.0_wp
-      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
-      !
-      mskisf_par(:,:) = 0
-      !
-      ! initialisation useful variable
-      r1_Lfusisf =  1._wp / rLfusisf
-      !
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            ikt = mikt(ji,jj)
-            ikb = mbkt(ji,jj)
-            bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1)
-            risfdep(ji,jj) = gdepw_0(ji,jj,ikt  )
-         END DO
-      END DO
-      !
-      IF ( ln_isfcav_mlt ) THEN
-         !
-         ! initialisation  of cav variable
-         CALL isf_cav_init()
-         !
-         ! read cav variable from restart
-         IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
-         !
-      END IF
-      !
-      IF ( ln_isfpar_mlt ) THEN
-         !
-         ! initialisation  of par variable
-         CALL isf_par_init()
-         !
-         ! read par variable from restart
-         IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
-         !
-      END IF
-      !
-  END SUBROUTINE isf_mlt_init
-   !!======================================================================
-END MODULE isfmlt
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfnxt.F90	(revision 11403)
@@ -1,22 +1,38 @@
 MODULE isfnxt
+   !!=========================================================================
+   !!                       ***  MODULE  isfnxt  ***
+   !! Ice shelf update: compute the dynnxt ice shelf contribution
+   !!=========================================================================
+   !! History :  OPA  !  2019-09  (P. Mathiot)  Original code
+   !!-------------------------------------------------------------------------
+  
+   !!-------------------------------------------------------------------------
+   !!   dyn_nxt       : obtain the next (after) horizontal velocity
+   !!-------------------------------------------------------------------------
 
-USE isf
-USE isfutils
+   USE isf
+   USE dom_oce
+   USE in_out_manager
 
-USE dom_oce
-USE in_out_manager
+   IMPLICIT NONE
 
-IMPLICIT NONE
+   PRIVATE
 
-PRIVATE
-
-PUBLIC isf_dynnxt !, isf_tranxt_mlt, isf_dynnxt_cpl, isf_tranxt_cpl
+   PUBLIC isf_dynnxt ! isf_tranxt
 
 CONTAINS
 
    SUBROUTINE isf_dynnxt ( pcoef )
-      !!----------------------------------------------------------------------
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_dynnxt  ***
+      !!
+      !! ** Purpose : compute the ice shelf volume filter correction for cavity, param, ice sheet coupling case
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
       REAL(wp),                     INTENT(in   ) :: pcoef           ! atfp * rdt * r1_rau0
-      !!----------------------------------------------------------------------
+      !!-------------------------- IN  -------------------------------------
+      !!--------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !
       ! ice shelf cavity
@@ -32,6 +48,12 @@
 
    SUBROUTINE isf_dynnxt_mlt ( ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, pcoef )
-   !atfp * rdt * r1_rau0
-      !!----------------------------------------------------------------------
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_dynnxt_mlt  ***
+      !!
+      !! ** Purpose : compute the ice shelf volume filter correction for cavity or param
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      !!-------------------------- IN  -------------------------------------
       INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot     ! top and bottom level of tbl
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl    ! fraction of bottom cell included in tbl, tbl thickness
@@ -46,5 +68,5 @@
       zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:)
       !
-      ! add the in depth increment
+      ! add the increment in the tbl
       DO jk = 1, jpkm1
          DO jj = 1, jpj
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 11403)
@@ -58,5 +58,4 @@
       INTEGER, INTENT(in) ::   kt                                           ! ocean time step
       !!---------------------------------------------------------------------
-      !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh
       !!---------------------------------------------------------------------
@@ -106,4 +105,10 @@
       CALL isf_alloc_par()
       !
+      ! par
+      misfkt_par(:,:)     = 1         ; misfkb_par(:,:)       = 1         
+      rhisf_tbl_par(:,:)  = 1e-20     ; rfrac_tbl_par(:,:)    = 0.0_wp
+      !
+      mskisf_par(:,:) = 0
+      !
       ! define isf tbl tickness, top and bottom indice
       CALL read_2dcstdta(TRIM(sn_isfpar_zmax%clname), TRIM(sn_isfpar_zmax%clvar), ztblmax)
@@ -126,5 +131,5 @@
       !
       ! compute misfkb_par, rhisf_tbl
-      rhisf_tbl(:,:) = rhisf0_tbl_par(:,:)
+      rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
       CALL isf_tbl_lvl( ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
       !
Index: MO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar_mlt.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar_mlt.F90	(revision 11402)
+++ 	(revision )
@@ -1,159 +1,0 @@
-MODULE isfparmlt
-   !!======================================================================
-   !!                       ***  MODULE  sbcisf  ***
-   !! Surface module :  update surface ocean boundary condition under ice
-   !!                   shelf
-   !!======================================================================
-   !! History :  4.0  ! original code
-   !!----------------------------------------------------------------------
-
-   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 in_out_manager ! I/O manager
-   USE iom            ! I/O library
-   USE fldread
-
-   IMPLICIT NONE
-
-   PRIVATE
-
-   PUBLIC  isfpar_mlt 
-   
-   !!----------------------------------------------------------------------
-   !! 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 SUBROUTINE ----------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
-
-  SUBROUTINE isfpar_mlt( kt, pqfwf, pqoce, pqhc )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_isf  ***
-      !!
-      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
-      !!              melting and freezing 
-      !!
-      !! ** Method  :  2 parameterizations are available according to XXXXX
-      !!                        2 : Beckmann & Goose parameterization
-      !!                        3 : Specified runoff in deptht (Mathiot & al. 2017)
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqfwf,pqoce, pqhc  ! fresh water, ice-ocean heat and heat content fluxes
-      !!-------------------------- IN  -------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !!---------------------------------------------------------------------
-      !
-      ! Choose among the available ice shelf parametrisation
-      SELECT CASE ( cn_isfpar_mlt )
-      CASE ( 'spe' )    ! specified runoff in depth (Mathiot et al., 2017 in preparation)
-         CALL isfpar_mlt_spe(kt, pqhc, pqoce, pqfwf)
-      CASE ( 'bg03' )    ! Beckmann and Goosse parametrisation 
-         CALL isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
-      CASE ( 'oasis' )
-         !CALL isfpar_mlt_oasis
-      CASE DEFAULT
-         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfpar (should not see this)')
-      END SELECT
-      !
-   END SUBROUTINE isfpar_mlt
-
-! -------------------------------------------------------------------------------------------------------
-! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
-
-   SUBROUTINE isfpar_mlt_spe(kt, pqhc, pqfwf, pqoce)
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_isf_bg03  ***
-      !!
-      !! ** Purpose : prescribed ice shelf melting in case ice shelf cavities are closed.
-      !!              data read into a forcing files.
-      !!
-      !! ** Reference : Mathiot et al. (2017)
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf,pqoce  ! fresh water and ice-ocean heat fluxes
-      !!-------------------------- IN  -------------------------------------
-      INTEGER,  INTENT(in) :: kt
-      !!--------------------------------------------------------------------
-      INTEGER :: jk
-      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: ztfrz3d
-      REAL(wp), DIMENSION(jpi,jpj)      :: ztfrz
-      !!--------------------------------------------------------------------
-      !
-      ! specified runoff in depth (Mathiot et al., 2017)
-      CALL fld_read ( kt, nn_fsbc, sf_isfpar_fwf   )
-      !
-      ! compute ptfrz
-      ! 0. ------------Mean freezing point
-      DO jk = 1,jpk
-         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
-      END DO
-      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
-      !
-      pqfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)      ! fresh water flux from the isf (fwfisf <0 mean melting) 
-      pqoce(:,:) =   pqfwf(:,:) * rLfusisf             ! ocean/ice shelf flux assume to be equal to latent heat flux
-      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
-      !
-   END SUBROUTINE isfpar_mlt_spe
-
-   SUBROUTINE isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_isf_bg03  ***
-      !!
-      !! ** Purpose : compute an estimate of ice shelf melting in case cavities are closed
-      !!              based on the far fields T and S properties. 
-      !!
-      !! ** Method  :   See reference
-      !!
-      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
-      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
-      !!         (hereafter BG)
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf,pqoce  ! fresh water and ice-ocean heat fluxes
-      !!-------------------------- IN  -------------------------------------
-      INTEGER,  INTENT(in) :: kt
-      !!--------------------------------------------------------------------
-      INTEGER :: jk
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d        ! freezing point
-      REAL(wp), DIMENSION(jpi,jpj)     :: ztfrz          ! freezing point
-      REAL(wp), DIMENSION(jpi,jpj)     :: ztavg          ! temperature avg
-      !!----------------------------------------------------------------------
-      !
-      ! 0. ------------Mean freezing point
-      DO jk = 1,jpk
-         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
-      END DO
-      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
-      !
-      ! 1. ------------Mean temperature
-      CALL isf_tbl(tsn(:,:,jk,jp_tem), ztavg, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
-      !
-      ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
-      pqoce(:,:) = - rau0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztfrz(:,:) - ztavg(:,:) ) * r1_e1e2t(:,:)
-      pqfwf(:,:) =   pqoce(:,:) / rLfusisf
-      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
-      !
-      ! output ttbl
-      CALL iom_put('ttbl_par', ztavg(:,:) )
-      !
-      ! output thermal driving
-      CALL iom_put('isfthermald_par',( ztfrz(:,:) - ztavg(:,:) ))
-      !
-      !
-   END SUBROUTINE isfpar_mlt_bg03
-
-   SUBROUTINE isfpar_mlt_oasis
-   !TODO
-   END SUBROUTINE isfpar_mlt_oasis
-
-END MODULE isfparmlt
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 11403)
@@ -0,0 +1,210 @@
+MODULE isfparmlt
+   !!======================================================================
+   !!                       ***  MODULE  isfparmlt  ***
+   !! Surface module :  update surface ocean boundary condition under ice
+   !!                   shelf
+   !!======================================================================
+   !! History :  4.0  ! original code
+   !!----------------------------------------------------------------------
+
+   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 in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE fldread
+   USE lib_fortran
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC  isfpar_mlt 
+   
+   !!----------------------------------------------------------------------
+   !! 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 SUBROUTINE ----------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+  SUBROUTINE isfpar_mlt( kt, pqfwf, pqoce, pqhc )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf  ***
+      !!
+      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
+      !!              melting and freezing 
+      !!
+      !! ** Method  :  2 parameterizations are available according
+      !!                        1 : Specified melt flux
+      !!                        2 : Beckmann & Goose parameterization
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqfwf,pqoce, pqhc  ! fresh water, ice-ocean heat and heat content fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      ! Choose among the available ice shelf parametrisation
+      SELECT CASE ( cn_isfpar_mlt )
+      CASE ( 'spe' )    ! specified runoff in depth (Mathiot et al., 2017 in preparation)
+         CALL isfpar_mlt_spe(kt, pqhc, pqoce, pqfwf)
+      CASE ( 'bg03' )    ! Beckmann and Goosse parametrisation 
+         CALL isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
+      CASE ( 'oasis' )
+         CALL isfpar_mlt_oasis( kt, pqhc, pqoce, pqfwf)
+      CASE DEFAULT
+         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfpar (should not see this)')
+      END SELECT
+      !
+   END SUBROUTINE isfpar_mlt
+
+! -------------------------------------------------------------------------------------------------------
+! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE isfpar_mlt_spe(kt, pqhc, pqfwf, pqoce)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfpar_mlt_spe  ***
+      !!
+      !! ** Purpose : prescribed ice shelf melting in case ice shelf cavities are closed.
+      !!              data read into a forcing files.
+      !!
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf,pqoce  ! fresh water and ice-ocean heat fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  INTENT(in) :: kt
+      !!--------------------------------------------------------------------
+      INTEGER :: jk
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: ztfrz3d
+      REAL(wp), DIMENSION(jpi,jpj)      :: ztfrz
+      !!--------------------------------------------------------------------
+      !
+      ! 0. ------------Read specified runoff
+      CALL fld_read ( kt, nn_fsbc, sf_isfpar_fwf   )
+      !
+      ! compute ptfrz
+      ! 1. ------------Mean freezing point
+      DO jk = 1,jpk
+         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
+      END DO
+      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+      !
+      pqfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)      ! fresh water flux from the isf (fwfisf <0 mean melting) 
+      pqoce(:,:) =   pqfwf(:,:) * rLfusisf             ! ocean/ice shelf flux assume to be equal to latent heat flux
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
+      !
+   END SUBROUTINE isfpar_mlt_spe
+
+   SUBROUTINE isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfpar_mlt_bg03  ***
+      !!
+      !! ** Purpose : compute an estimate of ice shelf melting and 
+      !!              latent, ocean-ice and heat content heat fluxes
+      !!              in case cavities are closed based on the far fields T and S properties. 
+      !!
+      !! ** Method  : The ice shelf melt is computed as proportional to the differences between the 
+      !!              mean temperature and mean freezing point in front of the ice shelf averaged 
+      !!              over the ice shelf min ice shelf draft and max ice shelf draft and the freezing point
+      !!
+      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!                interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf,pqoce  ! fresh water and ice-ocean heat fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  INTENT(in) :: kt
+      !!--------------------------------------------------------------------
+      INTEGER :: jk
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d        ! freezing point
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztfrz          ! freezing point
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztavg          ! temperature avg
+      !!----------------------------------------------------------------------
+      !
+      ! 0. ------------Mean freezing point
+      DO jk = 1,jpk
+         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
+      END DO
+      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+      !
+      ! 1. ------------Mean temperature
+      CALL isf_tbl(tsn(:,:,jk,jp_tem), ztavg, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+      !
+      ! 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
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
+      !
+      ! 3. ------------BG03 output
+      ! output ttbl
+      CALL iom_put('ttbl_par', ztavg(:,:) )
+      !
+      ! output thermal driving
+      CALL iom_put('isfthermald_par',( ztfrz(:,:) - ztavg(:,:) ))
+      !
+      !
+   END SUBROUTINE isfpar_mlt_bg03
+
+   SUBROUTINE isfpar_mlt_oasis(kt, pqhc , pqoce, pqfwf )
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
+      !!
+      !! ** Purpose    : - read ice shelf melt from forcing file => pattern
+      !!                 - total amount of fwf is given by sbccpl (fwfisf_cpl)
+      !!                 - scale fwf and compute heat fluxes
+      !!
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
+      !!--------------------------------------------------------------------
+      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), DIMENSION(jpi,jpj)      :: ztfrz                         ! tbl freezing temperature
+      REAL(wp), DIMENSION(jpi,jpj)      :: zfwf                          ! 2d fwf map after scaling
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: ztfrz3d
+      !!--------------------------------------------------------------------
+      !
+      ! 0. ------------Read specified runoff
+      CALL fld_read ( kt, nn_fsbc, sf_isfpar_fwf   )
+      !
+      ! compute ptfrz
+      ! 1. ------------Mean freezing point
+      DO jk = 1,jpk
+         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
+      END DO
+      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+      !
+      ! ice shelf 2d map
+      zfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)
+      !
+      ! compute glob sum from input file
+      zfwf_fld = glob_sum('isfcav_mlt', zfwf(:,:))
+      !
+      ! compute glob sum from atm->oce ice shelf fwf
+      zfwf_cpl = glob_sum('isfcav_mlt', fwfisf_cpl(:,:))
+      !
+      ! scale fwf
+      zfwf(:,:) = zfwf(:,:) * zfwf_cpl / zfwf_fld
+      ! 
+      ! define fwf and qoce
+      ! ocean heat flux is assume to be equal to the latent heat
+      pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out )
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out ) (assumed to be the latent heat flux)
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out )
+      !
+   END SUBROUTINE isfpar_mlt_oasis
+
+END MODULE isfparmlt
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90	(revision 11403)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE  isfrst  ***
-   !! iceshelf restart module :read and write iceshelf variables in/from restart
+   !! iceshelf restart module :read/write iceshelf variables from/in restart
    !!======================================================================
    !! History :  4.1  !  2019-07  (P. Mathiot) Original code
@@ -30,12 +30,14 @@
    ! 
    SUBROUTINE isfrst_read(cdisf, ptsc, pfwf, ptsc_b, pfwf_b )
-      !!----------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       !!   isfrst_read : read iceshelf variables from restart
-      !!----------------------------------------------------------------------
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT --------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) :: pfwf_b
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(  out) :: ptsc_b
+      !!-------------------------- IN  --------------------------------------
       CHARACTER(LEN=256)               , INTENT(in   ) :: cdisf
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: pfwf
-      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(inout) :: pfwf_b
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc
-      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc_b
       !!----------------------------------------------------------------------
       CHARACTER(LEN=256) :: cfwf_b, chc_b, csc_b
@@ -67,14 +69,16 @@
    ! 
    SUBROUTINE isfrst_write(kt, cdisf, ptsc, pfwf )
-      !!----------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       !!   isfrst_write : write iceshelf variables in restart
-      !!----------------------------------------------------------------------
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT --------------------------------------
+      !!-------------------------- IN  --------------------------------------
       INTEGER                     , INTENT(in   ) :: kt
       CHARACTER(LEN=256)          , INTENT(in   ) :: cdisf
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: pfwf
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc
-      !!----------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       CHARACTER(LEN=256) :: cfwf_b, chc_b, csc_b
-      !!----------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       !
       ! ocean output print
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 11403)
@@ -0,0 +1,238 @@
+MODULE isfstp
+   !!======================================================================
+   !!                       ***  MODULE  isfstp  ***
+   !! Surface module :  compute iceshelf load, melt and heat flux
+   !!======================================================================
+   !! 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 param/explicit ice shelf and re-organisation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfstp       : compute iceshelf melt and heat flux
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state
+   USE sbc_oce        ! surface boundary condition: ocean fields
+   USE zdfdrg         ! vertical physics: top/bottom drag coef.
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE fldread        ! read input field at current time step
+   USE lbclnk         !
+   USE lib_fortran    ! glob_sum
+   !
+   USE isfrst         ! iceshelf restart
+   USE isftbl         ! ice shelf boundary layer
+   USE isfpar         ! ice shelf parametrisation
+   USE isfcav         ! ice shelf cavity
+   USE isfload        ! ice shelf load
+   USE isf            ! isf variables
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC   isf_stp, isf_stp_init  ! routine called in sbcmod and divhor
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+ 
+  SUBROUTINE isf_stp( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_stp  ***
+      !!
+      !! ** Purpose : compute total heat flux and total fwf due to ice shelf melt
+      !!
+      !! ** Method  : For each case (parametrisation or explicity cavity) :
+      !!              - define the before fields
+      !!              - compute top boundary layer properties
+      !!                (in case of parametrisation, this is the 
+      !!                 depth range model array used to compute mean far fields properties)
+      !!              - compute fluxes
+      !!              - write restart variables
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      IF ( ln_isfcav_mlt ) THEN
+         !
+         ! before time step 
+         IF ( kt /= nit000 ) THEN 
+            risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
+            fwfisf_cav_b(:,:)      = fwfisf_cav(:,:)
+         END IF
+         !
+         ! compute tbl lvl/h
+         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         !
+         ! compute ice shelf melt
+         CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav)
+         !
+         ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+         IF (lrst_oce) CALL isfrst_write(kt, 'cav', risf_cav_tsc, fwfisf_cav)
+         !
+      END IF
+      ! 
+      IF ( ln_isfpar_mlt ) THEN
+         !
+         ! before time step 
+         IF ( kt /= nit000 ) THEN 
+            risf_par_tsc_b (:,:) = risf_par_tsc (:,:)
+            fwfisf_par_b(:,:)    = fwfisf_par(:,:)
+         END IF
+         !
+         ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         !
+         ! compute ice shelf melt
+         CALL isf_par( kt, risf_par_tsc, fwfisf_par)
+         !
+         ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+         IF (lrst_oce) CALL isfrst_write(kt, 'par', risf_par_tsc, fwfisf_par)
+         !
+      END IF
+      !
+   END SUBROUTINE isf_stp
+
+   SUBROUTINE isf_stp_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfstp_init  ***
+      !!
+      !! ** Purpose :   Initialisation of the ice shelf public variables
+      !!
+      !! ** Method  :   Read the namsbc namelist and set derived parameters
+      !!                Call init routines for all other SBC modules that have one
+      !!
+      !! ** Action  : - read namsbc parameters
+      !!              - allocate memory
+      !!              - call cav/param init routine
+      !!----------------------------------------------------------------------
+      INTEGER               :: inum, ierror
+      INTEGER               :: ios                  ! Local integer output status for namelist read
+      INTEGER               :: ikt, ikb
+      INTEGER               :: ji, jj
+      !!----------------------------------------------------------------------
+      NAMELIST/namisf/ ln_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
+      !!----------------------------------------------------------------------
+      !
+      ! Allocate public array
+      CALL isf_alloc()
+      !
+      riceload(:,:)       = 0.0_wp
+      fwfisf_cpl(:,:)     = 0.0_wp
+      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
+      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
+      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
+      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
+      !
+      ! terminate routine now if no ice shelf melt formulation specify
+      IF ( .NOT. ln_isf ) RETURN
+      !
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namisf )
+      !
+      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,*) ''
+      !
+      ! sanity check
+      ! melt in the cavity without cavity
+      IF ( ln_isfcav_mlt .AND. (.NOT. ln_isfcav) ) &
+         &   CALL ctl_stop('ice shelf melt in the cavity activated (ln_isfcav_mlt) but no cavity detected in domcfg (ln_isfcav), STOP' )
+      !
+      IF ( ln_cpl ) 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' )
+         !
+         ! oasis melt computation not tested (coded but not tested)
+         IF ( ln_isfcav_mlt .OR. ln_isfpar_mlt ) THEN
+            IF ( TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis not tested' )
+            IF ( TRIM(cn_isfpar_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis not tested' )
+         END IF
+         !
+         ! oasis melt computation with cavity open and cavity parametrised (not coded)
+         IF ( ln_isfcav_mlt .AND. ln_isfpar_mlt ) THEN
+            IF ( TRIM(cn_isfpar_mlt) == 'oasis' .AND. TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis and cn_isfcav_mlt = oasis not coded' )
+         END IF
+      END IF
+      !
+      ! initialisation ice shelf load
+      IF ( ln_isfcav ) THEN
+         !
+         ! compute ice shelf mask
+         mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
+         !
+         ! compute ice shelf load
+         CALL isf_load( risfload )
+         !
+      END IF
+      !
+      ! initialisation useful variable
+      r1_Lfusisf =  1._wp / rLfusisf
+      !
+      ! initialisation melt in the cavity
+      IF ( ln_isfcav_mlt ) THEN
+         !
+         ! initialisation  of cav variable
+         CALL isf_cav_init()
+         !
+         ! read cav variable from restart
+         IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
+         !
+      END IF
+      !
+      ! initialisation parametrised melt
+      IF ( ln_isfpar_mlt ) THEN
+         !
+         ! initialisation  of par variable
+         CALL isf_par_init()
+         !
+         ! read par variable from restart
+         IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
+         !
+      END IF
+      !
+  END SUBROUTINE isf_stp_init
+   !!======================================================================
+END MODULE isfstp
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 11403)
@@ -1,32 +1,43 @@
 MODULE isftbl
-
-USE in_out_manager
-USE dom_oce
-USE oce
-USE isf
-USE lbclnk
-
-IMPLICIT NONE
-
-PRIVATE
-
-PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isftbl_ktop, isftbl_kbot
+   !!======================================================================
+   !!                       ***  MODULE  isftbl  ***
+   !! isftbl module :  compute properties of top boundary layer
+   !!======================================================================
+   !! History :  4.1  !  2019-09  (P. Mathiot) original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isftbl       : routine to compute : 
+   !!                  - geometry of the ice shelf tbl (isf_tbl_lvl, isftbl_ktop, isftbl_kbot)
+   !!                    (top and bottom level, thickness and fraction of deepest level affected)
+   !!                  - tbl averaged properties (isf_tbl, isf_tbl_avg)
+   !!----------------------------------------------------------------------
+
+   USE dom_oce ! vertical scale factor
+   USE lbclnk  ! lbc_lnk subroutine
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isftbl_ktop, isftbl_kbot
 
 CONTAINS
 
    SUBROUTINE isf_tbl( pvarin, pvarout, cd_ptin, ktop, kbot, phtbl, pfrac )
-      !!----------------------------------------------------------------------
-      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
+      !!--------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_tbl  ***
       !!
       !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
       !!
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pvarout ! 2d average of pvarin
-      !!-------------------------- IN  -------------------------------------
-      CHARACTER(len=1),                       INTENT(in   ) :: cd_ptin ! point of variable in/out
-      INTEGER,  DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: ktop , kbot
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: phtbl, pfrac
-      REAL(wp), DIMENSION(jpi,jpj,jpk),       INTENT(in   ) :: pvarin  ! 3d variable to average over the tbl
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) :: pvarout ! 2d average of pvarin
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(len=1)                      , INTENT(in   ) :: cd_ptin       ! point of variable in/out
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      , INTENT(in   ) :: pvarin        ! 3d variable to average over the tbl
+      !!-------------------------- IN OPTIONAL -----------------------------
+      INTEGER,  DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: ktop , kbot   ! top and bottom level
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: phtbl, pfrac  ! tbl thickness and fraction of bottom cell affected
       !!--------------------------------------------------------------------
       INTEGER ::   ji, jj                   ! loop index
@@ -34,5 +45,5 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zhtbl ! thickness of the tbl
       REAL(wp), DIMENSION(jpi,jpj) :: zfrac ! thickness of the tbl
-      !!----------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       ! 
       SELECT CASE ( cd_ptin )
@@ -86,21 +97,24 @@
 
    SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !!                  ***  ROUTINE isf_tbl_lvl  ***
       !!
       !! ** Purpose : compute mean property in the boundary layer
       !!
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
-      !!-------------------------- IN  -------------------------------------
-      INTEGER,  DIMENSION(jpi,jpj), INTENT(in) :: ktop, kbot   ! top and bottom level of the top boundary layer
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac !! fraction of bottom level of the tbl to be affected by the tbl
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pe3      ! vertical scale factor
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pvarin   ! tbl property to average between ktop, kbot over phtbl
-      !!---------------------------------------------------------------------
+      !! ** Method  : Depth average is made between the top level ktop and the bottom level kbot
+      !!              over a thickness phtbl. The bottom level is partially counted (pfrac).
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pvarout      ! tbl property averaged over phtbl between level ktop and kbot
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop, kbot   ! top and bottom level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3          ! vertical scale factor
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvarin       ! tbl property to average between ktop, kbot over phtbl
+      !!--------------------------------------------------------------------
       INTEGER  :: ji,jj,jk                    ! loop indices
       INTEGER  :: ikt, ikb                    ! top and bottom levels
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !
       ! compute tbl top.bottom level and thickness
@@ -123,5 +137,5 @@
 
    SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !!                  ***  ROUTINE isf_tbl_lvl  ***
       !!
@@ -129,13 +143,13 @@
       !!              - thickness of the top boundary layer
       !!
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      INTEGER,  DIMENSION(jpi,jpj), INTENT(inout) :: kbot         ! bottom level of the top boundary layer
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtbl        ! top boundary layer thickness
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pfrac        ! top boundary layer thickness
-      !!-------------------------- IN  -------------------------------------
-      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in) :: ktop       ! top level of the top boundary layer
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in) :: phw        ! water column thickness
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pe3         ! vertical scale factor
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: phtbl  ! top boundary layer thickness
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pfrac  ! top boundary layer thickness
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phw    ! water column thickness
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
       !!---------------------------------------------------------------------
       INTEGER :: ji,jj,jk
@@ -178,18 +192,25 @@
    !
    SUBROUTINE isftbl_kbot(ktop, phtbl, pe3, kbot)
-      !!-------------------------- OUT -------------------------------------
-      INTEGER,  DIMENSION(jpi,jpj), INTENT(inout) :: kbot         ! bottom level of the top boundary layer
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl           ! top boundary layer thickness
-      INTEGER,  DIMENSION(jpi,jpj), INTENT(in) :: ktop            ! top level of the top boundary layer
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pe3         ! vertical scale factor
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_lvl  ***
+      !!
+      !! ** Purpose : compute bottom level of the isf top boundary layer
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl  ! top boundary layer thickness
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
+      !!--------------------------------------------------------------------
       INTEGER :: ji, jj
       INTEGER :: ikt, ikb
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !
       ! phtbl need to be bounded by water column thickness before
       ! test: if phtbl = water column thickness, should return mbathy
       ! test: if phtbl = 0 should return ktop
+      !
       ! get ktbl
       DO jj = 1,jpj
@@ -208,22 +229,24 @@
       !
    SUBROUTINE isftbl_ktop(pdep, ktop)
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !!                  ***  ROUTINE isf_tbl_lvl  ***
       !!
       !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
       !!
-      !!----------------------------------------------------------------------
-      !!-------------------------- OUT -------------------------------------
-      INTEGER,  DIMENSION(jpi,jpj), INTENT(inout) :: ktop         ! top level affected by the ice shelf parametrisation
-      !!-------------------------- IN  -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdep            ! top depth of the parametrisation influence
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj), INTENT(  out) :: ktop        ! top level affected by the ice shelf parametrisation
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pdep        ! top depth of the parametrisation influence
+      !!--------------------------------------------------------------------
       INTEGER :: ji,jj
       INTEGER :: ikt
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
+      !
       ! compute top level (need to be recomputed each time (z*, z~) 
       ! be sure pdep is already correctly bounded
       ! test: this routine run on isfdraft should return mikt
       ! test: this routine run with pdep = 0 should return 1
+      !
       DO ji = 1, jpi
          DO jj = 1, jpj
@@ -234,5 +257,5 @@
       END DO
       !
-   END SUBROUTINE
+   END SUBROUTINE isftbl_ktop
 
 END MODULE isftbl
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 11403)
@@ -1,35 +1,45 @@
 MODULE isfutils
+   !!======================================================================
+   !!                       ***  MODULE  isfutils  ***
+   !! istutils module : miscelenious useful routines
+   !!======================================================================
+   !! History :  4.1  !  2019-09  (P. Mathiot) original code
+   !!----------------------------------------------------------------------
 
-USE oce
-USE in_out_manager
-USE iom
-USE lib_fortran
-USE lib_mpp
+   !!----------------------------------------------------------------------
+   !!   isfutils       : - read_2dcstdta to read a constant input file with iom_get
+   !!                    - debug to print array sum, min, max in ocean.output
+   !!----------------------------------------------------------------------
 
-IMPLICIT NONE
+   USE iom
+   USE lib_fortran
+   USE lib_mpp
 
-PRIVATE
+   IMPLICIT NONE
 
-INTERFACE isf_debug
-   MODULE PROCEDURE isf_debug2d, isf_debug3d
-END INTERFACE isf_debug
+   PRIVATE
 
-PUBLIC read_2dcstdta, isf_debug
+   INTERFACE debug
+      MODULE PROCEDURE debug2d, debug3d
+   END INTERFACE debug
+
+   PUBLIC read_2dcstdta, debug
 
 CONTAINS
 
    SUBROUTINE read_2dcstdta(cdfile, cdvar, pvar)
-      !!---------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !!                  ***  ROUTINE read_2dcstdta  ***
       !!
       !! ** Purpose : read input file
       !!
-      !!----------------------------------------------------------------------
+      !!--------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pvar          ! output variable
       !!-------------------------- IN  -------------------------------------
-      CHARACTER(len=256), INTENT(in) :: cdfile, cdvar
+      CHARACTER(len=256)          , INTENT(in   ) :: cdfile, cdvar ! input file name and variable name
       !!--------------------------------------------------------------------
       INTEGER :: inum
+      !!--------------------------------------------------------------------
 
       CALL iom_open( cdfile, inum )
@@ -39,25 +49,52 @@
    END SUBROUTINE read_2dcstdta
 
-   SUBROUTINE isf_debug2d(cdtxt,pvar)
-      CHARACTER(LEN=256), INTENT(in) :: cdtxt
+   SUBROUTINE debug2d(cdtxt,pvar)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_debug2d  ***
+      !!
+      !! ** Purpose : add debug print
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(LEN=256)          , INTENT(in   ) :: cdtxt
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pvar
+      !!--------------------------------------------------------------------
       REAL(wp) :: zmin, zmax, zsum
-      zsum = glob_sum( 'isf_debug', pvar(:,:) )
-      zmin = MINVAL( pvar(:,:) ) ; CALL mpp_min( 'isf_debug', zmin ) ! min over the global domain
-      zmax = MAXVAL( pvar(:,:) ) ; CALL mpp_max( 'isf_debug', zmax ) ! max over the global domain
+      !!--------------------------------------------------------------------
+      !
+      zsum = glob_sum( 'debug', pvar(:,:) )
+      zmin = MINVAL( pvar(:,:) ) ; CALL mpp_min( 'debug', zmin ) ! min over the global domain
+      zmax = MAXVAL( pvar(:,:) ) ; CALL mpp_max( 'debug', zmax ) ! max over the global domain
+      !
       WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
+      !
       FLUSH(numout)
-   END SUBROUTINE isf_debug2d
+   END SUBROUTINE debug2d
 
-   SUBROUTINE isf_debug3d(cdtxt,pvar)
-      CHARACTER(LEN=256), INTENT(in) :: cdtxt
+   SUBROUTINE debug3d(cdtxt,pvar)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_debug3d  ***
+      !!
+      !! ** Purpose : add debug print
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(LEN=256)              , INTENT(in   ) :: cdtxt
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvar
+      !!--------------------------------------------------------------------
       REAL(wp) :: zmin, zmax, zsum
-      zsum = glob_sum( 'isf_debug', pvar(:,:) )
-      zmin = MINVAL( pvar(:,:) ) ; CALL mpp_min( 'isf_debug', zmin ) ! min over the global domain
-      zmax = MAXVAL( pvar(:,:) ) ; CALL mpp_max( 'isf_debug', zmax ) ! max over the global domain
+      !!--------------------------------------------------------------------
+      !
+      zsum = glob_sum( 'debug', pvar(:,:) )
+      zmin = MINVAL( pvar(:,:) ) ; CALL mpp_min( 'debug', zmin ) ! min over the global domain
+      zmax = MAXVAL( pvar(:,:) ) ; CALL mpp_max( 'debug', zmax ) ! max over the global domain
+      !
       WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
+      !
       FLUSH(numout)
-   END SUBROUTINE isf_debug3d
+      !
+   END SUBROUTINE debug3d
 
 END MODULE isfutils
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90	(revision 11403)
@@ -122,7 +122,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff        [Kg/m2/s]  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting [Kg/m2/s]  
-
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff                                 [Kg/m2/s]  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting                              [Kg/m2/s]  
    !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 11403)
@@ -36,5 +36,5 @@
    USE eosbn2         ! 
    USE sbcrnf  , ONLY : l_rnfcpl
-   USE isf     , ONLY : l_isfcpl, fwfisf_cav, fwfisf_par
+   USE isf     , ONLY : l_isfcpl, fwfisf_cpl
 #if defined key_cice
    USE ice_domain_size, only: ncat
@@ -1406,5 +1406,9 @@
              rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runfofs
          ENDIF
-         IF( srcv(jpr_isf)%laction )  fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)  
+         !
+         ! 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)  
+         END IF
         
          IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
@@ -1709,5 +1713,5 @@
       ENDIF
       IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
-        fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1)  
+        fwfisf_cpl(:,:) = - frcv(jpr_isf)%z3(:,:,1)  
       ENDIF
 
@@ -1748,5 +1752,5 @@
       ENDIF
       IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
-        fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1)
+        fwfisf_cpl(:,:) = - 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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traisf.F90	(revision 11403)
@@ -1,6 +1,6 @@
 MODULE traisf
    !!==============================================================================
-   !!                       ***  MODULE  trasbc  ***
-   !! Ocean active tracers:  surface boundary condition
+   !!                       ***  MODULE  traisf  ***
+   !! Ocean active tracers:  ice shelf boundary condition
    !!==============================================================================
    !! History :    4.0  !  2019-09  (P. Mathiot) original file
@@ -38,22 +38,9 @@
    SUBROUTINE tra_isf ( kt )
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_sbc  ***
+      !!                  ***  ROUTINE tra_isf  ***
       !!                   
-      !! ** Purpose :   Compute the tracer surface boundary condition trend of
-      !!      (flux through the interface, concentration/dilution effect)
-      !!      and add it to the general trend of tracer equations.
+      !! ** Purpose :  Compute the temperature trend due to the ice shelf melting (qhoce + qhc)
       !!
-      !! ** Method :   The (air+ice)-sea flux has two components: 
-      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 
-      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice. 
-      !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe,
-      !!             they are simply added to the tracer trend (tsa).
-      !!               In linear free surface case (ln_linssh=T), the volume of the
-      !!             ocean does not change with the water exchanges at the (air+ice)-sea
-      !!             interface. Therefore another term has to be added, to mimic the
-      !!             concentration/dilution effect associated with water exchanges.
-      !!
-      !! ** Action  : - Update tsa with the surface boundary condition trend 
-      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
+      !! ** Action  : - update tsa for cav, par and cpl case
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
@@ -62,11 +49,11 @@
       IF( ln_timing )   CALL timing_start('tra_sbc')
       !
-      !----------------------------------------
-      !       Ice Shelf effects (ISF)
-      !     tbl treated as in Losh (2008) JGR
-      !----------------------------------------
+      ! 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
       !CALL tra_isf_cpl()
       !
@@ -77,4 +64,11 @@
    SUBROUTINE tra_isf_mlt(ktop, kbot, phtbl, pfrac, ptsc, ptsc_b, pts)
       !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_isf_mlt  ***
+      !!
+      !! *** Purpose :  Compute the temperature trend due to the ice shelf melting (qhoce + qhc) for cav or par case
+      !!
+      !! *** Action :: Update tsa with the surface boundary condition trend 
+      !!
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts
       !!----------------------------------------------------------------------
@@ -82,17 +76,11 @@
       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   
+      INTEGER                      :: ikt, ikb  ! top and bottom level of the tbl
+      REAL(wp), DIMENSION(jpi,jpj) :: ztc       ! total ice shelf tracer trend
       !!----------------------------------------------------------------------
       !
-      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices  
-      INTEGER  ::   ikt, ikb                    ! local integers
-      REAL(wp) :: zeps = 1.e-20
-      REAL(wp), DIMENSION(jpi,jpj) :: ztc
-      !!----------------------------------------------------------------------
-      !
-      !----------------------------------------
-      !       Ice Shelf effects (ISF)
-      !     tbl treated as in Losh (2008) JGR
-      !----------------------------------------
-      ! 
       ! compute 2d total trend due to isf
       ztc(:,:) = 0.5_wp * ( ptsc(:,:,jp_tem) + ptsc_b(:,:,jp_tem) ) / phtbl(:,:)
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 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90	(revision 11403)
@@ -316,4 +316,6 @@
                      ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
                   ENDIF
+                  !
+                  ! river runoff
                   IF( ln_rnf_depth ) THEN
                      ! Rivers are not just at the surface must go down to nk_rnf(ji,jj)
@@ -327,15 +329,13 @@
                      ENDIF
                   ENDIF
-
+                  !
+                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
+                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
+                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
                   !
                   ! solar penetration (temperature only)
                   IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
                      &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
-                     !
-                  ! river runoff
-                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
-                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
-                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
-                     !
+                  !
                   ! ice shelf
                   IF( ll_isf ) THEN
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90	(revision 11403)
@@ -60,5 +60,5 @@
    USE diacfl         ! CFL diagnostics               (dia_cfl_init routine)
    USE step           ! NEMO time-stepping                 (stp     routine)
-   USE isfmlt         ! ice shelf                     (isf_mlt_init routine)
+   USE isfstp         ! ice shelf                     (isf_stp_init routine)
    USE icbini         ! handle bergs, initialisation
    USE icbstp         ! handle bergs, calving, themodynamics and transport
@@ -453,5 +453,5 @@
       !                                      ! Active tracers
       IF( ln_traqsr    )   CALL tra_qsr_init      ! penetrative solar radiation qsr
-      IF (ln_isf       )   CALL isf_mlt_init      ! ice shelf
+                           CALL isf_stp_init      ! ice shelf
                            CALL tra_bbc_init      ! bottom heat flux
                            CALL tra_bbl_init      ! advective (and/or diffusive) bottom boundary layer scheme
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90	(revision 11403)
@@ -113,5 +113,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     )   CALL isf_mlt ( kstp )
+      IF( ln_isf     )   CALL isf_stp ( kstp )
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
 
Index: /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step_oce.F90	(revision 11402)
+++ /NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step_oce.F90	(revision 11403)
@@ -22,5 +22,5 @@
    USE sbcwave         ! Wave intialisation
 
-   USE isfmlt          ! ice shelf boundary condition     (isf_mlt routine)
+   USE isfstp          ! ice shelf boundary condition     (isf_stp routine)
 
    USE traqsr          ! solar radiation penetration      (tra_qsr routine)
