  SUBROUTINE dyn3dy_tls_coef(vz,boz,uz,u,v,xec,  &
#ifdef key_lagoon
                xeuv,                           &
#endif
                   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---------------------------------------------------------------------

#include "toolcpp.h"
# define ipx !ipx
# define ipy !ipy

   !! * 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
#ifdef key_linear_bottom_stress
   USE comvars2d,  ONLY : frot_lin
#endif
   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
#ifdef key_lagoon
   USE comvars2d,  ONLY : nbgrau, igrau,jgrau
#endif
#if key_dyn_pg_djcs
   USE comcolig3d, ONLY : zr,pr
#endif
   USE comvars3d,  ONLY : sig,dsigu,sigw,wz,dsigw,rhoref       ! dyn3dy_ipg_pj.h(sig,dsigw)
#ifdef key_profil_back
   USE comvars3d,  ONLY : boz2
#endif
#ifdef key_wave_soulsby_bottom_stress
   USE wave,       ONLY : uorb,fwst,wave_uzk1,wave_vzk1
#endif
#ifdef key_wave_qeuler_dynamics
   USE comvars2d,  ONLY : dy_f                     ! dyn3dy_wave.h
   USE wave,       ONLY : usz_zint,wave_usz,wave_wsz,vsz_zint,dsigdy_v,   &
                          dsigdx_u,duzdsig,dvzdsig,wave_botdiss_v,hrms
   USE wave,       ONLY : wave_bhd,wave_tauwy,wave_tauoy
#endif
#ifdef key_wave_qeuler_diag
   USE wavediag,   ONLY : diag_wf_vortexy,diag_wf_smb1y
   USE wavediag,   ONLY : diag_wf_advy,diag_wf_cory
   USE wavediag,   ONLY : diag_wf_grdpy,diag_wf_difhy
   USE wavediag,   ONLY : diag_wf_difvy
#endif
#if defined key_siggen || defined key_gencoord
   USE comsiggen,  ONLY : hzuty,hzwty,hzex,hzey,hc_sig,csu_sig
#endif
   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)
#ifdef key_MPI_2D
   USE toolmpi,    ONLY :  ex_j_rsh
#endif
#ifdef key_diagdyn
   USE diagdyn,   ONLY : mask_diagy,diagy_cor,diagy_corcor,diagy_ipg,diagy_nlt,diagy_hdif,diagy_vdif,diagy_epg
#ifdef key_diagdyn_details
   USE diagdyn,   ONLY : diagy_nlt_udvdx,diagy_nlt_vdvdy,diagy_nlt_wdvdz,diagy_hdif_dnudvdx2,diagy_hdif_dnudvdy2
#endif
#endif
#ifdef key_obstructions
   USE obstructions,      ONLY : obst_fvz_e,obst_fvz_i,obst_fv_e,obst_fv_i
#endif
#ifdef key_adapt_imp_vert_adv
   USE comvars3d, ONLY : wmo_exp_y, theta_wz
   USE toolmath,  ONLY : f_alpha
#ifdef key_adapt_imp_vert_adv_suivi
   USE comvars3d, ONLY : adapt_imp_vz
#endif
#endif

   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
#ifdef key_lagoon
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)     ,INTENT(in) :: xeuv
#endif

   !! * Local declarations
   INTEGER                          :: i,j,k,j1,j2
   REAL(KIND=rsh),DIMENSION(kmax)   :: uz_v                     !dyn3dy_corio.h
#ifdef key_lagoon
   INTEGER                          :: i0,j0,ng
#endif
#ifdef key_dyn_pg_djcs
   REAL(KIND=rsh), PARAMETER        :: demi   = 0.5_rsh         !dyn3dy_ipg_cs.h
   REAL(KIND=rsh), PARAMETER        :: uncinq = 0.2_rsh
   REAL(KIND=rsh), PARAMETER        :: undouz = 1.0_rsh/12.0_rsh
   REAL(KIND=rsh)                   :: cff,tbc,dyzr,dyboz, &
                                       dyzrp1,dybozp1,     &
                                       dyzrp2,dybozp2
   REAL(KIND=rsh),DIMENSION(0:kmax) :: dzr,dboz
#elif key_dyn_pg_pint
   REAL(KIND=rsh),DIMENSION(kmax)   :: p_primu_j,p_primu_jp1,sec_gp
   REAL(KIND=rsh),DIMENSION(0:kmax) :: p_prim_interface_j,  &
                                       p_prim_interface_jp1
   REAL(KIND=rsh)                   :: hz_j_kp1,hz_j,hz_jp1_kp1,hz_jp1, &
                                       zk_j,zk_jp1,htxes,htxen,dyb
