Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/zpshde.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/zpshde.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/zpshde.F90	(revision 11949)
@@ -39,5 +39,5 @@
 CONTAINS
 
-   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
+   SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv,   &
       &                          prd, pgru, pgrv    )
       !!----------------------------------------------------------------------
@@ -85,4 +85,5 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  Kmm         ! ocean time level index
       INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
@@ -109,11 +110,11 @@
                iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
                ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-!!gm BUG ? when applied to before fields, e3w_b should be used....
-               ze3wu = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
-               ze3wv = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
+!!gm BUG ? when applied to before fields, e3w(:,:,:,Kbb) should be used....
+               ze3wu = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
+               ze3wv = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
                !
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku)
+                  zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
@@ -121,5 +122,5 @@
                   pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
-                  zmaxu = -ze3wu / e3w_n(ji,jj,iku)
+                  zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
@@ -130,5 +131,5 @@
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv)
+                  zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
@@ -136,5 +137,5 @@
                   pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
-                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv)
+                  zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
@@ -155,11 +156,11 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
-               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)     ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)     ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)     ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)     ! -     -      case 2
+               ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
+               ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)     ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)     ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)     ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)     ! -     -      case 2
                ENDIF
             END DO
@@ -173,6 +174,6 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
-               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
+               ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
+               ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
                ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
@@ -192,5 +193,5 @@
 
 
-   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
+   SUBROUTINE zps_hde_isf( kt, Kmm, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
       &                          prd, pgru, pgrv, pgrui, pgrvi )
       !!----------------------------------------------------------------------
@@ -241,4 +242,5 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  Kmm          ! ocean time level index
       INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields
@@ -270,10 +272,10 @@
                iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
                ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
-               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
+               ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
+               ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
                !
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku)
+                  zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
@@ -281,5 +283,5 @@
                   pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
-                  zmaxu = -ze3wu / e3w_n(ji,jj,iku)
+                  zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
                   ! interpolated values of tracers
                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
@@ -290,5 +292,5 @@
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv)
+                  zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
@@ -296,5 +298,5 @@
                   pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
-                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv)
+                  zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
                   ! interpolated values of tracers
                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
@@ -318,12 +320,12 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
-               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
-               !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
-               ENDIF
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
+               ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
+               ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
+               !
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
+               ENDIF
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
                ENDIF
 
@@ -340,6 +342,6 @@
                iku = mbku(ji,jj)
                ikv = mbkv(ji,jj)
-               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
-               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
+               ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
+               ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
 
                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
@@ -369,10 +371,10 @@
                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
-               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
+               ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
+               ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
 
                ! i- direction
                IF( ze3wu >= 0._wp ) THEN      ! case 1
-                  zmaxu = ze3wu / e3w_n(ji+1,jj,ikup1)
+                  zmaxu = ze3wu / e3w(ji+1,jj,ikup1,Kmm)
                   ! interpolated values of tracers
                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
@@ -380,5 +382,5 @@
                   pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                           ! case 2
-                  zmaxu = - ze3wu / e3w_n(ji,jj,ikup1)
+                  zmaxu = - ze3wu / e3w(ji,jj,ikup1,Kmm)
                   ! interpolated values of tracers
                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
@@ -389,5 +391,5 @@
                ! j- direction
                IF( ze3wv >= 0._wp ) THEN      ! case 1
-                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikvp1)
+                  zmaxv =  ze3wv / e3w(ji,jj+1,ikvp1,Kmm)
                   ! interpolated values of tracers
                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
@@ -395,5 +397,5 @@
                   pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                           ! case 2
-                  zmaxv =  - ze3wv / e3w_n(ji,jj,ikvp1)
+                  zmaxv =  - ze3wv / e3w(ji,jj,ikvp1,Kmm)
                   ! interpolated values of tracers
                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
@@ -416,13 +418,13 @@
                iku = miku(ji,jj)
                ikv = mikv(ji,jj)
-               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
-               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
-               !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
-               ENDIF
-
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
+               ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
+               ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
+               !
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
                ENDIF
 
@@ -437,6 +439,6 @@
                iku = miku(ji,jj) 
                ikv = mikv(ji,jj) 
-               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
-               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
+               ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
+               ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
 
                IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
