Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/domain.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 4812)
@@ -335,5 +335,5 @@
       INTEGER ::   ios
       !!
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
       !!----------------------------------------------------------------------
 
@@ -352,7 +352,8 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
-         WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco
-         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps
-         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco = ', ln_sco
+         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
+         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
+         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
+         WRITE(numout,*) '             ice shelf cavity               ln_isfcav = ', ln_isfcav
       ENDIF
 
@@ -361,5 +362,7 @@
       IF( ln_zps ) ioptio = ioptio + 1
       IF( ln_sco ) ioptio = ioptio + 1
-      IF ( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
+      IF( ln_isfcav ) ioptio = 33
+      IF ( ioptio /= 1  )   CALL ctl_stop( ' none or several vertical coordinate options used' )
+      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
 
    END SUBROUTINE dom_zgr
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 4812)
@@ -164,5 +164,6 @@
 
          CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
-         mbathy(:,:) = INT( zmbk(:,:) )
+         mbathy (:,:) = INT( zmbk(:,:) )
+         misfdep(:,:) = 1                                               ! ice shelf case not yet done
          
          CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
@@ -371,4 +372,5 @@
       !
       mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
+      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
       !                                     ! bottom k-index of W-level = mbkt+1
       DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 4812)
@@ -537,6 +537,11 @@
       CALL eos( pts, rhd, rhop, gdept_0(:,:,:) )   ! Time-filtered in situ density 
       CALL bn2( pts, rn2 )         ! before Brunt-Vaisala frequency
-      IF( ln_zps )   &
-         &  CALL zps_hde( kt, jpts, pts, gtsu, gtsv, rhd, gru, grv )  ! Partial steps: before Horizontal DErivative
+      IF( ln_zps )   & ! Partial steps: before Horizontal DErivative
+        &    CALL zps_hde( kt, jpts, pts, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient
+        &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             !
+        &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )  ! of t, s, rd at the last ocean level
+        ! only gtsu, gtsv, rhd, gru , grv are used 
+
+
          !                                                            ! of t, s, rd at the bottom ocean level
       CALL zdf_mxl( kt )            ! mixed layer depth
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 4812)
@@ -17,4 +17,5 @@
    USE timing          ! Timing
    USE oce             ! ocean dynamics and tracers 
+   USE sbcisf          ! ice shelf
    USE dom_oce         ! ocean space and time domain 
    USE phycst          ! physical constants
@@ -90,5 +91,5 @@
       ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
       ! -----------------------------------------------------------------------
-      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
+      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+rdivisf*fwfisf(:,:) ) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
       IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4812)
@@ -114,6 +114,6 @@
                ikbu = miku(ji,jj)       ! ocean top level at u- and v-points 
                ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
-               IF (ikbu .GE. 2) avmu(ji,jj,ikbu-1) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu-1)
-               IF (ikbv .GE. 2) avmv(ji,jj,ikbv-1) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv-1)
+               IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
+               IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
             END DO
          END DO
@@ -176,5 +176,5 @@
                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
-               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
+               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
@@ -184,6 +184,6 @@
       DO jj = 2, jpjm1        ! Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwi(ji,jj,miku(ji,jj)) = 0._wp
-            zwd(ji,jj,miku(ji,jj)) = 1._wp - zws(ji,jj,miku(ji,jj))
+            zwi(ji,jj,1) = 0._wp
+            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
          END DO
       END DO
@@ -275,6 +275,6 @@
       DO jj = 2, jpjm1        ! Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwi(ji,jj,mikv(ji,jj)) = 0._wp
-            zwd(ji,jj,mikv(ji,jj)) = 1._wp - zws(ji,jj,mikv(ji,jj))
+            zwi(ji,jj,1) = 0._wp
+            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
          END DO
       END DO
@@ -306,5 +306,5 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
+            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl   * fse3v_a(ji,jj,mikv(ji,jj)) 
 #if defined key_dynspg_ts            
             va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
@@ -349,20 +349,20 @@
       IF( ln_bfrimp ) THEN
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+        DO jj = 1, 1
+           DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
 # endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            avmu(ji,jj,ikbu+1) = 0.e0
-            avmv(ji,jj,ikbv+1) = 0.e0
-            ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
-            ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
-            avmu(ji,jj,ikbu-1) = 0.e0
-            avmv(ji,jj,ikbv-1) = 0.e0
-         END DO
-      END DO
+              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+              avmu(ji,jj,ikbu+1) = 0.e0
+              avmv(ji,jj,ikbv+1) = 0.e0
+              ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+              ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
+              IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
+              IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
+           END DO
+        END DO
       ENDIF
       !
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4812)
@@ -108,7 +108,9 @@
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
       REAL(wp) ::   zck, zfk,      zbw             !   -      -
+      REAL(wp) ::   zdepv, zdepu         !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zgru, zgrv
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhmlpu, zhmlpv
       !!----------------------------------------------------------------------
       !
