Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5462)
@@ -142,7 +142,12 @@
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-! IF should be useless check zpshde (PM)
-               IF ( mbku(ji,jj) > 1 ) zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
-               IF ( mbkv(ji,jj) > 1 ) zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
+                  zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
+                  zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
                IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 
                IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj)
@@ -151,6 +156,7 @@
          ENDIF
          !
-         zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
-         DO jk = 1, jpkm1
+         !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
+         ! interior value
+         DO jk = 2, jpkm1
             !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
             !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
@@ -162,9 +168,13 @@
          END DO
          ! surface initialisation 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-              zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp
-            END DO
-         END DO
+         zdzr(:,:,1) = 0._wp 
+         IF ( ln_isfcav ) THEN
+            ! if isf need to overwrite the interior value at at the first ocean point
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp
+               END DO
+            END DO
+         END IF
          !
          !                          !==   Slopes just below the mixed layer   ==!
@@ -175,14 +185,23 @@
          ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
          !
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji  ,jj)
-               IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji+1,jj)
-               IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj), hmlpt(ji+1,jj))
-               IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji  ,jj)
-               IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji,jj+1)
-               IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji,jj+1))
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
+                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp)
+                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp)
+                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
+                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp)
+                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp)
+               ENDDO
             ENDDO
-         ENDDO
+         ELSE
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
+                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
+               ENDDO
+            ENDDO
+         END IF
          DO jk = 2, jpkm1                            !* Slopes at u and v points
             DO jj = 2, jpjm1
@@ -198,21 +217,18 @@
                   zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  )
                   !                                      ! uslp and vslp output in zwz and zww, resp.
-                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
-                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
+                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj  ,jk) )
+                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) )
                   ! thickness of water column between surface and level k at u/v point
-                  zdepu = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                   &
-                             - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) )  &
-                             - fse3u(ji,jj,miku(ji,jj))                                         )
-                  zdepv = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                   &
-                             - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) &
-                             - fse3v(ji,jj,mikv(ji,jj))                                         )
-                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                              &
-                     &                 + zfi  * uslpml(ji,jj)                                                     &
-                     &                        * zdepu / MAX( zhmlpu(ji,jj), 5._wp )
-                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * umask(ji,jj,jk) * umask(ji,jj,jk-1)
-                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                              &
-                     &                 + zfj  * vslpml(ji,jj)                                                     &
-                     &                        * zdepv / MAX( zhmlpv(ji,jj), 5._wp ) 
-                  zww(ji,jj,jk) = zww(ji,jj,jk) * vmask(ji,jj,jk) * vmask(ji,jj,jk-1)
+                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              &
+                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) )
+                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              &
+                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) )
+                  !
+                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          &
+                     &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj)
+                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk)
+                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          &
+                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj) 
+                  zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk)
                   
                  
@@ -266,8 +282,8 @@
                   uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
                      &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   &
-                     &                            *   umask(ji,jj,jk-1) !* umask(ji,jj,jk) * umask(ji,jj,jk+1)
+                     &                            *   umask(ji,jj,jk-1)
                   vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   &
                      &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   &
-                     &                            *   vmask(ji,jj,jk-1) !* vmask(ji,jj,jk) * vmask(ji,jj,jk+1)
+                     &                            *   vmask(ji,jj,jk-1)
                END DO
             END DO
@@ -282,5 +298,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   !                                  !* Local vertical density gradient evaluated from N^2
-                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk)
                   !                                  !* Slopes at w point
                   !                                        ! i- & j-gradient of density at w-points
@@ -298,10 +314,10 @@
                   zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  )
                   !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
-                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )    ! zfk=1 in the ML otherwise zfk=0
+                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
                   zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp )
                   zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) &
-                     &            + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+                     &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
                   zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) &
-                     &            + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+                     &            + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
 
 !!gm  modif to suppress omlmask....  (as in Griffies operator)
@@ -356,6 +372,6 @@
                   zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   &
                      &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
-                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk)
-                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk)
+                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk)
                END DO
             END DO
@@ -423,7 +439,7 @@
                   vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) &
-                    &                              * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5 
+                    &                              * wmask(ji,jj,jk) * 0.5 
                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) &
-                    &                              * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5 
+                    &                              * wmask(ji,jj,jk) * 0.5 
                END DO 
             END DO 
@@ -736,6 +752,8 @@
             DO ji = 1, jpi
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
-               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
+               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN
+                  omlmask(ji,jj,jk) = 1._wp
+               ELSE
+                  omlmask(ji,jj,jk) = 0._wp
                ENDIF
             END DO
@@ -794,6 +812,6 @@
             zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  )
             !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
-            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
-            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
+            wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)
+            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)
          END DO
       END DO
