Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90	(revision 10883)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90	(revision 10884)
@@ -45,5 +45,5 @@
 CONTAINS
    
-   SUBROUTINE dyn_zdf( kt )
+   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_zdf  ***
@@ -54,16 +54,18 @@
       !!
       !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
-      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
-      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
+      !!         puu(:,:,:,Kaa) =         puu(:,:,:,Kbb) + 2*dt *       puu(:,:,:,Krhs)             vector form or linear free surf.
+      !!         puu(:,:,:,Kaa) = ( e3u_b*puu(:,:,:,Kbb) + 2*dt * e3u_n*puu(:,:,:,Krhs) ) / e3u(:,:,:,Kaa)   otherwise
       !!               - update the after velocity with the implicit vertical mixing.
       !!      This requires to solver the following system: 
-      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
+      !!         puu(:,:,:,Kaa) = puu(:,:,:,Kaa) + 1/e3u(:,:,:,Kaa) dk+1[ mi(avm) / e3uw(:,:,:,Kaa) dk[ua] ]
       !!      with the following surface/top/bottom boundary condition:
       !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
       !!
-      !! ** Action :   (ua,va)   after velocity 
+      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity 
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
       !
       INTEGER  ::   ji, jj, jk         ! dummy loop indices
@@ -96,27 +98,27 @@
       !
       !                             !* explicit top/bottom drag case
-      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
+      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
       !
       !
       IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
-         ztrdu(:,:,:) = ua(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:)
-      ENDIF
-      !
-      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
+      ENDIF
+      !
+      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
       !
       !                    ! time stepping except vertical diffusion
       IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
          DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
-            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)
          END DO
       ELSE                                      ! applied on thickness weighted velocity
          DO jk = 1, jpkm1
-            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
-               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
-               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
+            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  &
+               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  &
+               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk)
          END DO
       ENDIF
@@ -125,9 +127,9 @@
       !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 
       !                not lead to the effective stress seen over the whole barotropic loop. 
-      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
+      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
       IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
          DO jk = 1, jpkm1        ! remove barotropic velocities
-            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - ua_b(:,:) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - va_b(:,:) ) * vmask(:,:,jk)
          END DO
          DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
@@ -135,8 +137,8 @@
                iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
                ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
+               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
+               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -146,8 +148,8 @@
                   iku = miku(ji,jj)         ! top ocean level at u- and v-points 
                   ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
-                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
+                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
+                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
+                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
                END DO
             END DO
@@ -165,9 +167,9 @@
                DO jj = 2, jpjm1 
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                      zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
                      zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
@@ -182,7 +184,7 @@
                DO jj = 2, jpjm1 
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
-                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                      zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
                      zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
@@ -197,6 +199,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zwi(ji,jj,1) = 0._wp
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
-               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
+               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
                zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) )
                zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
@@ -210,9 +212,9 @@
                DO jj = 2, jpjm1 
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                      zwi(ji,jj,jk) = zzwi
                      zws(ji,jj,jk) = zzws
@@ -225,7 +227,7 @@
                DO jj = 2, jpjm1 
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
-                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                      zwi(ji,jj,jk) = zzwi
                      zws(ji,jj,jk) = zzws
@@ -254,5 +256,5 @@
             DO ji = 2, jpim1
                iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
+               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
                zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
             END DO
@@ -263,5 +265,5 @@
                   !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
                   iku = miku(ji,jj)       ! ocean top level at u- and v-points 
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
                   zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
                END DO
@@ -282,5 +284,5 @@
       !   m is decomposed in the product of an upper and a lower triangular matrix
       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
-      !   The solution (the after velocity) is in ua
+      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
       !-----------------------------------------------------------------------
       !
@@ -295,6 +297,6 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
-            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
+            puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
                &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
          END DO
@@ -303,5 +305,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
+               puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
             END DO
          END DO
@@ -310,5 +312,5 @@
       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+            puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
          END DO
       END DO
@@ -316,5 +318,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+               puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
             END DO
          END DO
@@ -331,9 +333,9 @@
                DO jj = 2, jpjm1 
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
+                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                      zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
                      zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
@@ -348,7 +350,7 @@
                DO jj = 2, jpjm1 
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
-                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
+                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                      zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
                      zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
@@ -363,6 +365,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zwi(ji,jj,1) = 0._wp
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
-               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
+               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
                zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) )
                zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
@@ -376,9 +378,9 @@
                DO jj = 2, jpjm1   
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
+                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                      zwi(ji,jj,jk) = zzwi
                      zws(ji,jj,jk) = zzws
@@ -391,7 +393,7 @@
                DO jj = 2, jpjm1   
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
-                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
+                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                      zwi(ji,jj,jk) = zzwi
                      zws(ji,jj,jk) = zzws
@@ -419,5 +421,5 @@
             DO ji = 2, jpim1
                ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
+               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
                zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
             END DO
@@ -427,5 +429,5 @@
                DO ji = 2, jpim1
                   ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
+                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
                   zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
                END DO
@@ -446,5 +448,5 @@
       !   m is decomposed in the product of an upper and lower triangular matrix
       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
-      !   The solution (after velocity) is in 2d array va
+      !   The solution (after velocity) is in 2d array pvv(:,:,:,Kaa)
       !-----------------------------------------------------------------------
       !
@@ -459,6 +461,6 @@
       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 ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
-            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
+            pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
                &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
          END DO
@@ -467,5 +469,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
+               pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
             END DO
          END DO
@@ -474,5 +476,5 @@
       DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+            pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
          END DO
       END DO
