#include "cppdefs.h"
      MODULE rp_step3d_uv_mod
#if defined TL_IOMS && defined SOLVE3D
!
!svn $Id: rp_step3d_uv.F 995 2020-01-10 04:01:28Z arango $
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group       Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  This routine time-steps the representers tangent linear horizontal  !
!  momentum equations.  The vertical viscosity terms are time-stepped  !
!  using an implicit algorithm.                                        !
!                                                                      !
!  BASIC STATE variables needed:  u, v, ru, rv, Akv,                   !
!                                 Hz, Huon, Hvom, z_r, z_w,            !
!                                 DU_avg1, DU_avg2, DV_avg1, DV_avg2,  !
!                                 Qsrc                                 !
!                                                                      !
!=======================================================================
!
      implicit none
!
      PRIVATE
      PUBLIC  :: rp_step3d_uv
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE rp_step3d_uv (ng, tile)
!***********************************************************************
!
      USE mod_param
      USE mod_coupling
# ifdef DIAGNOSTICS_UV
!!    USE mod_diags
# endif
      USE mod_forces
      USE mod_grid
      USE mod_mixing
      USE mod_ocean
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, iRPM, 34, __LINE__, __FILE__)
# endif
      CALL rp_step3d_uv_tile (ng, tile,                                 &
     &                        LBi, UBi, LBj, UBj,                       &
     &                        IminS, ImaxS, JminS, JmaxS,               &
     &                        nrhs(ng), nstp(ng), nnew(ng),             &
# ifdef MASKING
     &                        GRID(ng) % umask,                         &
     &                        GRID(ng) % vmask,                         &
# endif
# ifdef WET_DRY_NOT_YET
     &                        GRID(ng) % umask_wet,                     &
     &                        GRID(ng) % vmask_wet,                     &
# endif
     &                        GRID(ng) % om_v,                          &
     &                        GRID(ng) % on_u,                          &
     &                        GRID(ng) % pm,                            &
     &                        GRID(ng) % pn,                            &
     &                        GRID(ng) % Hz,                            &
     &                        GRID(ng) % tl_Hz,                         &
     &                        GRID(ng) % z_r,                           &
     &                        GRID(ng) % tl_z_r,                        &
     &                        GRID(ng) % z_w,                           &
     &                        GRID(ng) % tl_z_w,                        &
     &                        MIXING(ng) % Akv,                         &
     &                        MIXING(ng) % tl_Akv,                      &
     &                        COUPLING(ng) % DU_avg1,                   &
     &                        COUPLING(ng) % tl_DU_avg1,                &
     &                        COUPLING(ng) % DV_avg1,                   &
     &                        COUPLING(ng) % tl_DV_avg1,                &
     &                        COUPLING(ng) % DU_avg2,                   &
     &                        COUPLING(ng) % tl_DU_avg2,                &
     &                        COUPLING(ng) % DV_avg2,                   &
     &                        COUPLING(ng) % tl_DV_avg2,                &
     &                        OCEAN(ng) % tl_ru,                        &
     &                        OCEAN(ng) % tl_rv,                        &
# ifdef DIAGNOSTICS_UV
!!   &                        DIAGS(ng) % DiaU2wrk,                     &
!!   &                        DIAGS(ng) % DiaV2wrk,                     &
!!   &                        DIAGS(ng) % DiaU2int,                     &
!!   &                        DIAGS(ng) % DiaV2int,                     &
!!   &                        DIAGS(ng) % DiaU3wrk,                     &
!!   &                        DIAGS(ng) % DiaV3wrk,                     &
!!   &                        DIAGS(ng) % DiaRU,                        &
!!   &                        DIAGS(ng) % DiaRV,                        &
# endif
     &                        OCEAN(ng) % u,                            &
     &                        OCEAN(ng) % tl_u,                         &
     &                        OCEAN(ng) % v,                            &
     &                        OCEAN(ng) % tl_v,                         &
     &                        OCEAN(ng) % tl_ubar,                      &
     &                        OCEAN(ng) % tl_vbar,                      &
# ifdef NEARSHORE_MELLOR
     &                        OCEAN(ng) % ubar_stokes,                  &
     &                        OCEAN(ng) % tl_ubar_stokes,               &
     &                        OCEAN(ng) % vbar_stokes,                  &
     &                        OCEAN(ng) % tl_vbar_stokes,               &
     &                        OCEAN(ng) % u_stokes,                     &
     &                        OCEAN(ng) % tl_u_stokes,                  &
     &                        OCEAN(ng) % v_stokes,                     &
     &                        OCEAN(ng) % tl_v_stokes,                  &
# endif
     &                        GRID(ng) % Huon,                          &
     &                        GRID(ng) % tl_Huon,                       &
     &                        GRID(ng) % Hvom,                          &
     &                        GRID(ng) % tl_Hvom)
# ifdef PROFILE
      CALL wclock_off (ng, iRPM, 34, __LINE__, __FILE__)
# endif

      RETURN
      END SUBROUTINE rp_step3d_uv
!
!***********************************************************************
      SUBROUTINE rp_step3d_uv_tile (ng, tile,                           &
     &                              LBi, UBi, LBj, UBj,                 &
     &                              IminS, ImaxS, JminS, JmaxS,         &
     &                              nrhs, nstp, nnew,                   &
# ifdef MASKING
     &                              umask, vmask,                       &
# endif
# ifdef WET_DRY_NOT_YET
     &                              umask_wet, vmask_wet,               &
# endif
     &                              om_v, on_u, pm, pn,                 &
     &                              Hz, tl_Hz,                          &
     &                              z_r, tl_z_r,                        &
     &                              z_w, tl_z_w,                        &
     &                              Akv, tl_Akv,                        &
     &                              DU_avg1, tl_DU_avg1,                &
     &                              DV_avg1, tl_DV_avg1,                &
     &                              DU_avg2, tl_DU_avg2,                &
     &                              DV_avg2, tl_DV_avg2,                &
     &                              tl_ru, tl_rv,                       &
# ifdef DIAGNOSTICS_UV
!!   &                              DiaU2wrk, DiaV2wrk,                 &
!!   &                              DiaU2int, DiaV2int,                 &
!!   &                              DiaU3wrk, DiaV3wrk,                 &
!!   &                              DiaRU,    DiaRV,                    &
# endif
     &                              u, tl_u,                            &
     &                              v, tl_v,                            &
     &                              tl_ubar, tl_vbar,                   &
# ifdef NEARSHORE_MELLOR
     &                              ubar_stokes, tl_ubar_stokes,        &
     &                              vbar_stokes, tl_vbar_stokes,        &
     &                              u_stokes, tl_u_stokes,              &
     &                              v_stokes, tl_v_stokes,              &
# endif
     &                              Huon, tl_Huon,                      &
     &                              Hvom, tl_Hvom)
!***********************************************************************
!
      USE mod_param
      USE mod_ncparam
      USE mod_scalars
      USE mod_sources
!
      USE exchange_2d_mod
      USE exchange_3d_mod
# ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d
# endif
      USE rp_u3dbc_mod, ONLY : rp_u3dbc_tile
      USE rp_v3dbc_mod, ONLY : rp_v3dbc_tile
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: nrhs, nstp, nnew
!
# ifdef ASSUMED_SHAPE
#  ifdef MASKING
      real(r8), intent(in) :: umask(LBi:,LBj:)
      real(r8), intent(in) :: vmask(LBi:,LBj:)
#  endif
#  ifdef WET_DRY_NOT_YET
      real(r8), intent(in) :: umask_wet(LBi:,LBj:)
      real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
#  endif
      real(r8), intent(in) :: om_v(LBi:,LBj:)
      real(r8), intent(in) :: on_u(LBi:,LBj:)
      real(r8), intent(in) :: pm(LBi:,LBj:)
      real(r8), intent(in) :: pn(LBi:,LBj:)
      real(r8), intent(in) :: Hz(LBi:,LBj:,:)
      real(r8), intent(in) :: z_r(LBi:,LBj:,:)
      real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
      real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
      real(r8), intent(in) :: DU_avg1(LBi:,LBj:)
      real(r8), intent(in) :: DV_avg1(LBi:,LBj:)
      real(r8), intent(in) :: DU_avg2(LBi:,LBj:)
      real(r8), intent(in) :: DV_avg2(LBi:,LBj:)
      real(r8), intent(in) :: u(LBi:,LBj:,:,:)
      real(r8), intent(in) :: v(LBi:,LBj:,:,:)
#  ifdef NEARSHORE_MELLOR
      real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
      real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
      real(r8), intent(in) :: tl_ubar_stokes(LBi:,LBj:)
      real(r8), intent(in) :: tl_vbar_stokes(LBi:,LBj:)
#  endif
      real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
      real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
      real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
      real(r8), intent(in) :: tl_Akv(LBi:,LBj:,0:)
      real(r8), intent(in) :: tl_DU_avg1(LBi:,LBj:)
      real(r8), intent(in) :: tl_DV_avg1(LBi:,LBj:)
      real(r8), intent(in) :: tl_DU_avg2(LBi:,LBj:)
      real(r8), intent(in) :: tl_DV_avg2(LBi:,LBj:)
      real(r8), intent(in) :: tl_ru(LBi:,LBj:,0:,:)
      real(r8), intent(in) :: tl_rv(LBi:,LBj:,0:,:)

#  ifdef DIAGNOSTICS_UV
!!    real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
!!    real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
!!    real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
!!    real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
!!    real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
!!    real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
!!    real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
!!    real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
#  endif
      real(r8), intent(inout) :: Huon(LBi:,LBj:,:)
      real(r8), intent(inout) :: Hvom(LBi:,LBj:,:)
      real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
#  ifdef NEARSHORE_MELLOR
      real(r8), intent(inout) :: u_stokes(LBi:,LBj:,:)
      real(r8), intent(inout) :: v_stokes(LBi:,LBj:,:)
      real(r8), intent(inout) :: tl_u_stokes(LBi:,LBj:,:)
      real(r8), intent(inout) :: tl_v_stokes(LBi:,LBj:,:)
#  endif
      real(r8), intent(out) :: tl_ubar(LBi:,LBj:,:)
      real(r8), intent(out) :: tl_vbar(LBi:,LBj:,:)
      real(r8), intent(out) :: tl_Huon(LBi:,LBj:,:)
      real(r8), intent(out) :: tl_Hvom(LBi:,LBj:,:)

# else

#  ifdef MASKING
      real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
