Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trddyn.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trddyn.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trddyn.F90	(revision 11949)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE trd_dyn( putrd, pvtrd, ktrd, kt )
+   SUBROUTINE trd_dyn( putrd, pvtrd, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_mod  ***
@@ -55,4 +55,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
       !!----------------------------------------------------------------------
       !
@@ -66,20 +67,20 @@
       !   3D output of momentum and/or tracers trends using IOM interface
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+      IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt, Kmm )
          
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  Integral Constraints Properties for momentum and/or tracers trends
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt )
+      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt, Kmm )
 
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  Kinetic Energy trends
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt )
+      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt, Kmm )
 
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  Vorticity trends
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt )
+      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt, Kmm )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -91,5 +92,5 @@
 
 
-   SUBROUTINE trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+   SUBROUTINE trd_dyn_iom( putrd, pvtrd, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_dyn_iom  ***
@@ -100,4 +101,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -121,9 +123,9 @@
                               z3dx(:,:,:) = 0._wp                  ! U.dxU & V.dyV (approximation)
                               z3dy(:,:,:) = 0._wp
-                              DO jk = 1, jpkm1   ! no mask as un,vn are masked
+                              DO jk = 1, jpkm1   ! no mask as uu, vv are masked
                                  DO jj = 2, jpjm1
                                     DO ji = 2, jpim1
-                                       z3dx(ji,jj,jk) = un(ji,jj,jk) * ( un(ji+1,jj,jk) - un(ji-1,jj,jk) ) / ( 2._wp * e1u(ji,jj) )
-                                       z3dy(ji,jj,jk) = vn(ji,jj,jk) * ( vn(ji,jj+1,jk) - vn(ji,jj-1,jk) ) / ( 2._wp * e2v(ji,jj) )
+                                       z3dx(ji,jj,jk) = uu(ji,jj,jk,Kmm) * ( uu(ji+1,jj,jk,Kmm) - uu(ji-1,jj,jk,Kmm) ) / ( 2._wp * e1u(ji,jj) )
+                                       z3dy(ji,jj,jk) = vv(ji,jj,jk,Kmm) * ( vv(ji,jj+1,jk,Kmm) - vv(ji,jj-1,jk,Kmm) ) / ( 2._wp * e2v(ji,jj) )
                                     END DO
                                  END DO
@@ -142,6 +144,6 @@
                               !                                    ! wind stress trends
                               ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
-                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u_n(:,:,1) * rau0 )
-                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v_n(:,:,1) * rau0 )
+                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rau0 )
+                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rau0 )
                               CALL iom_put( "utrd_tau", z2dx )
                               CALL iom_put( "vtrd_tau", z2dy )
@@ -159,7 +161,7 @@
 !                                          ikbv = mbkv(ji,jj)
 !                                          z3dx(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) ) & 
-!                                               &         * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
+!                                               &         * uu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm)
 !                                          z3dy(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) ) &
-!                                               &         * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
+!                                               &         * vv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm)
 !                                    END DO
 !                                 END DO
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdglo.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdglo.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdglo.F90	(revision 11949)
@@ -59,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
+   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_glo  ***
@@ -72,4 +72,5 @@
       CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       !!
       INTEGER ::   ji, jj, jk      ! dummy loop indices
@@ -87,11 +88,11 @@
       	      DO jj = 1, jpj
       	         DO ji = 1, jpi        
-      	            zvm = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+      	            zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
       	            zvt = ptrdx(ji,jj,jk) * zvm
       	            zvs = ptrdy(ji,jj,jk) * zvm
          	         tmo(ktrd) = tmo(ktrd) + zvt   
          	         smo(ktrd) = smo(ktrd) + zvs
-         	         t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
-         	         s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
+         	         t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
+         	         s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
                   END DO
                END DO
@@ -99,8 +100,8 @@
             !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
             IF( ln_linssh .AND. ktrd == jptra_zad ) THEN  
-               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) * tmask_i(:,:) )
-               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:) * tmask_i(:,:)  )
-               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:) * tmask_i(:,:)  )
-               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:) * tmask_i(:,:)  )
+               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
+               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
+               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
+               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
             ENDIF
             !
@@ -121,10 +122,10 @@
                   DO ji = 1, jpim1
                      zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
-                        &                                     * e1e2u  (ji,jj) * e3u_n(ji,jj,jk)
+                        &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
                      zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
-                        &                                     * e1e2v  (ji,jj) * e3u_n(ji,jj,jk)
+                        &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
                      umo(ktrd) = umo(ktrd) + zvt
                      vmo(ktrd) = vmo(ktrd) + zvs
