  SUBROUTINE dyn3dx_tls_coef(vz,boz,uz,u,v,xec, &

                vstar)

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_tls_coef.F90  ***
!&E
!&E ** Purpose : estimate all the coefficient of the tridiagonal
!&E              linear system representative of the zonal momentum
!&E              and continuity equations (2D and 3D parts)
!&E
!&E ** Description :
!&E
!&E ** Called by : dyn3dxy.F90, dyn3dyx.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : h0,hc_sig,xec(xeu,xev),zr,boz,igdu,fu,vstar
!&E                     fu,vz    (dyn3dx_corio.h)
!&E                     xec,h0,boz   (dyn3dx_ipg_*.h)
!&E                     uz,wz,whx,ms,mn,tfiusr,dx,dy_f,dx_u,dy_u (dyn3dx_adv_*.h)
!&E                     uz,vish_xe,vish_phi,mn,ms,dx_u,dy_f  (dyn3dx_hdiff.h)
!&E                     uz,htx,z0b (dyn3dx_botstress.h)
!&E                     pa,tosx,dx_u,htx (dyn3dx_atmforc.h)
!&E                     potgene,dx_u (dyn3dx_genpot.h)
!&E                     u,v,xec(xeu,xev),igdu,hx,tfiusr,htx,hty,vis2d,fu,
!&E                     pa,potgene,intz_exp,dx,dx_u,dy_u,dy_f (dyn3dx_mtm.h)
!&E                     hex,hey,surf,dx_v,dy_u,vstar,phieau (dyn3dx_cont.h)
!&E
!&E
!&E ** Modified variables : zr,pr,yz_exp,intz_exp
!&E                         av,bv,cv,y2dv     (dyn3dx_mtm.h)
!&E                         axe,bxe,cxe,y2dxe (dyn3dx_cont.h)
!&E
!&E
!&E ** Reference :
!&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
!&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-04  (A. Mano) MPI2D
!&E       !  2010-03  (A.-C. Bennis)  Introduction of waves radiative stress
!&E       !  2011-03  (F. Dumas) bottom stress parametrisation + remove case 2 and 3 for tidal flat
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Review of the implementation of wave-current quasi-eulerian formalism
!&E       !  2014-07  (F. Ganthy) Obstruction-induced force on momentum (explicit part)
!&E       !  2015-04  (B. Thouvenin) 1DV case forcing by barotropic current or free surface gradient
!&E       !  2016-12  (B. Thouvenin, A. Fiandrino) introduction key_lagoon
!&E       !  2017-12  (F. Ganthy) Bug corrected for obstruction-induced force
!&E
!&E---------------------------------------------------------------------































































!! * Modules used
   USE parameters                              ! imin imax rsh etc
! cte(dyn3dx_adv_centered.h)
   USE comaree,    ONLY : l_genpot, potgene        ! dyn3dx_mtm.h dyn3dx_genpot.h
   USE comvars2d,  ONLY : ig,id                    ! dyn3dx_ipg_cs.h
   USE comvars2d,  ONLY : fu                       !dyn3dx_corio.h,dyn3dx_mtm.h
   USE comvars2d,  ONLY : icon,h0fond
   USE comvars2d,  ONLY : h0,fwet,fwetp              !dyn3dx_ipg_pj.h,dyn3dx_ipg_cs.h
   USE comvars2d,  ONLY : vish_phi,vish_xe         !diffh_u.h(vish_xe),dyn3dx_hdiff.h(vish_xe and vish_phi)
   USE comvars2d,  ONLY : goulet,govlet,surf       !dyn3dx_cont
   USE comvars2d,  ONLY : hex,hey,hty              !dyn3dx_cont
   USE comvars2d,  ONLY : dx_v,dy_u                !dyn3dx_cont.h,dyn3dx_adv_quickest.h
   USE comvars2d,  ONLY : htx                      !dyn3dx_cont.h, dyn3dx_atmforc.h, dyn3dx_mtm.h, tofd_x1.h, dyn3dx_botstress.h
   USE comvars2d,  ONLY : hx                       !dyn3dx_cont.h, dyn3dx_mtm.h
   USE comvars2d,  ONLY : dy_f                     !dyn3dx_adv_centered.h, dyn3dx_adv_quickest.h, dyn3dx_cont.h, dyn3dx_hdiff.h, dyn3dx_mtm.h
   USE comvars2d,  ONLY : dt,halfdt                !dyn3dx_cont.h, dyn3dx_mtm.h
   USE comvars2d,  ONLY : dx_u                     !grdpsl_xexp.h, dyn3dx_mtm.h
   USE comvars2d,  ONLY : dx                       !adv_u.h, dyn3dx_adv_centered.h, dyn3dx_adv_quickest.h, adv_uq.h,