#  endif
#  ifdef WET_DRY_NOT_YET
      real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: DU_avg1(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: DV_avg1(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: DU_avg2(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: DV_avg2(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
      real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
#  ifdef NEARSHORE_MELLOR
      real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: tl_ubar_stokes(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: tl_vbar_stokes(LBi:UBi,LBj:UBj)
#  endif

      real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: tl_Akv(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: tl_DU_avg1(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: tl_DV_avg1(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: tl_DU_avg2(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: tl_DV_avg2(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: tl_ru(LBi:UBi,LBj:UBj,0:N(ng),2)
      real(r8), intent(in) :: tl_rv(LBi:UBi,LBj:UBj,0:N(ng),2)
#  ifdef DIAGNOSTICS_UV
!!    real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
!!    real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
!!    real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
!!    real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
!!    real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
!!    real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
!!    real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
!!    real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
#  endif
      real(r8), intent(inout) :: Huon(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
      real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
#  ifdef NEARSHORE_MELLOR
      real(r8), intent(inout) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: tl_u_stokes(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: tl_v_stokes(LBi:UBi,LBj:UBj,N(ng))
#  endif
      real(r8), intent(out) :: tl_ubar(LBi:UBi,LBj:UBj,3)
      real(r8), intent(out) :: tl_vbar(LBi:UBi,LBj:UBj,3)
      real(r8), intent(out) :: tl_Huon(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(out) :: tl_Hvom(LBi:UBi,LBj:UBj,N(ng))
# endif
!
!  Local variable declarations.
!
      integer :: i, idiag, is, j, k

      real(r8) :: cff, cff1, cff2
      real(r8) :: tl_cff, tl_cff1, tl_cff2

      real(r8), dimension(IminS:ImaxS) :: CF1
      real(r8), dimension(IminS:ImaxS) :: FC1
# ifdef NEARSHORE_MELLOR
      real(r8), dimension(IminS:ImaxS) :: CFs1
      real(r8), dimension(IminS:ImaxS) :: DCs1
# endif
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: AK
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC1
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
# ifdef NEARSHORE_MELLOR
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CFs
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DCs
# endif
      real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
      real(r8), dimension(IminS:ImaxS,N(ng)) :: oHz
# ifdef DIAGNOSTICS_UV
!!    real(r8), dimension(IminS:ImaxS,0:N(ng)) :: wrk
!!    real(r8), dimension(IminS:ImaxS,1:NDM2d-1) :: Dwrk
# endif
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_AK
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_BC
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_CF
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DC
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FC
# ifdef NEARSHORE_MELLOR
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_CFs
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DCs
# endif
      real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hzk
      real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_oHz
!
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Time step momentum equation in the XI-direction.
!-----------------------------------------------------------------------
!
      DO j=Jstr,Jend
        DO i=IstrU,Iend
          AK(i,0)=0.5_r8*(Akv(i-1,j,0)+                                 &
     &                    Akv(i  ,j,0))
          tl_AK(i,0)=0.5_r8*(tl_Akv(i-1,j,0)+                           &
     &                       tl_Akv(i  ,j,0))
          DO k=1,N(ng)
            AK(i,k)=0.5_r8*(Akv(i-1,j,k)+                               &
     &                      Akv(i  ,j,k))
            tl_AK(i,k)=0.5_r8*(tl_Akv(i-1,j,k)+                         &
     &                         tl_Akv(i  ,j,k))
            Hzk(i,k)=0.5_r8*(Hz(i-1,j,k)+                               &
     &                       Hz(i  ,j,k))
            tl_Hzk(i,k)=0.5_r8*(tl_Hz(i-1,j,k)+                         &
     &                          tl_Hz(i  ,j,k))
# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
            oHz(i,k)=1.0_r8/Hzk(i,k)
            tl_oHz(i,k)=-oHz(i,k)*oHz(i,k)*tl_Hzk(i,k)+                 &
#  ifdef TL_IOMS
     &                  2.0_r8*oHz(i,k)
#  endif
# endif
          END DO
        END DO
!
!  Time step right-hand-side terms.
!
        IF (iic(ng).eq.ntfirst(ng)) THEN
          cff=0.25_r8*dt(ng)
        ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
          cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
        ELSE
          cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
        END IF
        DO i=IstrU,Iend
          DC(i,0)=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
        END DO
!
!  The BASIC STATE "u" used below must be in transport units, but "u"
!  retrived is in m/s so we multiply by "Hzk".
!
        DO k=1,N(ng)
          DO i=IstrU,Iend
!>          u(i,j,k,nnew)=u(i,j,k,nnew)+                                &
!>   &                    DC(i,0)*ru(i,j,k,nrhs)
!>
            tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)+                          &
     &                       DC(i,0)*tl_ru(i,j,k,nrhs)
# ifdef SPLINES_VVISC
!>          u(i,j,k,nnew)=u(i,j,k,nnew)*oHz(i,k)
!>
            tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*oHz(i,k)+                 &
     &                       (u(i,j,k,nnew)*Hzk(i,k))*tl_oHz(i,k)-      &
#  ifdef TL_IOMS
     &                       u(i,j,k,nnew)*Hzk(i,k)*oHz(i,k)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!          DO idiag=1,M3pgrd
!!            DiaU3wrk(i,j,k,idiag)=(DiaU3wrk(i,j,k,idiag)+             &
!!   &                               DC(i,0)*DiaRU(i,j,k,nrhs,idiag))*  &
!!   &                              oHz(i,k)
!!          END DO
#  if defined UV_VIS2 || defined UV_VIS4
!!          DiaU3wrk(i,j,k,M3xvis)=DiaU3wrk(i,j,k,M3xvis)*oHz(i,k)
!!          DiaU3wrk(i,j,k,M3yvis)=DiaU3wrk(i,j,k,M3yvis)*oHz(i,k)
!!          DiaU3wrk(i,j,k,M3hvis)=DiaU3wrk(i,j,k,M3hvis)*oHz(i,k)
#  endif
!!          DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)*oHz(i,k)
#  ifdef BODYFORCE
!!          DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+              &
!!   &                             DC(i,0)*DiaRU(i,j,k,nrhs,M3vvis)*    &
!!   &                             oHz(i,k)
#  endif
!!          DiaU3wrk(i,j,k,M3rate)=DiaU3wrk(i,j,k,M3rate)*oHz(i,k)
# endif
          END DO
        END DO

# ifdef SPLINES_VVISC
!
!  Apply spline algorithim to BASIC STATE "u" which should be
!  in units of velocity (m/s). DC will be used in the tangent
!  linear spline algorithm below.  Solve BASIC STATE tridiagonal
!  system.
!
        cff1=1.0_r8/6.0_r8
        DO k=1,N(ng)-1
          DO i=IstrU,Iend
            FC(i,k)=cff1*Hzk(i,k  )-dt(ng)*AK(i,k-1)*oHz(i,k  )
            CF(i,k)=cff1*Hzk(i,k+1)-dt(ng)*AK(i,k+1)*oHz(i,k+1)
          END DO
        END DO
        DO i=IstrU,Iend
          CF(i,0)=0.0_r8
          DC(i,0)=0.0_r8
        END DO
!
!  LU decomposition and forward substitution.
!
        cff1=1.0_r8/3.0_r8
        DO k=1,N(ng)-1
          DO i=IstrU,Iend
            BC(i,k)=cff1*(Hzk(i,k)+Hzk(i,k+1))+                         &
     &              dt(ng)*AK(i,k)*(oHz(i,k)+oHz(i,k+1))
            cff=1.0_r8/(BC(i,k)-FC(i,k)*CF(i,k-1))
            CF(i,k)=cff*CF(i,k)
            DC(i,k)=cff*(u(i,j,k+1,nnew)-u(i,j,k,nnew)-                 &
     &                   FC(i,k)*DC(i,k-1))
          END DO
        END DO
!
!  Backward substitution. Save DC for the tangent linearization.
!  DC is scaled later by AK.
!
        DO i=IstrU,Iend
          DC(i,N(ng))=0.0_r8
        END DO
        DO k=N(ng)-1,1,-1
          DO i=IstrU,Iend
            DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
          END DO
        END DO
!
!  Use conservative, parabolic spline reconstruction of tangent
!  linear vertical viscosity derivatives.  Then, time step vertical
!  viscosity term implicitly by solving a tridiagonal system.
#  ifdef TL_IOMS
!  If tl_AK is computed, then tl_FC=tl_FC+dt*AK*oHz and
!                             tl_CF=tl_CF+dt*AK*oHz below.
#  endif
!
        cff1=1.0_r8/6.0_r8
        DO k=1,N(ng)-1
          DO i=IstrU,Iend
            FC(i,k)=cff1*Hzk(i,k  )-                                    &
     &              dt(ng)*AK(i,k-1)*oHz(i,k  )
            tl_FC(i,k)=cff1*tl_Hzk(i,k  )-                              &
     &                 dt(ng)*(tl_AK(i,k-1)*oHz(i,k  )+                 &
     &                         AK(i,k-1)*tl_oHz(i,k  ))
            CF(i,k)=cff1*Hzk(i,k+1)-                                    &
     &               dt(ng)*AK(i,k+1)*oHz(i,k+1)
            tl_CF(i,k)=cff1*tl_Hzk(i,k+1)-                              &
     &                 dt(ng)*(tl_AK(i,k+1)*oHz(i,k+1)+                 &
     &                         AK(i,k+1)*tl_oHz(i,k+1))
          END DO
        END DO
        DO i=IstrU,Iend
          CF(i,0)=0.0_r8
          tl_CF(i,0)=0.0_r8
          tl_DC(i,0)=0.0_r8
        END DO
!
!  Tangent linear LU decomposition and forward substitution.
#  ifdef TL_IOMS
!  If tl_AK is computed, then tl_BC=tl_FC-dt*AK*oHz below.
#  endif
!
        cff1=1.0_r8/3.0_r8
        DO k=1,N(ng)-1
          DO i=IstrU,Iend
            BC(i,k)=cff1*(Hzk(i,k)+Hzk(i,k+1))+                         &
     &              dt(ng)*AK(i,k)*(oHz(i,k)+oHz(i,k+1))
            tl_BC(i,k)=cff1*(tl_Hzk(i,k)+tl_Hzk(i,k+1))+                &
     &                 dt(ng)*(tl_AK(i,k)*(oHz(i,k)+oHz(i,k+1))+        &
     &                         AK(i,k)*(tl_oHz(i,k)+tl_oHz(i,k+1)))
            cff=1.0_r8/(BC(i,k)-FC(i,k)*CF(i,k-1))
#  ifdef TL_IOMS
            tl_DC(i,k)=cff*(tl_u(i,j,k+1,nnew)-                         &
     &                      tl_u(i,j,k  ,nnew)-                         &
     &                     ((tl_FC(i,k)-FC(i,k))*DC(i,k-1)+             &
     &                      (tl_BC(i,k)-BC(i,k))*DC(i,k  )+             &
     &                      (tl_CF(i,k)-CF(i,k))*DC(i,k+1))-            &
     &                     FC(i,k)*tl_DC(i,k-1))
            CF(i,k)=cff*CF(i,k)
#  else
            CF(i,k)=cff*CF(i,k)
            tl_DC(i,k)=cff*(tl_u(i,j,k+1,nnew)-                         &
     &                      tl_u(i,j,k  ,nnew)-                         &
     &                      (tl_FC(i,k)*DC(i,k-1)+                      &
     &                       tl_BC(i,k)*DC(i,k  )+                      &
     &                       tl_CF(i,k)*DC(i,k+1))-                     &
     &                      FC(i,k)*tl_DC(i,k-1))
#  endif
          END DO
        END DO
!
!  Tangent linear backward substitution.
!
        DO i=IstrU,Iend
          tl_DC(i,N(ng))=0.0_r8
        END DO
        DO k=N(ng)-1,1,-1
          DO i=IstrU,Iend
            tl_DC(i,k)=tl_DC(i,k)-CF(i,k)*tl_DC(i,k+1)
          END DO
        END DO
!
!  Compute tl_DC before multiplying BASIC STATE spline gradients
!  DC by AK.
!
        DO k=1,N(ng)
          DO i=IstrU,Iend
            tl_DC(i,k)=tl_DC(i,k)*AK(i,k)+DC(i,k)*tl_AK(i,k)
#  ifdef TL_IOMS
!           tl_DC(i,k)=tl_DC(i,k)-                                      &
!     &                DC(i,k)*AK(i,k)    ! needed if tl_AK is computed
#  endif
            DC(i,k)=DC(i,k)*AK(i,k)
!>          cff=dt(ng)*oHz(i,k)*(DC(i,k)-DC(i,k-1))
!>
            tl_cff=dt(ng)*(tl_oHz(i,k)*(DC(i,k)-DC(i,k-1))+             &
     &                     oHz(i,k)*(tl_DC(i,k)-tl_DC(i,k-1)))-         &
#  ifdef TL_IOMS
     &             dt(ng)*oHz(i,k)*(DC(i,k)-DC(i,k-1))
#  endif
!>          u(i,j,k,nnew)=u(i,j,k,nnew)+cff
!>
            tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)+tl_cff
#  ifdef DIAGNOSTICS_UV
!!          DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+cff
#  endif
          END DO
        END DO
# else
!
!  Compute off-diagonal coefficients [lambda*dt*Akv/Hz] for the
!  implicit vertical viscosity term at horizontal U-points and
!  vertical W-points.
!
!  NOTE: The original code solves the tridiagonal system A*u=r where
!        A is a matrix and u and r are vectors. We need to solve the
!        tangent linear of this system which is A*tl_u+tl_A*u=tl_r.
!        Here, tl_A*u and tl_r are known, so we must solve for tl_u
!        by inverting A*tl_u=tl_r-tl_A*u.
!
        cff=-lambda*dt(ng)/0.5_r8
        DO k=1,N(ng)-1
          DO i=IstrU,Iend
            cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i-1,j,k+1)-                   &
     &                   z_r(i,j,k  )-z_r(i-1,j,k  ))
            tl_cff1=-cff1*cff1*(tl_z_r(i,j,k+1)+tl_z_r(i-1,j,k+1)-      &
     &                          tl_z_r(i,j,k  )-tl_z_r(i-1,j,k  ))+     &
#  ifdef TL_IOMS
     &              2.0_r8*cff1
#  endif
            FC(i,k)=cff*cff1*AK(i,k)
            tl_FC(i,k)=cff*(tl_cff1*AK(i,k)+cff1*tl_AK(i,k))
#  ifdef TL_IOMS
!
!  Uncomment if tl_Akv is computed.
!
!           tl_FC(i,k)=tl_FC(i,k)-FC(i,k)
#  endif
          END DO
        END DO
        DO i=IstrU,Iend
          FC(i,0)=0.0_r8
          tl_FC(i,0)=0.0_r8
          FC(i,N(ng))=0.0_r8
          tl_FC(i,N(ng))=0.0_r8
        END DO
!
!  Solve the tangent linear tridiagonal system.
!  (DC is a tangent linear variable here).
!
        DO k=1,N(ng)
          DO i=IstrU,Iend
            BC(i,k)=Hzk(i,k)-FC(i,k)-FC(i,k-1)
            tl_BC(i,k)=tl_Hzk(i,k)-tl_FC(i,k)-tl_FC(i,k-1)
          END DO
        END DO
        DO k=2,N(ng)-1
          DO i=IstrU,Iend
#  ifdef TL_IOMS
            DC(i,k)=tl_u(i,j,k,nnew)-                                   &
     &              ((tl_FC(i,k-1)-FC(i,k-1))*u(i,j,k-1,nnew)+          &
     &               (tl_BC(i,k  )-BC(i,k  ))*u(i,j,k  ,nnew)+          &
     &               (tl_FC(i,k  )-FC(i,k  ))*u(i,j,k+1,nnew))
#  else
            DC(i,k)=tl_u(i,j,k,nnew)-                                   &
     &              (tl_FC(i,k-1)*u(i,j,k-1,nnew)+                      &
     &               tl_BC(i,k  )*u(i,j,k  ,nnew)+                      &
     &               tl_FC(i,k  )*u(i,j,k+1,nnew))
#  endif
          END DO
        END DO
        DO i=IstrU,Iend
#  ifdef TL_IOMS
          DC(i,1)=tl_u(i,j,1,nnew)-                                     &
     &            ((tl_BC(i,1)-BC(i,1))*u(i,j,1,nnew)+                  &
     &             (tl_FC(i,1)-FC(i,1))*u(i,j,2,nnew))
          DC(i,N(ng))=tl_u(i,j,N(ng),nnew)-                             &
     &                ((tl_FC(i,N(ng)-1)-FC(i,N(ng)-1))*                &
     &                 u(i,j,N(ng)-1,nnew)+                             &
     &                 (tl_BC(i,N(ng)  )-BC(i,N(ng)  ))*                &
     &                 u(i,j,N(ng)  ,nnew))
#  else
          DC(i,1)=tl_u(i,j,1,nnew)-                                     &
     &            (tl_BC(i,1)*u(i,j,1,nnew)+                            &
     &             tl_FC(i,1)*u(i,j,2,nnew))
          DC(i,N(ng))=tl_u(i,j,N(ng),nnew)-                             &
     &                (tl_FC(i,N(ng)-1)*u(i,j,N(ng)-1,nnew)+            &
     &                 tl_BC(i,N(ng)  )*u(i,j,N(ng)  ,nnew))
#  endif
        END DO
        DO i=IstrU,Iend
          cff=1.0_r8/BC(i,1)
          CF(i,1)=cff*FC(i,1)
          DC(i,1)=cff*DC(i,1)
        END DO
        DO k=2,N(ng)-1
          DO i=IstrU,Iend
            cff=1.0_r8/(BC(i,k)-FC(i,k-1)*CF(i,k-1))
            CF(i,k)=cff*FC(i,k)
            DC(i,k)=cff*(DC(i,k)-FC(i,k-1)*DC(i,k-1))
          END DO
        END DO
!
!  Compute new solution by back substitution.
!  (DC is a tangent linear variable here).
!
        DO i=IstrU,Iend
#  ifdef DIAGNOSTICS_UV
!!        wrk(i,N(ng))=u(i,j,N(ng),nnew)*oHz(i,N(ng))
#  endif
          DC(i,N(ng))=(DC(i,N(ng))-FC(i,N(ng)-1)*DC(i,N(ng)-1))/        &
     &                (BC(i,N(ng))-FC(i,N(ng)-1)*CF(i,N(ng)-1))
!>        u(i,j,N(ng),nnew)=DC(i,N(ng))
!>
          tl_u(i,j,N(ng),nnew)=DC(i,N(ng))
#  ifdef DIAGNOSTICS_UV
!!        DiaU3wrk(i,j,N(ng),M3vvis)=DiaU3wrk(i,j,N(ng),M3vvis)+        &
!!   &                               u(i,j,N(ng),nnew)-wrk(i,N(ng))
#  endif
        END DO
        DO k=N(ng)-1,1,-1
          DO i=IstrU,Iend
#  ifdef DIAGNOSTICS_UV
!!          wrk(i,k)=u(i,j,k,nnew)*oHz(i,k)
#  endif
            DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
!>          u(i,j,k,nnew)=DC(i,k)
!>
            tl_u(i,j,k,nnew)=DC(i,k)
#  ifdef DIAGNOSTICS_UV
!!          DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+              &
!!   &                             u(i,j,k,nnew)-wrk(i,k)
#  endif
          END DO
        END DO
# endif
!
!  Replace INTERIOR POINTS incorrect vertical mean with more accurate
!  barotropic component, ubar=DU_avg1/(D*on_u). Recall that, D=CF(:,0).
!
        DO i=IstrU,Iend
          CF(i,0)=Hzk(i,1)
          tl_CF(i,0)=tl_Hzk(i,1)
          DC(i,0)=u(i,j,1,nnew)*Hzk(i,1)
          tl_DC(i,0)=tl_u(i,j,1,nnew)*Hzk(i,1)+                         &
     &               u(i,j,1,nnew)*tl_Hzk(i,1)-                         &
# ifdef TL_IOMS
     &               DC(i,0)
# endif
# ifdef NEARSHORE_MELLOR
          DCs(i,0)=u_stokes(i,j,1)*Hzk(i,1)
          tl_DCs(i,0)=tl_u_stokes(i,j,1)*Hzk(i,1)+                      &
     &                u_stokes(i,j,1)*tl_Hzk(i,1)-                      &
#  ifdef TL_IOMS
     &               DCs(i,0)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!        Dwrk(i,M2pgrd)=DiaU3wrk(i,j,1,M3pgrd)*Hzk(i,1)
!!        Dwrk(i,M2bstr)=DiaU3wrk(i,j,1,M3vvis)*Hzk(i,1)
#  ifdef UV_COR
!!        Dwrk(i,M2fcor)=DiaU3wrk(i,j,1,M3fcor)*Hzk(i,1)
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
!!        Dwrk(i,M2xvis)=DiaU3wrk(i,j,1,M3xvis)*Hzk(i,1)
!!        Dwrk(i,M2yvis)=DiaU3wrk(i,j,1,M3yvis)*Hzk(i,1)
!!        Dwrk(i,M2hvis)=DiaU3wrk(i,j,1,M3hvis)*Hzk(i,1)
#  endif
#  ifdef UV_ADV
!!        Dwrk(i,M2xadv)=DiaU3wrk(i,j,1,M3xadv)*Hzk(i,1)
!!        Dwrk(i,M2yadv)=DiaU3wrk(i,j,1,M3yadv)*Hzk(i,1)
!!        Dwrk(i,M2hadv)=DiaU3wrk(i,j,1,M3hadv)*Hzk(i,1)
#  endif
#  ifdef NEARSHORE_MELLOR
!!        Dwrk(i,M2hrad)=DiaU3wrk(i,j,1,M3hrad)*Hzk(i,1)
#  endif
# endif
        END DO
        DO k=2,N(ng)
          DO i=IstrU,Iend
            CF(i,0)=CF(i,0)+Hzk(i,k)
            tl_CF(i,0)=tl_CF(i,0)+tl_Hzk(i,k)
            DC(i,0)=DC(i,0)+u(i,j,k,nnew)*Hzk(i,k)
            tl_DC(i,0)=tl_DC(i,0)+                                      &
     &                 tl_u(i,j,k,nnew)*Hzk(i,k)+                       &
     &                 u(i,j,k,nnew)*tl_Hzk(i,k)-                       &
# ifdef TL_IOMS
     &                 u(i,j,k,nnew)*Hzk(i,k)
# endif
# ifdef NEARSHORE_MELLOR
            DCs(i,0)=DCs(i,0)+u_stokes(i,j,k)*Hzk(i,k)
            tl_DCs(i,0)=tl_DCs(i,0)+                                    &
     &                  tl_u_stokes(i,j,k)*Hzk(i,k)+                    &
     &                  u_stokes(i,j,k)*tl_Hzk(i,k)-                    &
#  ifdef TL_IOMS
     &                  u_stokes(i,j,k)*Hzk(i,k)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!          Dwrk(i,M2pgrd)=Dwrk(i,M2pgrd)+                              &
!!   &                     DiaU3wrk(i,j,k,M3pgrd)*Hzk(i,k)
!!          Dwrk(i,M2bstr)=Dwrk(i,M2bstr)+                              &
!!   &                     DiaU3wrk(i,j,k,M3vvis)*Hzk(i,k)
#  ifdef UV_COR
!!          Dwrk(i,M2fcor)=Dwrk(i,M2fcor)+                              &
!!   &                     DiaU3wrk(i,j,k,M3fcor)*Hzk(i,k)
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
!!          Dwrk(i,M2xvis)=Dwrk(i,M2xvis)+                              &
!!   &                     DiaU3wrk(i,j,k,M3xvis)*Hzk(i,k)
!!          Dwrk(i,M2yvis)=Dwrk(i,M2yvis)+                              &
!!   &                     DiaU3wrk(i,j,k,M3yvis)*Hzk(i,k)
!!          Dwrk(i,M2hvis)=Dwrk(i,M2hvis)+                              &
!!   &                     DiaU3wrk(i,j,k,M3hvis)*Hzk(i,k)
#  endif
#  ifdef UV_ADV
!!          Dwrk(i,M2xadv)=Dwrk(i,M2xadv)+                              &
!!   &                     DiaU3wrk(i,j,k,M3xadv)*Hzk(i,k)
!!          Dwrk(i,M2yadv)=Dwrk(i,M2yadv)+                              &
!!   &                     DiaU3wrk(i,j,k,M3yadv)*Hzk(i,k)
!!          Dwrk(i,M2hadv)=Dwrk(i,M2hadv)+                              &
!!   &                     DiaU3wrk(i,j,k,M3hadv)*Hzk(i,k)
#  endif
#  ifdef NEARSHORE_MELLOR
!!          Dwrk(i,M2hrad)=Dwrk(i,M2hrad)+                              &
!!   &                     DiaU3wrk(i,j,k,M3hrad)*Hzk(i,k)
#  endif
# endif
          END DO
        END DO
        DO i=IstrU,Iend
          DC1(i,0)=DC(i,0)                                ! intermediate
          cff1=1.0_r8/(CF(i,0)*on_u(i,j))
          tl_cff1=-cff1*cff1*tl_CF(i,0)*on_u(i,j)+                      &
# ifdef TL_IOMS
     &            2.0_r8*cff1
# endif
          DC(i,0)=(DC(i,0)*on_u(i,j)-DU_avg1(i,j))*cff1      ! recursive
          tl_DC(i,0)=(tl_DC(i,0)*on_u(i,j)-tl_DU_avg1(i,j))*cff1+       &
     &               (DC1(i,0)*on_u(i,j)-DU_avg1(i,j))*tl_cff1-         &
# ifdef TL_IOMS
     &               DC(i,0)
# endif
# ifdef NEARSHORE_MELLOR
          DCs1(i)=DCs(i,0)                                ! intermediate
          cff2=1.0_r8/CF(i,0)
          tl_cff2=-cff2*cff2*tl_CF(i,0)+                                &
#  ifdef TL_IOMS
     &            2.0_r8*cff2
#  endif
          DCs(i,0)=DCs(i,0)*cff2-ubar_stokes(i,j)            ! recursive
          tl_DCs(i,0)=tl_DCs(i,0)*cff2+                                 &
     &                DCs1(i)*tl_cff2-tl_ubar_stokes(i,j)-              &
#  ifdef TL_IOMS
     &                DCs(i,0)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!        DO idiag=1,M2pgrd
!!          Dwrk(i,idiag)=(Dwrk(i,idiag)*on_u(i,j)-                     &
!!   &                     DiaU2wrk(i,j,idiag))*cff1
!!        END DO
!!        Dwrk(i,M2bstr)=(Dwrk(i,M2bstr)*on_u(i,j)-                     &
!!   &                    DiaU2wrk(i,j,M2bstr)-DiaU2wrk(i,j,M2sstr))*   &
!!   &                   cff1
# endif
        END DO
!
!  Couple and update new solution.
!
        DO k=1,N(ng)
          DO i=IstrU,Iend
!>          u(i,j,k,nnew)=u(i,j,k,nnew)-DC(i,0)
!>
            tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-tl_DC(i,0)
# ifdef MASKING
!>          u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
!>
            tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*umask(i,j)
# endif
# ifdef WET_DRY_NOT_YET
!>          u(i,j,k,nnew)=u(i,j,k,nnew)*umask_wet(i,j)
!>
            tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*umask_wet(i,j)
# endif
# ifdef NEARSHORE_MELLOR
!>          u_stokes(i,j,k)=u_stokes(i,j,k)-DCs(i,0)
!>
            tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)-tl_DCs(i,0)
#  ifdef MASKING
!>          u_stokes(i,j,k)=u_stokes(i,j,k)*umask(i,j)
!>
            tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*umask(i,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>          u_stokes(i,j,k)=u_stokes(i,j,k)*umask_wet(i,j)
!>
            tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*umask_wet(i,j)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!          DiaU3wrk(i,j,k,M3pgrd)=DiaU3wrk(i,j,k,M3pgrd)-              &
!!   &                             Dwrk(i,M2pgrd)
!!          DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)-              &
!!   &                             Dwrk(i,M2bstr)
#  ifdef UV_COR
!!          DiaU3wrk(i,j,k,M3fcor)=DiaU3wrk(i,j,k,M3fcor)-              &
!!   &                             Dwrk(i,M2fcor)
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
!!          DiaU3wrk(i,j,k,M3xvis)=DiaU3wrk(i,j,k,M3xvis)-              &
!!   &                             Dwrk(i,M2xvis)
!!          DiaU3wrk(i,j,k,M3yvis)=DiaU3wrk(i,j,k,M3yvis)-              &
!!   &                             Dwrk(i,M2yvis)
!!          DiaU3wrk(i,j,k,M3hvis)=DiaU3wrk(i,j,k,M3hvis)-              &
!!   &                             Dwrk(i,M2hvis)
#  endif
#  ifdef UV_ADV
!!          DiaU3wrk(i,j,k,M3xadv)=DiaU3wrk(i,j,k,M3xadv)-              &
!!   &                             Dwrk(i,M2xadv)
!!          DiaU3wrk(i,j,k,M3yadv)=DiaU3wrk(i,j,k,M3yadv)-              &
!!   &                             Dwrk(i,M2yadv)
!!          DiaU3wrk(i,j,k,M3hadv)=DiaU3wrk(i,j,k,M3hadv)-              &
!!   &                             Dwrk(i,M2hadv)
#  endif
#  ifdef NEARSHORE_MELLOR
!!          DiaU3wrk(i,j,k,M3hrad)=DiaU3wrk(i,j,k,M3hrad)-              &
!!   &                             Dwrk(i,M2hrad)
#  endif
# endif
          END DO
        END DO

# if defined DIAGNOSTICS_UV && defined MASKING
!!      DO k=1,N(ng)
!!        DO i=IstrU,Iend
!!          DO idiag=1,NDM3d
!!            DiaU3wrk(i,j,k,idiag)=DiaU3wrk(i,j,k,idiag)*umask(i,j)
!!          END DO
!!        END DO
!!      END DO
# endif
!
!-----------------------------------------------------------------------
!  Time step momentum equation in the ETA-direction.
!-----------------------------------------------------------------------
!
        IF (j.ge.JstrV) THEN
          DO i=Istr,Iend
            AK(i,0)=0.5_r8*(Akv(i,j-1,0)+                               &
     &                      Akv(i,j  ,0))
            tl_AK(i,0)=0.5_r8*(tl_Akv(i,j-1,0)+                         &
     &                         tl_Akv(i,j  ,0))
            DO k=1,N(ng)
              AK(i,k)=0.5_r8*(Akv(i,j-1,k)+                             &
     &                        Akv(i,j  ,k))
              tl_AK(i,k)=0.5_r8*(tl_Akv(i,j-1,k)+                       &
     &                           tl_Akv(i,j  ,k))
              Hzk(i,k)=0.5_r8*(Hz(i,j-1,k)+                             &
     &                         Hz(i,j  ,k))
              tl_Hzk(i,k)=0.5_r8*(tl_Hz(i,j-1,k)+                       &
     &                            tl_Hz(i,j  ,k))
# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
              oHz(i,k)=1.0_r8/Hzk(i,k)
              tl_oHz(i,k)=-oHz(i,k)*oHz(i,k)*tl_Hzk(i,k)+               &
#  ifdef TL_IOMS
     &                    2.0_r8*oHz(i,k)
#  endif
# endif
            END DO
          END DO
!
!  Time step right-hand-side terms.
!
          IF (iic(ng).eq.ntfirst(ng)) THEN
            cff=0.25_r8*dt(ng)
          ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
            cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
          ELSE
            cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
          END IF
          DO i=Istr,Iend
            DC(i,0)=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
          END DO
!
!  The BASIC STATE "v" used below must be in transport units, but "v"
!  retrived is in m/s so we multiply by "Hzk".
!
          DO k=1,N(ng)
            DO i=Istr,Iend
!>            v(i,j,k,nnew)=v(i,j,k,nnew)+DC(i,0)*rv(i,j,k,nrhs)
!>
              tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)+                        &
     &                         DC(i,0)*tl_rv(i,j,k,nrhs)
# ifdef SPLINES_VVISC
!>            v(i,j,k,nnew)=v(i,j,k,nnew)*oHz(i,k)
!>
              tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*oHz(i,k)+               &
     &                         (v(i,j,k,nnew)*Hzk(i,k))*tl_oHz(i,k)-    &
#  ifdef TL_IOMS
     &                         v(i,j,k,nnew)*Hzk(i,k)*oHz(i,k)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!            DO idiag=1,M3pgrd
!!              DiaV3wrk(i,j,k,idiag)=(DiaV3wrk(i,j,k,idiag)+           &
!!   &                                 DC(i,0)*                         &
!!   &                                 DiaRV(i,j,k,nrhs,idiag))*        &
!!   &                                oHz(i,k)
!!            END DO
#  if defined UV_VIS2 || defined UV_VIS4
!!            DiaV3wrk(i,j,k,M3xvis)=DiaV3wrk(i,j,k,M3xvis)*oHz(i,k)
!!            DiaV3wrk(i,j,k,M3yvis)=DiaV3wrk(i,j,k,M3yvis)*oHz(i,k)
!!            DiaV3wrk(i,j,k,M3hvis)=DiaV3wrk(i,j,k,M3hvis)*oHz(i,k)
#  endif
!!            DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)*oHz(i,k)
#  ifdef BODYFORCE
!!            DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+            &
!!   &                               DC(i,0)*DiaRV(i,j,k,nrhs,M3vvis)*  &
!!   &                               oHz(i,k)
#  endif
!!            DiaV3wrk(i,j,k,M3rate)=DiaV3wrk(i,j,k,M3rate)*oHz(i,k)
# endif
            END DO
          END DO

# ifdef SPLINES_VVISC
!
!  Apply spline algorithim to BASIC STATE "v" which should be
!  in units of velocity (m/s). DC will be used in the tangent
!  linear spline algorithm below. Solve BASIC STATE tridiagonal
!  system.
!
          cff1=1.0_r8/6.0_r8
          DO k=1,N(ng)-1
            DO i=Istr,Iend
              FC(i,k)=cff1*Hzk(i,k  )-dt(ng)*AK(i,k-1)*oHz(i,k  )
              CF(i,k)=cff1*Hzk(i,k+1)-dt(ng)*AK(i,k+1)*oHz(i,k+1)
            END DO
          END DO
          DO i=Istr,Iend
            CF(i,0)=0.0_r8
            DC(i,0)=0.0_r8
          END DO
!
!  LU decomposition and forward substitution.
!
          cff1=1.0_r8/3.0_r8
          DO k=1,N(ng)-1
            DO i=Istr,Iend
              BC(i,k)=cff1*(Hzk(i,k)+Hzk(i,k+1))+                       &
     &                dt(ng)*AK(i,k)*(oHz(i,k)+oHz(i,k+1))
              cff=1.0_r8/(BC(i,k)-FC(i,k)*CF(i,k-1))
              CF(i,k)=cff*CF(i,k)
              DC(i,k)=cff*(v(i,j,k+1,nnew)-v(i,j,k,nnew)-               &
     &                     FC(i,k)*DC(i,k-1))
            END DO
          END DO
!
!  Backward substitution. Save DC for the tangent linearization.
!  DC is scaled later by AK.
!
          DO i=Istr,Iend
            DC(i,N(ng))=0.0_r8
          END DO
          DO k=N(ng)-1,1,-1
            DO i=Istr,Iend
              DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
            END DO
          END DO
!
!  Use conservative, parabolic spline reconstruction of tangent
!  linear vertical viscosity derivatives.  Then, time step vertical
!  viscosity term implicitly by solving a tridiagonal system.
#  ifdef TL_IOMS
!  If tl_AK is computed, then tl_FC=tl_FC+dt*AK*oHz and
!                             tl_CF=tl_CF+dt*AK*oHz below.
#  endif
!
          cff1=1.0_r8/6.0_r8
          DO k=1,N(ng)-1
            DO i=Istr,Iend
              FC(i,k)=cff1*Hzk(i,k  )-                                  &
     &                dt(ng)*AK(i,k-1)*oHz(i,k  )
              tl_FC(i,k)=cff1*tl_Hzk(i,k  )-                            &
     &                   dt(ng)*(tl_AK(i,k-1)*oHz(i,k  )+               &
     &                           AK(i,k-1)*tl_oHz(i,k  ))
              CF(i,k)=cff1*Hzk(i,k+1)-                                  &
     &                dt(ng)*AK(i,k+1)*oHz(i,k+1)
              tl_CF(i,k)=cff1*tl_Hzk(i,k+1)-                            &
     &                   dt(ng)*(tl_AK(i,k+1)*oHz(i,k+1)+               &
     &                           AK(i,k+1)*tl_oHz(i,k+1))
            END DO
          END DO
          DO i=Istr,Iend
            CF(i,0)=0.0_r8
            tl_CF(i,0)=0.0_r8
            tl_DC(i,0)=0.0_r8
          END DO
!
!  Tangent linear LU decomposition and forward substitution.
#  ifdef TL_IOMS
!  If tl_AK is computed, then tl_BC=tl_BC-dt*AK*oHz below.
#  endif
!
          cff1=1.0_r8/3.0_r8
          DO k=1,N(ng)-1
            DO i=Istr,Iend
              BC(i,k)=cff1*(Hzk(i,k)+Hzk(i,k+1))+                       &
     &                dt(ng)*AK(i,k)*(oHz(i,k)+oHz(i,k+1))
              tl_BC(i,k)=cff1*(tl_Hzk(i,k)+tl_Hzk(i,k+1))+              &
     &                   dt(ng)*(tl_AK(i,k)*(oHz(i,k)+oHz(i,k+1))+      &
     &                           AK(i,k)*(tl_oHz(i,k)+tl_oHz(i,k+1)))
              cff=1.0_r8/(BC(i,k)-FC(i,k)*CF(i,k-1))
#  ifdef TL_IOMS
              tl_DC(i,k)=cff*(tl_v(i,j,k+1,nnew)-                       &
     &                        tl_v(i,j,k  ,nnew)-                       &
     &                       ((tl_FC(i,k)-FC(i,k))*DC(i,k-1)+           &
     &                        (tl_BC(i,k)-BC(i,k))*DC(i,k  )+           &
     &                        (tl_CF(i,k)-CF(i,k))*DC(i,k+1))-          &
     &                       FC(i,k)*tl_DC(i,k-1))
              CF(i,k)=cff*CF(i,k)
#  else
              CF(i,k)=cff*CF(i,k)
              tl_DC(i,k)=cff*(tl_v(i,j,k+1,nnew)-                       &
     &                        tl_v(i,j,k  ,nnew)-                       &
     &                        (tl_FC(i,k)*DC(i,k-1)+                    &
     &                         tl_BC(i,k)*DC(i,k  )+                    &
     &                         tl_CF(i,k)*DC(i,k+1))-                   &
     &                        FC(i,k)*tl_DC(i,k-1))
#  endif
            END DO
          END DO
!
!  Tangent linear backward substitution.
!
          DO i=Istr,Iend
            tl_DC(i,N(ng))=0.0_r8
          END DO
          DO k=N(ng)-1,1,-1
            DO i=Istr,Iend
              tl_DC(i,k)=tl_DC(i,k)-CF(i,k)*tl_DC(i,k+1)
            END DO
          END DO
!
!  Compute tl_DC before multiplying BASIC STATE spline gradients
!  DC by AK.
!
          DO k=1,N(ng)
            DO i=Istr,Iend
              tl_DC(i,k)=tl_DC(i,k)*AK(i,k)+DC(i,k)*tl_AK(i,k)
#  ifdef TL_IOMS
!             tl_DC(i,k)=tl_DC(i,k)-                                    &
!    &                   DC(i,k)*AK(i,k)  ! needed if tl_AK is computed
#  endif
              DC(i,k)=DC(i,k)*AK(i,k)
!>            cff=dt(ng)*oHz(i,k)*(DC(i,k)-DC(i,k-1))
!>
              tl_cff=dt(ng)*(tl_oHz(i,k)*(DC(i,k)-DC(i,k-1))+           &
     &                       oHz(i,k)*(tl_DC(i,k)-tl_DC(i,k-1)))-       &
#  ifdef TL_IOMS
     &               dt(ng)*oHz(i,k)*(DC(i,k)-DC(i,k-1))
#  endif
!>            v(i,j,k,nnew)=v(i,j,k,nnew)+cff
!>
              tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)+tl_cff
#  ifdef DIAGNOSTICS_UV
!!            DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+cff
#  endif
            END DO
          END DO
# else
!
!  Compute off-diagonal coefficients [lambda*dt*Akv/Hz] for the
!  implicit vertical viscosity term at horizontal V-points and
!  vertical W-points.
!
!  NOTE: The original code solves the tridiagonal system A*v=r where
!        A is a matrix and v and r are vectors. We need to solve the
!        tangent linear of this system which is A*tl_v+tl_A*v=tl_r.
!        Here, tl_A*v and tl_r are known, so we must solve for tl_v
!        by inverting A*tl_v=tl_r-tl_A*v.
!
          cff=-lambda*dt(ng)/0.5_r8
          DO k=1,N(ng)-1
            DO i=Istr,Iend
              cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i,j-1,k+1)-                 &
     &                     z_r(i,j,k  )-z_r(i,j-1,k  ))
              tl_cff1=-cff1*cff1*(tl_z_r(i,j,k+1)+tl_z_r(i,j-1,k+1)-    &
     &                            tl_z_r(i,j,k  )-tl_z_r(i,j-1,k  ))+   &
#  ifdef TL_IOMS
     &                2.0_r8*cff1
#  endif
              FC(i,k)=cff*cff1*AK(i,k)
              tl_FC(i,k)=cff*(tl_cff1*AK(i,k)+cff1*tl_AK(i,k))
#  ifdef TL_IOMS
!
!  Uncomment if tl_Akv is computed.
!
!             tl_FC(i,k)=tl_FC(i,k)-FC(i,k)
#  endif
            END DO
          END DO
          DO i=Istr,Iend
            FC(i,0)=0.0_r8
            tl_FC(i,0)=0.0_r8
            FC(i,N(ng))=0.0_r8
            tl_FC(i,N(ng))=0.0_r8
          END DO
!
!  Solve the tangent linear tridiagonal system.
!
          DO k=1,N(ng)
            DO i=Istr,Iend
              BC(i,k)=Hzk(i,k)-FC(i,k)-FC(i,k-1)
              tl_BC(i,k)=tl_Hzk(i,k)-tl_FC(i,k)-tl_FC(i,k-1)
            END DO
          END DO
          DO k=2,N(ng)-1
            DO i=Istr,Iend
#  ifdef TL_IOMS
            DC(i,k)=tl_v(i,j,k,nnew)-                                   &
     &              ((tl_FC(i,k-1)-FC(i,k-1))*v(i,j,k-1,nnew)+          &
     &               (tl_BC(i,k  )-BC(i,k  ))*v(i,j,k  ,nnew)+          &
     &               (tl_FC(i,k  )-FC(i,k  ))*v(i,j,k+1,nnew))
#  else
              DC(i,k)=tl_v(i,j,k,nnew)-                                 &
     &                (tl_FC(i,k-1)*v(i,j,k-1,nnew)+                    &
     &                 tl_BC(i,k  )*v(i,j,k  ,nnew)+                    &
     &                 tl_FC(i,k  )*v(i,j,k+1,nnew))
#  endif
            END DO
          END DO
          DO i=Istr,Iend
#  ifdef TL_IOMS
            DC(i,1)=tl_v(i,j,1,nnew)-                                   &
     &              ((tl_BC(i,1)-BC(i,1))*v(i,j,1,nnew)+                &
     &               (tl_FC(i,1)-FC(i,1))*v(i,j,2,nnew))
            DC(i,N(ng))=tl_v(i,j,N(ng),nnew)-                           &
     &                  ((tl_FC(i,N(ng)-1)-FC(i,N(ng)-1))*              &
     &                   v(i,j,N(ng)-1,nnew)+                           &
     &                   (tl_BC(i,N(ng)  )-BC(i,N(ng)  ))*              &
     &                   v(i,j,N(ng)  ,nnew))
#  else
            DC(i,1)=tl_v(i,j,1,nnew)-                                   &
     &              (tl_BC(i,1)*v(i,j,1,nnew)+                          &
     &               tl_FC(i,1)*v(i,j,2,nnew))
            DC(i,N(ng))=tl_v(i,j,N(ng),nnew)-                           &
     &                  (tl_FC(i,N(ng)-1)*v(i,j,N(ng)-1,nnew)+          &
     &                   tl_BC(i,N(ng)  )*v(i,j,N(ng)  ,nnew))
#  endif
          END DO
          DO i=Istr,Iend
            cff=1.0_r8/BC(i,1)
            CF(i,1)=cff*FC(i,1)
            DC(i,1)=cff*DC(i,1)
          END DO
          DO k=2,N(ng)-1
            DO i=Istr,Iend
              cff=1.0_r8/(BC(i,k)-FC(i,k-1)*CF(i,k-1))
              CF(i,k)=cff*FC(i,k)
              DC(i,k)=cff*(DC(i,k)-FC(i,k-1)*DC(i,k-1))
            END DO
          END DO
!
!  Compute new solution by back substitution.
!  (DC is a tangent linear variable here).
!
          DO i=Istr,Iend
#  ifdef DIAGNOSTICS_UV
!!          wrk(i,N(ng))=v(i,j,N(ng),nnew)*oHz(i,N(ng))
#  endif
            DC(i,N(ng))=(DC(i,N(ng))-FC(i,N(ng)-1)*DC(i,N(ng)-1))/      &
     &                  (BC(i,N(ng))-FC(i,N(ng)-1)*CF(i,N(ng)-1))
!>          v(i,j,N(ng),nnew)=DC(i,N(ng))
!>
            tl_v(i,j,N(ng),nnew)=DC(i,N(ng))
#  ifdef DIAGNOSTICS_UV
!!          DiaV3wrk(i,j,N(ng),M3vvis)=DiaV3wrk(i,j,N(ng),M3vvis)+      &
!!   &                                 v(i,j,N(ng),nnew)-wrk(i,N(ng))
#  endif
          END DO
          DO k=N(ng)-1,1,-1
            DO i=Istr,Iend
#  ifdef DIAGNOSTICS_UV
!!            wrk(i,k)=v(i,j,k,nnew)*oHz(i,k)
#  endif
              DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
!>            v(i,j,k,nnew)=DC(i,k)
!>
              tl_v(i,j,k,nnew)=DC(i,k)
#  ifdef DIAGNOSTICS_UV
!!            DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+            &
!!   &                               v(i,j,k,nnew)-wrk(i,k)
#  endif
            END DO
          END DO
# endif
!
!  Replace INTERIOR POINTS incorrect vertical mean with more accurate
!  barotropic component, vbar=DV_avg1/(D*om_v). Recall that, D=CF(:,0).
!
          DO i=Istr,Iend
            CF(i,0)=Hzk(i,1)
            tl_CF(i,0)=tl_Hzk(i,1)
            DC(i,0)=v(i,j,1,nnew)*Hzk(i,1)
            tl_DC(i,0)=tl_v(i,j,1,nnew)*Hzk(i,1)+                       &
     &                 v(i,j,1,nnew)*tl_Hzk(i,1)-                       &
# ifdef TL_IOMS
     &                 DC(i,0)
# endif
# ifdef NEARSHORE_MELLOR
            DCs(i,0)=v_stokes(i,j,1)*Hzk(i,1)
            tl_DCs(i,0)=tl_v_stokes(i,j,1)*Hzk(i,1)+                    &
     &                  v_stokes(i,j,1)*tl_Hzk(i,1)-                    &
#  ifdef TL_IOMS
     &                  DCs(i,0)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!          Dwrk(i,M2pgrd)=DiaV3wrk(i,j,1,M3pgrd)*Hzk(i,1)
!!          Dwrk(i,M2bstr)=DiaV3wrk(i,j,1,M3vvis)*Hzk(i,1)
#  ifdef UV_COR
!!          Dwrk(i,M2fcor)=DiaV3wrk(i,j,1,M3fcor)*Hzk(i,1)
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
!!          Dwrk(i,M2xvis)=DiaV3wrk(i,j,1,M3xvis)*Hzk(i,1)
!!          Dwrk(i,M2yvis)=DiaV3wrk(i,j,1,M3yvis)*Hzk(i,1)
!!          Dwrk(i,M2hvis)=DiaV3wrk(i,j,1,M3hvis)*Hzk(i,1)
#  endif
#  ifdef UV_ADV
!!          Dwrk(i,M2xadv)=DiaV3wrk(i,j,1,M3xadv)*Hzk(i,1)
!!          Dwrk(i,M2yadv)=DiaV3wrk(i,j,1,M3yadv)*Hzk(i,1)
!!          Dwrk(i,M2hadv)=DiaV3wrk(i,j,1,M3hadv)*Hzk(i,1)
#  endif
#  ifdef NEARSHORE_MELLOR
!!          Dwrk(i,M2hrad)=DiaV3wrk(i,j,1,M3hrad)*Hzk(i,1)
#  endif
# endif
          END DO
          DO k=2,N(ng)
            DO i=Istr,Iend
              CF(i,0)=CF(i,0)+Hzk(i,k)
              tl_CF(i,0)=tl_CF(i,0)+tl_Hzk(i,k)
              DC(i,0)=DC(i,0)+v(i,j,k,nnew)*Hzk(i,k)
              tl_DC(i,0)=tl_DC(i,0)+                                    &
     &                   tl_v(i,j,k,nnew)*Hzk(i,k)+                     &
     &                   v(i,j,k,nnew)*tl_Hzk(i,k)-                     &
# ifdef TL_IOMS
     &                   v(i,j,k,nnew)*Hzk(i,k)
# endif
# ifdef NEARSHORE_MELLOR
              DCs(i,0)=DCs(i,0)+v_stokes(i,j,k)*Hzk(i,k)
              tl_DCs(i,0)=tl_DCs(i,0)+                                  &
     &                    tl_v_stokes(i,j,k)*Hzk(i,k)+                  &
     &                    v_stokes(i,j,k)*tl_Hzk(i,k)-                  &
#  ifdef TL_IOMS
     &                    v_stokes(i,j,k)*Hzk(i,k)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!            Dwrk(i,M2pgrd)=Dwrk(i,M2pgrd)+                            &
!!   &                       DiaV3wrk(i,j,k,M3pgrd)*Hzk(i,k)
!!            Dwrk(i,M2bstr)=Dwrk(i,M2bstr)+                            &
!!   &                       DiaV3wrk(i,j,k,M3vvis)*Hzk(i,k)
#  ifdef UV_COR
!!            Dwrk(i,M2fcor)=Dwrk(i,M2fcor)+                            &
!!   &                       DiaV3wrk(i,j,k,M3fcor)*Hzk(i,k)
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
!!            Dwrk(i,M2xvis)=Dwrk(i,M2xvis)+                            &
!!   &                       DiaV3wrk(i,j,k,M3xvis)*Hzk(i,k)
!!            Dwrk(i,M2yvis)=Dwrk(i,M2yvis)+                            &
!!   &                       DiaV3wrk(i,j,k,M3yvis)*Hzk(i,k)
!!            Dwrk(i,M2hvis)=Dwrk(i,M2hvis)+                            &
!!   &                       DiaV3wrk(i,j,k,M3hvis)*Hzk(i,k)
#  endif
#  ifdef UV_ADV
!!            Dwrk(i,M2xadv)=Dwrk(i,M2xadv)+                            &
!!   &                       DiaV3wrk(i,j,k,M3xadv)*Hzk(i,k)
!!            Dwrk(i,M2yadv)=Dwrk(i,M2yadv)+                            &
!!   &                       DiaV3wrk(i,j,k,M3yadv)*Hzk(i,k)
!!            Dwrk(i,M2hadv)=Dwrk(i,M2hadv)+                            &
!!   &                       DiaV3wrk(i,j,k,M3hadv)*Hzk(i,k)
#  endif
#  ifdef NEARSHORE_MELLOR
!!            Dwrk(i,M2hrad)=Dwrk(i,M2hrad)+                            &
!!   &                       DiaV3wrk(i,j,k,M3hrad)*Hzk(i,k)
#  endif
# endif
            END DO
          END DO
          DO i=Istr,Iend
            DC1(i,0)=DC(i,0)                              ! intermediate
            cff1=1.0_r8/(CF(i,0)*om_v(i,j))
            tl_cff1=-cff1*cff1*tl_CF(i,0)*om_v(i,j)+                    &
# ifdef TL_IOMS
     &              2.0_r8*cff1
# endif
            DC(i,0)=(DC(i,0)*om_v(i,j)-DV_avg1(i,j))*cff1    ! recursive
            tl_DC(i,0)=(tl_DC(i,0)*om_v(i,j)-tl_DV_avg1(i,j))*cff1+     &
     &                 (DC1(i,0)*om_v(i,j)-DV_avg1(i,j))*tl_cff1-       &
# ifdef TL_IOMS
     &                 DC(i,0)
# endif
# ifdef NEARSHORE_MELLOR
            DCs1(i)=DCs(i,0)                              ! intermediate
            cff2=1.0_r8/CF(i,0)
            tl_cff2=-cff2*cff2*tl_CF(i,0)+                              &
#  ifdef TL_IOMS
     &              2.0_r8*cff2
#  endif
            DCs(i,0)=DCs(i,0)*cff2-vbar_stokes(i,j)          ! recursive
            tl_DCs(i,0)=tl_DCs(i,0)*cff2+                               &
     &                  DCs1(i,0)*tl_cff2-tl_vbar_stokes(i,j)-          &
#  ifdef TL_IOMS
     &                  DCs(i,0)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!          DO idiag=1,M2pgrd
!!            Dwrk(i,idiag)=(Dwrk(i,idiag)*om_v(i,j)-                   &
!!   &                       DiaV2wrk(i,j,idiag))*cff1
!!          END DO
!!          Dwrk(i,M2bstr)=(Dwrk(i,M2bstr)*om_v(i,j)-                   &
!!   &                      DiaV2wrk(i,j,M2bstr)-DiaV2wrk(i,j,M2sstr))* &
!!   &                     cff1
# endif
          END DO
!
!  Couple and update new solution.
!
          DO k=1,N(ng)
            DO i=Istr,Iend
!>            v(i,j,k,nnew)=v(i,j,k,nnew)-DC(i,0)
!>
              tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)-tl_DC(i,0)
# ifdef MASKING
!>            v(i,j,k,nnew)=v(i,j,k,nnew)*vmask(i,j)
!>
              tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*vmask(i,j)
# endif
# ifdef WET_DRY_NOT_YET
!>            v(i,j,k,nnew)=v(i,j,k,nnew)*vmask_wet(i,j)
!>
              tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*vmask_wet(i,j)
# endif
# ifdef NEARSHORE_MELLOR
!>            v_stokes(i,j,k)=v_stokes(i,j,k)-DCs(i,0)
!>
              tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)-tl_DCs(i,0)