-                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
+                     hke(ktrd) = hke(ktrd) + uu(ji,jj,jk,Kmm) * zvt + vv(ji,jj,jk,Kmm) * zvs
                   END DO
                END DO
@@ -141,5 +142,5 @@
                      umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
                      vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
-                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
+                     hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
                   END DO
                END DO
@@ -155,9 +156,9 @@
 !                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
 !                        ikbv = mbkv(ji,jj)
-!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * un(ji,jj,ikbu) * e1e2u(ji,jj)
-!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vn(ji,jj,ikbv) * e1e2v(ji,jj)
+!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu(ji,jj,ikbu,Kmm) * e1e2u(ji,jj)
+!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv(ji,jj,ikbv,Kmm) * e1e2v(ji,jj)
 !                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
 !                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
-!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
+!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + uu(ji,jj,ikbu,Kmm) * zvt + vv(ji,jj,ikbv,Kmm) * zvs
 !                     END DO
 !                  END DO
@@ -183,5 +184,5 @@
 
 
-   SUBROUTINE glo_dyn_wri( kt )
+   SUBROUTINE glo_dyn_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE glo_dyn_wri  ***
@@ -190,4 +191,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! time level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -209,10 +211,10 @@
          zkepe(:,:,:) = 0._wp
    
-         CALL eos( tsn, rhd, rhop )       ! now potential density
+         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
 
          zcof = 0.5_wp / rau0             ! Density flux at w-point
          zkz(:,:,1) = 0._wp
          DO jk = 2, jpk
-            zkz(:,:,jk) = zcof * e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
+            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
          END DO
          
@@ -221,6 +223,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
-                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
+                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
+                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
                END DO
             END DO
@@ -233,5 +235,5 @@
                      &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
                      &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
-                     &              / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                     &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
@@ -242,5 +244,5 @@
          peke = 0._wp
          DO jk = 1, jpkm1
-            peke = peke + SUM( zkepe(:,:,jk) * gdept_n(:,:,jk) * e1e2t(:,:) * e3t_n(:,:,jk) )
+            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) * e3t(:,:,jk,Kmm) )
          END DO
          peke = grav * peke
