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

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

   !! * 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
#ifdef key_lagoon
   USE comvars2d,  ONLY : nbgrau, igrau,jgrau
#endif

   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
#if key_dyn_pg_djcs
   USE comcolig3d, ONLY : zr,pr
#endif
   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
#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 : dx_f                     ! dyn3dx_wave.h
   USE wave,       ONLY : wave_vsz,wave_wsz,duzdsig, &
                          dvzdsig,dsigdx_u,dsigdy_v,wave_botdiss_u, &
                          vsz_zint,usz_zint,hrms
   USE wave,       ONLY : wave_bhd,wave_tauwx,wave_tauox
#endif
#ifdef key_wave_qeuler_diag
   USE wavediag,   ONLY : diag_wf_vortexx,diag_wf_smb1x
   USE wavediag,   ONLY : diag_wf_advx,diag_wf_corx
   USE wavediag,   ONLY : diag_wf_grdpx,diag_wf_difhx
   USE wavediag,   ONLY : diag_wf_difvx
#endif
#if defined key_siggen || defined key_gencoord
   USE comsiggen,  ONLY : hzutx,hzwtx,hzex,hzey,hc_sig,csu_sig
#endif
   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
#ifdef key_MPI_2D
   USE toolmpi,    ONLY :  ex_i_rsh
#endif
#ifdef key_diagdyn
   USE diagdyn,   ONLY : mask_diagx,diagx_cor,diagx_corcor,diagx_ipg,diagx_nlt,diagx_hdif,diagx_vdif,diagx_epg,diag_boz
#ifdef key_diagdyn_details
   USE diagdyn,   ONLY : diagx_nlt_ududx,diagx_nlt_vdudy,diagx_nlt_wdudz,diagx_hdif_dnududx2,diagx_hdif_dnududy2
#endif
#endif
#ifdef key_obstructions
   USE obstructions,      ONLY : obst_fuz_e,obst_fuz_i,obst_fu_e,obst_fu_i
#endif
#ifdef key_casobstflume_cycle
   USE comvars2d, ONLY : uforcing
#endif
#if defined key_1DV_force_ubarotrope || defined key_1DV_force_ufile
   USE cometdyn,  ONLY : ubaro_force_1DV
#elif key_1DV_force_gradssh
   USE cometdyn,  ONLY : gradxe_force_1DV
#endif 
#ifdef key_adapt_imp_vert_adv
   USE comvars3d, ONLY : wmo_exp_x, theta_wz
   USE toolmath,  ONLY : f_alpha
#ifdef key_adapt_imp_vert_adv_suivi
   USE comvars3d, ONLY : adapt_imp_uz
#endif
#endif

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

   !! * Local declarations
   INTEGER                          :: i,j,k
   REAL(KIND=rsh),DIMENSION(kmax)   :: vz_u                     !dyn3dx_corio.h
#ifdef key_lagoon
   INTEGER                          :: i0,j0,ng
#endif
#ifdef key_dyn_pg_djcs
   REAL(KIND=rsh), PARAMETER        :: demi   = 0.5_rsh         !dyn3dx_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,dxzr,dxboz, &
                                       dxzrp1,dxbozp1,     &
                                       dxzrp2,dxbozp2
   REAL(KIND=rsh),DIMENSION(0:kmax) :: dzr,dboz
#elif key_dyn_pg_pint
   REAL(KIND=rsh),DIMENSION(kmax)   :: p_primu_i,p_primu_ip1,sec_gp
   REAL(KIND=rsh),DIMENSION(0:kmax) :: p_prim_interface_i,  &
                                       p_prim_interface_ip1
   REAL(KIND=rsh)                   :: hz_i_kp1,hz_i,hz_ip1_kp1,hz_ip1, &
                                       zk_i,zk_ip1,htxee,htxeo,dxb
#elif key_dyn_pg_dj
   REAL(KIND=rsh)    :: dxb,jacobian,dxbs,dszx,dsbx,dxzs,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,dxb,dhxx,dhcx,dxex,bhze,bhzo       !dyn3dx_ipg_pj_sig.h
#else
   REAL(KIND=rsh)    :: bms,bm,pri,tbc,dxb,dhxx,dhtx,htxee,htxeo   !dyn3dx_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,wdudz,uc,up,um,ududx,vp,vm,vdudy,uv !dyn3dx_adv_quickest.h, adv_uq.h
#else
   REAL(KIND=rsh)    :: wdus,wdui                                  !dyn3dx_adv_centered.h, adv_u.h
#endif
#endif
#endif
   REAL(KIND=rsh)    :: vitf,cdf,tofxexp                           !dyn3dx_botstress.h
#ifdef key_wave_soulsby_bottom_stress
   REAL(KIND=rsh)    :: uorb_u,shw,shc,shwc,fwst_u,vzk1_u          !dyn3dx_botstress.h
#endif
#ifdef key_wave_qeuler_dynamics
   REAL(KIND=rsh)    :: USY_u,dyuz_u,dxvz_u,USZ_u,dsigvz_u,dsigy_u !dyn3dx_wave.h
#endif
   REAL(KIND=rsh)    :: ft,gdt,vit,factft,dusdy,v_utmp             !dyn3dx_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                  ! dyn3dx_wave.h
#ifdef key_adapt_imp_vert_adv
   REAL(kind=rsh)    :: inv_f_alpha 
#endif

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

   PRINT_DBG*, 'enter DYN3DX_TLS_COEF routine'