@@ -116,4 +118,5 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
+      CALL wrk_alloc( jpi,jpj, zhmlpu, zhmlpv )
 
       IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN 
@@ -153,5 +156,5 @@
          !
          zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
-         DO jk = 2, jpkm1
+         DO jk = 1, jpkm1
             !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
             !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
@@ -162,4 +165,10 @@
                &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
          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
          !
          !                          !==   Slopes just below the mixed layer   ==!
@@ -170,4 +179,14 @@
          ! ===========================      | 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))
+            ENDDO
+         ENDDO
          DO jk = 2, jpkm1                            !* Slopes at u and v points
             DO jj = 2, jpjm1
@@ -183,14 +202,23 @@
                   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) )
-                  zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
-                     &                   + zfi  * uslpml(ji,jj)                                                     &
-                     &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
-                     &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk)
-                  zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
-                     &                   + zfj  * vslpml(ji,jj)                                                     &
-                     &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
-                     &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,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)
+                  
+                 
 !!gm  modif to suppress omlmask.... (as in Griffies case)
 !                  !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
@@ -242,8 +270,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   &
-                     &                            *   tmask(ji,jj,jk-1)
+                     &                            *   umask(ji,jj,jk-1) !* umask(ji,jj,jk) * 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   &
-                     &                            *   tmask(ji,jj,jk-1)
+                     &                            *   vmask(ji,jj,jk-1) !* vmask(ji,jj,jk) * vmask(ji,jj,jk+1)
                END DO
             END DO
@@ -258,5 +286,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. )
+                  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)
                   !                                  !* Slopes at w point
                   !                                        ! i- & j-gradient of density at w-points
@@ -266,7 +294,7 @@
                      &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj)
                   zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
-                     &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk)
+                     &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci
                   zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
-                     &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk)
+                     &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj
                   !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
                   !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
@@ -274,8 +302,8 @@
                   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
-                  zck = fsdepw(ji,jj,jk) / 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)
-                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
+                  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)
+                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
 
 !!gm  modif to suppress omlmask....  (as in Griffies operator)
@@ -330,6 +358,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)
-                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * tmask(ji,jj,jk-1)
+                  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)
                END DO
             END DO
@@ -387,6 +415,6 @@
                   uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
                   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) * 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) * 0.5 
+                  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 
+                  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 
                END DO 
             END DO 
@@ -410,4 +438,5 @@
       
       CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
+      CALL wrk_dealloc( jpi,jpj,     zhmlpu, zhmlpv)
       !
       IF( nn_timing == 1 )  CALL timing_stop('ldf_slp')
@@ -714,5 +743,5 @@
 # endif
                ik = nmln(ji,jj) - 1
-               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
+               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
                ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
                ENDIF
@@ -742,6 +771,6 @@
             !
             !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
-            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
-            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
+            iku = MIN(  MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
+            ikv = MIN(  MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
             zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
             zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 4812)
@@ -116,5 +116,5 @@
       CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj           ) 
       !
-      pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
+      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -144,5 +144,4 @@
                   ! gradient of  tracers
                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
-                  pgzu(ji,jj)    = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku)
                ELSE                           ! case 2
                   zmaxu = -ze3wu / fse3w(ji,jj,iku)
@@ -151,5 +150,4 @@
                   ! gradient of tracers
                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
-                  pgzu(ji,jj)    = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu)
                ENDIF
                ENDIF
@@ -163,5 +161,4 @@
                   ! gradient of tracers
                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
-                  pgzv(ji,jj)    = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 
                ELSE                           ! case 2
                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
@@ -170,5 +167,4 @@
                   ! gradient of tracers
                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
-                  pgzv(ji,jj)    = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv)
                ENDIF
                ENDIF
@@ -183,4 +179,5 @@
       ! horizontal derivative of density anomalies (rd)
       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
+         pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp 
 # if defined key_vectopt_loop
          jj = 1
@@ -224,4 +221,5 @@
                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
                IF( ze3wu >= 0._wp ) THEN 
+                  pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku)
                   pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
                   pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1 
@@ -230,4 +228,5 @@
                                    - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
                ELSE  
+                  pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu)
                   pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2
                   pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2
@@ -237,4 +236,5 @@
                ENDIF
                IF( ze3wv >= 0._wp ) THEN
+                  pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 
                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1
                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1
@@ -243,4 +243,5 @@
                                    - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
                ELSE 
+                  pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv)
                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2
                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2
@@ -284,12 +285,10 @@
                   ! gradient of tracers
                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
-                  sgzu(ji,jj)    = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)
                ELSE                           ! case 2
                   zmaxu = - ze3wu / fse3w(ji,jj,iku+1)
                   ! interpolated values of tracers
                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) )
+                  ! gradient of  tracers
                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
-                  ! gradient of  tracers
-                  sgzu(ji,jj)    = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
                ENDIF
                !
