Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA/diahth.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA/diahth.F90	(revision 15214)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA/diahth.F90	(revision 15231)
@@ -365,11 +365,11 @@
       !
       ilevel(:,:) = 1
-      DO jk = 2, jpkm1
+      DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( ( gdept_n(ji,jj,jk) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
-                   ilevel(ji,jj) = jk
-                   zthick(ji,jj) = zthick(ji,jj) + e3t_n(ji,jj,jk)
-                   phtc  (ji,jj) = phtc  (ji,jj) + e3t_n(ji,jj,jk) * ptn(ji,jj,jk)
+               IF( ( gdepw_n(ji,jj,jk+1) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
+                  ilevel(ji,jj) = jk+1
+                  zthick(ji,jj) = zthick(ji,jj) + e3t_n(ji,jj,jk)
+                  phtc  (ji,jj) = phtc  (ji,jj) + e3t_n(ji,jj,jk) * ptn(ji,jj,jk)
                ENDIF
             ENDDO
@@ -380,7 +380,8 @@
          DO ji = 1, jpi
             ik = ilevel(ji,jj)
-            zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep
-            phtc(ji,jj)   = phtc(ji,jj) + ptn(ji,jj,ik+1) * MIN( e3t_n(ji,jj,ik+1), zthick(ji,jj) ) &
-                                                          * tmask(ji,jj,ik+1)
+            IF( tmask(ji,jj,ik) == 1 ) THEN
+               zthick(ji,jj) = MIN ( gdepw_n(ji,jj,ik+1), pdep ) - zthick(ji,jj)   ! remaining thickness to reach dephw pdep
+               phtc(ji,jj)   = phtc(ji,jj) + ptn(ji,jj,ik) * zthick(ji,jj)
+            ENDIF
          END DO
       ENDDO
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/FLO/flodom.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/FLO/flodom.F90	(revision 15214)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/FLO/flodom.F90	(revision 15231)
@@ -366,11 +366,4 @@
       REAL(wp) ::   zabt, zbct, zcdt, zdat, zabpt, zbcpt, zcdpt, zdapt
       !!---------------------------------------------------------------------
-      !! Statement function
-      REAL(wp) ::   fsline
-      REAL(wp) ::   psax, psay, psbx, psby, psx, psy
-      fsline( psax, psay, psbx, psby, psx, psy ) = psy  * ( psbx - psax )   &
-         &                                       - psx  * ( psby - psay )   &
-         &                                       + psax *   psby - psay * psbx
-      !!---------------------------------------------------------------------
       
       ! 4 semi plane defined by the 4 points and including the T point
@@ -407,4 +400,19 @@
    END SUBROUTINE flo_findmesh
 
+   FUNCTION fsline( psax, psay, psbx, psby, psx, psy )
+      !! ---------------------------------------------------------------------
+      !!                 ***  Function fsline  ***
+      !!          
+      !! ** Purpose :
+      !! ** Method  : 
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   fsline
+      REAL(wp), INTENT(in) ::   psax, psay, psbx, psby, psx, psy
+      !!---------------------------------------------------------------------
+      fsline = psy  * ( psbx - psax )   &
+         &   - psx  * ( psby - psay )   &
+         &   + psax *   psby - psay * psbx
+      !
+   END FUNCTION fsline
 
    FUNCTION flo_dstnce( pla1, phi1, pla2, phi2 )
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcisf.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcisf.F90	(revision 15214)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcisf.F90	(revision 15231)
@@ -354,6 +354,5 @@
             DO jj = 1, jpj
                 ik = 2
-!!gm potential bug: use gdepw_0 not _n
-                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
+                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_0(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
                 misfkt(ji,jj) = ik-1
             END DO