@@ -506,5 +508,5 @@
 
 
-   SUBROUTINE trd_glo_init
+   SUBROUTINE trd_glo_init( Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_glo_init  ***
@@ -512,4 +514,5 @@
       !! ** Purpose :   Read the namtrd namelist
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kmm   ! time level index
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -524,5 +527,5 @@
       tvolt = 0._wp
       DO jk = 1, jpkm1
-         tvolt = tvolt + SUM( e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
+         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
       END DO
       CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
@@ -541,6 +544,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u_n(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v_n(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
+               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
+               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
             END DO
          END DO
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdini.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdini.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdini.F90	(revision 11949)
@@ -34,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE trd_init
+   SUBROUTINE trd_init( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_init  ***
@@ -40,4 +40,5 @@
       !! ** Purpose :   Initialization of trend diagnostics
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kmm  ! time level index
       INTEGER ::   ios   ! local integer
       !!
@@ -96,5 +97,5 @@
 
       !                             ! diagnostic initialization  
-      IF( ln_glo_trd )   CALL trd_glo_init      ! global domain averaged trends
+      IF( ln_glo_trd )   CALL trd_glo_init( Kmm )      ! global domain averaged trends
       IF( ln_tra_mxl )   CALL trd_mxl_init      ! mixed-layer          trends  
       IF( ln_vor_trd )   CALL trd_vor_init      ! barotropic vorticity trends
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdken.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdken.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdken.F90	(revision 11949)
@@ -59,5 +59,5 @@
 
 
-   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt )
+   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_ken  ***
@@ -67,5 +67,5 @@
       !! ** Method  : - apply lbc to the input masked velocity trends 
       !!              - compute the associated KE trend:
-      !!          zke = 0.5 * (  mi-1[ un * putrd * bu ] + mj-1[ vn * pvtrd * bv]  ) / bt
+      !!          zke = 0.5 * (  mi-1[ uu(Kmm) * putrd * bu ] + mj-1[ vv(Kmm) * pvtrd * bv]  ) / bt
       !!      where bu, bv, bt are the volume of u-, v- and t-boxes. 
       !!              - vertical diffusion case (jpdyn_zdf): 
@@ -80,4 +80,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
       !
       INTEGER ::   ji, jj, jk       ! dummy loop indices
@@ -92,7 +93,7 @@
       nkstp = kt
       DO jk = 1, jpkm1
-         bu   (:,:,jk) =    e1e2u(:,:) * e3u_n(:,:,jk)
-         bv   (:,:,jk) =    e1e2v(:,:) * e3v_n(:,:,jk)
-         r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t_n(:,:,jk) * tmask(:,:,jk)
+         bu   (:,:,jk) =    e1e2u(:,:) * e3u(:,:,jk,Kmm)
+         bv   (:,:,jk) =    e1e2v(:,:) * e3v(:,:,jk,Kmm)
+         r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
       !
@@ -103,8 +104,8 @@
          DO jj = 2, jpj
             DO ji = 2, jpi
-               zke(ji,jj,jk) = 0.5_wp * rau0 *( un(ji  ,jj,jk) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
-                  &                           + un(ji-1,jj,jk) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
-                  &                           + vn(ji,jj  ,jk) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
-                  &                           + vn(ji,jj-1,jk) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
+               zke(ji,jj,jk) = 0.5_wp * rau0 *( uu(ji  ,jj,jk,Kmm) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
+                  &                           + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
+                  &                           + vv(ji,jj  ,jk,Kmm) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
+                  &                           + vv(ji,jj-1,jk,Kmm) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
             END DO
          END DO
@@ -122,6 +123,6 @@
          !                   !                                          ! wind stress trends
                                  ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
-                           z2dx(:,:) = un(:,:,1) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
-                           z2dy(:,:) = vn(:,:,1) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
+                           z2dx(:,:) = uu(:,:,1,Kmm) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
+                           z2dy(:,:) = vv(:,:,1,Kmm) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
                            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
                            DO jj = 2, jpj
@@ -141,6 +142,6 @@
 !                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
 !                  ikbv = mbkv(ji,jj)   
-!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu)
-!                  z2dy(ji,jj) = vn(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv)
+!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
+!                  z2dy(ji,jj) = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
 !               END DO
 !            END DO
@@ -157,5 +158,5 @@
          CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf"   , zke )    ! asselin filter trends 
 !! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdyn_bfr avant le swap !!!
-!! reflechir a une possible sauvegarde du "vrai" un,vn pour le calcul de atf....
+!! reflechir a une possible sauvegarde du "vrai" uu(Kmm),vv(Kmm) pour le calcul de atf....
 !
 !         IF( ln_drgimp ) THEN                                          ! bottom friction (implicit case)
@@ -164,6 +165,6 @@
 !                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
 !                  ikbv = mbkv(ji,jj)
-!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
-!                  z2dy(ji,jj) = un(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
+!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm)
+!                  z2dy(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm)
 !               END DO
 !            END DO
@@ -179,9 +180,9 @@
         CASE( jpdyn_ken )   ;   ! kinetic energy
                     ! called in dynnxt.F90 before asselin time filter
-                    ! with putrd=ua and pvtrd=va
+                    ! with putrd=uu(Krhs) and pvtrd=vv(Krhs)
                     zke(:,:,:) = 0.5_wp * zke(:,:,:)
                     CALL iom_put( "KE", zke )
                     !
-                    CALL ken_p2k( kt , zke )
+                    CALL ken_p2k( kt , zke, Kmm )
                       CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
          !
@@ -191,5 +192,5 @@
 
 
-   SUBROUTINE ken_p2k( kt , pconv )
+   SUBROUTINE ken_p2k( kt , pconv, Kmm )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE ken_p2k  ***
@@ -202,4 +203,5 @@
       !!---------------------------------------------------------------------- 
       INTEGER                   , INTENT(in   ) ::   kt      ! ocean time-step index
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pconv   ! 
       !
@@ -214,9 +216,9 @@
       
       !  Surface value (also valid in partial step case)
-      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * wn(:,:,1) * e3w_n(:,:,1)
+      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * ww(:,:,1) * e3w(:,:,1,Kmm)
 
       ! interior value (2=<jk=<jpkm1)
       DO jk = 2, jpk
-         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * wn(:,:,jk) * e3w_n(:,:,jk)
+         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * ww(:,:,jk) * e3w(:,:,jk,Kmm)
       END DO
 
@@ -225,5 +227,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcoef = 0.5_wp / e3t_n(ji,jj,jk)
+               zcoef = 0.5_wp / e3t(ji,jj,jk,Kmm)
                pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
             END DO
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdmxl.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdmxl.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdmxl.F90	(revision 11949)
@@ -86,5 +86,5 @@
 
 
-   SUBROUTINE trd_tra_mxl( ptrdx, ptrdy, ktrd, kt, p2dt, kmxln )
+   SUBROUTINE trd_tra_mxl( ptrdx, ptrdy, ktrd, kt, p2dt, kmxln, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_mng  ***
@@ -98,4 +98,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
       INTEGER                   , INTENT(in   ) ::   kt      ! time step index
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       REAL(wp)                  , INTENT(in   ) ::   p2dt    ! time step  [s]
       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   kmxln   ! number of t-box for the vertical average 
@@ -119,5 +120,5 @@
             DO jj = 1,jpj
                DO ji = 1,jpi
-                  IF( jk - kmxln(ji,jj) < 0 )   wkx(ji,jj,jk) = e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  IF( jk - kmxln(ji,jj) < 0 )   wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -136,6 +137,6 @@
          tml(:,:) = 0._wp   ;   sml(:,:) = 0._wp
          DO jk = 1, jpktrd
-            tml(:,:) = tml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_tem)
-            sml(:,:) = sml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_sal)
+            tml(:,:) = tml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_tem,Kmm)
+            sml(:,:) = sml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_sal,Kmm)
          END DO
          !
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdpen.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdpen.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdpen.F90	(revision 11949)
@@ -55,5 +55,5 @@
 
 
-   SUBROUTINE trd_pen( ptrdx, ptrdy, ktrd, kt, pdt )
+   SUBROUTINE trd_pen( ptrdx, ptrdy, ktrd, kt, pdt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_mng  ***
@@ -66,4 +66,5 @@
       INTEGER                   , INTENT(in) ::   ktrd           ! tracer trend index
       INTEGER                   , INTENT(in) ::   kt             ! time step index
+      INTEGER                   , INTENT(in) ::   Kmm            ! time level index
       REAL(wp)                  , INTENT(in) ::   pdt            ! time step [s]
       !
@@ -77,5 +78,5 @@
       IF( kt /= nkstp ) THEN     ! full eos: set partial derivatives at the 1st call of kt time step
          nkstp = kt
-         CALL eos_pen( tsn, rab_PE, zpe )
+         CALL eos_pen( ts(:,:,:,:,Kmm), rab_PE, zpe, Kmm )
          CALL iom_put( "alphaPE", rab_pe(:,:,:,jp_tem) )
          CALL iom_put( "betaPE" , rab_pe(:,:,:,jp_sal) )
@@ -95,8 +96,8 @@
                                 IF( ln_linssh ) THEN                   ! cst volume : adv flux through z=0 surface
                                    ALLOCATE( z2d(jpi,jpj) )
-                                   z2d(:,:) = wn(:,:,1) * ( &
-                                     &   - ( rab_n(:,:,1,jp_tem) + rab_pe(:,:,1,jp_tem) ) * tsn(:,:,1,jp_tem)    &
-                                     &   + ( rab_n(:,:,1,jp_sal) + rab_pe(:,:,1,jp_sal) ) * tsn(:,:,1,jp_sal)    &
-                                     & ) / e3t_n(:,:,1)
+                                   z2d(:,:) = ww(:,:,1) * ( &
+                                     &   - ( rab_n(:,:,1,jp_tem) + rab_pe(:,:,1,jp_tem) ) * ts(:,:,1,jp_tem,Kmm)    &
+                                     &   + ( rab_n(:,:,1,jp_sal) + rab_pe(:,:,1,jp_sal) ) * ts(:,:,1,jp_sal,Kmm)    &
+                                     & ) / e3t(:,:,1,Kmm)
                                    CALL iom_put( "petrd_sad" , z2d )
                                    DEALLOCATE( z2d )
@@ -112,12 +113,4 @@
       CASE ( jptra_bbc  )   ;   CALL iom_put( "petrd_bbc" , zpe )   ! bottom bound cond (geoth flux)
       CASE ( jptra_atf  )   ;   CALL iom_put( "petrd_atf" , zpe )   ! asselin time filter (last trend)
-                                !IF( ln_linssh ) THEN                   ! cst volume : ssh term (otherwise include in e3t variation)
-                                !   ALLOCATE( z2d(jpi,jpj) )
-                                !   z2d(:,:) = ( ssha(:,:) - sshb(:,:) )                 &
-                                !      &     * (   dPE_dt(:,:,1) * tsn(:,:,1,jp_tem)    &
-                                !      &         + dPE_ds(:,:,1) * tsn(:,:,1,jp_sal)  ) / ( e3t_n(:,:,1) * pdt )
-                                !   CALL iom_put( "petrd_sad" , z2d )
-                                !   DEALLOCATE( z2d )
-                                !ENDIF
          !
       END SELECT
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdtra.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdtra.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdtra.F90	(revision 11949)
@@ -60,5 +60,5 @@
 
 
-   SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pun, ptra )
+   SUBROUTINE trd_tra( kt, Kmm, Krhs, ctype, ktra, ktrd, ptrd, pu, ptra )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra  ***
@@ -77,6 +77,7 @@
       INTEGER                         , INTENT(in)           ::   ktra    ! tracer index
       INTEGER                         , INTENT(in)           ::   ktrd    ! tracer trend index
+      INTEGER                         , INTENT(in)           ::   Kmm, Krhs ! time level indices
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::   ptrd    ! tracer trend  or flux
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pun     ! now velocity 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pu      ! now velocity 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
       !
@@ -94,11 +95,11 @@
          SELECT CASE( ktrd )
          !                            ! advection: transform the advective flux into a trend
-         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx ) 
-         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty ) 
-         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt  ) 
+         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'X', trdtx, Kmm ) 
+         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'Y', trdty, Kmm ) 
+         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'Z', trdt, Kmm )
          CASE( jptra_bbc,    &        ! qsr, bbc: on temperature only, send to trd_tra_mng
             &  jptra_qsr )   ;   trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
                                  ztrds(:,:,:) = 0._wp
