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

  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---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d, ONLY: ig,id,hx,htx,halfdt,nbouc,z0b
   USE comvars3d, ONLY: nz,dsigu,dsigw,sig

   USE comcolig2d,ONLY : igdu
   USE comcolig3d,ONLY : cdvit,yz_exp
   USE toolgeom,  ONLY : f_zu_x




  USE comvars3d, ONLY : wmo_exp_x,theta_wz


   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

   REAL(kind=rsh)                                     :: d

   REAL(kind=rlg),DIMENSION(kmax)                     :: cz,bz,az,pz0,qz0
   REAL(kind=rlg),DIMENSION(kmax)                     :: u0


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

   !PRINT*, 'enter DYN3DX_ZTLS_SOLVE routine'

!$OMP Do schedule(runtime) private(i,j,k,nzms          &

!$OMP&   ,d                                            &


!$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

! 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

       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
!-----------------------------------------------------



         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
!----------------------------------------------------------


         d=htx(i,j)


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

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


         az(k)=-nzms*halfdt/(d*d*dsigu(k)*dsigw(k))
         bz(k)=1.0_rlg+cdvit(i,j)*halfdt/(dsigu(k)*d)-az(k)


         cz(k)=0.0_rlg


         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


         nzmi=nzms

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

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


           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)




           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


           nzmi=nzms

         END DO

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


         cz(k)=-nzmi*halfdt/(d*d*dsigu(k)*dsigw(k-1))

         bz(k)=1.0_rlg-cz(k)



         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


!-------------------------------------------------
! 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
 



       END IF

     END DO
   END DO
!$OMP end do

   !PRINT*, 'exit DYN3DX_ZTLS_SOLVE routine'

  END SUBROUTINE dyn3dx_ztls_solve