#elif key_dyn_pg_dj
   REAL(KIND=rsh)    :: dyb,jacobian,dybs,dszy,dsby,dyzs,z1,z2,z3,z4,dzb,b2surf,b4surf
   REAL(KIND=rsh)    :: f31,f43,f42,f21
#else
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rlg)    :: bms,bm,pri,tbc,dyb,dhyy,dhcy,dxey,bhzn,bhzs       !dyn3dy_ipg_pj_sig.h
#else
   REAL(KIND=rsh)    :: bms,bm,pri,tbc,dyb,dhyy,dhty,htxen,htxes !dyn3dy_ipg_pj.h
#endif
#endif
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
#if defined key_dyn_adv_quick || defined key_dyn_adv_quickest
   REAL(KIND=rsh)    :: wmo,wm,wp,wdvdz,vp,um,vc,vdvdy,up,vm,udvdx,uv !dyn3dy_adv_quickest.h
#else
   REAL(KIND=rsh)    :: wdvi,wdvs                                     !dyn3dy_adv_centered.h
#endif
#endif
#endif
   REAL(KIND=rsh)    :: vitf,cdf,tofyexp                              !dyn3dy_botstress.h
#ifdef key_wave_soulsby_bottom_stress
   REAL(KIND=rsh)    :: uorb_v,shw,shc,shwc,fwst_v,uzk1_v             !dyn3dx_botstress.h
#endif
#ifdef key_wave_qeuler_dynamics
   REAL(KIND=rsh)    :: USX_v,dxvz_v,dyuz_v,USZ_v,dsiguz_v,dsigx_v    !dyn3dy_wave.h
#endif
   REAL(KIND=rsh)    :: ft,gdt,vit,factft,dvsdx,u_vtmp                !dyn3dy_mtm.h
#ifdef key_diagdyn  
   REAL(KIND=rsh)    :: zp,effet_pression ! zp :save previous value of yz_exp(kmax) before adding wind stress
   REAL(kind=rsh),DIMENSION(1:kmax)  :: z_prev
#endif
   REAL(KIND=rsh)    :: immer                  ! dyn3dy_wave.h
#ifdef key_adapt_imp_vert_adv
   REAL(kind=rsh)    :: inv_f_alpha
#endif

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

   PRINT_DBG*, 'enter DYN3DY_TLS_COEF routine'

#if key_dyn_pg_djcs
   !************************
   !     Preliminary step
   !************************

#ifdef key_MPI_2D
!$OMP DO schedule(runtime) private(i,j,k,cff,dboz,dzr)
#else
!$OMP DO schedule(dynamic,5) private(i,j,k,cff,dboz,dzr)
#endif
   DO j = ljmin, ljmax
     DO i = liminp2,limaxm1
       IF (j >= jb(i) .AND. j<=jh(i) ) THEN
!         IF (h0(i,j) > -valmanq) THEN
! if dzr too small over wet-drying cells
         IF (h0(i,j)+xec(i,j) > h0fond) THEN

           zr(kmax,i,j) = f_zu(h0(i,j),xec(i,j),kmax,i,j)
           DO k = kmax-1, 1, -1
             zr(k,i,j) = f_zu(h0(i,j),xec(i,j),k,i,j)
             dzr(k) = zr(k+1,i,j)-zr(k,i,j)
#ifdef key_profil_back
             dboz(k) = boz2(k+1,i,j)-boz2(k,i,j)
#else
             dboz(k) = boz(k+1,i,j)-boz(k,i,j)
#endif
           END DO
           dzr(kmax) = dzr(kmax-1)
           dboz(kmax) = dboz(kmax-1)
           dzr(0) = dzr(1)
           dboz(0) = dboz(1)

           DO k = kmax,1,-1
             cff = 2.0_rsh*dboz(k)*dboz(k-1)
             IF (cff > epsilon) THEN
               dboz(k) = cff / (dboz(k)+dboz(k-1))
             ELSE
               dboz(k) = 0.0_rsh
             ENDIF
             dzr(k) = 2.0_rsh*dzr(k)*dzr(k-1) / (dzr(k)+dzr(k-1))
           END DO

! pressure vertically interpolated
#ifdef key_profil_back
           pr(kmax,i,j) = (  boz2(kmax,i,j) + &
                             0.5_rsh * (xec(i,j)-zr(kmax,i,j)) * (boz2(kmax,i,j)-boz2(kmax-1,i,j)) / (zr(kmax,i,j)-zr(kmax-1,i,j)) ) &
                                     * (xec(i,j)-zr(kmax,i,j))
#else
           pr(kmax,i,j) = (  boz(kmax,i,j) + &
                             0.5_rsh * (xec(i,j)-zr(kmax,i,j)) * (boz(kmax,i,j)-boz(kmax-1,i,j)) / (zr(kmax,i,j)-zr(kmax-1,i,j)) ) &
                                     * (xec(i,j)-zr(kmax,i,j))