-                                 CALL trd_tra_mng( trdt, ztrds, ktrd, kt )
+                                 CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
  !!gm Gurvan, verify the jptra_evd trend please !
          CASE( jptra_evd )   ;   avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
@@ -114,10 +115,10 @@
          !                            ! advection: transform the advective flux into a trend
          !                            !            and send T & S trends to trd_tra_mng
-         CASE( jptra_xad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'X'  , ztrds ) 
-                                  CALL trd_tra_mng( trdtx, ztrds, ktrd, kt   )
-         CASE( jptra_yad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Y'  , ztrds ) 
-                                  CALL trd_tra_mng( trdty, ztrds, ktrd, kt   )
-         CASE( jptra_zad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Z'  , ztrds ) 
-                                  CALL trd_tra_mng( trdt , ztrds, ktrd, kt   )
+         CASE( jptra_xad  )   ;   CALL trd_tra_adv( ptrd , pu  , ptra, 'X'  , ztrds, Kmm ) 
+                                  CALL trd_tra_mng( trdtx, ztrds, ktrd, kt, Kmm   )
+         CASE( jptra_yad  )   ;   CALL trd_tra_adv( ptrd , pu  , ptra, 'Y'  , ztrds, Kmm ) 
+                                  CALL trd_tra_mng( trdty, ztrds, ktrd, kt, Kmm   )
+         CASE( jptra_zad  )   ;   CALL trd_tra_adv( ptrd , pu  , ptra, 'Z'  , ztrds, Kmm ) 
+                                  CALL trd_tra_mng( trdt , ztrds, ktrd, kt, Kmm   )
          CASE( jptra_zdfp )           ! diagnose the "PURE" Kz trend (here: just before the swap)
             !                         ! iso-neutral diffusion case otherwise jptra_zdf is "PURE"
@@ -127,28 +128,28 @@
             zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
             DO jk = 2, jpk
-               zwt(:,:,jk) = avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
-               zws(:,:,jk) = avs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
             END DO
             !
             ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t_n(:,:,jk)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t_n(:,:,jk) 
-            END DO
-            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )  
+               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
+               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm) 
+            END DO
+            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt, Kmm )  
             !
             !                         ! Also calculate EVD trend at this point. 
             zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes
             DO jk = 2, jpk
