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

                   ustar)

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_tls_coef.F90  ***
!&E
!&E ** Purpose : estimate all the coefficient of the tridiagonal
!&E              linear system representative of the meridional 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,igdv,fv,ustar
!&E                     fv,uz    (dyn3dy_corio.h)
!&E                     xec,h0,boz   (dyn3dy_ipg_*.h)
!&E                     vz,wz,why,mo,me,tfivsr,dx_v,dy_v,dx_f,dy (dyn3dy_adv_*.h)
!&E                     vz,vish_xe,vish_phi,me,mo,dy,dx_f  (dyn3dy_hdiff.h)
!&E                     vz,hty,z0b (dyn3dy_botstress.h)
!&E                     pa,tosy,dy_v,hty (dyn3dy_atmforc.h)
!&E                     potgene,dy_v (dyn3dy_genpot.h)
!&E                     u,v,xec(xeu,xev),igdv,hy,tfivsr,htx,hty,vis2d,fv,
!&E                     pa,potgene,intz_exp,dx_f,dx_v,dy,dy_v (dyn3dy_mtm.h)
!&E                     hex,hey,surf,dx_v,dy_u,ustar,phieau (dyn3dy_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-08  (T. Odaka) MPI
!&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       !  2017-08  (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(dyn3dy_adv_centered.h)
   USE comaree,    ONLY : l_genpot,potgene        ! l_genpot, potgene(dyn3dy_genpot)
   USE comvars2d,  ONLY : ig,id,icon,h0fond
   USE comvars2d,  ONLY : h0,hx,hy                 ! dyn3dy_ipg_pj.h(h0)
   USE comvars2d,  ONLY : fwet,fwetp
   USE comvars2d,  ONLY : dx_v,dy_u                ! dyn3dy_cont.h
   USE comvars2d,  ONLY : dy_v,dx_f                ! dyn3dy_mtm.h
   USE comvars2d,  ONLY : dy,tfivsr,fv             ! dyn3dy_ipg_pj.h(h0,hx)
   USE comvars2d,  ONLY : vish_xe,vish_phi         ! dyn3dy_adv_centered.h ( dy,tfiusr)
   USE comvars2d,  ONLY : htx,surf,govlet
   USE comvars2d,  ONLY : jb,jh,vis2d,goulet,hey,hex,z0b,hminfrot,hty,botstressmax,dt,halfdt

   USE comcolig2d, ONLY : igdv                     ! fv(dyn3dy_corio)
   USE comcolig2d, ONLY : y2dxe,y2dv,axe,bxe,cxe,cv,bv,av,y2d_expv    ! igdv
   USE comcolig2d, ONLY : me,mo,why                ! dyn3dy_adv_scheme.h
   USE comcolig3d, ONLY : yz_exp,intz_exp          ! yz_exp(dyn3dy_corio,dyn3dy_ipg_pj.h,dyn3dx_adv_centered.h)
   USE comcolig3d, ONLY : cdvit


   USE comvars3d,  ONLY : sig,dsigu,sigw,wz,dsigw,rhoref       ! dyn3dy_ipg_pj.h(sig,dsigw)





   USE toolgeom,   ONLY : f_hzw,f_zu,f_hzu,f_zu_y
   USE cometdyn,   ONLY : pa,tosy                   ! forcmet
   USE comfluxdo,  ONLY : phieau                   ! phieau(dyn3dy_cont)

   USE toolmpi,    ONLY :  ex_j_rsh




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



   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) :: ustar
   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,j1,j2
   REAL(KIND=rsh),DIMENSION(kmax)   :: uz_v                     !dyn3dy_corio.h



   REAL(KIND=rsh)    :: bms,bm,pri,tbc,dyb,dhyy,dhty,htxen,htxes !dyn3dy_ipg_pj.h





   REAL(KIND=rsh)    :: wmo,wm,wp,wdvdz,vp,um,vc,vdvdy,up,vm,udvdx,uv !dyn3dy_adv_quickest.h



   REAL(KIND=rsh)    :: vitf,cdf,tofyexp                              !dyn3dy_botstress.h


   REAL(KIND=rsh)    :: ft,gdt,vit,factft,dvsdx,u_vtmp                !dyn3dy_mtm.h

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

   REAL(kind=rsh)    :: inv_f_alpha


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

   !PRINT*, 'enter DYN3DY_TLS_COEF routine'



!$OMP DO schedule(dynamic,5)                                    &
!$OMP& private(i,j,k,uz_v                                       &


!$OMP& ,bms,bm,pri,tbc,dyb,dhyy,dhty,htxen,htxes                &





!$OMP& ,vm,vp,vdvdy,up,udvdx,wm,wp,wdvdz,wmo,um,vc,uv           &




!$OMP& ,vitf,cdf,tofyexp                                        &
!$OMP& ,ft,gdt,vit,factft,dvsdx,u_vtmp                          &
!$OMP& ,j1,j2                                                   &


!$OMP& ,inv_f_alpha                                             &


!$OMP& )
!do j = jmin+1, jmax-1
   DO j=ljminp1,ljmaxm1
!     DO i=liminp2,limaxm1
     DO i=MAX0(limin,2+ishift_obc_west), MIN0(limax,imax-1-ishift_obc_east)
       IF(igdv(i,j) == 1) THEN
!
!  debut des calculs des termes de tendances
!-------------------------------------------
! 1- coriolis a l instant n-1/2
!************    including dyn3dy_corio.h       ***********************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_corio.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the Coriolis
!&E              parameter for the meridional 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 : uz,fv
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E
!&E---------------------------------------------------------------------


   uz_v(:)=0.25_rsh*(uz(:,i,j)+uz(:,i,j+1)+uz(:,i-1,j)+uz(:,i-1,j+1))
   yz_exp(:,i,j)=-fv(i,j)*uz_v(:)




!*************    end including dyn3dy_corio.h     ************



!****************** including dyn3dy_intz  *********************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_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 meridional equations
!&E
!&E ** Description :
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : yz_exp,hzuty,hty
!&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 dyn3dy_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 dyn3dy_ipg_pj.h      *******************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_ipg_pj.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the internal
!&E              pressure gradient for the meridional momentum forcing
!&E              (3D equation)
!&E
!&E ** Description : explicit
!&E                  vertical integration of the pressure jacobian
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xec,h0,boz,dx_v
!&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---------------------------------------------------------------------

    dyb=dy_v(i,j)

!    dhyy=(h0(i,j+1)-h0(i,j))/dyb
    dhyy=(xec(i,j+1)-xec(i,j))/dyb
    htxen=h0(i,j+1)+xec(i,j+1)
    htxes=h0(i,j)+xec(i,j)
    dhty=(htxen-htxes)/dyb
    pri=0.0_rlg

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

    bms=(htxen*boz(k,i,j+1)-htxes*boz(k,i,j))/dyb
    pri=pri-bms*sig(kmax)
    tbc=pri+(boz(k,i,j)+boz(k,i,j+1))*0.5_rlg*(sig(k)*dhty+dhyy)


    yz_exp(k,i,j)=yz_exp(k,i,j)+tbc
!
    DO k=kmax-1,1,-1

      bm=(htxen*boz(k,i,j+1)-htxes*boz(k,i,j))/dyb
      pri=pri+(bm+bms)*0.5_rlg*dsigw(k)
      tbc=pri+(boz(k,i,j)+boz(k,i,j+1))*0.5_rlg*(sig(k)*dhty+dhyy)

      bms=bm

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

!********* end  including dyn3dy_ipg_pj.h      ***********************





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



!**************    including dyn3dy_adv_quick.h      ***********************

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn3dy_adv_quick.h  ***
!&E
!&E ** Purpose : estimate of advection terms (3d) of the momentum
!&E              equation along the y-direction
!&E
!&E ** Description : explicit - QUICK scheme (decentre amont)
!&E
!&E ** Called by : dyn3dy_tls_coef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : uz_v,vz,tfivsr,wz,yz_exp,dx_v,dy,dx_f,dy_v
!&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+3 .AND. i <= imax-2 .AND. j <= jmax-2 .AND. j >= jmin+2) THEN 


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

     DO k=1,kmax