@@ -480,5 +482,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+               pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
             END DO
          END DO
@@ -486,6 +488,6 @@
       !
       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
-         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
+         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:)
+         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
          DEALLOCATE( ztrdu, ztrdv ) 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10883)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10884)
@@ -92,5 +92,5 @@
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean transport components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90	(revision 10883)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90	(revision 10884)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE tra_zdf( kt )
+   SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf  ***
@@ -50,5 +50,7 @@
       !! ** Purpose :   compute the vertical ocean tracer physics.
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                                  , INTENT(in)    :: kt                  ! ocean time-step index
+      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs, Kaa ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts                 ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   jk   ! Dummy loop indices
@@ -70,10 +72,10 @@
       IF( l_trdtra )   THEN                  !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
       ENDIF
       !
       !                                      !* compute lateral mixing trend and add it to the general trend
-      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) 
+      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, pts, Kaa, jpts ) 
 
 !!gm WHY here !   and I don't like that !
@@ -81,13 +83,13 @@
       ! JMM avoid negative salinities near river outlet ! Ugly fix
       ! JMM : restore negative salinities to small salinities:
-      WHERE( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp
+      WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp
 !!gm
 
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &
-               &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &
-              &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)
+            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) &
+               &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
+              &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk)
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
@@ -107,5 +109,5 @@
 
  
-   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 
+   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, pt, Kaa, kjpt ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf_imp  ***
@@ -125,13 +127,13 @@
       !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
       !!
-      !! ** Action  : - pta  becomes the after tracer
-      !!---------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
-      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field
+      !! ** Action  : - pt(:,:,:,:,Kaa)  becomes the after tracer
+      !!---------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
+      INTEGER                                  , INTENT(in   ) ::   kit000   ! first time step index
+      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt     ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt       ! tracers and RHS of tracer equation
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -181,7 +183,7 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
-                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
-                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
-                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   &
+                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm)
+                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   &
                            &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
                         zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
@@ -194,7 +196,7 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk)
-                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
-                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm)
+                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
                     END DO
                   END DO
@@ -218,5 +220,5 @@
             !   The solution will be in the 4d array pta.
             !   The 3d array zwt is used as a work space array.
-            !   En route to the solution pta is used a to evaluate the rhs and then 
+            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 
             !   used as a work space array: its value is modified.
             !
@@ -238,5 +240,5 @@
          DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
             DO ji = fs_2, fs_jpim1
-               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)
+               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Kaa)
             END DO
          END DO
@@ -244,6 +246,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side
-                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
+                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Kaa)   ! zrhs=right hand side
+                  pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
                END DO
             END DO
@@ -252,5 +254,5 @@
          DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
             DO ji = fs_2, fs_jpim1
-               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+               pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
             END DO
          END DO
@@ -258,5 +260,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
+                  pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   &
                      &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10883)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10884)
@@ -208,5 +208,5 @@
       ENDIF
       
-                         CALL dyn_zdf       ( kstp )  ! vertical diffusion
+                         CALL dyn_zdf( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion  	==> after
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -261,5 +261,5 @@
       IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
 !!gm
-                         CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields
+                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer	==> after
       IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10883)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90	(revision 10884)
@@ -77,5 +77,5 @@
          IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc       ! tracers sponge
 #endif
-                                CALL trc_zdf    ( kt )      ! vertical mixing and after tracer fields
+                                CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
                                 CALL trc_nxt    ( kt )      ! tracer fields at next time step     
          IF( ln_trcrad )        CALL trc_rad    ( kt )      ! Correct artificial negative concentrations
@@ -86,5 +86,5 @@
                                 CALL trc_sbc( kt )            ! surface boundary condition
          IF( ln_trcdmp )        CALL trc_dmp( kt )            ! internal damping trends
-                                CALL trc_zdf( kt )            ! vertical mixing and after tracer fields
+                                CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
                                 CALL trc_nxt( kt )            ! tracer fields at next time step     
           IF( ln_trcrad )       CALL trc_rad( kt )            ! Correct artificial negative concentrations
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90	(revision 10883)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90	(revision 10884)
@@ -36,5 +36,5 @@
 CONTAINS
 
-   SUBROUTINE trc_zdf( kt )
+   SUBROUTINE trc_zdf( kt, Kbb, Kmm, Krhs, ptr, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_zdf  ***
@@ -43,5 +43,7 @@
       !!              an implicit time-stepping scheme.
       !!---------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::  kt      ! ocean time-step index
+      INTEGER                                   , INTENT(in   ) ::   kt                   ! ocean time-step index
+      INTEGER                                   , INTENT(in   ) ::   Kbb, Kmm, Krhs, Kaa  ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::   ptr                  ! passive tracers and RHS of tracer equation
       !
       INTEGER               ::  jk, jn
@@ -52,12 +54,12 @@
       IF( ln_timing )   CALL timing_start('trc_zdf')
       !
-      IF( l_trdtrc )   ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+      IF( l_trdtrc )   ztrtrd(:,:,:,:)  = ptr(:,:,:,:,Krhs)
       !
-      CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, trb, tra, jptra )    !   implicit scheme          
+      CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, Kbb, Kmm, ptr, Kaa, jptra )    !   implicit scheme          
       !
       IF( l_trdtrc )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jn = 1, jptra
             DO jk = 1, jpkm1
-               ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dttrc ) - ztrtrd(:,:,jk,jn)
+               ztrtrd(:,:,jk,jn) = ( ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kbb) ) / r2dttrc ) - ztrtrd(:,:,jk,jn)
             END DO
             CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) )