-               zwt(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
-               zws(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk)
+               zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
+               zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
             END DO
             !
             ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t_n(:,:,jk)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t_n(:,:,jk) 
-            END DO
-            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt )  
+               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
+               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm) 
+            END DO
+            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt, Kmm )  
             !
             DEALLOCATE( zwt, zws, ztrdt )
@@ -156,5 +157,5 @@
          CASE DEFAULT                 ! other trends: mask and send T & S trends to trd_tra_mng
             ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
-            CALL trd_tra_mng( trdt, ztrds, ktrd, kt )  
+            CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )  
          END SELECT
       ENDIF
@@ -164,12 +165,12 @@
          SELECT CASE( ktrd )
          !                            ! advection: transform the advective flux into a masked trend
-         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'X', ztrds ) 
-         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Y', ztrds ) 
-         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Z', ztrds ) 
+         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'X', ztrds, Kmm ) 
+         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'Y', ztrds, Kmm ) 
+         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'Z', ztrds, Kmm ) 
          CASE DEFAULT                 ! other trends: just masked 
                                  ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
          END SELECT
          !                            ! send trend to trd_trc
-         CALL trd_trc( ztrds, ktra, ktrd, kt ) 
+         CALL trd_trc( ztrds, ktra, ktrd, kt, Kmm ) 
          !
       ENDIF