! v * dv/dy

       vc=vz(k,i,j)
       vp=0.5_rsh*(vc+ABS(vc))
       vm=0.5_rsh*(vc-ABS(vc))
       vdvdy=(vp*vz(k,i,j-2)*why(i,j-1)+(-7.0_rsh*vp-3.0_rsh*vm)*vz(k,i,j-1)+  &
             (3.0_rsh*vp-3.0_rsh*vm)*vc+(3.0_rsh*vp+7.0_rsh*vm)*vz(k,i,j+1) &
             -vm*vz(k,i,j+2)*why(i,j+1))/(8.0_rsh*dy_v(i,j))

! u * dv/dx

       up=0.5_rsh*(uz_v(k)+ABS(uz_v(k)))
       um=0.5_rsh*(uz_v(k)-ABS(uz_v(k)))
       udvdx=(up*vz(k,i-2,j)*mo(i-1,j)*mo(i,j)+(-7.0_rsh*up-3.0_rsh*um)*vz(k,i-1,j)+  &
             (3.0_rsh*up-3.0_rsh*um)*vc+(3.0_rsh*up+7.0_rsh*um)*vz(k,i+1,j) &
             -um*vz(k,i+2,j)*me(i+1,j)*me(i,j))/(8.0_rsh*dx_v(i,j))