#endif
           DO k = kmax-1, 1, -1
#ifdef key_profil_back
             pr(k,i,j) = pr(k+1,i,j) +                                               &
                         demi*(  (boz2(k+1,i,j)+boz2(k,i,j))*(zr(k+1,i,j)-zr(k,i,j)) &
                         - uncinq*( (dboz(k+1)-dboz(k))*(zr(k+1,i,j)-zr(k,i,j)   &
                                             -undouz*(dzr(k+1)+dzr(k)))          &
                                   -(dzr(k+1)-dzr(k))*(boz2(k+1,i,j)-boz2(k,i,j) &
                                             -undouz*(dboz(k+1)+dboz(k))) )  )
#else
             pr(k,i,j) = pr(k+1,i,j) +                                               &
                         demi*(  (boz(k+1,i,j)+boz(k,i,j))*(zr(k+1,i,j)-zr(k,i,j))   &
                         - uncinq*( (dboz(k+1)-dboz(k))*(zr(k+1,i,j)-zr(k,i,j)   &
                                             -undouz*(dzr(k+1)+dzr(k)))          &
                                   -(dzr(k+1)-dzr(k))*(boz(k+1,i,j)-boz(k,i,j)   &
                                             -undouz*(dboz(k+1)+dboz(k))) )  )
#endif
           END DO

         END IF
       END IF
     END DO
   END DO
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,zr(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh( 0,1,kmax,limin  ,limaxp1,ljmin  ,ljmaxp1,pr(:,limin  :limaxp1,ljmin  :ljmaxp1))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(zr,pr)

#endif

!$OMP DO schedule(dynamic,5)                                    &
!$OMP& private(i,j,k,uz_v                                       &
#ifdef key_dyn_pg_djcs
!$OMP& ,cff,tbc,dyzr,dyboz,dyzrp1,dybozp1,dyzrp2,dybozp2        &
#elif key_dyn_pg_pint
!$OMP& ,p_primu_j,p_primu_jp1,p_prim_interface_j                &
!$OMP& ,p_prim_interface_jp1,sec_gp,zk_j,zk_jp1,htxes,htxen,dyb &
!$OMP& ,hz_j_kp1,hz_j,hz_jp1_kp1,hz_jp1                         &
#elif key_dyn_pg_dj
!$OMP& ,dyb,jacobian,dybs,dszy,dsby,dyzs,z1,z2,z3,z4            &
!$OMP& ,dzb,b2surf,b4surf,f31,f43,f42,f21                       &
#else
#if defined key_siggen || defined key_gencoord
!$OMP& ,bms,bm,pri,tbc,dyb,dhyy,dhcy,dxey,bhzn,bhzs             &
#else
!$OMP& ,bms,bm,pri,tbc,dyb,dhyy,dhty,htxen,htxes                &
#endif
#endif
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
#if defined key_dyn_adv_quick || defined key_dyn_adv_quickest
!$OMP& ,vm,vp,vdvdy,up,udvdx,wm,wp,wdvdz,wmo,um,vc,uv           &
# else
!$OMP& ,wdvi,wdvs                                               &
#endif
#endif
#endif
#if defined key_diagdyn
!$OMP& ,z_prev                                                  &
#endif
!$OMP& ,vitf,cdf,tofyexp                                        &
!$OMP& ,ft,gdt,vit,factft,dvsdx,u_vtmp                          &
!$OMP& ,j1,j2                                                   &
#ifdef key_wave_qeuler_dynamics
!$OMP& ,immer,dyuz_v,dxvz_v,USX_v,USZ_v,dsiguz_v,dsigx_v,uzk1_v &
#endif
#ifdef key_adapt_imp_vert_adv
!$OMP& ,inv_f_alpha                                             &
#ifdef key_adapt_imp_vert_adv_suivi
!$OMP& ,adapt_imp_vz                                            &
#endif
#endif
!$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
#        include "dyn3dy_corio.h"

#ifdef key_diagdyn
         IF (mask_diagy(i,j)) diagy_cor(:,i,j)=fv(i,j)*uz_v(:)
#endif

#        include "dyn3dy_intz.h"

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

#ifdef key_diagdyn
         ! premiere correction pour Coriolis
         IF (mask_diagy(i,j)) diagy_corcor(:,i,j)=yz_exp(:,i,j)
         IF (icon<=1 .AND. mask_diagy(i,j)) z_prev(:)=yz_exp(:,i,j)
#endif

                                         ! W1- coupling with WW3 (terms depending to the depth)
#ifdef key_wave_qeuler_dynamics
#     include "dyn3dy_wave_3d.h"
#endif

         IF(icon==2) THEN                ! 2- gradient de pression interne a n-1/2
 
#ifdef key_diagdyn
           IF (mask_diagy(i,j)) z_prev(:)=yz_exp(:,i,j)
#endif	

#ifdef key_dyn_pg_djcs
#          include "dyn3dy_ipg_cs.h"
#elif key_dyn_pg_pint
#          include "dyn3dy_ipg_pint.h"
#elif key_dyn_pg_dj
#          include "dyn3dy_ipg_dj.h"
#else
#if defined key_siggen || defined key_gencoord
#          include "dyn3dy_ipg_pj_sig.h"
#else
#          include "dyn3dy_ipg_pj.h"
#endif
#endif

#ifdef key_diagdyn
           IF (mask_diagy(i,j)) THEN
             diagy_ipg(:,i,j)=yz_exp(:,i,j)-z_prev(:)
             z_prev(:)=yz_exp(:,i,j)
           END IF
#endif

         END IF
                                     ! 3- termes d advection a l instant n-1
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
#ifdef key_dyn_adv_quick
#        include "dyn3dy_adv_quick.h"
#elif key_dyn_adv_quickest
#        include "dyn3dy_adv_quickest.h"
#else
#        include "dyn3dy_adv_centered.h"
#endif
#endif
#endif

#if defined key_diagdyn && ! defined key_diagdyn_details
         IF (mask_diagy(i,j)) THEN
           diagy_nlt(:,i,j)=-(yz_exp(:,i,j)-z_prev(:))
           z_prev(:)=yz_exp(:,i,j)
         END IF
#endif

                                     ! 4- terme de diffusion horizontale a n-1
#        include "dyn3dy_hdiff.h" 

#if defined key_diagdyn && ! defined key_diagdyn_details
         IF (mask_diagy(i,j)) diagy_hdif(:,i,j)=yz_exp(:,i,j)-z_prev(:)
#endif	
                                     ! 5- integration des termes precedents sur
                                     !    la verticale a destination du modele
                                     !    barotrope
#        include "dyn3dy_intz.h"

                                     ! 6- tension sur le fond a l instant n-1
#if ! defined key_no_bottom_stress
#        include "dyn3dy_botstress.h"
#endif

                                         ! W2 - coupling with WW3 (surface terms)
#ifdef key_wave_qeuler_dynamics
#     include "dyn3dy_wave_surf.h"
#endif

                                     ! 8- forcage meteo tension de surface et
                                     !    pression a n-1/2
#ifdef key_diagdyn
         IF (mask_diagy(i,j)) THEN
           zp=yz_exp(kmax,i,j)  ! on sauve en surface pour recuperer ensuite  le vent
           z_prev(:)=yz_exp(:,i,j) ! sur la verticale on va ensuite ajouter la pression
         END IF
#endif

#        include "dyn3dy_atmforc.h"

#ifdef key_diagdyn
         IF (mask_diagy(i,j)) THEN
           effet_pression=yz_exp(1,i,j)-z_prev(1)
           diagy_vdif(kmax,i,j)=yz_exp(kmax,i,j)-zp-effet_pression
           diagy_epg(:,i,j)=effet_pression ! on aura pu creer un tableau en plus...
         ENDIF
#endif

         IF (l_genpot) THEN             ! 9- potentiel generateur de maree a l instant n-1/2
#          include "dyn3dy_genpot.h" 
         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))