@@ -178,5 +179,5 @@
 
 
-   SUBROUTINE trd_tra_adv( pf, pun, ptn, cdir, ptrd )
+   SUBROUTINE trd_tra_adv( pf, pu, pt, cdir, ptrd, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_adv  ***
@@ -191,8 +192,9 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pf      ! advective flux in one direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pun     ! now velocity   in one direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptn     ! now or before tracer 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu      ! now velocity   in one direction
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt      ! now or before tracer 
       CHARACTER(len=1)                , INTENT(in   ) ::   cdir    ! X/Y/Z direction
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   ptrd    ! advective trend in one direction
+      INTEGER,  INTENT(in)                            ::   Kmm     ! time level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -215,6 +217,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
-                 &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   &
-                 &              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                 &                  - ( pu(ji,jj,jk) - pu(ji-ii,jj-ij,jk-ik) ) * pt(ji,jj,jk)  )   &
+                 &              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
             END DO
          END DO
@@ -224,5 +226,5 @@
 
 
-   SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt )
+   SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_mng  ***
@@ -236,4 +238,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       !!----------------------------------------------------------------------
 
@@ -243,11 +246,11 @@
 
       !                   ! 3D output of tracers trends using IOM interface
-      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt )
+      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt, Kmm )
 
       !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt )
+      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt, Kmm )
 
       !                   ! Potential ENergy trends
-      IF( ln_PE_trd  )   CALL trd_pen( ptrdx, ptrdy, ktrd, kt, r2dt )
+      IF( ln_PE_trd  )   CALL trd_pen( ptrdx, ptrdy, ktrd, kt, r2dt, Kmm )
 
       !                   ! Mixed layer trends for active tracers
@@ -290,5 +293,5 @@
 
 
-   SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt )
+   SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_tra_iom  ***
@@ -300,4 +303,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       !!
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -326,6 +330,6 @@
                                   IF( ln_linssh ) THEN                   ! cst volume : adv flux through z=0 surface
                                      ALLOCATE( z2dx(jpi,jpj), z2dy(jpi,jpj) )
-                                     z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / e3t_n(:,:,1)
-                                     z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / e3t_n(:,:,1)
+                                     z2dx(:,:) = ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) / e3t(:,:,1,Kmm)
+                                     z2dy(:,:) = ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) / e3t(:,:,1,Kmm)
                                      CALL iom_put( "ttrd_sad", z2dx )
                                      CALL iom_put( "strd_sad", z2dy )
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdtrc.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdtrc.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdtrc.F90	(revision 11949)
@@ -9,6 +9,7 @@
 CONTAINS
 
-   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt )
+   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt, Kmm )
       INTEGER ::   kt, kjn, ktrd   