#if key_dyn_pg_djcs

   ! pressure : vertical interpolation (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 = ljminp2, ljmaxm1
     DO i = MAX0(limin,ig(j)), MIN0(limax,id(j))
!       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
#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 DO
   END DO
!$OMP END DO

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

#endif

#ifdef key_MPI_2D
!$OMP DO schedule(runtime)                                     &
#else
!$OMP DO schedule(dynamic,5)                                   &
#endif
!$OMP& private(i,j,k,vz_u                                      &
#ifdef key_dyn_pg_djcs
!$OMP& ,cff,tbc,dxzr,dxboz,dxzrp1,dxbozp1,dxzrp2,dxbozp2       &
#elif key_dyn_pg_pint
!$OMP& ,p_primu_i,p_primu_ip1,p_prim_interface_i               &
!$OMP& ,p_prim_interface_ip1,sec_gp,zk_i,zk_ip1,htxee,htxeo    &
!$OMP& ,hz_i_kp1,hz_i,hz_ip1_kp1,hz_ip1,dxb                    &
#elif key_dyn_pg_dj
!$OMP& ,dxb,jacobian,dxbs,dszx,dsbx,dxzs,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,dxb,dhxx,dhcx,dxex,bhze,bhzo            &
#else
!$OMP& ,bms,bm,pri,tbc,dxb,dhxx,dhtx,htxee,htxeo               &
#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& ,wmo,wm,wp,wdudz,uc,up,um,ududx,vp,vm,vdudy,uv          &
#else
!$OMP& ,wdus,wdui                                              &
#endif
#endif
#endif
#if defined key_diagdyn
!$OMP& ,z_prev                                                 &
#endif
!$OMP& ,vitf,cdf,tofxexp                                       &
!$OMP& ,gdt,v_utmp,vit,factft,ft,dusdy                         &
#ifdef key_wave_qeuler_dynamics
!$OMP& ,immer,dxvz_u,dyuz_u,USY_u,USZ_u,dsigvz_u,dsigy_u,vzk1_u  &
#endif
#ifdef key_adapt_imp_vert_adv
!$OMP& ,inv_f_alpha                                            &
#ifdef key_adapt_imp_vert_adv_suivi
!$OMP& ,adapt_imp_uz                                           &
#endif
#endif
!$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
#        include "dyn3dx_corio.h"

#ifdef key_diagdyn
         IF (mask_diagx(i,j)) THEN
           diagx_cor(:,i,j)=-fu(i,j)*vz_u(:)
           diag_boz(:,i,j)=boz(:,i,j)
         END IF
#endif

#        include "dyn3dx_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_diagx(i,j)) diagx_corcor(:,i,j)=yz_exp(:,i,j)
         IF (icon<=1 .AND. mask_diagx(i,j)) z_prev(:)=yz_exp(:,i,j)
#endif

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

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

#ifdef key_diagdyn
           IF (mask_diagx(i,j)) z_prev(:)=yz_exp(:,i,j)
#endif	 
 
#ifdef key_dyn_pg_djcs
#          include "dyn3dx_ipg_cs.h"
#elif key_dyn_pg_pint
#          include "dyn3dx_ipg_pint.h"
#elif key_dyn_pg_dj
#          include "dyn3dx_ipg_dj.h"
#else
#if defined key_siggen || defined key_gencoord
#          include "dyn3dx_ipg_pj_sig.h"
#else
#          include "dyn3dx_ipg_pj.h"
#endif
#endif

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

         ENDIF

                                      ! 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 "dyn3dx_adv_quick.h"
#elif defined key_dyn_adv_quickest
#        include "dyn3dx_adv_quickest.h"
#else
#        include "dyn3dx_adv_centered.h"
#endif
#endif
#endif

#if defined key_diagdyn && ! defined key_diagdyn_details
         IF (mask_diagx(i,j)) THEN
           diagx_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 "dyn3dx_hdiff.h"

#if defined key_diagdyn && ! defined key_diagdyn_details
         IF (mask_diagx(i,j)) diagx_hdif(:,i,j)=yz_exp(:,i,j)-z_prev(:)
#endif	 

                                      ! 5- integration des termes precedents 1 a 4
                                      !    sur la verticale a destination du mode
                                      !    barotrope.
#        include "dyn3dx_intz.h" 

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

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

                                      ! 8- forcing meteo tension de surface et
                                      !    pression a l instant n-1/2
#ifdef key_diagdyn
         IF (mask_diagx(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 "dyn3dx_atmforc.h"

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

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

#ifdef key_casobstflume_cycle
         yz_exp(:,i,j) = yz_exp(:,i,j) + 0.0025_rlg*uforcing
#endif
#ifdef key_obstructions
         ! obstruction force : 3D explicit Part
         yz_exp(:,i,j) = yz_exp(:,i,j) - obst_fuz_e(:,i,j)
#endif

#ifdef key_diagdyn
         ! seconde correction pour Coriolis
         IF (mask_diagx(i,j)) diagx_corcor(:,i,j)= &
                         -(diagx_corcor(:,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	  

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

#ifdef key_casobstflume_cycle
         intz_exp(i,j) = intz_exp(i,j) + 0.0025_rlg*uforcing
#endif
       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

#      include "dyn3dx_mtm.h"
#      include "dyn3dx_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 DYN3DX_TLS_COEF routine'

  END SUBROUTINE dyn3dx_tls_coef
