Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynzdf.F90	(revision 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynzdf.F90	(revision 11844)
@@ -428,5 +428,5 @@
                   ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
-                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
+                  zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
                END DO
             END DO
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 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isf.F90	(revision 11844)
@@ -31,4 +31,5 @@
    CHARACTER(LEN=256), PUBLIC :: cn_isfdir
    LOGICAL           , PUBLIC :: ln_isf
+   LOGICAL           , PUBLIC :: ln_isfdebug=.TRUE.
    !
    ! 0.2 -------- ice shelf cavity opened namelist parameter -------------
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 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 11844)
@@ -15,4 +15,5 @@
    USE oce            ! ocean dynamics and tracers
    USE isf            ! ice shelf public variables
+   USE isfutils
    USE isftbl         ! ice shelf top boundary layer properties
    USE isfcavmlt      ! ice shelf melt formulation
@@ -139,4 +140,10 @@
       ptsc(:,:,jp_tem) = - zqh(:,:) * r1_rau0_rcp
       !
+      IF ( ln_isfdebug ) THEN
+         CALL debug('isf_cav: ptsc T',ptsc(:,:,1))
+         CALL debug('isf_cav: ptsc S',ptsc(:,:,2))
+         CALL debug('isf_cav: pqfwf fwf',pqfwf(:,:))
+      END IF
+      !
    END SUBROUTINE isf_cav
 
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 11844)
@@ -12,4 +12,5 @@
    USE oce            ! ocean dynamics and tracers
    USE isf
+   USE isfutils
    USE isftbl
    USE dom_oce        ! ocean space and time domain
@@ -91,4 +92,9 @@
       CALL iom_put('isfgammas', pgs(:,:))
       !
+      IF (ln_isfdebug) THEN
+         CALL debug( 'isfcav_gam pgt:', pgt(:,:) )
+         CALL debug( 'isfcav_gam pgs:', pgs(:,:) )
+      END IF
+      !
    END SUBROUTINE isfcav_gammats
    !
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 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 11844)
@@ -13,4 +13,5 @@
    USE oce            ! ocean dynamics and tracers
    USE isf            ! ice shelf public variables
+   USE isfutils       ! ice shelf debug subroutine
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -74,4 +75,10 @@
       END SELECT
       !
+      IF (ln_isfdebug) THEN
+         CALL debug( 'isfcav_mlt:', pqhc (:,:) )
+         CALL debug( 'isfcav_mlt:', pqoce(:,:) )
+         CALL debug( 'isfcav_mlt:', pqfwf(:,:) )
+      END IF
+      !
    END SUBROUTINE isfcav_mlt
 
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 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 11844)
@@ -16,4 +16,6 @@
    USE dom_oce ! vertical scale factor
    USE lbclnk  ! lbc_lnk subroutine
+   USE isfutils
+   USE isf
 
    IMPLICIT NONE
@@ -48,8 +50,9 @@
       REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: pfrac         ! fraction of bottom cell affected by tbl
       !!--------------------------------------------------------------------
-      INTEGER ::   ji, jj                   ! loop index
-      INTEGER , DIMENSION(jpi,jpj) :: ikbot ! bottom level of the tbl
-      REAL(wp), DIMENSION(jpi,jpj) :: zhtbl ! thickness of the tbl
-      REAL(wp), DIMENSION(jpi,jpj) :: zfrac ! thickness of the tbl
+      INTEGER ::   ji, jj                     ! loop index
+      INTEGER , DIMENSION(jpi,jpj) :: ikbot   ! bottom level of the tbl
+      REAL(wp), DIMENSION(jpi,jpj) :: zvarout ! 2d average of pvarin
+      REAL(wp), DIMENSION(jpi,jpj) :: zhtbl   ! thickness of the tbl
+      REAL(wp), DIMENSION(jpi,jpj) :: zfrac   ! thickness of the tbl
       !!--------------------------------------------------------------------
       ! 
@@ -64,15 +67,17 @@
          !
          ! compute tbl property at U point
-         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u_n, pvarin, pvarout )
+         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u_n, pvarin, zvarout )
+         !
+         ! check if needed (probably yes)
+         CALL lbc_lnk('sbcisf', pvarout,'U',-1.)
          !
          ! compute tbl property at T point
+         pvarout(1,:) = 0._wp
          DO jj = 1, jpj
             DO ji = 2, jpi
-               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
+               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
             END DO
          END DO
-         !
-         ! check if needed (probably yes)
-         CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
+         ! lbclnk not needed as a final communication is done after the computation of fwf
          ! 
       CASE ( 'V' )
@@ -85,15 +90,17 @@
          !
          ! compute tbl property at V point
-         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v_n, pvarin, pvarout )
+         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v_n, pvarin, zvarout )
+         !
+         ! check if needed (probably yes)
+         CALL lbc_lnk('sbcisf', pvarout,'V',-1.)
          !
          ! pvarout is an averaging of wet point
+         pvarout(:,1) = 0._wp
          DO jj = 2, jpj
             DO ji = 1, jpi
-               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
+               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
             END DO
          END DO
-         !
-         ! check if needed (probably yes)
-         CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
+         ! lbclnk not needed as a final communication is done after the computation of fwf
          !
       CASE ( 'T' )
@@ -103,4 +110,6 @@
          !
       END SELECT
+      !
+      IF (ln_isfdebug) CALL debug('isf_tbl pvarout:',pvarout)
       !
    END SUBROUTINE isf_tbl
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 11823)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 11844)
@@ -65,11 +65,12 @@
       !!--------------------------------------------------------------------
       !
-      zsum = glob_sum_full( '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
+      zsum = glob_sum( 'debug', pvar(:,:) )
+      zmin = glob_min( 'debug', pvar(:,:) )
+      zmax = glob_max( 'debug', pvar(:,:) )
       !
-      WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
+      IF (lwp) WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
       !
-      FLUSH(numout)
+      CALL FLUSH(numout)
+      !
    END SUBROUTINE debug2d
 
@@ -89,11 +90,11 @@
       !!--------------------------------------------------------------------
       !
-      zsum = glob_sum_full( '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
+      zsum = glob_sum( 'debug', pvar(:,:,:) )
+      zmin = glob_min( 'debug', pvar(:,:,:) )
+      zmax = glob_max( 'debug', pvar(:,:,:) )
       !
-      WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
+      IF (lwp) WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
       !
-      FLUSH(numout)
+      CALL FLUSH(numout)
       !
    END SUBROUTINE debug3d
