   !!======================================================================

  SUBROUTINE dyn3dx_ztls_solve(uz,uv)

   !&E---------------------------------------------------------------------
   !&E                 ***  dyn3dy_ztls_solve.F90  ***
   !&E
   !&E ** Purpose : estimate the 3D zonal velocity all along the vertical
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : dyn3dxy.F90, dyn3dyx.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : igdu,uv,uz,yz_exp,htx,nz,cdvit,hzutx,hzwtx 
   !&E 
   !&E ** Modified variables : uz
   !&E 
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  1987-02  (P. Lazure)  Original code
   !&E       !  1990-10  (F. Dumas ) new formulation
   !&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
   !&E       !  2009-04  (A. Mano) MPI2D
   !&E       !  2011-12  (S. Petton, F. Dumas) log profile for the velocity over small depths (ig=5)
   !&E       !  2011-12  (S. Petton) barotropic field for the velocity over small depths (igd=5 and hz(1)<z0fond)
   !&E       !  2014-07  (F. Ganthy) Obstruction-induced force on momentum (implicit part)
   !&E       !  2017-12  (F. Ganthy) Bug corrected for obstruction-induced force 
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY: ig,id,hx,htx,halfdt,nbouc,z0b
   USE comvars3d, ONLY: nz,dsigu,dsigw,sig
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hzutx,hzwtx
#endif
   USE comcolig2d,ONLY : igdu
   USE comcolig3d,ONLY : cdvit,yz_exp
   USE toolgeom,  ONLY : f_zu_x
#ifdef key_diagdyn
   USE diagdyn,   ONLY : mask_diagx,diagx_vdif,diagx_dudt
#endif  
#ifdef key_obstructions
  USE obstructions,      ONLY : obst_fuz_i,obst_fu_i
#endif
#ifdef key_wave_qeuler_diag
  USE wavediag,  ONLY :  diag_wf_botfricx,diag_wf_dudt
#endif
#ifdef key_adapt_imp_vert_adv
  USE comvars3d, ONLY : wmo_exp_x,theta_wz
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(     liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in)    :: uv
   REAL(kind=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(inout) :: uz

   !! * Local declarations
   INTEGER                                            :: i,j,k
   REAL(kind=rsh)                                     :: nzmi,nzms
#if ! defined key_siggen && ! defined key_gencoord
   REAL(kind=rsh)                                     :: d
#endif
   REAL(kind=rlg),DIMENSION(kmax)                     :: cz,bz,az,pz0,qz0
   REAL(kind=rlg),DIMENSION(kmax)                     :: u0
#if defined key_diagdyn
   REAL(kind=rsh),DIMENSION(1:kmax)  :: z_prev
#endif

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter DYN3DX_ZTLS_SOLVE routine'

!$OMP Do schedule(runtime) private(i,j,k,nzms          &
#if ! defined key_siggen && ! defined key_gencoord
!$OMP&   ,d                                            &
#endif
#if defined key_diagdyn
!$OMP& ,z_prev                                         &
#endif
!$OMP&   ,nzmi,cz,bz,az,pz0,qz0,u0)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)

       !---------------------------------------------
       !  pas d eau pas de vitesse
       !---------------------------------------------

       IF(igdu(i,j) == 4) uz(:,i,j)=0.0_rsh

       !---------------------------------------------
       !  peu d eau ou banc decouvrant :
       ! la structure 3D est un profil log issu de l ecoulement barotrope.
       !---------------------------------------------

!       IF(igdu(i,j) == 2.OR.igdu(i,j) == 3.OR.igdu(i,j) == 5) uz(:,i,j)=uv(i,j)
       IF(igdu(i,j) == 5) THEN
#ifdef key_no_bottom_stress
         uz(:,i,j)=uv(i,j)
#else
            ! rq : we use the old ssh
            ! depth from the bottom : hx + f_zu_x(hx,xe,k,i,j)
         IF (hx(i,j)  + f_zu_x(hx(i,j),htx(i,j)-hx(i,j),1,i,j) <= z0b(i,j)) THEN
           uz(:,i,j)=uv(i,j)
         ELSE
           DO k=1,kmax
              uz(k,i,j)=uv(i,j)* LOG( ( hx(i,j) + f_zu_x(hx(i,j),htx(i,j)-hx(i,j),k,i,j) ) /z0b(i,j) ) & 
                               / LOG( htx(i,j) / (exp1*z0b(i,j)) )
           END DO
         END IF
#endif
       END IF

       !-------------------------------------------------
       ! cas en pleine eau. Implicitation de la viscosite
       ! verticale => systeme tridiagonal a inverser
       !-------------------------------------------------

       IF(igdu(i,j) == 1) THEN

         ! intermediate estimate of the 3d zonal velocity
         ! at this step, the estimate is fully explicit,
         ! the vertical diffusivity (implicit) has not been taken
         ! into account
         !-----------------------------------------------------

#ifdef key_diagdyn
         IF (mask_diagx(i,j)) z_prev(:)=uz(:,i,j)
#endif 

         uz(:,i,j)=uz(:,i,j)+halfdt*yz_exp(:,i,j)

         ! Solve the vertical tridiagonal system taking the internal
         ! vertical diffusivity or mixing into account
         !----------------------------------------------------------

#if ! defined key_siggen && ! defined key_gencoord
         d=htx(i,j)
#endif

         !******************************
         k=1
         !******************************

         nzms=0.5*(nz(k,i,j)+nz(k,i+1,j))