#  ifdef MASKING
!>            v_stokes(i,j,k)=v_stokes(i,j,k)*vmask(i,j)
!>
              tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*vmask(i,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>            v_stokes(i,j,k)=v_stokes(i,j,k)*vmask_wet(i,j)
!>
              tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*vmask_wet(i,j)
#  endif
# endif
# ifdef DIAGNOSTICS_UV
!!            DiaV3wrk(i,j,k,M3pgrd)=DiaV3wrk(i,j,k,M3pgrd)-            &
!!   &                               Dwrk(i,M2pgrd)
!!            DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)-            &
!!   &                               Dwrk(i,M2bstr)
#  ifdef UV_COR
!!            DiaV3wrk(i,j,k,M3fcor)=DiaV3wrk(i,j,k,M3fcor)-            &
!!   &                               Dwrk(i,M2fcor)
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
!!            DiaV3wrk(i,j,k,M3xvis)=DiaV3wrk(i,j,k,M3xvis)-            &
!!   &                               Dwrk(i,M2xvis)
!!            DiaV3wrk(i,j,k,M3yvis)=DiaV3wrk(i,j,k,M3yvis)-            &
!!   &                               Dwrk(i,M2yvis)
!!            DiaV3wrk(i,j,k,M3hvis)=DiaV3wrk(i,j,k,M3hvis)-            &
!!   &                               Dwrk(i,M2hvis)
#  endif
#  ifdef UV_ADV
!!            DiaV3wrk(i,j,k,M3xadv)=DiaV3wrk(i,j,k,M3xadv)-            &
!!   &                               Dwrk(i,M2xadv)
!!            DiaV3wrk(i,j,k,M3yadv)=DiaV3wrk(i,j,k,M3yadv)-            &
!!   &                               Dwrk(i,M2yadv)
!!            DiaV3wrk(i,j,k,M3hadv)=DiaV3wrk(i,j,k,M3hadv)-            &
!!   &                               Dwrk(i,M2hadv)
#  endif
#  ifdef NEARSHORE_MELLOR
!!            DiaV3wrk(i,j,k,M3hrad)=DiaV3wrk(i,j,k,M3hrad)-            &
!!   &                               Dwrk(i,M2hrad)
#  endif
# endif
            END DO
          END DO

