MODULE dynzdf
   !!==============================================================================
   !!                 ***  MODULE  dynzdf  ***
   !! Ocean dynamics :  vertical component of the momentum mixing trend
   !!==============================================================================
   !! History :  1.0  !  2005-11  (G. Madec)  Original code
   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option + avm at t-point
   !!            4.5  !  2022-06  (S. Techene, G, Madec) refactorization to reduce local memory usage
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers variables
   USE phycst         ! physical constants
   USE dom_oce        ! ocean space and time domain variables 
   USE sbc_oce        ! surface boundary condition: ocean
   USE zdf_oce        ! ocean vertical physics variables
   USE zdfdrg         ! vertical physics: top/bottom drag coef.
   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
#if defined key_loop_fusion
   USE dynldf_iso_lf,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing 
#else
   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing 
#endif
   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator
   USE trd_oce        ! trends: ocean variables
   USE trddyn         ! trend manager: dynamics
   !
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! MPP library
   USE prtctl         ! Print control
   USE timing         ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dyn_zdf   !  routine called by step.F90

   !! * Substitutions
#  include "do_loop_substitute.h90"
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: dynzdf.F90 14547 2021-02-25 17:07:15Z techene $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS
   
   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_zdf  ***
      !!
      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
      !!              together with the Leap-Frog time stepping using an 
      !!              implicit scheme.
      !!
      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after   otherwise
      !!               - update the after velocity with the implicit vertical mixing.
      !!      This requires to solver the following system: 
      !!         u(after) = u(after) + 1/e3u_after  dk+1[ mi(avm) / e3uw_after 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 :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity 
      !!---------------------------------------------------------------------
      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
      INTEGER  ::   iku, ikv             ! local integers
      REAL(wp) ::   zzwi, ze3ua, zDt_2   ! local scalars
      REAL(wp) ::   zzws, ze3va          !   -      -
      REAL(wp) ::   z1_e3ua, z1_e3va     !   -      -
      REAL(wp) ::   zWu , zWv            !   -      -
      REAL(wp) ::   zWui, zWvi           !   -      -
      REAL(wp) ::   zWus, zWvs           !   -      -
      REAL(wp), DIMENSION(A1Di(0),jpk)        ::   zwi, zwd, zws  ! 2D workspace
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
      !!---------------------------------------------------------------------
      !
      IF( ln_timing )   CALL timing_start('dyn_zdf')
      !
      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
         IF( kt == nit000 ) THEN       !* initialization
            IF(lwp) WRITE(numout,*)
            IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
         ENDIF
      ENDIF
      !
      zDt_2 = rDt * 0.5_wp
      !
      !                             !* explicit top/bottom drag case
      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, 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(:,:,:) = puu(:,:,:,Krhs)
         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
      ENDIF
      !
      !                                               ! ================= !
      DO_1Dj( 0, 0 )                                  !  i-k slices loop  !
         !                                            ! ================= !
         !
         !              !==  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_2Dik( 0, 0,    1, jpkm1, 1 )
               puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
               pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
            END_2D
         ELSE                                      ! applied on thickness weighted velocity
            DO_2Dik( 0, 0,    1, jpkm1, 1 )
               puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb )  &
                  &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) &
                  &                        / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
               pvv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb )  &
                  &                  + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs)  ) &
                  &                        / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
            END_2D
         ENDIF
         !                    ! add top/bottom friction 
         !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
         !     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(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
         IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
            DO_2Dik( 0, 0,     1, jpkm1, 1 )      ! remove barotropic velocities
               puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
               pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
            END_2D
            DO_1Di( 0, 0 )      ! Add bottom/top stress due to barotropic component only
               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa)   &
                  &                                            / e3u(ji,jj,iku,Kaa)
               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa)   &
                  &                                            / e3v(ji,jj,ikv,Kaa)
            END_1D
            IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
               DO_1Di( 0, 0 )
                  iku = miku(ji,jj)         ! top ocean level at u- and v-points 
                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa)   &
                     &                                            / e3u(ji,jj,iku,Kaa)
                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa)   &
                     &                                            / e3v(ji,jj,ikv,Kaa)
               END_1D
            END IF
         ENDIF
         !
         !              !==  Vertical diffusion on u  ==!
         !
         !
         !                    !* Matrix construction
         IF( ln_zad_Aimp ) THEN      !- including terms associated with partly implicit vertical advection 
            SELECT CASE( nldf_dyn )
            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
               DO_2Dik( 0, 0,   1, jpkm1, 1 )
                  z1_e3ua =  1._wp  / e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
                  zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
                     &         / e3uw(ji,jj,jk  ,Kmm) * z1_e3ua * wumask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
                     &         / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
                  zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) * z1_e3ua
                  zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
                  zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) 
                  zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
                  zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
               END_2D
            CASE DEFAULT               ! iso-level lateral mixing
               DO_2Dik( 0, 0,   1, jpkm1, 1 )
                  z1_e3ua =  1._wp  / e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
                  zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
                     &         / e3uw(ji,jj,jk  ,Kmm) * z1_e3ua * wumask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
                     &         / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
                  zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) * z1_e3ua
                  zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
                  zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
                  zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
                  zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
               END_2D
            END SELECT
            !
            zwi(:,1) = 0._wp
            DO_1Di( 0, 0 )     !* Surface boundary conditions
               zwi(ji,1) = 0._wp
               zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
                  &         / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
               zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
               zws(ji,1) = zzws - zDt_2 * MAX( zWus, 0._wp )
               zwd(ji,1) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
            END_1D
         ELSE                       !- only vertical diffusive terms
            SELECT CASE( nldf_dyn )
            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
               DO_2Dik( 0, 0,   1, jpkm1, 1 )
                  zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
                     &         / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
                     &         / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                  zwi(ji,jk) = zzwi
                  zws(ji,jk) = zzws
                  zwd(ji,jk) = 1._wp - zzwi - zzws
               END_2D
            CASE DEFAULT               ! iso-level lateral mixing
               DO_2Dik( 0, 0,   1, jpkm1, 1 )
                  zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    &
                     &         / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    &
                     &         / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                  zwi(ji,jk) = zzwi
                  zws(ji,jk) = zzws
                  zwd(ji,jk) = 1._wp - zzwi - zzws
               END_2D
            END SELECT
            !
            zwi(:,1) = 0._wp
            DO_1Di( 0, 0 )     !* Surface boundary conditions
               zwd(ji,1) = 1._wp - zws(ji,1)
            END_1D
         ENDIF
         !
         !
         !              !==  Apply semi-implicit bottom friction  ==!
         !
         !     Only needed for semi-implicit bottom friction setup. The explicit
         !     bottom friction has been included in "u(v)a" which act as the R.H.S
         !     column vector of the tri-diagonal matrix equation
         !
         IF ( ln_drgimp ) THEN      ! implicit bottom friction
            DO_1Di( 0, 0 )
               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
               zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )   &
                  &                                    / e3u(ji,jj,iku,Kaa)
            END_1D
            IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
               DO_1Di( 0, 0 )
               !!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 
                  zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )    &
                     &                                    / e3u(ji,jj,iku,Kaa)
               END_1D
            ENDIF
         ENDIF
         !
         ! Matrix inversion starting from the first level
         !-----------------------------------------------------------------------
         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
         !
         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
         !        (        ...               )( ...  ) ( ...  )
         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
         !
         !   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 puu(:,:,:,Kaa)
         !-----------------------------------------------------------------------
         !
         DO_2Dik( 0, 0,    2, jpkm1, 1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
            zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
         END_2D
         !
         DO_1Di( 0, 0 )                    !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
#if defined key_RK3
            !                                  ! RK3: use only utau (not utau_b)
            puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj)   &
                 &                                 / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else
            !                                  ! MLF: average of utau and utau_b
            puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) )   &
                 &                                 / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#endif
         END_1D
         DO_2Dik( 0, 0,     2, jpkm1, 1 )
            puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * puu(ji,jj,jk-1,Kaa)
         END_2D
         !
         DO_1Di( 0, 0 )                    !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
            puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
         END_1D
         DO_2Dik( 0, 0,    jpk-2, 1, -1 )
            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
         END_2D
         !
         !
         !              !==  Vertical diffusion on v  ==!
         !
         !                       !* Matrix construction
         IF( ln_zad_Aimp ) THEN   !!
            SELECT CASE( nldf_dyn )
            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
               DO_2Dik( 0, 0,    1, jpkm1, 1 )
                  z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
                  zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
                     &         / e3vw(ji,jj,jk  ,Kmm) * z1_e3va * wvmask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
                     &         / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
                  zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * z1_e3va
                  zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
                  zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
                  zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
                  zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
               END_2D
            CASE DEFAULT               ! iso-level lateral mixing
               DO_2Dik( 0, 0,    1, jpkm1, 1 )
                  z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
                  zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
                     &         / e3vw(ji,jj,jk  ,Kmm) * z1_e3va * wvmask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
                     &         / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
                  zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * z1_e3va
                  zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
                  zwi(ji,jk) = zzwi  + zDt_2 * MIN( zWvi, 0._wp )
                  zws(ji,jk) = zzws  - zDt_2 * MAX( zWvs, 0._wp )
                  zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
               END_2D
            END SELECT
            DO_1Di( 0, 0 )   !* Surface boundary conditions
               zwi(ji,1) = 0._wp
               zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
                  &         / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
               zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
               zws(ji,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
               zwd(ji,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
            END_1D
         ELSE
            SELECT CASE( nldf_dyn )
            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
               DO_2Dik( 0, 0,    1, jpkm1, 1 )
                  zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
                     &         / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
                     &         / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                  zwi(ji,jk) = zzwi
                  zws(ji,jk) = zzws
                  zwd(ji,jk) = 1._wp - zzwi - zzws
               END_2D
            CASE DEFAULT               ! iso-level lateral mixing
               DO_2Dik( 0, 0,    1, jpkm1, 1 )
                  zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
                     &         / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
                  zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
                     &         / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                  zwi(ji,jk) = zzwi
                  zws(ji,jk) = zzws
                  zwd(ji,jk) = 1._wp - zzwi - zzws
               END_2D
            END SELECT
            DO_1Di( 0, 0 )        !* Surface boundary conditions
               zwi(ji,1) = 0._wp
               zwd(ji,1) = 1._wp - zws(ji,1)
            END_1D
         ENDIF
         !
         !              !==  Apply semi-implicit top/bottom friction  ==!
         !
         !     Only needed for semi-implicit bottom friction setup. The explicit
         !     bottom friction has been included in "u(v)a" which act as the R.H.S
         !     column vector of the tri-diagonal matrix equation
         !
         IF( ln_drgimp ) THEN
            DO_1Di( 0, 0 )
               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
               zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )   &
                  &                                   / e3v(ji,jj,ikv,Kaa)
            END_1D
            IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
               DO_1Di( 0, 0 )
                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
                  zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )   &
                     &                                   / e3v(ji,jj,ikv,Kaa)
               END_1D
            ENDIF
         ENDIF
   
         ! Matrix inversion
         !-----------------------------------------------------------------------
         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
         !
         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
         !        (        ...               )( ...  ) ( ...  )
         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
         !
         !   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
         !-----------------------------------------------------------------------
         !
         DO_2Dik( 0, 0,    2, jpkm1, 1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
            zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
         END_2D
         !
         DO_1Di( 0, 0 )                    !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
#if defined key_RK3
            !                                  ! RK3: use only vtau (not vtau_b)
            pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj)   &
               &                                   / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else
            !                                  ! MLF: average of vtau and vtau_b
            pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtauV(ji,jj) )   &
                 &                                 / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#endif
         END_1D
         DO_2Dik( 0, 0,    2, jpkm1, 1 )
            pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * pvv(ji,jj,jk-1,Kaa)
         END_2D
         !
         DO_1Di( 0, 0 )                    !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
            pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
         END_1D
         DO_2Dik( 0, 0,   jpk-2, 1, -1 )
            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
         END_2D
         !                                            ! ================= !
      END_1D                                          !  i-k slices loop  !      
      !                                               ! ================= !
      !
      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) )*r1_Dt - ztrdv(:,:,:)
         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
         DEALLOCATE( ztrdu, ztrdv ) 
      ENDIF
      !                                          ! print mean trends (used for debugging)
      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
         !
      IF( ln_timing )   CALL timing_stop('dyn_zdf')
      !
   END SUBROUTINE dyn_zdf

   !!==============================================================================
END MODULE dynzdf