#ifdef key_obstructions
         ! obstruction force : 3D explicit Part 
         yz_exp(:,i,j) = yz_exp(:,i,j) - obst_fvz_e(:,i,j)
#endif

#ifdef key_diagdyn
         ! seconde correction pour Coriolis
         IF (mask_diagy(i,j)) diagy_corcor(:,i,j)= &
            -(diagy_corcor(:,i,j)-fv(i,j)*0.25_rlg*(ustar(i,j)+ustar(i,j+1)+ustar(i-1,j)+ustar(i-1,j+1)))
#endif	  

       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

#      include "dyn3dy_mtm.h"
#      include "dyn3dy_cont.h"

     END DO
   END DO
!$OMP end do

#ifdef key_lagoon
!$OMP SINGLE
   ! prise en compte des Graus : limite ouverte est situee a l interieur du domaine
   DO ng=1,nbgrau
     i0=igrau(ng)
     j0=jgrau(ng)
     IF_MPI ( (i0 >= limin) .AND. (i0 <= limaxp1) .AND.(j0 >= ljmin) .AND. (j0 <= ljmaxp1) ) THEN
       axe(i0,j0)=0.0_rsh
       cxe(i0,j0)=0.0_rsh
       y2dxe(i0,j0)=xeuv(i0,j0)*fwetp(i0,j0)
     ENDIF_MPI
   END DO
!$OMP END SINGLE
#endif

   PRINT_DBG*, 'exit DYN3DY_TLS_COEF routine'

  END SUBROUTINE dyn3dy_tls_coef