@@ -301,5 +300,4 @@
                   ! gradient of tracers
                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
-                  sgzv(ji,jj)    = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 
                ELSE                           ! case 2
                   zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1)
@@ -308,5 +306,4 @@
                   ! gradient of tracers
                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
-                  sgzv(ji,jj)    = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)
                ENDIF
 # if ! defined key_vectopt_loop
@@ -361,4 +358,5 @@
                ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
                IF( ze3wu >= 0._wp ) THEN
+                 sgzu  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)
                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1
                  smru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1 
@@ -367,4 +365,5 @@
                                    - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
                ELSE
+                 sgzu  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2
                  smru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
@@ -374,4 +373,5 @@
                ENDIF
                IF( ze3wv >= 0._wp ) THEN
+                 sgzv  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 
                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1
                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
@@ -381,4 +381,5 @@
                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco
                ELSE
+                 sgzv  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)
                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2
                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 4812)
@@ -138,5 +138,5 @@
                ENDIF
          ! mask zmsk in order to have avt and avs masked
-               zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
+               zmsks(ji,jj) = zmsks(ji,jj) * tmask(ji,jj,jk)
 
 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 4812)
@@ -99,6 +99,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rho_c )   nmln(ji,jj) = jk      ! Mixed layer
-               IF( avt (ji,jj,jk) < avt_c                     )   imld(ji,jj) = jk      ! Turbocline 
+               IF( rhop(ji,jj,jk) > rhop(ji,jj,MAX(mikt(ji,jj),nla10)) + rho_c )   nmln(ji,jj) = MAX(jk,mikt(ji,jj))      ! Mixed layer
+               IF( avt (ji,jj,jk) < avt_c                     )   imld(ji,jj) = MAX(mikt(ji,jj),jk)      ! Turbocline 
             END DO
          END DO
@@ -109,7 +109,7 @@
             iiki = imld(ji,jj)
             iikn = nmln(ji,jj)
-            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * tmask(ji,jj,1)    ! Turbocline depth 
-            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * tmask(ji,jj,1)    ! Mixed layer depth
-            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1)                     ! depth of the last T-point inside the mixed layer
+            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,mikt(ji,jj) )            ) * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX(mikt(ji,jj),nla10) ) ) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,mikt(ji,jj) )            ) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 4812)
@@ -73,16 +73,16 @@
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra        )  ! iso-level laplacian
       CASE ( 1 )                                                                                            ! rotated laplacian
                        IF( ln_traldf_grif ) THEN
                           CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
                        ELSE
-                          CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+                          CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra, rn_ahtb_0 )
                        ENDIF
-      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra            )  ! iso-level bilaplacian
       CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, nittrc000, 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
          !
       CASE ( -1 )                                     ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf0 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
@@ -90,9 +90,9 @@
             CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
          ELSE
-            CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+            CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra, rn_ahtb_0 )
          ENDIF
          WRITE(charout, FMT="('ldf1 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf2 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 4812)
@@ -82,5 +82,5 @@
       IF( .NOT. Agrif_Root())   CALL Agrif_Update_Trc( kstp )   ! Update tracer at AGRIF zoom boundaries : children only
 #endif
-         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
+         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi )  ! Partial steps: now horizontal gradient of passive
                                                                 ! tracers at the bottom ocean level
          !
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 4812)
@@ -42,4 +42,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)           ::  gtru           !: hor. gradient at u-points at bottom ocean level
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)           ::  gtrv           !: hor. gradient at v-points at bottom ocean level
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)           ::  gtrui          !: hor. gradient at u-points at top    ocean level
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)           ::  gtrvi          !: hor. gradient at v-points at top    ocean level
    
    !! passive tracers  (input and output)
@@ -183,5 +185,6 @@
       !
       ALLOCATE( trn(jpi,jpj,jpk,jptra), trb(jpi,jpj,jpk,jptra), tra(jpi,jpj,jpk,jptra),       &  
-         &      gtru(jpi,jpj,jpk)     , gtrv(jpi,jpj,jpk)                             ,       &
+         &      gtru (jpi,jpj,jptra)  , gtrv (jpi,jpj,jptra)                          ,       &
+         &      gtrui(jpi,jpj,jptra)  , gtrvi(jpi,jpj,jptra)                          ,       &
          &      cvol(jpi,jpj,jpk)     , rdttrc(jpk)           , trai(jptra)           ,       &
          &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       & 
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 4811)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 4812)
@@ -143,7 +143,6 @@
  
       tra(:,:,:,:) = 0._wp
-      
       IF( ln_zps .AND. .NOT. lk_c1d )   &              ! Partial steps: before horizontal gradient of passive
-        &    CALL zps_hde( nit000, jptra, trn, gtru, gtrv )       ! tracers at the bottom ocean level
+        &    CALL zps_hde( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi )       ! tracers at the bottom ocean level
 
       !