#if defined key_siggen || defined key_gencoord
         az(k)=-nzms*halfdt/(hzutx(k,i,j)*dsigu(k)*hzwtx(k,i,j)*dsigw(k))
         bz(k)=1.0_rlg+cdvit(i,j)*halfdt/(hzutx(k,i,j)*dsigu(k))-az(k)
#else
         az(k)=-nzms*halfdt/(d*d*dsigu(k)*dsigw(k))
         bz(k)=1.0_rlg+cdvit(i,j)*halfdt/(dsigu(k)*d)-az(k)
#endif
#ifdef key_obstructions
         ! obstruction force : 3D implicit Part
         bz(k) = bz(k) + obst_fuz_i(k,i,j)*halfdt
#endif
         cz(k)=0.0_rlg

#ifdef key_adapt_imp_vert_adv
         az(k)=az(k)+halfdt*wmo_exp_x(k,i,j)*theta_wz
         bz(k)=bz(k)-halfdt*wmo_exp_x(k,i,j)*theta_wz
#endif

         nzmi=nzms

         !******************************
         DO k=2,kmax-1
         !******************************

           nzms=0.5_rsh*(nz(k,i,j)+nz(k,i+1,j))

#if defined key_siggen || defined key_gencoord
           cz(k)=-nzmi*halfdt/(hzutx(k,i,j)*dsigu(k)*hzwtx(k-1,i,j)*dsigw(k-1))
           az(k)=-nzms*halfdt/(hzutx(k,i,j)*dsigu(k)*hzwtx(k,i,j)*dsigw(k))
           bz(k)=1.0_rlg-az(k)-cz(k)
#else
           cz(k)=-nzmi*halfdt/(d*d*dsigu(k)*dsigw(k-1))
           az(k)=-nzms*halfdt/(d*d*dsigu(k)*dsigw(k))
           bz(k)=1.0_rlg-az(k)-cz(k)
#endif
#ifdef key_obstructions
           ! obstruction force : 3D implicit Part
           bz(k) = bz(k) + obst_fuz_i(k,i,j)*halfdt
#endif

#ifdef key_adapt_imp_vert_adv
           az(k)=az(k)+halfdt*wmo_exp_x(k,i,j)*theta_wz
           cz(k)=cz(k)-halfdt*wmo_exp_x(k,i,j)*theta_wz
#endif

           nzmi=nzms

         END DO

         !************************
         k=kmax
         !************************

#if defined key_siggen || defined key_gencoord
         cz(k)=-nzmi*halfdt/(hzutx(k,i,j)*dsigu(k)*hzwtx(k-1,i,j)*dsigw(k-1))
#else
         cz(k)=-nzmi*halfdt/(d*d*dsigu(k)*dsigw(k-1))
#endif
         bz(k)=1.0_rlg-cz(k)
#ifdef key_obstructions
         ! obstruction force : 3D explicit Part
         bz(k) = bz(k) + obst_fuz_i(k,i,j)*halfdt
#endif

#ifdef key_adapt_imp_vert_adv
         bz(k)=bz(k)+halfdt*wmo_exp_x(k,i,j)*theta_wz
         cz(k)=cz(k)-halfdt*wmo_exp_x(k,i,j)*theta_wz
#endif

         !-------------------------------------------------
         ! inversion  du syteme tridiagonal obtenu.
         !-------------------------------------------------

         qz0(kmax)=-cz(kmax)/bz(kmax)
         pz0(kmax)= bz(kmax)
         DO k=kmax-1,1,-1
           pz0(k)=az(k)*qz0(k+1)+bz(k)
           qz0(k)=-cz(k)/pz0(k)
         END DO

         u0(kmax)=uz(kmax,i,j)/pz0(kmax)
         DO k=kmax-1,1,-1
           u0(k)=(uz(k,i,j)-az(k)*u0(k+1))/pz0(k)
         END DO

         uz(1,i,j)=u0(1)
         DO k=2,kmax
           uz(k,i,j)=qz0(k)*uz(k-1,i,j)+u0(k)
         END DO
 
#ifdef key_diagdyn
         IF (mask_diagx(i,j)) THEN 
           k=1
           diagx_vdif(k,i,j)=-(az(k)*uz(k+1,i,j)+(bz(k)-1.)*uz(k,i,j))/halfdt ! implicit bottom stress (in bz)
           diagx_dudt(k,i,j)=(uz(k,i,j)-z_prev(k))/halfdt
           DO k=2,kmax-1
             diagx_vdif(k,i,j)=-(az(k)*uz(k+1,i,j)+(bz(k)-1.)*uz(k,i,j)+cz(k)*uz(k-1,i,j))/halfdt
             diagx_dudt(k,i,j)=(uz(k,i,j)-z_prev(k))/halfdt
           END DO
           k=kmax
           ! cumul diffusion verticale explicite (dyn3dx_tls_coef) et implicit (ci dessous)
           diagx_vdif(k,i,j)=diagx_vdif(k,i,j)-((bz(k)-1.)*uz(k,i,j)+cz(k)*uz(k-1,i,j))/halfdt
           diagx_dudt(k,i,j)=(uz(k,i,j)-z_prev(k))/halfdt
         END IF
#endif  
#ifdef key_wave_qeuler_diag
         diag_wf_botfricx(i,j)=cdvit(i,j)*uz(1,i,j)
#endif

       END IF

     END DO
   END DO
!$OMP end do

   PRINT_DBG*, 'exit DYN3DX_ZTLS_SOLVE routine'

  END SUBROUTINE dyn3dx_ztls_solve