+      INTEGER ::   Kmm            ! time level index
       REAL    ::   ptrtrd(:,:,:)  
       WRITE(*,*) 'trd_trc : You should not have seen this print! error?', ptrtrd(1,1,1)
Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdvor.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdvor.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdvor.F90	(revision 11949)
@@ -78,5 +78,5 @@
 
 
-   SUBROUTINE trd_vor( putrd, pvtrd, ktrd, kt )
+   SUBROUTINE trd_vor( putrd, pvtrd, ktrd, kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_vor  ***
@@ -88,4 +88,5 @@
       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
       !
       INTEGER ::   ji, jj   ! dummy loop indices
@@ -94,27 +95,27 @@
 
       SELECT CASE( ktrd ) 
-      CASE( jpdyn_hpg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_prg )   ! Hydrostatique Pressure Gradient 
-      CASE( jpdyn_keg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_keg )   ! KE Gradient 
-      CASE( jpdyn_rvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_rvo )   ! Relative Vorticity 
-      CASE( jpdyn_pvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_pvo )   ! Planetary Vorticity Term 
-      CASE( jpdyn_ldf )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_ldf )   ! Horizontal Diffusion 
-      CASE( jpdyn_zad )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_zad )   ! Vertical Advection 
-      CASE( jpdyn_spg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_spg )   ! Surface Pressure Grad. 
+      CASE( jpdyn_hpg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_prg, Kmm )   ! Hydrostatique Pressure Gradient 
+      CASE( jpdyn_keg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_keg, Kmm )   ! KE Gradient 
+      CASE( jpdyn_rvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_rvo, Kmm )   ! Relative Vorticity 
+      CASE( jpdyn_pvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_pvo, Kmm )   ! Planetary Vorticity Term 
+      CASE( jpdyn_ldf )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_ldf, Kmm )   ! Horizontal Diffusion 
+      CASE( jpdyn_zad )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_zad, Kmm )   ! Vertical Advection 
+      CASE( jpdyn_spg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_spg, Kmm )   ! Surface Pressure Grad. 
       CASE( jpdyn_zdf )                                                      ! Vertical Diffusion 
          ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
          DO jj = 2, jpjm1                                                             ! wind stress trends
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u_n(ji,jj,1) * rau0 )
-               ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v_n(ji,jj,1) * rau0 )
+               ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rau0 )
+               ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rau0 )
             END DO
          END DO
          !
-         CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf )                             ! zdf trend including surf./bot. stresses 
-         CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                             ! surface wind stress 
+         CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf, Kmm )                             ! zdf trend including surf./bot. stresses 
+         CALL trd_vor_zint( ztswu, ztswv, jpvor_swf, Kmm )                             ! surface wind stress 
       CASE( jpdyn_bfr )
-         CALL trd_vor_zint( putrd, pvtrd, jpvor_bfr )                             ! Bottom stress
+         CALL trd_vor_zint( putrd, pvtrd, jpvor_bfr, Kmm )                             ! Bottom stress
          !
       CASE( jpdyn_atf )       ! last trends: perform the output of 2D vorticity trends
-         CALL trd_vor_iom( kt )
+         CALL trd_vor_iom( kt, Kmm )
       END SELECT
       !
@@ -122,5 +123,5 @@
 
 
-   SUBROUTINE trd_vor_zint_2d( putrdvor, pvtrdvor, ktrd )
+   SUBROUTINE trd_vor_zint_2d( putrdvor, pvtrdvor, ktrd, Kmm )
       !!----------------------------------------------------------------------------
       !!                  ***  ROUTINE trd_vor_zint  ***
@@ -150,4 +151,5 @@
       !!----------------------------------------------------------------------
       INTEGER                     , INTENT(in   ) ::   ktrd       ! ocean trend index
+      INTEGER                     , INTENT(in   ) ::   Kmm        ! time level index
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   putrdvor   ! u vorticity trend 
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pvtrdvor   ! v vorticity trend
@@ -175,18 +177,18 @@
                ikbu = mbkv(ji,jj)
                ikbv = mbkv(ji,jj)            
-               zudpvor(ji,jj) = putrdvor(ji,jj) * e3u_n(ji,jj,ikbu) * e1u(ji,jj) * umask(ji,jj,ikbu)
-               zvdpvor(ji,jj) = pvtrdvor(ji,jj) * e3v_n(ji,jj,ikbv) * e2v(ji,jj) * vmask(ji,jj,ikbv)
+               zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,ikbu,Kmm) * e1u(ji,jj) * umask(ji,jj,ikbu)
+               zvdpvor(ji,jj) = pvtrdvor(ji,jj) * e3v(ji,jj,ikbv,Kmm) * e2v(ji,jj) * vmask(ji,jj,ikbv)
             END DO
          END DO
          !
       CASE( jpvor_swf )        ! wind stress
-         zudpvor(:,:) = putrdvor(:,:) * e3u_n(:,:,1) * e1u(:,:) * umask(:,:,1)
-         zvdpvor(:,:) = pvtrdvor(:,:) * e3v_n(:,:,1) * e2v(:,:) * vmask(:,:,1)
+         zudpvor(:,:) = putrdvor(:,:) * e3u(:,:,1,Kmm) * e1u(:,:) * umask(:,:,1)
+         zvdpvor(:,:) = pvtrdvor(:,:) * e3v(:,:,1,Kmm) * e2v(:,:) * vmask(:,:,1)
          !
       END SELECT
 
       ! Average except for Beta.V
-      zudpvor(:,:) = zudpvor(:,:) * r1_hu_n(:,:)
-      zvdpvor(:,:) = zvdpvor(:,:) * r1_hv_n(:,:)
+      zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
+      zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
    
       ! Curl
