Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traldf_triad.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traldf_triad.F90	(revision 11949)
@@ -48,7 +48,7 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
-      &                                                     pgui, pgvi,   &
-      &                                         ptb , ptbb, pta , kjpt, kpass )
+  SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv,               &
+      &                                              pgu , pgv  , pgui, pgvi , &
+      &                                         pt , pt2, pt_rhs, kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_triad  ***
@@ -66,5 +66,5 @@
       !!      see documentation for the desciption
       !!
-      !! ** Action :   pta   updated with the before rotated diffusion
+      !! ** Action :   pt_rhs   updated with the before rotated diffusion
       !!               ah_wslp2 ....
       !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp)
@@ -75,10 +75,11 @@
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
       INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      INTEGER                              , INTENT(in)    ::   Kmm        ! ocean time level indices
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -142,10 +143,10 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
-                        zbu   = e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
+                        zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
                         zah   = 0.25_wp * pahu(ji,jj,jk)
                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
                         zslope2 = zslope2 *zslope2
                         ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
@@ -166,11 +167,11 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp)
-                        zbv   = e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                        ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,Kmm)
+                        zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                         zah   = 0.25_wp * pahv(ji,jj,jk)
                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                         ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
                         !    (do this by *adding* gradient of depth)
-                        zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                        zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
                         zslope2 = zslope2 * zslope2
                         ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
@@ -193,5 +194,5 @@
                      DO ji = 1, fs_jpim1
                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                      END DO
                   END DO
@@ -201,5 +202,5 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, fs_jpim1
-                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
+                        ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
@@ -213,5 +214,5 @@
          ENDIF
          !
-         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
+         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
          !
       ENDIF                                  !==  end 1st pass only  ==!
@@ -229,6 +230,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -257,9 +258,9 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
+            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
             !
             !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
-            ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
+            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
             ENDIF
             !
@@ -273,10 +274,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
                            !
-                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
                            zah = pahu(ji,jj,jk)
@@ -296,9 +297,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
                            zah = pahv(ji,jj,jk)
@@ -320,10 +321,10 @@
                            ze1ur = r1_e1u(ji,jj)
                            zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp)
+                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
                            zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
                            !
-                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk)
+                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
@@ -343,9 +344,9 @@
                            ze2vr = r1_e2v(ji,jj)
                            zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp)
+                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm)
                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                            ! ln_botmix_triad is .F. mask zah for bottom half cells
                            zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
@@ -362,7 +363,7 @@
             DO jj = 2 , jpjm1
                DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
                      &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
-                     &                                        / (  e1e2t(ji,jj) * e3t_n(ji,jj,jk)  )
+                     &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
                END DO
             END DO
@@ -375,7 +376,7 @@
                DO jj = 1, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)   &
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
-                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                   END DO
                END DO
@@ -387,16 +388,16 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)             &
-                           &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
+                           &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                      END DO
                   END DO
                END DO 
-            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)                      &
-                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
-                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
+                           &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
                      END DO
                   END DO
@@ -405,9 +406,9 @@
          ENDIF
          !
-         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-                     &                                        / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
+                     &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
                END DO
             END DO