! u**2 * tg(phi) / R

       uv=uz_v(k)*uz_v(k)*tfivsr(i,j)

! add all terms

       yz_exp(k,i,j) = yz_exp(k,i,j) - vdvdy - udvdx - uv




     END DO

   ELSE

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

     DO k=1,kmax

! v * dv/dy

       vc=vz(k,i,j)
       vp=0.5_rsh*(vc+ABS(vc))
       vm=0.5_rsh*(vc-ABS(vc))
       vdvdy=vp*(vz(k,i,j)-vz(k,i,j-1))/dy(i,j)+vm*(vz(k,i,j+1)-vz(k,i,j))/dy(i,j+1)

! u * dv/dx

       up=0.5_rsh*(uz_v(k)+ABS(uz_v(k)))
       um=0.5_rsh*(uz_v(k)-ABS(uz_v(k)))
       udvdx=up*(vz(k,i,j)-vz(k,i-1,j))/dx_f(i-1,j)+um*(vz(k,i+1,j)-vz(k,i,j))/dx_f(i,j)

! u**2 * tg(phi) / R

       uv=uz_v(k)*uz_v(k)*tfivsr(i,j)

! add all terms

       yz_exp(k,i,j) = yz_exp(k,i,j) - vdvdy - udvdx - uv



     END DO

   END IF


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

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

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

     wmo_exp_y(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / dsigw(k)
     wdvdz = wdvdz * inv_f_alpha + wmo_exp_y(k,i,j) * (1.0_rsh-theta_wz) * (vz(k+1,i,j)-vz(k,i,j))

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



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

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

     wmo_exp_y(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / (dsigw(k)+dsigw(k-1))
     wdvdz = wdvdz * inv_f_alpha + wmo_exp_y(k,i,j) * (1.0_rsh-theta_wz) * (vz(k+1,i,j)-vz(k-1,i,j))

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



   DO k=3,kmax-2
     wmo=0.25_rsh*(wz(k,i,j)+wz(k,i,j+1)+wz(k-1,i,j)+wz(k-1,i,j+1))
     wp=0.5_rsh*(wmo+ABS(wmo))
     wm=0.5_rsh*(wmo-ABS(wmo))
     wdvdz=(wp*vz(k-2,i,j)+(-7.0_rsh*wp-3.0_rsh*wm)*vz(k-1,i,j)+         &
           (3.0_rsh*wp-3.0_rsh*wm)*vz(k,i,j)+(3.0_rsh*wp+7.0_rsh*wm)*vz(k+1,i,j) &
     	   -wm*vz(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_y(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / (dsigw(k)+dsigw(k-1))
     wdvdz = wdvdz * inv_f_alpha + wmo_exp_y(k,i,j) * (1.0_rsh-theta_wz) * (vz(k+1,i,j)-vz(k-1,i,j))

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



   END DO

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

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

     wmo_exp_y(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / (dsigw(k)+dsigw(k-1))
     wdvdz = wdvdz * inv_f_alpha + wmo_exp_y(k,i,j) * (1.0_rsh-theta_wz) * (vz(k+1,i,j)-vz(k-1,i,j))

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



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

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

     wmo_exp_y(k,i,j)=(1.0_rsh-inv_f_alpha) * wmo / dsigw(k-1)
     wdvdz = wdvdz * inv_f_alpha + wmo_exp_y(k,i,j) * (1.0_rsh-theta_wz) * (vz(k,i,j)-vz(k-1,i,j))

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



!*****************     end  including dyn3dy_adv_quick.h      ******************






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

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_hdiff.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to horizontal viscosity
!&E              terms for the meridional momentum equation
!&E              (3D equation)
!&E
!&E ** Description : explicit part
!&E                  d(nu*d(vz)/dx)/dx + d(nu*d(vz)/dy)/dy
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : yz_exp,me,mo,dy,dy_v,dx_f,vz,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---------------------------------------------------------------------

     yz_exp(:,i,j)=yz_exp(:,i,j)+ &
!-----------------------------------
!    d2v/dy2
!-----------------------------------
        (vish_xe(:,i,j+1)*(vz(:,i,j+1)-vz(:,i,j  ))/dy(i,j+1) &
        -vish_xe(:,i,j  )*(vz(:,i,j  )-vz(:,i,j-1))/dy(i,j  ))&
         /dy_v(i,j) &
!-----------------------------------
!    d2v/dx2
!-----------------------------------
       +(vish_phi(:,i  ,j)*me(i,j)*(vz(:,i+1,j)-vz(:,i  ,j))/dx_f(i  ,j)  &
        -vish_phi(:,i-1,j)*mo(i,j)*(vz(:,i  ,j)-vz(:,i-1,j))/dx_f(i-1,j)) &
        /dx_v(i,j)



!**************** end including dyn3dy_hdiff  *********************


! 5- integration des termes precedents sur
!    la verticale a destination du modele
!    barotrope
!****************** including dyn3dy_intz  *********************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_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 meridional equations
!&E
!&E ** Description :
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : yz_exp,hzuty,hty
!&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 dyn3dy_intz  **************

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

!*************include  dyn3dy_botstress   **********************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_botstress.h  ***
!&E
!&E ** Purpose : estimate the drag coefficient and the trend term
!&E              relative to the bottom stress
!&E              for the meridional momentum forcing (2D and 3D equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : uz_v,vz,hzuty,hty,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(uz_v(1)*uz_v(1)+vz(1,i,j)*vz(1,i,j))


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



   cdvit(i,j)=cdf*vitf





    tofyexp=cdvit(i,j)*vz(1,i,j)/hty(i,j)*cfroe
    intz_exp(i,j)=intz_exp(i,j)-tofyexp

!*************   end include dyn3dy_botstress   *******************


! W2 - coupling with WW3 (surface terms)


! 8- forcage meteo tension de surface et
!    pression a n-1/2


!**************** include dyn3dy_atmforc.h *********************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_atmforc.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the met. forcing
!&E              (atmospherical pressure gradient and wind stress)
!&E              for the meridional momentum forcing (2D and 3D equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : pa,tosy,hzuty,hty,dy_v
!&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,j+1)-pa(i,j))/(dy_v(i,j)*rhoref)
    intz_exp(i,j) = intz_exp(i,j)-(pa(i,j+1)-pa(i,j))/(dy_v(i,j)*rhoref)

! wind stress

    yz_exp(kmax,i,j) = yz_exp(kmax,i,j)+tosy(i,j)/MAX(hty(i,j),hminfrot)/rhoref/dsigu(kmax)

    intz_exp(i,j)  =  intz_exp(i,j)+tosy(i,j)/MAX(hty(i,j),hminfrot)/rhoref




!**************** end include dyn3dy_atmforc ********************



         IF (l_genpot) THEN             ! 9- potentiel generateur de maree a l instant n-1/2
!**************      include dyn3dy_genpot.h    ************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_genpot.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the tidal potential
!&E              for the meridional momentum forcing (2D and 3D equation)
!&E
!&E ** Description : explicit
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : potgene,dy_v
!&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,j+1)-potgene(i,j))/ &
                                 (dy_v(i,j)*rhoref)
     intz_exp(i,j)=intz_exp(i,j)+(potgene(i,j+1)-potgene(i,j))/ &
                                 (dy_v(i,j)*rhoref)

!************** end include dyn3dy_genpot  **********************
         ENDIF

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





       ENDIF                             ! fin test igdv

! 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 de
!     quantite de mouvement dyn3dy_mtm ) et des termes de divergence
!     du flux Nord-Sud (equation de conservation de la masse dyn3dy_cont)
!     et des vecteurs second-membre

!&E   ***************   include dyn3dy_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              de Saint-Venant (2D) dans la direction Nord-Sud
!&E              sur l ensemble de la grille.
!&E
!&E ** Description :
!&E
!&E ** Called by : dyn3dy_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       !  2017-12  (F. Ganthy) Bug corrected for obstruction-induced force
!&E
!&E-------------------------------------------------------------------------

     gdt=grav*dt

     j1=jb(i)
     j2=jh(i)
     IF (j >= j1 .AND. j <= j2-1) THEN

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

      IF (igdv(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 (igdv(i,j) == 1) THEN

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

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


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


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

! bottom stress : implicit part
                +ft*cfroi)



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


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

        u_vtmp=0.25_rsh*(u(i,j)+u(i,j+1)+u(i-1,j)+u(i-1,j+1))

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



! bottom stress (implicit part)
        vit=SQRT(u_vtmp*u_vtmp+v(i,j)*v(i,j))
        ft=0.5_rsh*dt*vit/hty(i,j)*(0.4_rsh/LOG(hty(i,j)/(exp1*z0b(i,j))+1.0_rsh))**2


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

! bottom stress : implicit part
                +ft*cfroi)



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


! NL advection : u*u*tan(phi)/R (spherical coordinates)
                  -halfdt*tfivsr(i,j)*u_vtmp**2          &
! non-linear (vdv/dy) : implicit part
                  -halfdt*v(i,j)/(dy(i,j)+dy(i,j+1))*(v(i,j+1)-v(i,j-1)) &


! viscosity (along y)
                  +halfdt*( vis2d(i,j+1)*(v(i,j+1)-v(i,j  ))/dy(i,j+1)   &
                           -vis2d(i,j  )*(v(i,j  )-v(i,j-1))/dy(i,j  ) )/dy_v(i,j) 



! viscosity and non-linear term
        IF ((htx(i-1,j)+htx(i-1,j+1)) >= 0.0_rsh) THEN
          dvsdx=(v(i,j)-v(i-1,j))/dx_f(i-1,j)
          y2dv(i,j)=y2dv(i,j)                            &


! non-linear (udv/dx) : explicit part
                    -0.25_rlg*dt*u_vtmp*dvsdx             &


! viscosity (along x) : explicit
                    -halfdt*vis2d(i-1,j)*dvsdx/dx_v(i,j)
        END IF
        IF ((htx(i,j)+htx(i,j+1)) >= 0.0_rsh) THEN
          dvsdx=(v(i+1,j)-v(i,j))/dx_f(i,j)
          y2dv(i,j)=y2dv(i,j)                            &


! non-linear (udv/dx) : explicit part
                    -0.25_rlg*dt*u_vtmp*dvsdx             &


! viscosity (along x) : explicit
                    +halfdt*vis2d(i,j)*dvsdx/dx_v(i,j)
        END IF

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

! wind stress in y-direction : explicit
                  +halfdt*tosy(i,j)/MAX(hty(i,j),hminfrot)/rhoref &

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



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

      END IF
     END IF

!&E   ***************   end include dyn3dy_mtm   ******************************
!&E   ***************   include dyn3dy_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 Nord-Sud.
!&E
!&E
!&E ** Called by : dyn2dy_tls_coef
!&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-03
!&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-------------------------------------------------------------------------

       j1=jb(i)
       j2=jh(i)
       IF (j >= j1+1.AND.j <= j2-1) THEN

! 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 meridienne
!    (flux au sud  du point xe(i,j) considere). Remplissage
!    du tableau hey qui ressert dans les flux de traceurs
!------------------------------------------------------------


      axe(i,j)=-halfdt*dx_v(i,j-1)/surf(i,j)*hey(i,j-1)

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

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

       cxe(i,j)=halfdt*dx_v(i,j)/surf(i,j)*hey(i,j)


! 4- terme du vecteur second membre du systeme lineaire
!    issu de la divergence de flux zonale
!    (flux a l est et a l ouest du point xe(i,j) considere)
!    Prise en compte de tous les apports d eau dans
!    la colonne d eau. Remplissage
!    du tableau hex qui ressert dans les flux de traceurs
!------------------------------------------------------------

       y2dxe(i,j)=xec(i,j)*fwetp(i,j)-halfdt/surf(i,j)*(              &
                              ( hex(i  ,j)*ustar(i  ,j)*dy_u(i  ,j)   &
                              - hex(i-1,j)*ustar(i-1,j)*dy_u(i-1,j) ) &
                              - SUM(phieau(:,i,j))                    &

                                )

       END IF

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

     END DO
   END DO
!$OMP end do



   !PRINT*, 'exit DYN3DY_TLS_COEF routine'

  END SUBROUTINE dyn3dy_tls_coef