# if defined DIAGNOSTICS_UV && defined MASKING
!!        DO k=1,N(ng)
!!          DO i=Istr,Iend
!!            DO idiag=1,NDM3d
!!              DiaV3wrk(i,j,k,idiag)=DiaV3wrk(i,j,k,idiag)*vmask(i,j)
!!            END DO
!!          END DO
!!        END DO
# endif
        END IF
      END DO
!
!-----------------------------------------------------------------------
!  Set lateral boundary conditions.
!-----------------------------------------------------------------------
!
!>    CALL u3dbc_tile (ng, tile,                                        &
!>   &                 LBi, UBi, LBj, UBj, N(ng),                       &
!>   &                 IminS, ImaxS, JminS, JmaxS,                      &
!>   &                 nstp, nnew,                                      &
!>   &                 u)
!>
      CALL rp_u3dbc_tile (ng, tile,                                     &
     &                    LBi, UBi, LBj, UBj, N(ng),                    &
     &                    IminS, ImaxS, JminS, JmaxS,                   &
     &                    nstp, nnew,                                   &
     &                    tl_u)
!>    CALL v3dbc_tile (ng, tile,                                        &
!>   &                 LBi, UBi, LBj, UBj, N(ng),                       &
!>   &                 IminS, ImaxS, JminS, JmaxS,                      &
!>   &                 nstp, nnew,                                      &
!>   &                 v)
!>
      CALL rp_v3dbc_tile (ng, tile,                                     &
     &                    LBi, UBi, LBj, UBj, N(ng),                    &
     &                    IminS, ImaxS, JminS, JmaxS,                   &
     &                    nstp, nnew,                                   &
     &                    tl_v)
