Index: NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90	(revision 14356)
+++ NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90	(revision 14357)
@@ -148,5 +148,5 @@
          END_2D
       ENDIF
-      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
+      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the top ocean level
          DO_2D( 1, 0, 1, 0 )
             IF( miku(ji,jj) > 1 )   zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 
@@ -356,4 +356,15 @@
             END_2D
          ENDIF
+         IF( ln_zps .AND. ln_isfcav .AND. l_grad_zps) THEN           ! partial steps correction at the top ocean level
+            DO_2D( 1, 0, 1, 0 )
+               iku  = miku(ji,jj)          ;   ikv  = mikv(ji,jj)        ! first wet ocean level (u- & v-points)
+               zdit = gtui(ji,jj,jp_tem)   ;   zdjt = gtvi(ji,jj,jp_tem) ! i- & j-gradient of Temperature
+               zdis = gtui(ji,jj,jp_sal)   ;   zdjs = gtvi(ji,jj,jp_sal) ! i- & j-gradient of Salinity
+               zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj)
+               zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj)
+               zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
+               zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
+            END_2D
+         ENDIF
          !
       END DO
@@ -368,4 +379,11 @@
                zdks = 0._wp
             ENDIF
+            IF ( ln_isfcav ) THEN
+               IF( mikt(ji,jj) > 1 .AND.  jk+kp <= mikt(ji,jj) )THEN
+                  zdkt = 0._wp                                          ! 1st wet level gradient set to zero
+                  zdks = 0._wp
+               ENDIF
+            ENDIF
+
             zdzrho_raw = ( - rab_b(ji,jj,jk   ,jp_tem) * zdkt & 
                        &   + rab_b(ji,jj,jk   ,jp_sal) * zdks &
@@ -382,10 +400,7 @@
       !                                       !==  intialisations to zero  ==!
       !
-      wslp2  (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field set to zero
-      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
-      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp
-      !!gm _iso set to zero missing
-      triadi  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
-      triadj  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp
+      wslp2   (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field
+      triadi_g(:,:,:,:,:) = 0._wp ; triadj_g(:,:,:,:,:) = 0._wp
+      triadi  (:,:,:,:,:) = 0._wp ; triadj  (:,:,:,:,:) = 0._wp
 
       !-------------------------------------!
@@ -445,6 +460,7 @@
                   !
                   ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
-                  zti_coord = znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj)
-                  ztj_coord = znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked
+                  zti_coord =  wmask(ji,jj,jk+kp) * znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj)
+                  ztj_coord =  wmask(ji,jj,jk+kp) * znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked
+
                   zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces
                   ztj_g_raw = ztj_raw - ztj_coord
@@ -464,8 +480,8 @@
                   zti_g_lim =          ( zfacti   * zti_g_lim                       &
                      &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
-                     &                           * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
+                     &                           * (gdepw(ji+ip,jj,jk+kp,Kmm)-risfdep(ji+ip,jj)) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
                   ztj_g_lim =          ( zfactj   * ztj_g_lim                       &
                      &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
-                     &                           * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
+                     &                           * (gdepw(ji,jj+jp,jk+kp,Kmm)-risfdep(ji,jj+jp)) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
                   !
                   triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim
@@ -478,4 +494,7 @@
                   ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms
                   ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2
+                  !
+                  zti_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS(zti_lim ) ), zti_lim )
+                  ztj_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS(ztj_lim ) ), ztj_lim )
                   !
                   zti_lim  = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp)    ! remove coordinate slope => relative to coordinate surfaces
Index: NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldftra.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldftra.F90	(revision 14356)
+++ NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldftra.F90	(revision 14357)
@@ -248,6 +248,4 @@
       ENDIF
       !
-      IF( ln_isfcav .AND. ln_traldf_triad )   CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' )
-           !
       IF(  nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. &
          & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required
Index: NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/TRA/traldf_triad.F90	(revision 14356)
+++ NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/TRA/traldf_triad.F90	(revision 14357)
@@ -277,5 +277,5 @@
             !                    !==  Vertical tracer gradient at level jk and jk+1
             DO_2D( 1, 1, 1, 1 )
-               zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
+               zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)
             END_2D
             !
@@ -284,5 +284,5 @@
             ELSE
                DO_2D( 1, 1, 1, 1 )
-                  zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                  zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
                END_2D
             ENDIF
@@ -307,5 +307,5 @@
                         IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew
                         zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr
+                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) -                (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr
                      END_2D
                   END DO
@@ -327,5 +327,5 @@
                         IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew
                         zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr
+                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) -                (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr
                      END_2D
                   END DO
@@ -350,5 +350,5 @@
                         IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew
                         zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
+                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) -                (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr
                      END_2D
                   END DO
@@ -369,6 +369,6 @@
                         zah_slp = zah * zslope_iso
                         IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew
-                        zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
+                        zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
+                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) -                (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr
                      END_2D
                   END DO
@@ -388,5 +388,5 @@
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
             DO_3D( 0, 0, 1, 0, 2, jpkm1 )
-               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
+               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
                   &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
                   &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
@@ -396,10 +396,10 @@
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
                DO_3D( 0, 0, 1, 0, 2, jpkm1 )
-                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
+                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)             &
                      &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                END_3D
             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
                DO_3D( 0, 0, 1, 0, 2, jpkm1 )
-                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
+                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                      &
                      &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
                      &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