!dyn3dx_ipg_pj.h, dyn3dx_ipg_cs.h
!diffh_u.h, dyn3dx_hdiff.h, dyn3dx_atmforc.h,
!grdpsl_xexp.h, dyn3dx_mtm.h, dyn3dx_genpot.h,
   USE comvars2d,  ONLY : frot_lin,z0b             !dyn3dx_botstress.h, tofd_x1.h
   USE comvars2d,  ONLY : hminfrot
   USE comvars2d,  ONLY : vis2d,botstressmax     !dyn3dx_mtm.h
   USE comvars2d,  ONLY : tfiusr                   !dyn3dx_adv_centered.h, dyn3dx_adv_quickest.h


   USE comcolig2d, ONLY : igdu                    ! dyn3dx_mtm.h
   USE comcolig2d, ONLY : mn,ms,whx               ! dyn3dx_adv_scheme.h,dyn3dx_hdiff.h
   USE comcolig2d, ONLY : av,bv,cv,y2dv           ! dyn3dx_mtm.h
   USE comcolig2d, ONLY : axe,bxe,cxe,y2dxe       ! dyn3dx_cont.h

   USE comcolig3d, ONLY : yz_exp                  !  dyn3dx_adv_centered.h,dyn3dx_adv_quickest.h,dyn3dx_corio.h,dyn3dx_ipg_pj.h,dyn3dx_ipg_cs.h,
!dyn3dx_hdiff.h,dyn3dx_atmforc.h,grdpsl_xexp.h,
!dyn3dx_intz.h,dyn3dx_genpot.h,dyn3dx_botstress.h
   USE comcolig3d, ONLY : intz_exp                !dyn3dx_atmforc.h,grdpsl_xexp.h,dyn3dx_intz.h,
!dyn3dx_mtm.h,dyn3dx_genpot.h,dyn3dx_botstress.h
   USE comcolig3d, ONLY : cdvit                   !tofd_x1.h, dyn3dx_botstress.h
   USE comvars3d,  ONLY :  sig,sigw                     !dyn3dx_ipg_pj.h, tofd_x1.h, dyn3dx_botstress.h
   USE comvars3d,  ONLY :  dsigw                     !adv_u.h,dyn3dx_adv_centered.h,dyn3dx_adv_quickest.h,adv_uq.h,dyn3dx_ipg_pj.h
   USE comvars3d,  ONLY :  dsigu                     !dyn3dx_atmforc.h, integration_x.h,
!dyn3dx_intz.h, dyn3dx_botstress.h
   USE comvars3d,  ONLY :  rhoref                  !dyn3dx_atmforc.h,dyn3dx_mtm.h,dyn3dx_genpot.h
   USE comvars3d,  ONLY :  wz                      !adv_u.h,dyn3dx_adv_centered.h,dyn3dx_adv_quickest.h,adv_uq.h





   USE toolgeom,   ONLY : f_hzw,f_zu,f_hzu,f_zu_x
   USE cometdyn,   ONLY : pa,tosx                  !dyn3dx_atmforc.h and dyn3dx_mtm.h
   USE comfluxdo,  ONLY : phieau                   !dyn3dx_cont.h

   USE toolmpi,    ONLY :  ex_i_rsh






   USE comvars3d, ONLY : wmo_exp_x, theta_wz
   USE toolmath,  ONLY : f_alpha



!!!FG(05/11/2018)
   USE comfluxdo, ONLY : phiriver
!!!FG(05/11/2018)

   IMPLICIT NONE

!! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),     INTENT(in) :: u,v
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),     INTENT(in) :: xec
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),     INTENT(in) :: vstar
   REAL(KIND=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: uz,vz
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: boz


!! * Local declarations
   INTEGER                          :: i,j,k
   REAL(KIND=rsh),DIMENSION(kmax)   :: vz_u                     !dyn3dx_corio.h



   REAL(KIND=rsh)    :: bms,bm,pri,tbc,dxb,dhxx,dhtx,htxee,htxeo   !dyn3dx_ipg_pj.h





   REAL(KIND=rsh)    :: wmo,wm,wp,wdudz,uc,up,um,ududx,vp,vm,vdudy,uv !dyn3dx_adv_quickest.h, adv_uq.h



   REAL(KIND=rsh)    :: vitf,cdf,tofxexp                           !dyn3dx_botstress.h


   REAL(KIND=rsh)    :: ft,gdt,vit,factft,dusdy,v_utmp             !dyn3dx_mtm.h

   REAL(KIND=rsh)    :: immer                  ! dyn3dx_wave.h

   REAL(kind=rsh)    :: inv_f_alpha 


!!!FG(05/11/2018)
   REAL(kind=rsh) :: phitot
   phitot=SUM(phiriver(2,:))   !k=2: Garonne
!!!FG(05/11/2018)
!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter DYN3DX_TLS_COEF routine'




!$OMP DO schedule(runtime)                                     &

!$OMP& private(i,j,k,vz_u                                      &


!$OMP& ,bms,bm,pri,tbc,dxb,dhxx,dhtx,htxee,htxeo               &





!$OMP& ,wmo,wm,wp,wdudz,uc,up,um,ududx,vp,vm,vdudy,uv          &




!$OMP& ,vitf,cdf,tofxexp                                       &
!$OMP& ,gdt,v_utmp,vit,factft,ft,dusdy                         &


!$OMP& ,inv_f_alpha                                            &


!$OMP& )

   DO j=ljminp2,ljmaxm1
     DO i=MAX0(ig(j),limin),MIN0(limax,id(j)-1)

       IF(igdu(i,j) == 1) THEN
!
!  debut des calculs des termes de tendances
!-------------------------------------------
! 1- coriolis
!*********    including dyn3dx_corio.h      ****************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_corio.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the Coriolis
!&E              parameter for the zonal momentum forcing (3D equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : vz,fu
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E
!&E---------------------------------------------------------------------


   vz_u(:)=0.25_rsh*(vz(:,i,j)+vz(:,i+1,j)+vz(:,i,j-1)+vz(:,i+1,j-1))
   yz_exp(:,i,j)=fu(i,j)*vz_u(:)


!********* end  including dyn3dx_corio.h      ****************************



!****************** including dyn3dx_intz  *********************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_intz.h  ***
!&E
!&E ** Purpose : vertical integration of 3D trend terms :
!&E              advection, diffusion, external and internal pressure
!&E              gradient,generator potential,coriolis, bottom stress,
!&E              wind stress...
!&E              The integrated results is sent to the 2D trilinear
!&E              linear system relative to the zonal equations
!&E
!&E ** Description :
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : yz_exp,hzutx,htx
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  1990-10  (F. Dumas ) new formulation
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E
!&E---------------------------------------------------------------------


   intz_exp(i,j)=SUM(yz_exp(:,i,j)*dsigu(:))


!*************** end including dyn3dx_intz  *********************

! removes barotropic part of coriolis term at time n (or n+1/2 depending on the order of calculation of u and v)
         yz_exp(:,i,j) = yz_exp(:,i,j)-intz_exp(i,j)



! W1 - coupling with WW3 (terms depending to the depth)


         IF(icon==2) THEN                ! 2- gradient de pression interne a n-1/2


 


!*********    including dyn3dx_ipg_pj.h      ****************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_ipg_pj.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the internal
!&E              pressure gradient for the zonal momentum forcing
!&E              (3D equation)
!&E
!&E ** Description : explicit
!&E                  vertical integration of the pressure jacobian
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xec,h0,boz,dx_u
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E
!&E---------------------------------------------------------------------

    dxb=dx_u(i,j)

!    dhxx=(h0(i+1,j)-h0(i,j))/dxb
    dhxx=(xec(i+1,j)-xec(i,j))/dxb
    htxee=h0(i+1,j)+xec(i+1,j)
    htxeo=h0(i,j)+xec(i,j)
    dhtx=(htxee-htxeo)/dxb
    pri=0.0_rlg

    k=kmax
! internal pressure gradient : dP/dx part

    bms=(htxee*boz(kmax,i+1,j)-htxeo*boz(kmax,i,j))/dxb
    pri=pri-bms*sig(kmax)
    tbc=pri+(boz(k,i,j)+boz(k,i+1,j))*0.5_rlg*(sig(k)*dhtx+dhxx)


    yz_exp(k,i,j) = yz_exp(k,i,j)+tbc

    DO k=kmax-1,1,-1

      bm=(htxee*boz(k,i+1,j)-htxeo*boz(k,i,j))/dxb
      pri=pri+(bm+bms)*0.5_rlg*dsigw(k)
      tbc=pri+(boz(k,i,j)+boz(k,i+1,j))*0.5_rlg*(sig(k)*dhtx+dhxx)


      yz_exp(k,i,j) = yz_exp(k,i,j)+tbc
      bms=bm
    END DO

!********* end  including dyn3dx_ipg_pj.h      ****************************





         ENDIF

! 3- termes d advection a l instant n-1



!**************    including dyn3dx_adv_quick.h      ***********************

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn3dx_adv_quick.h  ***
!&E
!&E ** Purpose : estimate of advection terms (3d) of the momentum
!&E              equation along the x-direction
!&E
!&E ** Description : explicit - QUICK scheme (decentre amont)
!&E
!&E ** Called by : dyn3dx_tls_coef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : uz,vz_u,tfiusr,wz,yz_exp,dx_u,dx,dy_f,dy_u
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference : Kowalik, Murty "Numerical Modeling of Ocean Dynamics"
!&E
!&E ** History :
!&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2009-02 (F. Dumas) Upwind scheme along open boundaries + QUICK at imin+2 and jmin+2
!&E       !  2009-03  (V. Garnier) land mask
!&E       !  2011-11  (F. Vandermeirsch) correct spatial extension on which we can estimate quick scheme
!&E
!&E---------------------------------------------------------------------

! ----------------
! Horizontal terms
! ----------------


   IF (i >= imin+2 .AND. i <= imax-2 .AND. j <= jmax-2 .AND. j >= jmin+3) THEN


! QUICK scheme
! ------------

     DO k=1,kmax

! u * du/dx

       uc=uz(k,i,j)
       up=0.5_rsh*(uc+ABS(uc))
       um=0.5_rsh*(uc-ABS(uc))
       ududx=(up*uz(k,i-2,j)*whx(i-1,j)+(-7.0_rsh*up-3.0_rsh*um)*uz(k,i-1,j)+  &
             (3.0_rsh*up-3.0_rsh*um)*uc+(3.0_rsh*up+7.0_rsh*um)*uz(k,i+1,j) &
             -um*uz(k,i+2,j)*whx(i+1,j))/(8.0_rsh*dx_u(i,j))

! v * du/dy

       vp=0.5_rsh*(vz_u(k)+ABS(vz_u(k)))
       vm=0.5_rsh*(vz_u(k)-ABS(vz_u(k)))
       vdudy=(vp*uz(k,i,j-2)*ms(i,j-1)*ms(i,j)+(-7.0_rsh*vp-3.0_rsh*vm)*uz(k,i,j-1)+  &
             (3.0_rsh*vp-3.0_rsh*vm)*uc+(3.0_rsh*vp+7.0_rsh*vm)*uz(k,i,j+1) &
             -vm*uz(k,i,j+2)*mn(i,j+1)*mn(i,j))/(8.0_rsh*dy_u(i,j))

! u * v * tg(phi) / R

       uv=uz(k,i,j)*vz_u(k)*tfiusr(i,j)

! add all terms

       yz_exp(k,i,j) = yz_exp(k,i,j) - ududx - vdudy + uv



     END DO

   ELSE

! Upwind scheme
! -------------

     DO k=1,kmax

! u * du/dx

       uc=uz(k,i,j)
       up=0.5_rsh*(uc+ABS(uc))
       um=0.5_rsh*(uc-ABS(uc))
       ududx=up*(uz(k,i,j)-uz(k,i-1,j))/dx(i,j)+um*(uz(k,i+1,j)-uz(k,i,j))/dx(i+1,j)

! v * du/dy

       vp=0.5_rsh*(vz_u(k)+ABS(vz_u(k)))
       vm=0.5_rsh*(vz_u(k)-ABS(vz_u(k)))
       vdudy=vp*(uz(k,i,j)-uz(k,i,j-1))/dy_f(i,j-1)+vm*(uz(k,i,j+1)-uz(k,i,j))/dy_f(i,j)

! u * v * tg(phi) / R

       uv=uz(k,i,j)*vz_u(k)*tfiusr(i,j)

! add all terms

       yz_exp(k,i,j) = yz_exp(k,i,j) - ududx - vdudy + uv



     END DO

   END IF


! -------------------------
! Vertical term : w * du/dz
! -------------------------

   k=1
     wmo=0.25_rsh*(wz(k,i,j)+wz(k,i+1,j))
     wm=0.5_rsh*(wmo-ABS(wmo))
     wdudz=wm*(uz(k+1,i,j)-uz(k,i,j))/dsigw(k)

     inv_f_alpha = f_alpha( abs(wm*halfdt/dsigw(k)) )

     wmo_exp_x(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / dsigw(k)
     wdudz = wdudz * inv_f_alpha + wmo_exp_x(k,i,j) * (1.0_rsh-theta_wz) * (uz(k+1,i,j)-uz(k,i,j))

     yz_exp(k,i,j) = yz_exp(k,i,j) - wdudz



   k=2
     wmo=0.25_rsh*(wz(k,i,j)+wz(k,i+1,j)+wz(k-1,i,j)+wz(k-1,i+1,j))
     wp=0.5_rsh*(wmo+ABS(wmo))
     wm=0.5_rsh*(wmo-ABS(wmo))
     wdudz=wm*(uz(k+1,i,j)-uz(k,i,j))/dsigw(k) &
          +wp*(uz(k,i,j)-uz(k-1,i,j))/dsigw(k-1)

     inv_f_alpha = f_alpha( abs(wm*halfdt/dsigw(k)+wp*halfdt/dsigw(k-1)) )
     wmo_exp_x(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / (dsigw(k)+dsigw(k-1))
     wdudz = wdudz * inv_f_alpha + wmo_exp_x(k,i,j) * (1.0_rsh-theta_wz) * (uz(k+1,i,j)-uz(k-1,i,j))

     yz_exp(k,i,j) = yz_exp(k,i,j) - wdudz



   DO k=3,kmax-2
     wmo=0.25_rsh*(wz(k,i,j)+wz(k,i+1,j)+wz(k-1,i,j)+wz(k-1,i+1,j))
     wp=0.5_rsh*(wmo+ABS(wmo))
     wm=0.5_rsh*(wmo-ABS(wmo))
     wdudz=(wp*uz(k-2,i,j)+(-7.0_rsh*wp-3.0_rsh*wm)*uz(k-1,i,j)+         &
           (3.0_rsh*wp-3.0_rsh*wm)*uz(k,i,j)+(3.0_rsh*wp+7.0_rsh*wm)*uz(k+1,i,j) &
           -wm*uz(k+2,i,j))/(8.0_rsh*dsigw(k))

     inv_f_alpha = f_alpha( abs(wm*halfdt/dsigw(k)+wp*halfdt/dsigw(k-1)) )

     wmo_exp_x(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / (dsigw(k)+dsigw(k-1))
     wdudz = wdudz * inv_f_alpha + wmo_exp_x(k,i,j) * (1.0_rsh-theta_wz) * (uz(k+1,i,j)-uz(k-1,i,j))

     yz_exp(k,i,j) = yz_exp(k,i,j) - wdudz



   END DO

   k=kmax-1
     wmo=0.25_rsh*(wz(k,i,j)+wz(k,i+1,j)+wz(k-1,i,j)+wz(k-1,i+1,j))
     wp=0.5_rsh*(wmo+ABS(wmo))
     wm=0.5_rsh*(wmo-ABS(wmo))
     wdudz=wm*(uz(k+1,i,j)-uz(k,i,j))/dsigw(k) &
          +wp*(uz(k,i,j)-uz(k-1,i,j))/dsigw(k-1)

     inv_f_alpha = f_alpha( abs(wm*halfdt/dsigw(k)+wp*halfdt/dsigw(k-1)) )

     wmo_exp_x(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / (dsigw(k)+dsigw(k-1))
     wdudz = wdudz * inv_f_alpha + wmo_exp_x(k,i,j) * (1.0_rsh-theta_wz) * (uz(k+1,i,j)-uz(k-1,i,j))

     yz_exp(k,i,j) = yz_exp(k,i,j) - wdudz



   k=kmax
     wmo=0.25_rsh*(wz(k-1,i+1,j)+wz(k-1,i,j))
     wp=0.5_rsh*(wmo+ABS(wmo))
     wdudz=wp*(uz(k,i,j)-uz(k-1,i,j))/dsigw(k-1)

     inv_f_alpha = f_alpha( abs(wp*halfdt/dsigw(k-1)) )

     wmo_exp_x(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / dsigw(k-1)
     wdudz = wdudz * inv_f_alpha + wmo_exp_x(k,i,j) * (1.0_rsh-theta_wz) * (uz(k,i,j)-uz(k-1,i,j))

     yz_exp(k,i,j) = yz_exp(k,i,j) - wdudz



!************** end  including dyn3dx_adv_quick.h    ****************************






! 4- terme de diffusion horizontale a n-1
!**************** including dyn3dx_hdiff  *****************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_hdiff.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to horizontal viscosity
!&E              terms for the zonal momentum equation
!&E              (3D equation)
!&E
!&E ** Description : explicit part
!&E                  (d (nu duz/dx)/dx+d (nu duz/dy)/dy)
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : yz_exp,mm,ms,dx,dx_u,dy_f,uz,vish_xe,vish_phi
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  1990-10  (F. Dumas ) new formulation
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E
!&E---------------------------------------------------------------------

!     included from    : dyn3dx_tls_coef
!                      :
!     variables in     : uz
!     variables in     : vish_phi, vish_xe(comvers2d)
!     variables in     : dx,dx_u,dy_u,dy_f (comvers2d)
!     variables in     : cne(parameters)
!     variables in     : ms, mn (comcolig2d)
!     variables inout  : yz_exp(comcolig3d)
!                      :

     yz_exp(:,i,j)=yz_exp(:,i,j)+ &
!-----------------------------------
!   d nu(du/dx) /dx
!-----------------------------------
         (vish_xe(:,i+1,j)*(uz(:,i+1,j)-uz(:,i  ,j))/dx(i+1,j)    &
         -vish_xe(:,i  ,j)*(uz(:,i  ,j)-uz(:,i-1,j))/dx(i  ,j) )  &
                 /dx_u(i,j) &
!
!-----------------------------------
!   d nu(du/dy) /dy
!-----------------------------------
!
        +(vish_phi(:,i,j  )*mn(i,j)*(uz(:,i,j+1)-uz(:,i,j  ))/dy_f(i,j  ) &
         -vish_phi(:,i,j-1)*ms(i,j)*(uz(:,i,j  )-uz(:,i,j-1))/dy_f(i,j-1) )/dy_u(i,j)


!
!**************** end including dyn3dx_hdiff  *****************************



! 5- integration des termes precedents 1 a 4
!    sur la verticale a destination du mode
!    barotrope.
!****************** including dyn3dx_intz  *********************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_intz.h  ***
!&E
!&E ** Purpose : vertical integration of 3D trend terms :
!&E              advection, diffusion, external and internal pressure
!&E              gradient,generator potential,coriolis, bottom stress,
!&E              wind stress...
!&E              The integrated results is sent to the 2D trilinear
!&E              linear system relative to the zonal equations
!&E
!&E ** Description :
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : yz_exp,hzutx,htx
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  1990-10  (F. Dumas ) new formulation
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E
!&E---------------------------------------------------------------------


   intz_exp(i,j)=SUM(yz_exp(:,i,j)*dsigu(:))


!*************** end including dyn3dx_intz  *********************

! 6- tension sur le fond a l instant n-1

!*************include  dyn3dx_botstress   **********************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_botstress.h  ***
!&E
!&E ** Purpose : estimate the drag coefficient and the trend term
!&E              relative to the bottom stress
!&E              for the zonal momentum forcing (2D and 3D equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : vz_u,uz,hzutx,htx,z0b,cdvit
!&E
!&E ** Modified variables : yz_exp,intz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original cod
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-07  (F. Dumas) Generalized sigmas
!&E       !  2013-08  (A-C. Bennis and F. Dumas) Soulsby (1995) and
!&E                      Yu and Slinn (2003) parameterizations
!&E       !  2011-12  (S. Petton, F. Dumas) Implicitation of the bottom stress (consistence between 2D en 3D)
!&E                                         removing of min,max
!&E       !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Review of the implementation of wave-current quasi-eulerian formalism
!&E       !  2016-09  (V. Garnier) Review of the formulation of cdvit for the wave-current quasi-eulerian formalism (key_wave_soulsby_bottom_stress)
!&E
!&E---------------------------------------------------------------------






   vitf=SQRT(vz_u(1)*vz_u(1)+uz(1,i,j)*uz(1,i,j))


   cdf=(0.4_rsh/LOG((sig(1)+1.0_rsh)*htx(i,j)/z0b(i,j)+1.0_rsh))**2



   cdvit(i,j)=cdf*vitf





   tofxexp=cdvit(i,j)*uz(1,i,j)/htx(i,j)*cfroe
   intz_exp(i,j)=intz_exp(i,j)-tofxexp

!************* end include  dyn3dx_botstress   **********************************


! W2 - coupling with WW3 (surface terms)


! 8- forcing meteo tension de surface et
!    pression a l instant n-1/2


!**************** include dyn3dx_atmforc *********************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_atmforc.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the met. forcing
!&E              (atmospherical pressure gradient and wind stress)
!&E              for the zonal momentum forcing (2D and 3d equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : pa,tosx,hzutx,htx,dx_u
!&E
!&E ** Modified variables : yz_exp,intz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E
!&E---------------------------------------------------------------------

! Atmospherical pressure
   yz_exp(:,i,j) = yz_exp(:,i,j)-(pa(i+1,j)-pa(i,j))/(dx_u(i,j)*rhoref)
   intz_exp(i,j) = intz_exp(i,j)-(pa(i+1,j)-pa(i,j))/(dx_u(i,j)*rhoref)

! wind stress

   yz_exp(kmax,i,j)= yz_exp(kmax,i,j)+tosx(i,j)/MAX(htx(i,j),hminfrot)/rhoref/dsigu(kmax)

   intz_exp(i,j)   = intz_exp(i,j)   +tosx(i,j)/MAX(htx(i,j),hminfrot)/rhoref




!**************** en d include dyn3dx_atmforc *********************************



         IF (l_genpot) THEN              ! 9- calcul du potentiel generateur de maree
!    a l instant n-1/2
!**************include dyn3dx_genpot  ***********************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_genpot.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the tidal potential
!&E              for the zonal momentum forcing (2D and 3D equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : potgene,dx_u
!&E
!&E ** Modified variables : yz_exp,intz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E
!&E---------------------------------------------------------------------


     yz_exp(:,i,j)=yz_exp(:,i,j)+(potgene(i+1,j)-potgene(i,j))/ &
                                 (dx_u(i,j)*rhoref)
     intz_exp(i,j)=intz_exp(i,j)+(potgene(i+1,j)-potgene(i,j))/ &
                                 (dx_u(i,j)*rhoref)

!************** end include dyn3dx_genpot  ***********************************
         ENDIF

! adds barotropic part of coriolis term estimated with velocity extrapolated in time (vstar)
         yz_exp(:,i,j)=yz_exp(:,i,j)+fu(i,j)*0.25_rlg*(vstar(i,j)+vstar(i+1,j)+vstar(i,j-1)+vstar(i+1,j-1))






         intz_exp(i,j)=intz_exp(i,j)+fu(i,j)*0.25_rlg*(vstar(i,j)+vstar(i+1,j)+vstar(i,j-1)+vstar(i+1,j-1))


       ENDIF                              ! fin test igdu

! 15- calcul des coefficients des matrices tridiagonales pour chaque colonne de la
!     grille C, issus de l implicitation des termes de gradient de pression externe
!     (equation quantite de mouvement dyn3dx_mtm) et des termes de divergence du flux
!     Nord-Sud (equation de conservation de la masse conty) et des vecteurs 2nd-membre

!&E   ***************   include dyn3dx_mtm   **********************************
!&E
!&E ** Purpose : remplissage des termes du systeme lineaire
!&E              provenant de la discretisation de l equation de conservation
!&E              de la quantite de mouvement du systeme d equations
!&E              relatif au mode barotrope du modele 3D dans la direction Est-Ouest
!&E              sur l ensemble de la grille.
!&E
!&E ** Description :
!&E
!&E ** Called by : dyn3dx_tls_coef.F90
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
!&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
!&E
!&E ** History : 2006-05
!&E              2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2010-07  (L. Debreu, F. Dumas) rewritting of ADI
!&E       !  2011-03  (F. Dumas) bottom stress parametrisation + remove case 2 and 3 for tidal flat processing
!&E       !  2011-10  (L. Pineau, F. Dumas, V. garnier) correction bottom stress parametrisation
!&E       !  2011-12  (S. Petton, F. Dumas) Implicitation of the bottom stress (consistence between 2D en 3D)
!&E                                         removing of min,max
!&E       !  2014-01  (F. Dumas, V. Garnier) Correction bottom stress parametrisation to take into account
!&E                                          the bottom velocity in the implicit part of the 2D momentum eq.
!&E       !  2015-04  (B. Thouvenin, F. Dumas) 1DV case : forcing dynamic
!&E                                          with barotropic courrent (key_1DV_force_ubarotrope or key_1DV_force_ufile )
!&E                                          or with free surface slope (key_1DV_force_gradssh)
!&E       !  2017-12  (F. Ganthy) Bug corrected for obstruction-induced force
!&E
!&E
!&E  RQ V Garnier vis2d est homogene spatialement donc OK
!&E     neanmoins, il faudrait plutot utiliser vish_xe(1,:,:) et vis
!&E     car ainsi le coef pourrait varier spatialement
!&E     Aujourd hui, faux si cela varie spatialement, il faut definir
!&E     coef de viscosite en u et v
!&E
!&E-------------------------------------------------------------------------

     gdt = grav*dt

! cas d un point u sans eau de part ni d autre
! drying u-point : no water on the right and left

     IF (igdu(i,j) == 4) THEN
       av(i,j)=0.0_rlg
       bv(i,j)=1.0_rlg
       cv(i,j)=0.0_rlg
       y2dv(i,j)=0.0_rlg


! cas d un point u avec de l eau de part et d autre
! u-point surrounded by water

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

! sea surface pressure gradient : implicit part
       av(i,j)=-alphan*0.5_rlg*gdt/dx_u(i,j)                ! i-1
       cv(i,j)= alphan*0.5_rlg*gdt/dx_u(i,j)                ! i

! time evolution (u at n+1/2)


! bottom stress (implicit part)
       ft=halfdt*cdvit(i,j)/htx(i,j)


       bv(i,j)=1.0_rlg/(                                &
! time evolution (u at n+1/2)
               1.0_rlg                                  &

! bottom stress : implicit part
               +ft*cfroi)


                 

! rhs : explicit terms
       y2dv(i,j)=                                       &
! time evolution (u at n-1/2)
                 u(i,j)                                 &
! other explicit terms integrated from the bottom to the sea surface
                 +halfdt*intz_exp(i,j)                  &
! explicit part of the sea surface pressure gradient
                 -0.5_rlg*gdt/dx_u(i,j)*alphap*(xec(i+1,j)-xec(i,j)) &
! bottom stress : taking into account the bottom velocity to correct the implicit part
                 -ft*(uz(1,i,j)-u(i,j))*cfroi


                 
     ELSE IF (igdu(i,j) == 5) THEN
! low depth

       v_utmp=0.25_rsh*(v(i,j)+v(i+1,j)+v(i,j-1)+v(i+1,j-1))

! sea surface pressure gradient : implicit part
       av(i,j)=-alphan*0.5_rlg*gdt/dx_u(i,j)                ! i-1
       cv(i,j)= alphan*0.5_rlg*gdt/dx_u(i,j)                ! i



! bottom stress (implicit part)
       vit=SQRT(u(i,j)*u(i,j)+v_utmp*v_utmp)
!       factft=1.0_rsh+vit-botstressmax+ABS(vit-botstressmax)
!       ft=factft*0.5_rsh*dt*vit/htx(i,j)*MIN((0.4_rsh/LOG(MAX(MIN(htx(i,j),hmaxfrot),hminfrot)/(exp1*z0b(i,j))))**2,0.1_rsh)
       ft=0.5_rsh*dt*vit/htx(i,j)*(0.4_rsh/LOG(htx(i,j)/(exp1*z0b(i,j))+1.0_rsh))**2


       bv(i,j)=1.0_rlg/(                                &
! time evolution (u at n+1/2)
               1.0_rlg                                  &

! bottom stress : implicit part
               +ft*cfroi)


! rhs : explicit terms
       y2dv(i,j)=                                       &
! time evolution (u at n-1/2)
                 u(i,j)                                 &
! bottom stress (u at n) explicit part
                 -ft*cfroe*u(i,j)                       &
! coriolis
                 +fu(i,j)*v_utmp*halfdt                 &
! explicit part of the sea surface pressure gradient
                 -0.5_rlg*gdt/dx_u(i,j)*alphap*(xec(i+1,j)-xec(i,j))      &


! NL advection : u*v*tan(phi)/R (spherical coordinates)
                 +halfdt*u(i,j)*v_utmp*tfiusr(i,j)      &
! non-linear (udu/dx) : explicit part
                 -halfdt*u(i,j)/(dx(i,j)+dx(i+1,j))*(u(i+1,j)-u(i-1,j))   &


! viscosity (along x)
                 +halfdt*( (vis2d(i+1,j)*(u(i+1,j)-u(i  ,j)))/dx(i+1,j)   &
                          -(vis2d(i  ,j)*(u(i  ,j)-u(i-1,j)))/dx(i  ,j) )/dx_u(i,j)



! viscosity and non-linear term
       IF ((hty(i,j-1)+hty(i+1,j-1)) >= 0.0_rsh) THEN
         dusdy=(u(i,j  )-u(i,j-1))/dy_f(i,j-1)
         y2dv(i,j)=y2dv(i,j)                            &


! non-linear (vdu/dy) : explicit part
                   -0.25_rlg*dt*v_utmp*dusdy             &


! viscosity (along y) : explicit
                   -halfdt*vis2d(i,j-1)*dusdy/dy_u(i,j)
       END IF
       IF ((hty(i,j  )+hty(i+1,j  )) >= 0.0_rsh) THEN
         dusdy=(u(i,j+1)-u(i,j  ))/dy_f(i,j)
         y2dv(i,j)=y2dv(i,j)                            &


! non-linear (vdu/dy) : explicit part
                   -0.25_rlg*dt*v_utmp*dusdy             &


! viscosity (along y) : explicit
                   +halfdt*vis2d(i,j  )*dusdy/dy_u(i,j)
       END IF

! wind stress and atmospherical pressure gradient
       y2dv(i,j)=y2dv(i,j)                              &

! wind stress in x-direction : explicit
                 +halfdt*tosx(i,j)/MAX(htx(i,j),hminfrot)/rhoref &

! atmospherical pressure gradient : explicit
                 -halfdt*(pa(i+1,j)-pa(i,j))/(dx_u(i,j)*rhoref)



! generator potential : explicit
       IF (l_genpot) &
         y2dv(i,j)=y2dv(i,j)+halfdt/dx_u(i,j)*(potgene(i+1,j)-potgene(i,j))

     END IF 

!&E   ***************   end include dyn3dx_mtm   ******************************
!&E   ***************   include dyn3dx_cont.h   **********************************
!&E
!&E ** Purpose : remplissage des termes du systeme lineaire
!&E              provenant de la discretisation de l equation de
!&E              continuite dans le systeme d equation de Saint-
!&E              Venant sur l ensemble de la grille de calcul.
!&E
!&E ** Description : Equation couplee a celle de la conservation de la
!&E                  quantite de mouvement dans la direction Ouest-Est.
!&E
!&E
!&E ** Called by : dyn2dx
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
!&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
!&E
!&E ** History : 2006-05
!&E              2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-08  (T. Duhaut) time scheme Rq : schema des traceurs
!&E                               semble instable si u est decomposee en
!&E                               alpha u et (alpha - 1) u
!&E       !  2013-08  (AC. Bennis) Add Stokes drift divergence
!&E       !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Review of the implementation of wave-current quasi-eulerian formalism
!&E       !  2015-04  (B. Thouvenin, F. Dumas) 1DV case => axe and cxe=0 (key_1DV)
!&E
!&E-------------------------------------------------------------------------

! 1- terme de diagonale de la matrice tridiagonale
!    du modele barotrope
!--------------------------------------------------

! 2- terme de sous-diagonale de la matrice tridiagonale
!    issu de l implicitation de la divergence de flux zonale
!    (flux a l ouest du point xe(i,j) considere). Remplissage
!    du tableau hex qui ressert dans les flux de traceurs
!------------------------------------------------------------


      axe(i,j)=-halfdt*dy_u(i-1,j)/surf(i,j)*hex(i-1,j)

      bxe(i,j)=fwet(i,j)

! 3- terme de sous-diagonale de la matrice tridiagonale
!    issu de l implicitation de la divergence de flux zonale
!    (flux a l est du point xe(i,j) considere). Remplissage
!    du tableau hex qui ressert dans les flux de traceurs
!-----------------------------------------------------------

       cxe(i,j)=halfdt*dy_u(i,j)/surf(i,j)*hex(i,j)


! 4- terme du vecteur second membre du systeme lineaire
!    issu de la divergence de flux meridienne
!    (flux au nord et sud  du point xe(i,j) considere)
!    Prise en compte de tous les apports d eau dans
!    la colonne d eau. Remplissage
!    du tableau hey qui ressert dans les flux de traceurs
!----------------------------------------------------------

       y2dxe(i,j)=xec(i,j)*fwetp(i,j)-halfdt/surf(i,j)*(              &
                           ( hey(i,j  )*vstar(i,j  )*dx_v(i,j  )      &
                           - hey(i,j-1)*vstar(i,j-1)*dx_v(i,j-1) )    &
                           - SUM(phieau(:,i,j))                       &

                                ) 

               

!&E   ***************   end include dyn3dx_cont.h   ******************************

     END DO
   END DO
!$OMP END DO



   !PRINT*, 'exit DYN3DX_TLS_COEF routine'

  END SUBROUTINE dyn3dx_tls_coef