!
!-----------------------------------------------------------------------
!  Apply momentum transport point sources (like river runoff), if any.
!-----------------------------------------------------------------------
!
      IF (LuvSrc(ng)) THEN
        DO is=1,Nsrc(ng)
          i=SOURCES(ng)%Isrc(is)
          j=SOURCES(ng)%Jsrc(is)
          IF (((IstrR.le.i).and.(i.le.IendR)).and.                      &
     &        ((JstrR.le.j).and.(j.le.JendR))) THEN
            IF (INT(SOURCES(ng)%Dsrc(is)).eq.0) THEN
              DO k=1,N(ng)
                cff1=1.0_r8/(on_u(i,j)*                                 &
     &                       0.5_r8*(z_w(i-1,j,k)-z_w(i-1,j,k-1)+       &
     &                               z_w(i  ,j,k)-z_w(i  ,j,k-1)))
                tl_cff1=-cff1*cff1*on_u(i,j)*                           &
     &                  0.5_r8*(tl_z_w(i-1,j,k)-tl_z_w(i-1,j,k-1)+      &
     &                          tl_z_w(i  ,j,k)-tl_z_w(i  ,j,k-1))+     &
# ifdef TL_IOMS
     &                  2.0_r8*cff1
# endif
!>              u(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*cff1
!>
                tl_u(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*tl_cff1
              END DO
            ELSE
              DO k=1,N(ng)
                cff1=1.0_r8/(om_v(i,j)*                                 &
     &                       0.5_r8*(z_w(i,j-1,k)-z_w(i,j-1,k-1)+       &
     &                               z_w(i,j  ,k)-z_w(i,j  ,k-1)))
                tl_cff1=-cff1*cff1*om_v(i,j)*                           &
     &                  0.5_r8*(tl_z_w(i,j-1,k)-tl_z_w(i,j-1,k-1)+      &
     &                          tl_z_w(i,j  ,k)-tl_z_w(i,j  ,k-1))+     &
# ifdef TL_IOMS
     &                  2.0_r8*cff1
# endif
!>              v(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*cff1
!>
                tl_v(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*tl_cff1
              END DO
            END IF
          END IF
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Couple 2D and 3D momentum equations.
!-----------------------------------------------------------------------
!
!  Couple velocity component in the XI-direction.
!
      DO j=JstrT,JendT
        DO i=IstrP,IendT
          DC(i,0)=0.0_r8
          tl_DC(i,0)=0.0_r8
          CF(i,0)=0.0_r8
          tl_CF(i,0)=0.0_r8
# ifdef NEARSHORE_MELLOR
          CFs(i,0)=0.0_r8
          tl_CFs(i,0)=0.0_r8
# endif
          FC(i,0)=0.0_r8
          tl_FC(i,0)=0.0_r8
        END DO
!
!  Compute thicknesses of U-boxes DC(i,1:N), total depth of the water
!  column DC(i,0), and incorrect vertical mean CF(i,0).  Notice that
!  barotropic component is replaced with its fast-time averaged
!  values.
!
        DO k=1,N(ng)
          DO i=IstrP,IendT
            cff=0.5_r8*on_u(i,j)
            DC(i,k)=cff*(Hz(i,j,k)+Hz(i-1,j,k))
            tl_DC(i,k)=cff*(tl_Hz(i,j,k)+tl_Hz(i-1,j,k))
            DC(i,0)=DC(i,0)+DC(i,k)
            tl_DC(i,0)=tl_DC(i,0)+tl_DC(i,k)
            CF(i,0)=CF(i,0)+                                            &
     &              DC(i,k)*u(i,j,k,nnew)
            tl_CF(i,0)=tl_CF(i,0)+                                      &
     &                 tl_DC(i,k)*u(i,j,k,nnew)+                        &
     &                 DC(i,k)*tl_u(i,j,k,nnew)-                        &
# ifdef TL_IOMS
     &                 DC(i,k)*u(i,j,k,nnew)
# endif
# ifdef NEARSHORE_MELLOR
            CFs(i,0)=CFs(i,0)+                                          &
     &               DC(i,k)*u_stokes(i,j,k)
            tl_CFs(i,0)=tl_CFs(i,0)+                                    &
     &                  tl_DC(i,k)*u_stokes(i,j,k)+                     &
     &                  DC(i,k)*tl_u_stokes(i,j,k)-                     &
#  ifdef TL_IOMS
     &                  DC(i,k)*u_stokes(i,j,k)
#  endif
# endif
          END DO
        END DO
        DO i=IstrP,IendT
          DC1(i,0)=DC(i,0)                              ! intermediate
# ifdef NEARSHORE_MELLOR
          cff2=DC(i,0)*ubar_stokes(i,j)
          tl_cff2=tl_DC(i,0)*ubar_stokes(i,j)+                          &
     &            DC(i,0)*tl_ubar_stokes(i,j)-                          &
#  ifdef TL_IOMS
     &            cff2
#  endif
# endif
          DC(i,0)=1.0_r8/DC(i,0)                        ! recursive
          tl_DC(i,0)=-tl_DC(i,0)/(DC1(i,0)*DC1(i,0))+                   &
# ifdef TL_IOMS
     &                2.0_r8*DC(i,0)
# endif
          CF1(i)=CF(i,0)                                ! intermediate
          CF(i,0)=DC(i,0)*(CF(i,0)-DU_avg1(i,j))        ! recursive
          tl_CF(i,0)=tl_DC(i,0)*(CF1(i)-DU_avg1(i,j))+                  &
     &               DC(i,0)*(tl_CF(i,0)-tl_DU_avg1(i,j))-              &
# ifdef TL_IOMS
     &               CF(i,0)
# endif
# ifdef NEARSHORE_MELLOR
          CFs1(i)=CFs(i,0)                              ! intermediate
          CFs(i,0)=DC(i,0)*(CFs(i,0)-cff2)              ! recursive
          tl_CFs(i,0)=tl_DC(i,0)*(CFs1(i)-cff2)+                        &
     &                DC(i,0)*(tl_CFs(i,0)-tl_cff2)-                    &
#  ifdef TL_IOMS
     &                CFs(i,0)
#  endif
# endif
!>        ubar(i,j,1)=DC(i,0)*DU_avg1(i,j)
!>
          tl_ubar(i,j,1)=tl_DC(i,0)*DU_avg1(i,j)+                       &
     &                   DC(i,0)*tl_DU_avg1(i,j)-                       &
# ifdef TL_IOMS
     &                   DC(i,0)*DU_avg1(i,j)
# endif
!>        ubar(i,j,2)=ubar(i,j,1)
!>
          tl_ubar(i,j,2)=tl_ubar(i,j,1)
# ifdef DIAGNOSTICS_UV
!!        DiaU2wrk(i,j,M2rate)=ubar(i,j,1)-DiaU2int(i,j,M2rate)*DC(i,0)
!!        DiaU2int(i,j,M2rate)=ubar(i,j,1)*DC1(i,0)
# endif
        END DO
# ifdef DIAGNOSTICS_UV
!!
!! Convert the units of the fast-time integrated change in mass flux
!! diagnostic values to change in velocity (m s-1).
!!
!!      DO idiag=1,NDM2d-1
!!        DO i=IstrP,IendT
!!          DiaU2wrk(i,j,idiag)=DC(i,0)*DiaU2wrk(i,j,idiag)
#  ifdef MASKING
!!          DiaU2wrk(i,j,idiag)=DiaU2wrk(i,j,idiag)*umask(i,j)
#  endif
!!        END DO
!!      END DO
# endif
!
!  Replace only BOUNDARY POINTS incorrect vertical mean with more
!  accurate barotropic component, ubar=DU_avg1/(D*on_u). Recall that,
!  D=CF(:,0).
!
!  NOTE:  Only the BOUNDARY POINTS need to be replaced. Avoid redundant
!         update in the interior again for computational purposes which
!         will not affect the nonlinear code.  However, the adjoint
!         code is wrong because the interior solution is corrected
!         twice. The replacement is avoided if the boundary edge is
!         periodic. The J-loop is pipelined, so we need to do a special
!         test for the southern and northern domain edges.
!
        IF (.not.(CompositeGrid(iwest,ng).or.EWperiodic(ng))) THEN
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            DO k=1,N(ng)
!>            u(Istr,j,k,nnew)=u(Istr,j,k,nnew)-CF(Istr,0)
!>
              tl_u(Istr,j,k,nnew)=tl_u(Istr,j,k,nnew)-                  &
     &                            tl_CF(Istr,0)
# ifdef MASKING
!>            u(Istr,j,k,nnew)=u(Istr,j,k,nnew)*                        &
!>   &                         umask(Istr,j)
!>
              tl_u(Istr,j,k,nnew)=tl_u(Istr,j,k,nnew)*                  &
     &                            umask(Istr,j)
# endif
# ifdef WET_DRY_NOT_YET
!>            u(Istr,j,k,nnew)=u(Istr,j,k,nnew)*                        &
!>   &                         umask_wet(Istr,j)
!>
              tl_u(Istr,j,k,nnew)=tl_u(Istr,j,k,nnew)*                  &
     &                            umask_wet(Istr,j)
# endif
# ifdef NEARSHORE_MELLOR
!>            u_stokes(Istr,j,k)=u_stokes(Istr,j,k)-CFs(Istr,0)
!>
              tl_u_stokes(Istr,j,k)=tl_u_stokes(Istr,j,k)-              &
     &                              tl_CFs(Istr,0)
#  ifdef MASKING
!>            u_stokes(Istr,j,k)=u_stokes(Istr,j,k)*                    &
!>   &                           umask(Istr,j)
!>
              tl_u_stokes(Istr,j,k)=tl_u_stokes(Istr,j,k)*              &
     &                              umask(Istr,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>            u_stokes(Istr,j,k)=u_stokes(Istr,j,k)*                    &
!>   &                           umask_wet(Istr,j)
!>
              tl_u_stokes(Istr,j,k)=tl_u_stokes(Istr,j,k)*              &
     &                              umask_wet(Istr,j)
#  endif
# endif
            END DO
          END IF
        END IF
!
        IF (.not.(CompositeGrid(ieast,ng).or.EWperiodic(ng))) THEN
          IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
            DO k=1,N(ng)
!>            u(Iend+1,j,k,nnew)=u(Iend+1,j,k,nnew)-CF(Iend+1,0)
!>
              tl_u(Iend+1,j,k,nnew)=tl_u(Iend+1,j,k,nnew)-              &
     &                              tl_CF(Iend+1,0)
# ifdef MASKING
!>            u(Iend+1,j,k,nnew)=u(Iend+1,j,k,nnew)*                    &
!>   &                           umask(Iend+1,j)
!>
              tl_u(Iend+1,j,k,nnew)=tl_u(Iend+1,j,k,nnew)*              &
     &                              umask(Iend+1,j)
# endif
# ifdef WET_DRY_NOT_YET
!>            u(Iend+1,j,k,nnew)=u(Iend+1,j,k,nnew)*                    &
!>   &                           umask_wet(Iend+1,j)
!>
              tl_u(Iend+1,j,k,nnew)=tl_u(Iend+1,j,k,nnew)*              &
     &                              umask_wet(Iend+1,j)
# endif
# ifdef NEARSHORE_MELLOR
!>            u_stokes(Iend+1,j,k)=u_stokes(Iend+1,j,k)-CFs(Iend+1,0)
!>
              tl_u_stokes(Iend+1,j,k)=tl_u_stokes(Iend+1,j,k)-          &
     &                                tl_CFs(Iend+1,0)
#  ifdef MASKING
!>            u_stokes(Iend+1,j,k)=u_stokes(Iend+1,j,k)*                &
!>   &                             umask(Iend+1,j)
!>
              tl_u_stokes(Iend+1,j,k)=tl_u_stokes(Iend+1,j,k)*          &
     &                                umask(Iend+1,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>            u_stokes(Iend+1,j,k)=u_stokes(Iend+1,j,k)*                &
!>   &                             umask_wet(Iend+1,j)
!>
              tl_u_stokes(Iend+1,j,k)=tl_u_stokes(Iend+1,j,k)*          &
     &                                umask_wet(Iend+1,j)
#  endif
# endif
            END DO
          END IF
        END IF
!
        IF (.not.(CompositeGrid(isouth,ng).or.NSperiodic(ng))) THEN
          IF (j.eq.0) THEN                           ! southern boundary
            DO k=1,N(ng)                             ! J-loop pipelined
              DO i=IstrU,Iend
!>              u(i,j,k,nnew)=u(i,j,k,nnew)-CF(i,0)
!>
                tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-                      &
     &                           tl_CF(i,0)
# ifdef MASKING
!>              u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
!>
                tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*                      &
     &                           umask(i,j)
# endif
# ifdef WET_DRY_NOT_YET
!>              u(i,j,k,nnew)=u(i,j,k,nnew)*                            &
!>   &                        umask_wet(i,j)
!>
                tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*                      &
     &                           umask_wet(i,j)
# endif
# ifdef NEARSHORE_MELLOR
!>              u_stokes(i,j,k)=u_stokes(i,j,k)-CFs(i,0)
!>
                tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)-                  &
     &                             tl_CFs(i,0)
#  ifdef MASKING
!>              u_stokes(i,j,k)=u_stokes(i,j,k)*                        &
!>   &                          umask(i,j)
!>
                tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*                  &
     &                             umask(i,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>              u_stokes(i,j,k)=u_stokes(i,j,k)*                        &
!>   &                          umask_wet(i,j)
!>
                tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*                  &
     &                             umask_wet(i,j)
#  endif
# endif
              END DO
            END DO
          END IF
        END IF
!
        IF (.not.(CompositeGrid(inorth,ng).or.NSperiodic(ng))) THEN
          IF (j.eq.Mm(ng)+1) THEN                    ! northern boundary
            DO k=1,N(ng)                             ! J-loop pipelined
              DO i=IstrU,Iend
!>              u(i,j,k,nnew)=u(i,j,k,nnew)-CF(i,0)
!>
                tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-                      &
     &                           tl_CF(i,0)
# ifdef MASKING
!>              u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
!>
                tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*                      &
     &                           umask(i,j)
# endif
# ifdef WET_DRY_NOT_YET
!>              u(i,j,k,nnew)=u(i,j,k,nnew)*umask_wet(i,j)
!>
                tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*                      &
     &                           umask_wet(i,j)
# endif
# ifdef NEARSHORE_MELLOR
!>              u_stokes(i,j,k)=u_stokes(i,j,k)-CFs(i,0)
!>
                tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)-                  &
     &                             tl_CFs(i,0)
#  ifdef MASKING
!>              u_stokes(i,j,k)=u_stokes(i,j,k)*                        &
!>   &                          umask(i,j)
!>
                tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*                  &
     &                             umask(i,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>              u_stokes(i,j,k)=u_stokes(i,j,k)*                        &
!>   &                          umask_wet(i,j)
!>
                tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*                  &
     &                             umask_wet(i,j)
#  endif
# endif
              END DO
            END DO
          END IF
        END IF
!
!  Compute correct mass flux, Hz*u/n.
!
        DO k=N(ng),1,-1
          DO i=IstrP,IendT
            Huon(i,j,k)=0.5_r8*(Huon(i,j,k)+u(i,j,k,nnew)*DC(i,k))
            tl_Huon(i,j,k)=0.5_r8*(tl_Huon(i,j,k)+                      &
     &                             tl_u(i,j,k,nnew)*DC(i,k)+            &
     &                             u(i,j,k,nnew)*tl_DC(i,k))-           &
# ifdef TL_IOMS
     &                     0.5_r8*u(i,j,k,nnew)*DC(i,k)
# endif
# ifdef NEARSHORE_MELLOR
            Huon(i,j,k)=Huon(i,j,k)+0.5_r8*u_stokes(i,j,k)*DC(i,k)
            tl_Huon(i,j,k)=tl_Huon(i,j,k)+                              &
     &                     0.5_r8*(tl_u_stokes(i,j,k)*DC(i,k)+          &
     &                             u_stokes(i,j,k)*tl_DC(i,k))-         &
#  ifdef TL_IOMS
     &                     0.5_r8*u_stokes(i,j,k)*DC(i,k)
#  endif
# endif
            FC(i,0)=FC(i,0)+Huon(i,j,k)
            tl_FC(i,0)=tl_FC(i,0)+tl_Huon(i,j,k)
# ifdef DIAGNOSTICS_UV
!!          DiaU3wrk(i,j,k,M3rate)=u(i,j,k,nnew)-DiaU3wrk(i,j,k,M3rate)
# endif
          END DO
        END DO
        DO i=IstrP,IendT
          FC1(i)=FC(i,0)                                ! intermediate
          FC(i,0)=DC(i,0)*(FC(i,0)-DU_avg2(i,j))        ! recursive
          tl_FC(i,0)=tl_DC(i,0)*(FC1(i)-DU_avg2(i,j))+                  &
     &               DC(i,0)*(tl_FC(i,0)-tl_DU_avg2(i,j))-              &
# ifdef TL_IOMS
     &               FC(i,0)
# endif
        END DO
        DO k=1,N(ng)
          DO i=IstrP,IendT
            Huon(i,j,k)=Huon(i,j,k)-DC(i,k)*FC(i,0)
            tl_Huon(i,j,k)=tl_Huon(i,j,k)-                              &
     &                     tl_DC(i,k)*FC(i,0)-                          &
     &                     DC(i,k)*tl_FC(i,0)+                          &
# ifdef TL_IOMS
     &                     DC(i,k)*FC(i,0)
# endif
          END DO
        END DO
!
!  Couple velocity component in the ETA-direction.
!
        IF (j.ge.Jstr) THEN
          DO i=IstrT,IendT
            DC(i,0)=0.0_r8
            tl_DC(i,0)=0.0_r8
            CF(i,0)=0.0_r8
            tl_CF(i,0)=0.0_r8
# ifdef NEARSHORE_MELLOR
            CFs(i,0)=0.0_r8
            tl_CFs(i,0)=0.0_r8
# endif
            FC(i,0)=0.0_r8
            tl_FC(i,0)=0.0_r8
          END DO
!
!  Compute thicknesses of V-boxes DC(i,1:N), total depth of the water
!  column DC(i,0), and incorrect vertical mean CF(i,0).  Notice that
!  barotropic component is replaced with its fast-time averaged
!  values.
!
          DO k=1,N(ng)
            DO i=IstrT,IendT
              cff=0.5_r8*om_v(i,j)
              DC(i,k)=cff*(Hz(i,j,k)+Hz(i,j-1,k))
              tl_DC(i,k)=cff*(tl_Hz(i,j,k)+tl_Hz(i,j-1,k))
              DC(i,0)=DC(i,0)+DC(i,k)
              tl_DC(i,0)=tl_DC(i,0)+tl_DC(i,k)
              CF(i,0)=CF(i,0)+                                          &
     &                DC(i,k)*v(i,j,k,nnew)
              tl_CF(i,0)=tl_CF(i,0)+                                    &
     &                   tl_DC(i,k)*v(i,j,k,nnew)+                      &
     &                   DC(i,k)*tl_v(i,j,k,nnew)-                      &
# ifdef TL_IOMS
     &                   DC(i,k)*v(i,j,k,nnew)
# endif
# ifdef NEARSHORE_MELLOR
              CFs(i,0)=CFs(i,0)+                                        &
     &                 DC(i,k)*v_stokes(i,j,k)
              tl_CFs(i,0)=tl_CFs(i,0)+                                  &
     &                    tl_DC(i,k)*v_stokes(i,j,k)+                   &
     &                    DC(i,k)*tl_v_stokes(i,j,k)-                   &
#  ifdef TL_IOMS
     &                    DC(i,k)*v_stokes(i,j,k)
#  endif
# endif
            END DO
          END DO
          DO i=IstrT,IendT
            DC1(i,0)=DC(i,0)                            ! intermediate
# ifdef NEARSHORE_MELLOR
            cff2=DC(i,0)*vbar_stokes(i,j)
            tl_cff2=tl_DC(i,0)*vbar_stokes(i,j)+                        &
     &              DC(i,0)*tl_vbar_stokes(i,j)-                        &
#  ifdef TL_IOMS
     &              cff2
#  endif
# endif
            DC(i,0)=1.0_r8/DC(i,0)                      ! recursive
            tl_DC(i,0)=-tl_DC(i,0)/(DC1(i,0)*DC1(i,0))+                 &
# ifdef TL_IOMS
     &                 2.0_r8*DC(i,0)
# endif
            CF1(i)=CF(i,0)                              ! intermediate
            CF(i,0)=DC(i,0)*(CF(i,0)-DV_avg1(i,j))      ! recursive
            tl_CF(i,0)=tl_DC(i,0)*(CF1(i)-DV_avg1(i,j))+                &
     &                 DC(i,0)*(tl_CF(i,0)-tl_DV_avg1(i,j))-            &
# ifdef TL_IOMS
     &                 CF(i,0)
# endif
# ifdef NEARSHORE_MELLOR
            CFs1(i)=CFs(i,0)
            CFs(i,0)=DC(i,0)*(CFs(i,0)-cff2)            ! recursive
            tl_CFs(i,0)=tl_DC(i,0)*(CFs1(i)-cff2)+                      &
     &                  DC(i,0)*(tl_CFs(i,0)-tl_cff2)-                  &
#  ifdef TL_IOMS
     &                  CF(i,0)
#  endif
# endif
!>          vbar(i,j,1)=DC(i,0)*DV_avg1(i,j)
!>
            tl_vbar(i,j,1)=tl_DC(i,0)*DV_avg1(i,j)+                     &
     &                     DC(i,0)*tl_DV_avg1(i,j)-                     &
# ifdef TL_IOMS
     &                     DC(i,0)*DV_avg1(i,j)
# endif
!>          vbar(i,j,2)=vbar(i,j,1)
!>
            tl_vbar(i,j,2)=tl_vbar(i,j,1)
# ifdef DIAGNOSTICS_UV
!!          DiaV2wrk(i,j,M2rate)=vbar(i,j,1)-                           &
!!   &                           DiaV2int(i,j,M2rate)*DC(i,0)
!!          DiaV2int(i,j,M2rate)=vbar(i,j,1)*DC1(i,0)
!!          DiaV2wrk(i,j,M2rate)=vbar(i,j,1)-DiaV2int(i,j,M2rate)
!!          DiaV2int(i,j,M2rate)=vbar(i,j,1)
# endif
          END DO
# ifdef DIAGNOSTICS_UV
!!
!! Convert the units of the fast-time integrated change in mass flux
!! diagnostic values to change in velocity (m s-1).
!!
!!        DO idiag=1,NDM2d-1
!!          DO i=IstrT,IendT
!!            DiaV2wrk(i,j,idiag)=DC(i,0)*DiaV2wrk(i,j,idiag)
#  ifdef MASKING
!!            DiaV2wrk(i,j,idiag)=DiaV2wrk(i,j,idiag)*vmask(i,j)
#  endif
!!          END DO
!!        END DO
# endif
!
!  Replace only BOUNDARY POINTS incorrect vertical mean with more
!  accurate barotropic component, vbar=DV_avg1/(D*om_v).  Recall that,
!  D=CF(:,0).
!
!  NOTE:  Only the BOUNDARY POINTS need to be replaced. Avoid redundant
!         update in the interior again for computational purposes which
!         will not affect the nonlinear code.  However, the adjoint
!         code is wrong because the interior solution is corrected
!         twice. The replacement is avoided if the boundary edge is
!         periodic. The J-loop is pipelined, so we need to do a special
!         test for the southern and northern domain edges.
!
          IF (.not.(CompositeGrid(iwest,ng).or.EWperiodic(ng))) THEN
            IF (DOMAIN(ng)%Western_Edge(tile)) THEN
              DO k=1,N(ng)
!>              v(Istr-1,j,k,nnew)=v(Istr-1,j,k,nnew)-CF(Istr-1,0)
!>
                tl_v(Istr-1,j,k,nnew)=tl_v(Istr-1,j,k,nnew)-            &
     &                                tl_CF(Istr-1,0)
# ifdef MASKING
!>              v(Istr-1,j,k,nnew)=v(Istr-1,j,k,nnew)*                  &
!>   &                             vmask(Istr-1,j)
!>
                tl_v(Istr-1,j,k,nnew)=tl_v(Istr-1,j,k,nnew)*            &
     &                                vmask(Istr-1,j)
# endif
# ifdef WET_DRY_NOT_YET
!>              v(Istr-1,j,k,nnew)=v(Istr-1,j,k,nnew)*                  &
!>   &                             vmask_wet(Istr-1,j)
!>
                tl_v(Istr-1,j,k,nnew)=tl_v(Istr-1,j,k,nnew)*            &
     &                                vmask_wet(Istr-1,j)
# endif
# ifdef NEARSHORE_MELLOR
!>              v_stokes(Istr-1,j,k)=v_stokes(Istr-1,j,k)-              &
!>   &                               CFs(Istr-1,0)
!>
                tl_v_stokes(Istr-1,j,k)=tl_v_stokes(Istr-1,j,k)-        &
     &                                  tl_CFs(Istr-1,0)
#  ifdef MASKING
!>              v_stokes(Istr-1,j,k)=v_stokes(Istr-1,j,k)*              &
!>   &                               vmask(Istr-1,j)
!>
                tl_v_stokes(Istr-1,j,k)=tl_v_stokes(Istr-1,j,k)*        &
     &                                  vmask(Istr-1,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>              v_stokes(Istr-1,j,k)=v_stokes(Istr-1,j,k)*              &
!>   &                               vmask_wet(Istr-1,j)
!>
                tl_v_stokes(Istr-1,j,k)=tl_v_stokes(Istr-1,j,k)*        &
     &                                  vmask_wet(Istr-1,j)
#  endif
# endif
              END DO
            END IF
          END IF
!
          IF (.not.(CompositeGrid(ieast,ng).or.EWperiodic(ng))) THEN
            IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
              DO k=1,N(ng)
!>              v(Iend+1,j,k,nnew)=v(Iend+1,j,k,nnew)-CF(Iend+1,0)
!>
                tl_v(Iend+1,j,k,nnew)=tl_v(Iend+1,j,k,nnew)-            &
     &                                tl_CF(Iend+1,0)
# ifdef MASKING
!>              v(Iend+1,j,k,nnew)=v(Iend+1,j,k,nnew)*                  &
!>   &                             vmask(Iend+1,j)
!>
                tl_v(Iend+1,j,k,nnew)=tl_v(Iend+1,j,k,nnew)*            &
     &                                vmask(Iend+1,j)
# endif
# ifdef WET_DRY_NOT_YET
!>              v(Iend+1,j,k,nnew)=v(Iend+1,j,k,nnew)*                  &
!>   &                             vmask_wet(Iend+1,j)
!>
                tl_v(Iend+1,j,k,nnew)=tl_v(Iend+1,j,k,nnew)*            &
     &                                vmask_wet(Iend+1,j)
# endif
# ifdef NEARSHORE_MELLOR
!>              v_stokes(Iend+1,j,k)=v_stokes(Iend+1,j,k)-              &
!>   &                               CFs(Iend+1,0)
!>
                tl_v_stokes(Iend+1,j,k)=tl_v_stokes(Iend+1,j,k)-        &
     &                                  tl_CFs(Iend+1,0)
#  ifdef MASKING
!>              v_stokes(Iend+1,j,k)=v_stokes(Iend+1,j,k)*              &
!>   &                               vmask(Iend+1,j)
!>
                tl_v_stokes(Iend+1,j,k)=tl_v_stokes(Iend+1,j,k)*        &
     &                                  vmask(Iend+1,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>              v_stokes(Iend+1,j,k)=v_stokes(Iend+1,j,k)*              &
!>   &                               vmask_wet(Iend+1,j)
!>
                tl_v_stokes(Iend+1,j,k)=tl_v_stokes(Iend+1,j,k)*        &
     &                                  vmask_wet(Iend+1,j)
#  endif
# endif
              END DO
            END IF
          END IF
!
          IF (.not.(CompositeGrid(isouth,ng).or.NSperiodic(ng))) THEN
            IF (j.eq.1) THEN                         ! southern boundary
              DO k=1,N(ng)                           ! J-loop pipelined
                DO i=Istr,Iend
!>                v(i,j,k,nnew)=v(i,j,k,nnew)-CF(i,0)
!>
                  tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)-                    &
     &                             tl_CF(i,0)
# ifdef MASKING
!>                v(i,j,k,nnew)=v(i,j,k,nnew)*
!>   &                          vmask(i,j)
!>
                  tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*                    &
     &                             vmask(i,j)
# endif
# ifdef WET_DRY_NOT_YET
!>                v(i,j,k,nnew)=v(i,j,k,nnew)*                          &
!>   &                          vmask_wet(i,j)
!>
                  tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*
     &                             vmask_wet(i,j)
# endif
# ifdef NEARSHORE_MELLOR
!>                v_stokes(i,j,k)=v_stokes(i,j,k)-CFs(i,0)
!>
                  tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)-                &
     &                               tl_CFs(i,0)
#  ifdef MASKING
!>                v_stokes(i,j,k)=v_stokes(i,j,k)*vmask(i,j)
!>
                  tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*                &
     &                               vmask(i,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>                v_stokes(i,j,k)=v_stokes(i,j,k)*                      &
!>   &                            vmask_wet(i,j)
!>
                  tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*                &
     &                               vmask_wet(i,j)
#  endif
# endif
                END DO
              END DO
            END IF
          END IF
!
          IF (.not.(CompositeGrid(inorth,ng).or.NSperiodic(ng))) THEN
            IF (j.eq.Mm(ng)+1) THEN                  ! northern boundary
              DO k=1,N(ng)                           ! J-loop pipelined
                DO i=Istr,Iend
!>                v(i,j,k,nnew)=v(i,j,k,nnew)-CF(i,0)
!>
                  tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)-                    &
     &                             tl_CF(i,0)
# ifdef MASKING
!>                v(i,j,k,nnew)=v(i,j,k,nnew)*                          &
!>   &                          vmask(i,j)
!>
                  tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*                    &
     &                             vmask(i,j)
# endif
# ifdef WET_DRY_NOT_YET
!>                v(i,j,k,nnew)=v(i,j,k,nnew)*                          &
!>   &                          vmask_wet(i,j)
!>
                  tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*                    &
     &                             vmask_wet(i,j)
# endif
# ifdef NEARSHORE_MELLOR
!>                v_stokes(i,j,k)=v_stokes(i,j,k)-CFs(i,0)
!>
                  tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)-                &
     &                               tl_CFs(i,0)
#  ifdef MASKING
!>                v_stokes(i,j,k)=v_stokes(i,j,k)*                      &
!>   &                            vmask(i,j)
!>
                  tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*                &
     &                               vmask(i,j)
#  endif
#  ifdef WET_DRY_NOT_YET
!>                v_stokes(i,j,k)=v_stokes(i,j,k)*                      &
!>   &                            vmask_wet(i,j)
!>
                  tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*                &
     &                               vmask_wet(i,j)
#  endif
# endif
                END DO
              END DO
            END IF
          END IF
!
!  Compute correct mass flux, Hz*v/m.
!
          DO k=N(ng),1,-1
            DO i=IstrT,IendT
              Hvom(i,j,k)=0.5_r8*(Hvom(i,j,k)+v(i,j,k,nnew)*DC(i,k))
              tl_Hvom(i,j,k)=0.5_r8*(tl_Hvom(i,j,k)+                    &
     &                               tl_v(i,j,k,nnew)*DC(i,k)+          &
     &                               v(i,j,k,nnew)*tl_DC(i,k))-         &
# ifdef TL_IOMS
     &                       0.5_r8*v(i,j,k,nnew)*DC(i,k)
# endif
# ifdef NEARSHORE_MELLOR
              Hvom(i,j,k)=Hvom(i,j,k)+0.5_r8*v_stokes(i,j,k)*DC(i,k)
              tl_Hvom(i,j,k)=tl_Hvom(i,j,k)+                            &
     &                       0.5_r8*(tl_v_stokes(i,j,k)*DC(i,k)+        &
     &                               v_stokes(i,j,k)*tl_DC(i,k))-       &
#  ifdef TL_IOMS
     &                       0.5_r8*v_stokes(i,j,k)*DC(i,k)
#  endif
# endif
              FC(i,0)=FC(i,0)+Hvom(i,j,k)
              tl_FC(i,0)=tl_FC(i,0)+tl_Hvom(i,j,k)
# ifdef DIAGNOSTICS_UV
!!            DiaV3wrk(i,j,k,M3rate)=v(i,j,k,nnew)-                     &
!!   &                               DiaV3wrk(i,j,k,M3rate)
# endif
            END DO
          END DO
          DO i=IstrT,IendT
            FC1(i)=FC(i,0)                              ! intermediate
            FC(i,0)=DC(i,0)*(FC(i,0)-DV_avg2(i,j))      ! recursive
            tl_FC(i,0)=tl_DC(i,0)*(FC1(i)-DV_avg2(i,j))+                &
     &                 DC(i,0)*(tl_FC(i,0)-tl_DV_avg2(i,j))-            &
# ifdef TL_IOMS
     &                FC(i,0)
# endif
          END DO
          DO k=1,N(ng)
            DO i=IstrT,IendT
              Hvom(i,j,k)=Hvom(i,j,k)-DC(i,k)*FC(i,0)
              tl_Hvom(i,j,k)=tl_Hvom(i,j,k)-                            &
     &                       tl_DC(i,k)*FC(i,0)-                        &
     &                       DC(i,k)*tl_FC(i,0)+                        &
# ifdef TL_IOMS
     &                       DC(i,k)*FC(i,0)
# endif
            END DO
          END DO
        END IF
      END DO
!
!-----------------------------------------------------------------------
!  Exchange boundary data.
!-----------------------------------------------------------------------
!
      IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
!>      CALL exchange_u3d_tile (ng, tile,                               &
!>   &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
!>   &                          u(:,:,:,nnew))
!>
        CALL exchange_u3d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
     &                          tl_u(:,:,:,nnew))
!>      CALL exchange_v3d_tile (ng, tile,                               &
!>   &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
!>   &                          v(:,:,:,nnew))
!>
        CALL exchange_v3d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
     &                          tl_v(:,:,:,nnew))

        CALL exchange_u3d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
     &                          Huon)
        CALL exchange_u3d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
     &                          tl_Huon)

        CALL exchange_v3d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
     &                          Hvom)
        CALL exchange_v3d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj, 1, N(ng),           &
     &                          tl_Hvom)
