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

  SUBROUTINE dyn3dy_ztls_solve(vz,uv)

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_ztls_solve.F90  ***
!&E
!&E ** Purpose : estimate the 3D meridional 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 : igdv,uv,yz_exp,hty,nz,hzuty,hzwty,cdvit,vz
!&E
!&E ** Modified variables : vz
!&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 (igd=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
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d, ONLY : hy,hty,halfdt,nbouc,z0b
   USE comvars3d, ONLY : dsigu,dsigw

   USE comcolig2d,ONLY : igdv
   USE comcolig3d,ONLY : cdvit,yz_exp
   USE comvars3d, ONLY : nz
   USE toolgeom,  ONLY : f_zu_y




  USE comvars3d, ONLY : wmo_exp_y,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) :: vz

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


!$OMP DO schedule(runtime),private(u0,cz,bz,az,pz0,qz0,i,j,k       &

!$OMP&   ,d                                                        &


!$OMP&   ,nzmi,nzms)
   DO j=ljminp1,ljmaxm1
   DO i=liminp2,limaxm1
!
!-----------------------------
! si pas d eau pas de vitesse
!-----------------------------
!
     IF(igdv(i,j)==4) vz(:,i,j)=0.0_rsh
!
!---------------------------------------------------
! si bancs decouvrants on ne fait que du 2D et le
! profil sur la verticale est un profil log deduit du courant barotrope
!---------------------------------------------------
!
!     IF(igdv(i,j)==2 .OR. igdv(i,j)==3 .OR. igdv(i,j)==5) vz(:,i,j)=uv(i,j)
       IF(igdv(i,j) == 5) THEN
! rq : we use the old ssh
! depth from the bottom : hy + f_zu_y(hy,xe,k,i,j)

         IF (hy(i,j) + f_zu_y(hy(i,j),hty(i,j)-hy(i,j),1,i,j)  <= z0b(i,j)) THEN
           vz(:,i,j)=uv(i,j)
         ELSE
           DO k=1,kmax
              vz(k,i,j)=uv(i,j)* LOG( ( hy(i,j) + f_zu_y(hy(i,j),hty(i,j)-hy(i,j),k,i,j) ) /z0b(i,j) ) &
                               / LOG( hty(i,j) / (exp1*z0b(i,j))  )
           END DO
         END IF

       END IF
!
     IF(igdv(i,j)==1) THEN

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



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

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

!---------------------------------------------------
! si de l eau traitement implicite sur la verticale
! de la viscosite => systeme tridiagonale
!---------------------------------------------------

       d=hty(i,j)


       k = 1
!----

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


       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_y(k,i,j)*theta_wz
       bz(k)=bz(k)-halfdt*wmo_exp_y(k,i,j)*theta_wz

       nzmi=nzms

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

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

         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_y(k,i,j)*theta_wz
         cz(k)=cz(k)-halfdt*wmo_exp_y(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_y(k,i,j)*theta_wz
         cz(k)=cz(k)-halfdt*wmo_exp_y(k,i,j)*theta_wz


!-----------------------------------
! inversion du systeme tridiagonal
!-----------------------------------
       pz0(kmax)= bz(kmax)
       qz0(kmax)=-cz(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)=vz(kmax,i,j)/pz0(kmax)
       DO k=kmax-1,1,-1
         u0(k)=(vz(k,i,j)-az(k)*u0(k+1))/pz0(k)
       END DO

       vz(1,i,j)=u0(1)
       DO k=2,kmax
         vz(k,i,j)=qz0(k)*vz(k-1,i,j)+u0(k)
       ENDDO




     END IF
   END DO
   END DO
!$OMP end do


END SUBROUTINE dyn3dy_ztls_solve