@@ -207,5 +209,5 @@
 
 
-   SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd )
+   SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd , Kmm )
       !!----------------------------------------------------------------------------
       !!                  ***  ROUTINE trd_vor_zint  ***
@@ -236,4 +238,5 @@
       !
       INTEGER                         , INTENT(in   ) ::   ktrd       ! ocean trend index
+      INTEGER                         , INTENT(in   ) ::   Kmm        ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   putrdvor   ! u vorticity trend 
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvtrdvor   ! v vorticity trend
@@ -257,6 +260,6 @@
       ! putrdvor and pvtrdvor terms
       DO jk = 1,jpk
-        zudpvor(:,:) = zudpvor(:,:) + putrdvor(:,:,jk) * e3u_n(:,:,jk) * e1u(:,:) * umask(:,:,jk)
-        zvdpvor(:,:) = zvdpvor(:,:) + pvtrdvor(:,:,jk) * e3v_n(:,:,jk) * e2v(:,:) * vmask(:,:,jk)
+        zudpvor(:,:) = zudpvor(:,:) + putrdvor(:,:,jk) * e3u(:,:,jk,Kmm) * e1u(:,:) * umask(:,:,jk)
+        zvdpvor(:,:) = zvdpvor(:,:) + pvtrdvor(:,:,jk) * e3v(:,:,jk,Kmm) * e2v(:,:) * vmask(:,:,jk)
       END DO
 
@@ -273,10 +276,10 @@
          END DO
          ! Average of the Curl and Surface mask
-         vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * r1_hu_n(:,:) * fmask(:,:,1)
+         vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * r1_hu(:,:,Kmm) * fmask(:,:,1)
       ENDIF
       !
       ! Average 
-      zudpvor(:,:) = zudpvor(:,:) * r1_hu_n(:,:)
-      zvdpvor(:,:) = zvdpvor(:,:) * r1_hv_n(:,:)
+      zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
+      zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
       !
       ! Curl
@@ -298,5 +301,5 @@
 
 
-   SUBROUTINE trd_vor_iom( kt )
+   SUBROUTINE trd_vor_iom( kt , Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_vor  ***
@@ -306,9 +309,10 @@
       !!----------------------------------------------------------------------
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
       !
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       INTEGER  ::   it, itmod        ! local integers
       REAL(wp) ::   zmean            ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) :: zun, zvn
+      REAL(wp), DIMENSION(jpi,jpj) :: zuu, zvv
       !!----------------------------------------------------------------------
 
@@ -327,6 +331,6 @@
 
       vor_avr   (:,:) = 0._wp
-      zun       (:,:) = 0._wp
-      zvn       (:,:) = 0._wp
+      zuu       (:,:) = 0._wp
+      zvv       (:,:) = 0._wp
       vor_avrtot(:,:) = 0._wp
       vor_avrres(:,:) = 0._wp
@@ -334,16 +338,16 @@
       ! Vertically averaged velocity
       DO jk = 1, jpk - 1
-         zun(:,:) = zun(:,:) + e1u(:,:) * un(:,:,jk) * e3u_n(:,:,jk)
-         zvn(:,:) = zvn(:,:) + e2v(:,:) * vn(:,:,jk) * e3v_n(:,:,jk)
+         zuu(:,:) = zuu(:,:) + e1u(:,:) * uu(:,:,jk,Kmm) * e3u(:,:,jk,Kmm)
+         zvv(:,:) = zvv(:,:) + e2v(:,:) * vv(:,:,jk,Kmm) * e3v(:,:,jk,Kmm)
       END DO
  
-      zun(:,:) = zun(:,:) * r1_hu_n(:,:)
-      zvn(:,:) = zvn(:,:) * r1_hv_n(:,:)
+      zuu(:,:) = zuu(:,:) * r1_hu(:,:,Kmm)
+      zvv(:,:) = zvv(:,:) * r1_hv(:,:,Kmm)
 
       ! Curl
       DO ji = 1, jpim1
          DO jj = 1, jpjm1
-            vor_avr(ji,jj) = (  ( zvn(ji+1,jj) - zvn(ji,jj) )    &
-               &              - ( zun(ji,jj+1) - zun(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1)
+            vor_avr(ji,jj) = (  ( zvv(ji+1,jj) - zvv(ji,jj) )    &
+               &              - ( zuu(ji,jj+1) - zuu(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1)
          END DO
       END DO