!
        DO k=1,2
!>        CALL exchange_u2d_tile (ng, tile,                             &
!>   &                            LBi, UBi, LBj, UBj,                   &
!>   &                            ubar(:,:,k))
!>
          CALL exchange_u2d_tile (ng, tile,                             &
     &                            LBi, UBi, LBj, UBj,                   &
     &                            tl_ubar(:,:,k))
!>        CALL exchange_v2d_tile (ng, tile,                             &
!>   &                            LBi, UBi, LBj, UBj,                   &
!>   &                            vbar(:,:,k))
!>
          CALL exchange_v2d_tile (ng, tile,                             &
     &                            LBi, UBi, LBj, UBj,                   &
     &                            tl_vbar(:,:,k))
        END DO
      END IF

# ifdef DISTRIBUTE
!>    CALL mp_exchange3d (ng, tile, iNLM, 4,                            &
!>   &                    LBi, UBi, LBj, UBj, 1, N(ng),                 &
!>   &                    NghostPoints,                                 &
!>   &                    EWperiodic(ng), NSperiodic(ng),               &
!>   &                    u(:,:,:,nnew), v(:,:,:,nnew),                 &
!>   &                    Huon, Hvom)
!>
      CALL mp_exchange3d (ng, tile, iRPM, 4,                            &
     &                    LBi, UBi, LBj, UBj, 1, N(ng),                 &
     &                    NghostPoints,                                 &
     &                    EWperiodic(ng), NSperiodic(ng),               &
     &                    tl_u(:,:,:,nnew), tl_v(:,:,:,nnew),           &
     &                    tl_Huon, tl_Hvom)
      CALL mp_exchange3d (ng, tile, iRPM, 2,                            &
     &                    LBi, UBi, LBj, UBj, 1, N(ng),                 &
     &                    NghostPoints,                                 &
     &                    EWperiodic(ng), NSperiodic(ng),               &
     &                    Huon, Hvom)
!
!>    CALL mp_exchange2d (ng, tile, iNLM, 4,                            &
!>   &                    LBi, UBi, LBj, UBj,                           &
!>   &                    NghostPoints,                                 &
!>   &                    EWperiodic(ng), NSperiodic(ng),               &
!>   &                    ubar(:,:,1), vbar(:,:,1),                     &
!>   &                    ubar(:,:,2), vbar(:,:,2))
!>
      CALL mp_exchange2d (ng, tile, iRPM, 4,                            &
     &                    LBi, UBi, LBj, UBj,                           &
     &                    NghostPoints,                                 &
     &                    EWperiodic(ng), NSperiodic(ng),               &
     &                    tl_ubar(:,:,1), tl_vbar(:,:,1),               &
     &                    tl_ubar(:,:,2), tl_vbar(:,:,2))
# endif

      RETURN
      END SUBROUTINE rp_step3d_uv_tile
#endif
      END MODULE rp_step3d_uv_mod
