Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90	(revision 10922)
@@ -382,5 +382,5 @@
 
 
-   SUBROUTINE ldf_tra( kt )
+   SUBROUTINE ldf_tra( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_tra  ***
@@ -403,4 +403,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -411,5 +412,5 @@
          !                                ! =F(growth rate of baroclinic instability)
          !                                ! max value aeiv_0 ; decreased to 0 within 20N-20S
-         CALL ldf_eiv( kt, aei0, aeiu, aeiv )
+         CALL ldf_eiv( kt, aei0, aeiu, aeiv, Kmm )
       ENDIF
       !
@@ -424,5 +425,5 @@
             ahtv(:,:,1) = aeiv(:,:,1)
          ELSE                                            ! compute aht. 
-            CALL ldf_eiv( kt, aht0, ahtu, ahtv )
+            CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm )
          ENDIF
          !
@@ -448,11 +449,11 @@
          IF( ln_traldf_lap     ) THEN          !   laplacian operator |u| e /12
             DO jk = 1, jpkm1
-               ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12   ! n.b. ub,vb are masked
-               ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12
+               ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12   ! n.b. uu,vv are masked
+               ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12
             END DO
          ELSEIF( ln_traldf_blp ) THEN      ! bilaplacian operator      sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e
             DO jk = 1, jpkm1
-               ahtu(:,:,jk) = SQRT(  ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12  ) * e1u(:,:)
-               ahtv(:,:,jk) = SQRT(  ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12  ) * e2v(:,:)
+               ahtu(:,:,jk) = SQRT(  ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12  ) * e1u(:,:)
+               ahtv(:,:,jk) = SQRT(  ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12  ) * e2v(:,:)
             END DO
          ENDIF
@@ -625,5 +626,5 @@
 
 
-   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )
+   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv  ***
@@ -637,4 +638,5 @@
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kmm            ! ocean time level indices
       REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
@@ -658,9 +660,9 @@
                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
-                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm)
                   ! Compute elements required for the inverse time scale of baroclinic
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
                   zhw(ji,jj) = zhw(ji,jj) + ze3w
@@ -676,9 +678,9 @@
                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
-                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm)
                   ! Compute elements required for the inverse time scale of baroclinic
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
                      &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
@@ -725,5 +727,5 @@
 
 
-   SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype )
+   SUBROUTINE ldf_eiv_trp( kt, kit000, pu, pv, pw, cdtype, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv_trp  ***
@@ -741,12 +743,13 @@
       !!                                    velocity and heat transport (call ldf_eiv_dia)
       !!
-      !! ** Action  : pun, pvn increased by the eiv transport
+      !! ** Action  : pu, pv increased by the eiv transport
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
       INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
+      INTEGER                         , INTENT(in   ) ::   Kmm      ! ocean time level indices
       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean transport components   [m3/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean transport components   [m3/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv                [m3/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu      ! in : 3 ocean transport components   [m3/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv      ! out: 3 ocean transport components   [m3/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw      ! increased by the eiv                [m3/s]
       !!
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
@@ -780,6 +783,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.               
-               pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
-               pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
+               pu(ji,jj,jk) = pu(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
+               pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
             END DO
          END DO
@@ -788,5 +791,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               pwn(ji,jj,jk) = pwn(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
+               pw(ji,jj,jk) = pw(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
                   &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) )
             END DO
@@ -795,10 +798,10 @@
       !
       !                              ! diagnose the eddy induced velocity and associated heat transport
-      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 )
       !
     END SUBROUTINE ldf_eiv_trp
 
 
-   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )
+   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv_dia  ***
@@ -811,4 +814,5 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   psi_uw, psi_vw   ! streamfunction   [m3/s]
+      INTEGER                         , INTENT(in   ) ::   Kmm   ! ocean time level indices
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
@@ -831,10 +835,10 @@
       !
       DO jk = 1, jpkm1                                         ! e2u e3u u_eiv = -dk[psi_uw]
-         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) )
+         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u(:,:,jk,Kmm) )
       END DO
       CALL iom_put( "uoce_eiv", zw3d )
       !
       DO jk = 1, jpkm1                                         ! e1v e3v v_eiv = -dk[psi_vw]
-         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) )
+         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v(:,:,jk,Kmm) )
       END DO
       CALL iom_put( "voce_eiv", zw3d )
@@ -859,6 +863,6 @@
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                    &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+                    &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) ) 
                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
               END DO
@@ -875,6 +879,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+                  &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) ) 
                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
             END DO
@@ -894,6 +898,6 @@
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
-                    &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
+                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+                    &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) ) 
                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
               END DO
@@ -903,5 +907,5 @@
         CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. )
         CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
-        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
+        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                ! salt transport in i-direction
       ENDIF
       zw2d(:,:) = 0._wp 
@@ -910,6 +914,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
-                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
+               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+                  &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) ) 
                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
             END DO
