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

  SUBROUTINE dyn3dxy(u,v,xev,uv,uv2,xeuv,uz,vz,boz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dyn3dxy  ***
   !&E
   !&E ** Purpose : resolution of Saint-Venant equations (barotropic flow)
   !&E              and of 3D momentum equations (barotropic + baroclinic)
   !&E
   !&E ** Description : estimate of uz,vz and xeuv. xeuv and vz are the last
   !&E                  estimates and balance each other.
   !&E
   !&E ** Called by : step3dxy.F90
   !&E
   !&E ** External calls : turb_smago,dyn3d_height_mask,
   !&E                     obc_apply_u,obc_apply_uz,obc_apply_v,obc_apply_vz
   !&E                     dyn2dx_tls_obc,dyn2dx_tls_solve,dyn2dx_cyc_solve
   !&E                     dyn2dy_tls_obc,dyn2dy_tls_solve,dyn2dy_cyc_solve
   !&E                     dyn3dx_tls_coef,dyn3dx_ztls_solve,dyn3dx_epg
   !&E                     dyn3dy_tls_coef,dyn3dy_ztls_solve,dyn3dy_epg
   !&E
   !&E ** Used ij-arrays : igdu,igdv,igdv_old,xev,u,v,boz
   !&E
   !&E ** Modified variables : uv,xeuv,uv2,uz,vz,vstar,igdu_old,xeint
   !&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       !  2007-09 (V. Garnier) homogeneneity of dx dy...
   !&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
   !&E       !  2008-07  (T. Odaka) correction MPI
   !&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       !  2009-04  (A. Mano V. Garnier) MPI2D
   !&E       !  2009-04  (F. Dumas) generalized sigma
   !&E       !  2009-05  (F. Dumas)  AGRIF to new time scheme
   !&E       !  2009-05  (F. Vandermeirsch) characteristics method
   !&E       !  2009-07  (F. Vandermeirsch)  MC 2nd order scheme
   !&E       !  2009-07  (V. Garnier) MPI2d for characteristics method
   !&E       !  2010-04  (T. Odaka) MPI2D cyclic
   !&E       !  2010-07  (L. Debreu, F. Dumas) rewritting of ADI
   !&E       !  2010-09  (L. Debreu, F. Dumas) FCT to avoid negative water depths if wet-drying
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2011-12  (S. Petton, F. Dumas) Implicitation of the bottom stress (consistence between 2D en 3D)
   !&E       !  2012-02  (F. Vandermeirsch) Complement of previous evolution for key_no_bottom_stress
   !&E       !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Review of the implementation of wave-current quasi-eulerian formalism
   !&E       !  2016-12  (A. Fiandrino, B. Thouvenin) introduction key_lagoon
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : hex,hey,h0,surf,dx_v,dy_v,dy_u,htx,z0b
   USE comvars2d,    ONLY : halfdt
   USE comvars2d,    ONLY : l_obc_cycl_y,l_obc_cycl_x,l_smagor
   USE comvars2d,    ONLY : jb,jh,ig,id,fwet,fwetp
   USE comvars2d,    ONLY : frot_lin
   USE comfluxdo,    ONLY : fxeau,fyeau,phieau
   USE comvars3d,    ONLY : dsigu,sig
#if defined key_siggen || defined key_gencoord
   USE comsiggen,    ONLY : hzex,hzey,hzutx
#endif
#ifdef key_wave_qeuler_dynamics
   USE wave,         ONLY : wave_usz,wave_vsz
#if ! defined key_palm_ww3
   USE waveqeulerdyn,ONLY : wave_qeuler
#endif
#ifdef key_wave_breaking_visco
   USE waveqeulerdyn,ONLY : wave_visc
#endif
#endif
   USE obcapplyu
   USE obcapplyv
   USE obcapplyuz
   USE obcapplyvz
   USE obccombine,   ONLY : l_obc_char,l_obc_north,l_obc_south,l_obc_east,l_obc_west
   USE comcolig2d,   ONLY : igdu_old,igdv_old,igdu,igdv,av,bv,cv,y2dv
   USE comcolig3d,   ONLY : intz_exp,cdvit,intz_exp_ustar,intz_exp_vstar
   USE_MPI toolmpi,  ONLY : ex_j_rsh,ex_j_rlg
   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_i_rlg,ex_peri_x_rsh,ex_peri_y_rsh
   USE_MPI toolmpi,  ONLY : ex_peri_x_avcv,ex_peri_y_avcv
   USE_MPI comcolig_mpi, ONLY : transpose_row,trback_row
   USE_MPI comcolig_mpi, ONLY : transpose_col,trback_col
   USE_MPI comcolig_mpi, ONLY : cy2dxe,caxe,cbxe,ccxe,cxeuv
   USE_MPI comcolig_mpi, ONLY : ry2dxe,raxe,rbxe,rcxe,rxeuv
#ifdef key_lagoon
   USE comvars2d,     ONLY : igrau,jgrau,nbgrau,type_grau,dx_u   
   USE obccombine,    ONLY : obc_sshs,obc_sshn,obc_sshe,obc_sshw
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)     ,INTENT(in   ) :: xev
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)     ,INTENT(in   ) :: u,v
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)     ,INTENT(inout) :: xeuv
   REAL(KIND=rsh),DIMENSION(     liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(inout) :: uv,uv2
   REAL(KIND=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(inout) :: vz,uz
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in   ) :: boz

   !! * Local declarations
   INTEGER                                  :: i,j,i1,i2,j1,j2,test
   LOGICAL                                  :: Agrif_Root_OMP,Agrif_GrandMother_OMP,bidon
   REAL(KIND=rsh)                           :: umoy
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2) :: xeint,vstar
   REAL(KIND=rsh),DIMENSION(limin:limax)                     :: tfh,tfb
   REAL(KIND=rsh),DIMENSION(ljmin:ljmax)                     :: tfg,tfd
#ifdef key_lagoon
   INTEGER                                  :: i0,j0,ng
#endif
#ifdef key_dyn_wetdry_fct
   REAL(KIND=rlg)                           :: A
   REAL(KIND=rlg),DIMENSION(limin:limaxp1,ljmin:ljmaxp1) :: phimax
#endif
   REAL(kind=rsh) vz_u,vitf,cdf

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

   PRINT_DBG*, 'enter DYN3DXY routine'

#ifdef key_MPI_2D
   ! To make sure AGRIF conv will take into account these required variables
   INST_AGRIF bidon = allocated(ry2dxe)
   INST_AGRIF bidon = allocated(raxe)
   INST_AGRIF bidon = allocated(rbxe)
   INST_AGRIF bidon = allocated(rcxe)
   INST_AGRIF bidon = allocated(cy2dxe)
   INST_AGRIF bidon = allocated(caxe)
   INST_AGRIF bidon = allocated(cbxe)
   INST_AGRIF bidon = allocated(ccxe)
#endif

! because the variables (pointer?) used in Agrif_Root function are not seen by all the cpus in OMP
   Agrif_Root_OMP = agrif_root()
   INST_AGRIF Agrif_GrandMother_OMP = agrif_grandmother()

!$OMP PARALLEL default(none) shared(xev,igdu,igdv,dsigu     &
!$OMP&,u,v,uz,vz,xeuv,uv,uv2,boz,xeint,halfdt               &
!$OMP&,l_smagor,tfg,tfd,tfh,tfb,fxeau,fyeau,hex,hey,htx     &
!$OMP&,vstar,dy_v,igdu_old,igdv_old,intz_exp                &
!$OMP&,limin,liminm1,liminm2,liminp1,liminp2                &
!$OMP&,ljmin,ljminm1,ljminm2,ljminp1,ljminp2                &
!$OMP&,limax,limaxm1,limaxp2                                &
!$OMP&,ljmax,ljmaxm1,ljmaxp2                                &
#ifdef key_MPI_2D
#  ifdef key_dyn_wetdry_fct
!$OMP&,ljmaxp1,limaxp1                                      &
#  endif
!$OMP&,ry2dxe,raxe,rbxe,rcxe,rxeuv                          &
!$OMP&,cy2dxe,caxe,cbxe,ccxe,cxeuv                          &
#endif
!$OMP&,ig,id,jb,jh,av,bv,cv,y2dv                            &
!$OMP&,z0b,cdvit,intz_exp_ustar,intz_exp_vstar              &
#ifdef key_wave_qeuler_dynamics
!$OMP&,wave_usz,wave_vsz                                    &
#endif
#if defined key_siggen || defined key_gencoord
!$OMP&,hzex,hzey,hzutx                                      &
#else
!$OMP&,sig                                                  &
#endif
#ifdef key_dyn_wetdry_fct
!$OMP&,phimax,h0,surf,dx_v,dy_u,phieau,fwet,fwetp           &
#endif
#if key_linear_bottom_stress
!$OMP&,frot_lin                                             &
#endif
#ifdef key_lagoon
!$OMP&, igrau,jgrau,type_grau,nbgrau,obc_sshe,obc_sshw      &
!$OMP&, obc_sshs,obc_sshn,dx_u                              &
#endif
!$OMP&,l_obc_char,l_obc_cycl_x,l_obc_cycl_y,test            &
!$OMP&,l_obc_north,l_obc_south,l_obc_east,l_obc_west)       &
#ifdef key_agrif
!$OMP& SHARED(imin,imax,jmin,jmax)                          &
!$OMP& FIRSTPRIVATE(Agrif_GrandMother_OMP)                  &
#endif
!$OMP& FIRSTPRIVATE(Agrif_Root_OMP)                         &
!$OMP& PRIVATE(umoy,vz_u,vitf,cdf,                          &
#ifdef key_dyn_wetdry_fct
!$OMP& A,                                                   &
#endif
#ifdef key_lagoon
!$OMP& i0,j0, ng,                                           &
#endif
!$OMP& i,j,i1,i2,j1,j2)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                             Part I                               !
   ! Estimate of intermediate variables :                             !
   !                         XE^n+1/2*, u_i^n+1/2 and uz_i^n+1/2      !
   !                                                                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   !-----------------------------------------------------------------------
   ! Calculation of heights hex and hey
   ! used in continuity equation (and for tracers advection)
   ! Calculation of heights htx and hty, hzutx and hzutx, hzwtx and hzwtx
   ! used in momentum equation
   ! Definition of the type of cell (igdu and igdv),
   ! the type of extremity of the row (tfg,tfd) used in dyn2dxtlsobc
   !-----------------------------------------------------------------------

   test = 1

   CALL dyn3d_height_mask(xev,tfg,tfd,tfb,tfh,test)

   IF_AGRIF ( .NOT. Agrif_GrandMother_OMP) THEN

   ! Estimate of viscosity coefficients
   ! (smagorinsky formulation)
   !-----------------------------------

!   IF ( l_smagor .AND. mod(nbouc,5) == 0 ) &  !!! ATTENTION MPI si utilise
   IF ( l_smagor ) THEN

     CALL turb_smago(uz,vz,igdu)

   ELSE

#if defined  key_wave_qeuler_dynamics && defined key_wave_breaking_visco
     CALL wave_visc
#endif

   END IF

   ! Estimate the terms (or intermediate fields) relative to the wave-current interactions
   ! quasi-eulerian formalism
   ! -------------------------------------------------------------------------------------

#if defined key_wave_qeuler_dynamics && ! defined key_palm_ww3
   CALL wave_qeuler
#endif


   ! Estimate of vstar :
   ! Add the external pressure gradient to Vstar to complete Vstar estimate.
   ! Finaly, Vstar contains the verticaly integrated right hand side
   ! explicitly evaluated in stepvu and involving the use
   ! of uz^n-1/2, vz^n-1/2 and v_i^n :
   ! v*^n+1/2 = v^n + RHS(uz^n-1/2,vz^n-1/2) - gdt/2*d(xe^n)/dx
   !---------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=limin,limax
! strange, the following lines should be sufficient
   !DO j=ljminp1,ljmaxm1
     !DO i = liminp2, limaxm1
       IF (igdv(i,j)==1 .AND. igdv_old(i,j)==1) THEN
!         vstar(i,j)=vstar(i,j)-grav*halfdt*(xeuv(i,j+1)-xeuv(i,j))/dy_v(i,j) mpext
         vstar(i,j) = v(i,j)+halfdt*(intz_exp_vstar(i,j)-grav*(xev(i,j+1)-xev(i,j))/dy_v(i,j) mpext)
       ELSE
         vstar(i,j) = v(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

   ! Apply open boundary conditions to vstar
   !----------------------------------------
OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_char) THEN
! Attention pas de ustar alors on utilise (u)n-.5
     CALL obc_apply_v_mc(u,u,v,vstar,xev,xeuv,.FALSE.)
   ELSE
     CALL obc_apply_v(vstar)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vstar)

   ! Construction of 2D and 3D matrixes for calculation along x
   ! 1st part : estimate of right hand side :
   ! RHS(uz^n,vz_i^n) for NL advection, internal pressure gradient,
   !                      horizontal viscosity, atmospherical forcing, 
   !                      bottom friction (2d part only)
   ! RHS(baroclinic vz_i^n,vstar) for coriolis term
   ! 2nd part : estimate of tridiagonal linear system coefficients :
   ! av,bv,cv,y2dv,axe,bxe,cxe,y2dxe
   !-------------------------------------------------------------------

   CALL dyn3dx_tls_coef(vz,boz,uz,u,v,xev, &
#ifdef key_lagoon
                        xeuv,              &
#endif
                        vstar)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Solve the barotropic system along x !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   !-----------------------------------------------------------------------
   ! Add the open boundary conditions to the right hand side
   !-----------------------------------------------------------------------

   IF ( .NOT.l_obc_cycl_x .AND. (l_obc_east .OR. l_obc_west) ) THEN
     IF_MPI ( imin==limin .OR. imax==limax ) THEN
#ifdef key_MPI_2D
!$OMP DO SCHEDULE(RUNTIME)
#else
!$OMP DO SCHEDULE(DYNAMIC,2)
#endif
     DO j=ljminp2,ljmaxm1
       !!     introduction des conditions aux limites en niveaux
       !!     liees a l implicitation du gradient de pression
       !!     Modification de coefficient
       i1=ig(j)
       i2=id(j)
       CALL dyn2dx_tls_obc(i1,i2,j,tfg(j),tfd(j),u,xev,xeuv)
     END DO
!$OMP END DO
#ifdef key_lagoon
!$OMP SINGLE
       DO ng=1,nbgrau
         i0=igrau(ng)
         j0=jgrau(ng)
         IF (type_grau(ng)==2 .AND.i0/=imax) THEN
           ! grau sur la frontiere Est
           IF_MPI ( (i0 >= limin) .AND. (i0 <= limaxp2) .AND.(j0 >= ljmin) .AND. (j0 <= ljmaxp2) ) &
           y2dv(i0-1,j0)=y2dv(i0-1,j0)-alphan*grav*halfdt/dx_u(i0-1,j0)*obc_sshe(j0)*tfd(j0)
         END IF
         IF (type_grau(ng)==4 .AND.i0/=imin+1) THEN
           ! grau sur la frontiere Ouest
           IF_MPI ( (i0 >= limin) .AND. (i0 <= limaxp2) .AND.(j0 >= ljmin) .AND. (j0 <= ljmaxp2) ) &
           y2dv(i0,j0)=y2dv(i0,j0)+alphan*grav*halfdt/dx_u(i0,j0)*obc_sshw(j0)*tfg(j0)
         END IF
       END DO
!$OMP END SINGLE
#endif

     ENDIF_MPI
   END IF

   !------------------------------------------------------------------
   ! Resolve the 2D tridiagonal system to evaluate
   ! an intermediate estimate of the zonal velocity : U_i^n+1/2
   ! and a first guess for the free surface: XE^n+1/2,*
   !------------------------------------------------------------------

   ! pre solve
   ! ---------

   CALL dyn2dx_tls_presolve

   ! ----------------------------
   ! solve the sea surface height
   ! ----------------------------

OMPMPI barrier
OMPMPI master
   ! echange throughout the lines for columns from rjminp2 to rjmaxm1
   CALL_MPI transpose_row
OMPMPI end master
OMPMPI barrier
OMPMPI flush(ry2dxe,raxe,rbxe,rcxe)

#ifdef key_MPI_2D
   IF (l_obc_cycl_x) THEN
     CALL dyn2dx_cyc_solve(rxeuv)
   ELSE
     CALL dyn2dx_tls_solve(rxeuv)
   END IF
#else
   IF (l_obc_cycl_x) THEN
     CALL dyn2dx_cyc_solve(xeuv)
   ELSE
     CALL dyn2dx_tls_solve(xeuv)
   END IF
#endif

OMPMPI barrier
OMPMPI master
   ! exchange back to each CPU (limin:limax and ljmi:ljmax)
   CALL_MPI trback_row(xeuv)
OMPMPI end master
OMPMPI barrier
OMPMPI flush(xeuv)

   !---------------------------------------------------------
   ! Apply open boundary conditions to the sea surface height
   !---------------------------------------------------------

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc_curv
     CALL_MPI ex_peri_x_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax  ! changed from DO j=jmin+2,jmax-1
       xeuv(ig(j)-1,j)=xeuv(id(j)-2,j)
       xeuv(ig(j)  ,j)=xeuv(id(j)-1,j)
       xeuv(id(j)  ,j)=xeuv(ig(j)+1,j)
#ifdef key_shiftobc_curv
       xeuv(ig(j)-2,j)=xeuv(id(j)-3,j)
       xeuv(id(j)+1,j)=xeuv(ig(j)+2,j)
       xeuv(id(j)+2,j)=xeuv(ig(j)+3,j)
#endif
     END DO
#endif
   END IF
   IF(l_obc_cycl_y)THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc_curv
     CALL_MPI ex_peri_y_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO i=limin,limax
       xeuv(i,jb(i)-1)=xeuv(i,jh(i)-2)
       xeuv(i,jb(i)  )=xeuv(i,jh(i)-1)
       xeuv(i,jh(i)  )=xeuv(i,jb(i)+1)
#ifdef key_shiftobc_curv
       xeuv(i,jb(i)-2)=xeuv(i,jh(i)-3)
       xeuv(i,jh(i)+1)=xeuv(i,jb(i)+2)
       xeuv(i,jh(i)+2)=xeuv(i,jb(i)+3)
#endif
     END DO
#endif
   END IF
!$OMP END SINGLE

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(xeuv)

   ! --------------------------------------------------
   ! post solve : estimate of the barotropic velocity U
   ! --------------------------------------------------

!$OMP DO SCHEDULE(runtime)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       IF (limin==imin) av(ig(j),j)=0.0_rlg
       IF (limax==imax) cv(id(j)-1,j)=0.0_rlg
       uv(i,j)=(y2dv(i,j)-(cv(i,j)*xeuv(i+1,j)+av(i,j)*xeuv(i,j)))*bv(i,j)
     END DO
   END DO
!$OMP END DO

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN

     CALL_MPI ex_peri_x_avcv(liminm1,limax,ljminm1,ljmax,av(liminm1:limax,ljminm1:ljmax),cv(liminm1:limax,ljminm1:ljmax))

     IF_MPI (limax == imax) THEN
     DO j=ljminp2,ljmaxm1
#ifdef key_MPI_2D
       !uv(imax-1,j) = uv(imax-1,j)-cv(imax-1,j)*xeuv(imax,j)*bv(imax-1,j)
       uv(id(j)-1,j) = uv(id(j)-1,j)-cv(id(j)-1,j)*xeuv(id(j),j)*bv(id(j)-1,j)
#else
       !uv(imax-1,j) = uv(imax-1,j)-cv(imin+1,j)*xeuv(imax,j)*bv(imax-1,j)
       uv(id(j)-1,j) = uv(id(j)-1,j)-cv(ig(j),j)*xeuv(id(j),j)*bv(id(j)-1,j)
#endif
     END DO
     ENDIF_MPI

     IF_MPI (limin == imin) THEN
     DO j=ljminp2,ljmaxm1
#ifdef key_MPI_2D
       !uv(imin+1,j) = uv(imin+1,j)-av(imin+1,j)*xeuv(imin+1,j)*bv(imin+1,j)
       uv(ig(j),j) = uv(ig(j),j)-av(ig(j),j)*xeuv(ig(j),j)*bv(ig(j),j)
#else
       !uv(imin+1,j) = uv(imin+1,j)-av(imax-1,j)*xeuv(imin+1,j)*bv(imin+1,j)
       uv(ig(j),j) = uv(ig(j),j)-av(id(j)-1,j)*xeuv(ig(j),j)*bv(ig(j),j)
#endif
     END DO
     ENDIF_MPI

   END IF
!$OMP END SINGLE

   !------------------------------------------------------
   ! apply open boundary conditions to barotropic velocity
   !------------------------------------------------------

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_char) THEN
     IF (.NOT.l_obc_cycl_x .AND. Agrif_Root_OMP)                  &
       CALL obc_apply_u_mc(u,uv,v,vstar,xev,xeuv,.TRUE.)
     IF_AGRIF (.NOT.l_obc_cycl_x .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_u(uv)
   ELSE
     IF_AGRIF (.NOT. Agrif_Root_OMP) CALL agrif_obc_u(uv)
     IF (Agrif_Root_OMP) CALL obc_apply_u(uv)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uv)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Intermediate estimate of the 3d zonal velocity)                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   ! Add the external pressure gradient (using XE^n+1/2,*)
   ! to right hand side of the zonal momentum equation
   !------------------------------------------------------

   ! 22- gradient de pression externe
   CALL dyn3dx_epg(xev,xeuv)

   ! complete estimate of the 3d zonal velocity : UZ_i^n+1/2 taking
   ! all the explicit terms
   ! and the implicit vertical diffusivity or mixing into account
   !----------------------------------------------------------

   CALL dyn3dx_ztls_solve(uz,uv)

   ! Apply open boundary conditions to intermediate estimate of UZ_i^n+1/2
   !----------------------------------------------------------------------

#if defined AGRIF_MPI && defined key_agrif_obc_sponge_uv
   IF_AGRIF (.NOT.Agrif_Root_OMP) THEN
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(uz)
   ENDIF_AGRIF
#endif


   IF (l_obc_char) THEN
OMPAGRIF BARRIER
OMPAGRIF MASTER
     IF (.NOT.l_obc_cycl_x .AND. Agrif_Root_OMP)                     &
       CALL obc_apply_uz_mc(uz,vz,boz)
     IF_AGRIF (.NOT.l_obc_cycl_x.AND..NOT.Agrif_Root_OMP) CALL agrif_obc_uz(uz,uv,xeuv)
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uz)
   ELSE
     IF_AGRIF (.NOT.Agrif_Root_OMP) CALL agrif_obc_uz(uz,uv,xeuv)
OMPAGRIF BARRIER
OMPAGRIF MASTER
     IF (Agrif_Root_OMP) CALL obc_apply_uz(uz)
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uz)
   END IF


   ! convergence of the 3d velocity towards the barotropic estimate
   !----------------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=liminp1,limaxm1
       IF (igdu(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         umoy=SUM(uz(:,i,j)*hzex(:,i,j)*dsigu(:))/hex(i,j)
#else
         umoy=SUM(uz(:,i,j)*dsigu(:))
#endif
         uz(:,i,j)=uz(:,i,j)-umoy+uv(i,j)
       ELSE
         uz(:,i,j)=uv(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_cycl_x .OR. l_obc_cycl_y) THEN
     IF_AGRIF (.NOT.Agrif_Root_OMP) CALL agrif_obc_uz(uz,uv,xeuv)
     IF (Agrif_Root_OMP) CALL obc_apply_uz(uz)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uz)

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
!part really usefull CALL ex_j_rsh(-2,2,kmax,liminm1,limax,ljminm2,ljmaxp2,uz(:,liminm1:limax,ljminm2:ljmaxp2))
!larger exchanges to make transfer in step2d* and step3d* easier
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(uz)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                             Part II                              !
   ! Estimate of final variables :                                    !
   !                         XE^n+1/2, v^n+1/2, and vz^n+1/2          !
   !                                                                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   ! Start building Ustar from the integrated right hand side
   ! At this step, Ustar does not contains the external
   ! pressure gradient term :
   ! u*^n+1/2 = u_i^n+1/2 + RHS(uz^n,vz_i^n)
   !---------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       IF (igdu(i,j)==1) THEN

         vz_u=0.25_rsh*(vz(1,i,j)+vz(1,i+1,j)+vz(1,i,j-1)+vz(1,i+1,j-1))
         vitf=SQRT(vz_u*vz_u+uz(1,i,j)*uz(1,i,j))
#if defined key_siggen || defined key_gencoord
         cdf=(0.4_rsh/LOG(0.5_rsh*dsigu(1)*hzutx(1,i,j)/z0b(i,j)+1.0_rsh))**2
#else
         cdf=(0.4_rsh/LOG((sig(1)+1.0_rsh)*htx(i,j)/z0b(i,j)+1.0_rsh))**2
#endif
#ifdef key_no_bottom_stress
         cdf=0.0_rsh
#endif
#if key_linear_bottom_stress
         cdvit(i,j)=frot_lin
#else
         cdvit(i,j)=cdf*vitf
#endif
         intz_exp(i,j)=intz_exp(i,j)-cdvit(i,j)*uz(1,i,j)/htx(i,j)*cfroi
         intz_exp_ustar(i,j) = intz_exp(i,j)
!         ustar(i,j)=uv(i,j)+halfdt*intz_exp(i,j)
       END IF
       igdu_old(i,j)=igdu(i,j)
     END DO
   END DO
!$OMP END DO

   !-----------------------------------------------------------------------
   ! Update of heights hex and hey
   ! used in continuity equation (and for tracers advection)
   ! Update of heights htx and hty, estimates of heights hzuty and hzuty,
   ! hzwty and hzwty, heights in momentum equation
   ! Definition of the type of cell (igdu and igdv),
   ! the type of extremity of the column (tfb,tfh) used in dyn2dytlsobc
   ! the mask used for the viscosity term (mo and me)
   !-----------------------------------------------------------------------

   test = 2

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminm2,ljmaxp2
     xeint(liminm2:limaxp2,j) = alphap*xev(liminm2:limaxp2,j) + alphan*xeuv(liminm2:limaxp2,j)
   END DO
!$OMP END DO

   CALL dyn3d_height_mask(xeint,tfg,tfd,tfb,tfh,test)

   ! Construction of 2D and 3D matrixes for calculation along y
   ! 1st part : estimate of right hand side :
   ! RHS(uz_i^n+1/2,vz_i^n) for NL advection, internal pressure gradient
   !                      horizontal viscosity, atmospherical forcing
   !                      bottom friction (2d part only)
   ! RHS(baroclinic vz^n,u_i^n+1/2) for coriolis term
   ! 2nd part : estimate of tridiagonal linear system coefficients :
   ! av,bv,cv,y2dv,axe,bxe,cxe,y2dxe
   !-----------------------------------------------------------------

   CALL dyn3dy_tls_coef(vz,boz,uz,u,v,xev,  &
#ifdef key_lagoon
                    xeuv,    &
#endif
                    uv)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Solve the barotropic system along y !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   !-----------------------------------------------------------------------
   ! Add the open boundary conditions to the right hand side
   !-----------------------------------------------------------------------

   IF ( .NOT.l_obc_cycl_y .AND. (l_obc_north .OR. l_obc_south) ) THEN
     IF_MPI ( jmin==ljmin .OR. jmax==ljmax ) THEN
!$OMP DO SCHEDULE(DYNAMIC,2)
     DO i=liminp2,limaxm1
       !!     introduction des conditions aux limites en niveaux
       !!     liees a l implicitation du gradient de pression
       !!     Modification de coefficient
       j1=jb(i)
       j2=jh(i)
       CALL dyn2dy_tls_obc(j1,j2,i,tfb(i),tfh(i),v,xev,xeuv)
     END DO
!$OMP END DO

#ifdef key_lagoon
!$OMP SINGLE
       DO ng=1,nbgrau
         i0=igrau(ng)
         j0=jgrau(ng)
         IF (type_grau(ng)==3 .AND. j0/=jmin+1) THEN
           ! grau sur la frontiere Sud
           IF_MPI ( (i0 >= limin) .AND. (i0 <= limaxp2) .AND.(j0 >= ljmin) .AND. (j0 <= ljmaxp2) ) &
           y2dv(i0,j0 )=y2dv(i0,j0  )+alphan*grav*halfdt/dy_v(i0,j0)*obc_sshs(i0)*tfb(i0)
         END IF
         IF (type_grau(ng)==1 .AND. j0/=jmax) THEN
           ! grau sur la frontiere Nord
           IF_MPI ( (i0 >= limin) .AND. (i0 <= limaxp2) .AND.(j0 >= ljmin) .AND. (j0 <= ljmaxp2) ) &
           y2dv(i0,j0-1)=y2dv(i0,j0-1)-alphan*grav*halfdt/dy_v(i0,j0-1)*obc_sshn(i0)*tfh(i0)
         END IF
       END DO
!$OMP END SINGLE
#endif

     ENDIF_MPI
   END IF

   !------------------------------------------------------------------
   ! Resolve the 2D tridiagonal system to evaluate
   ! an intermediate estimate of the meridional velocity : V^n+1/2
   ! and a final estimate of the free surface: XE^n+1/2
   !------------------------------------------------------------------

   ! pre solve
   ! ---------

   CALL dyn2dy_tls_presolve

   ! ----------------------------
   ! solve the sea surface height
   ! ----------------------------

OMPMPI barrier
OMPMPI master
   ! echange throughout the columns for lines from ciminp2 to cimaxm1
   CALL_MPI transpose_col
OMPMPI end master
OMPMPI barrier
OMPMPI flush(cy2dxe,caxe,cbxe,ccxe)

#ifdef key_MPI_2D
   IF (l_obc_cycl_y) THEN
     CALL dyn2dy_cyc_solve(cxeuv)
   ELSE
     CALL dyn2dy_tls_solve(cxeuv)
   END IF
#else
   IF (l_obc_cycl_y) THEN
     CALL dyn2dy_cyc_solve(xeuv)
   ELSE
     CALL dyn2dy_tls_solve(xeuv)
   END IF
#endif

OMPMPI barrier
OMPMPI master
   ! exchange back to each CPU (limin:limax and ljmi:ljmax)
   CALL_MPI trback_col(xeuv)
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(xeuv)

   !---------------------------------------------------------
   ! Apply open boundary conditions to the sea surface height
   !---------------------------------------------------------

!$OMP SINGLE
   IF (l_obc_cycl_y) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc_curv
     CALL_MPI ex_peri_y_rsh(3,1,liminm2,limaxp2,ljmin,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO i=limin,limax ! changed from DO i=imin+2,imax-1
       xeuv(i,jb(i)-1)=xeuv(i,jh(i)-2)
       xeuv(i,jb(i)  )=xeuv(i,jh(i)-1)
       xeuv(i,jh(i)  )=xeuv(i,jb(i)+1)
#ifdef key_shiftobc_curv
       xeuv(i,jb(i)-2)=xeuv(i,jh(i)-3)
       xeuv(i,jh(i)+1)=xeuv(i,jb(i)+2)
       xeuv(i,jh(i)+2)=xeuv(i,jb(i)+3)
#endif
     END DO
#endif
   END IF
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc_curv
     CALL_MPI ex_peri_x_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
       xeuv(ig(j)-1,j)=xeuv(id(j)-2,j)      ! modif N. Bruneau (val utile ?)
       xeuv(ig(j)  ,j)=xeuv(id(j)-1,j)
       xeuv(id(j)  ,j)=xeuv(ig(j)+1,j)
#ifdef key_shiftobc_curv
       xeuv(ig(j)-2,j)=xeuv(id(j)-3,j)
       xeuv(id(j)+1,j)=xeuv(ig(j)+2,j)
       xeuv(id(j)+2,j)=xeuv(ig(j)+3,j)
#endif
     END DO
#endif
   END IF
!$OMP END SINGLE

   ! --------------------------------------------------
   ! post solve : estimate of the barotropic velocity V
   ! --------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=liminp2,limaxm1
       IF (j.GE.jb(i) .AND. j.LE.jh(i)-1) THEN
         IF (ljmin==jmin) av(i,jb(i))=0.0_rlg
         IF (ljmax==jmax) cv(i,jh(i)-1)=0.0_rlg
         uv2(i,j)=(y2dv(i,j)-(cv(i,j)*xeuv(i,j+1)+av(i,j)*xeuv(i,j)))*bv(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

!$OMP SINGLE
   IF (l_obc_cycl_y) THEN

     CALL_MPI ex_peri_y_avcv(liminm1,limax,ljminm1,ljmax,av(liminm1:limax,ljminm1:ljmax),cv(liminm1:limax,ljminm1:ljmax))

     IF_MPI (ljmax == jmax) THEN
     DO i=liminp2,limaxm1
#ifdef key_MPI_2D
       !uv2(i,jmax-1) = uv2(i,jmax-1)-cv(i,jmax-1)*xeuv(i,jmax)*bv(i,jmax-1)
       uv2(i,jh(i)-1) = uv2(i,jh(i)-1)-cv(i,jh(i)-1)*xeuv(i,jh(i))*bv(i,jh(i)-1)
#else
       !uv2(i,jmax-1) = uv2(i,jmax-1)-cv(i,jmin+1)*xeuv(i,jmax)*bv(i,jmax-1)
       uv2(i,jh(i)-1) = uv2(i,jh(i)-1)-cv(i,jb(i))*xeuv(i,jh(i))*bv(i,jh(i)-1)
#endif
     END DO
     ENDIF_MPI

     IF_MPI (ljmin == jmin) THEN
     DO i=liminp2,limaxm1
#ifdef key_MPI_2D
       !uv2(i,jmin+1) = uv2(i,jmin+1)-av(i,jmin+1)*xeuv(i,jmin+1)*bv(i,jmin+1)
       uv2(i,jb(i)) = uv2(i,jb(i))-av(i,jb(i))*xeuv(i,jb(i))*bv(i,jb(i))
#else
       !uv2(i,jmin+1) = uv2(i,jmin+1)-av(i,jmax-1)*xeuv(i,jmin+1)*bv(i,jmin+1)
       uv2(i,jb(i)) = uv2(i,jb(i))-av(i,jh(i)-1)*xeuv(i,jb(i))*bv(i,jb(i))
#endif
     END DO
     ENDIF_MPI

   END IF
!$OMP END SINGLE


   !----------------------------------------------------------
   ! Apply open boundary conditions to the barotropic velocity
   !----------------------------------------------------------

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_char) THEN
     IF (.NOT.l_obc_cycl_y .AND. Agrif_Root_OMP)                   &
       CALL obc_apply_v_mc(u,uv,v,uv2,xev,xeuv,.TRUE.)
     IF_AGRIF (.NOT.l_obc_cycl_y .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_v(uv2)
   ELSE
     IF_AGRIF (.NOT. Agrif_Root_OMP) CALL agrif_obc_v(uv2)
     IF (Agrif_Root_OMP) CALL obc_apply_v(uv2)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uv2)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Solve the complete system (barotropic+baroclinic meridional velocity) !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   ! Add the external pressure gradient (using XE^n+1/2)
   ! to right hand side of the meridional momentum equation
   !-------------------------------------------------------

   CALL dyn3dy_epg(xev,xeuv)

   ! Estimate of the meridional velocity VZ^n+1/2 taking
   ! all the explicit terms
   ! and the implicit vertical diffusivity or mixing into account
   !----------------------------------------------------------

   CALL dyn3dy_ztls_solve(vz,uv2)

   ! apply open boundary conditions
   !-------------------------------

#if defined AGRIF_MPI && defined key_agrif_obc_sponge_uv
   IF_AGRIF (.NOT.Agrif_Root_OMP) THEN
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(vz)
   ENDIF_AGRIF
#endif


   IF (l_obc_char) THEN
OMPAGRIF BARRIER
OMPAGRIF MASTER
     IF (.NOT.l_obc_cycl_y .AND. Agrif_Root_OMP)                     &
       CALL obc_apply_vz_mc(uz,vz,boz)
     IF_AGRIF (.NOT.l_obc_cycl_y .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_vz(vz,uv2,xeuv)
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vz)
   ELSE
     IF_AGRIF (.NOT.Agrif_Root_OMP) CALL agrif_obc_vz(vz,uv2,xeuv)
OMPAGRIF BARRIER
OMPAGRIF MASTER
     IF (Agrif_Root_OMP) CALL obc_apply_vz(vz)
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vz)
   END IF


   ! convergence of the 3d velocity towards the barotropic estimate
   ! final estimate of the meridional velocity VZ^n+1/2
   !----------------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF (igdv(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         umoy=SUM(vz(:,i,j)*hzey(:,i,j)*dsigu(:))/hey(i,j)
#else
         umoy=SUM(vz(:,i,j)*dsigu(:))
#endif
         vz(:,i,j)=vz(:,i,j)-umoy+uv2(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_cycl_x .OR. l_obc_cycl_y) THEN
     IF_AGRIF (.NOT.Agrif_Root_OMP) CALL agrif_obc_vz(vz,uv2,xeuv)
     IF (Agrif_Root_OMP) CALL obc_apply_vz(vz)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vz)

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(vz)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                             Part III                             !
   ! Estimate of water fluxes used by the tracer advection part       !
   !                                                                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     IF ( j>=jmin+2 .AND. j<=jmax-1 ) THEN
       DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
#if defined key_siggen || defined key_gencoord
         fxeau(:,i,j)=hzex(:,i,j)*uz(:,i,j)*halfdt
#else
#ifdef key_wave_qeuler_dynamics
         fxeau(:,i,j)=hex(i,j)*(uz(:,i,j)+wave_usz(:,i,j))*halfdt
#else
         fxeau(:,i,j)=hex(i,j)*uz(:,i,j)*halfdt
#endif
#endif
       END DO
     END IF
     DO i=liminp2,limaxm1
#if defined key_siggen || defined key_gencoord
       IF ( j>=jb(i) .AND. j<=jh(i)-1 ) fyeau(:,i,j)=hzey(:,i,j)*vz(:,i,j)*halfdt
#else
#ifdef key_wave_qeuler_dynamics
       IF ( j>=jb(i) .AND. j<=jh(i)-1 ) fyeau(:,i,j)=hey(i,j)*(vz(:,i,j)+wave_vsz(:,i,j))*halfdt
#else
       IF ( j>=jb(i) .AND. j<=jh(i)-1 ) fyeau(:,i,j)=hey(i,j)*vz(:,i,j)*halfdt
#endif
#endif
     END DO
   END DO
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fxeau(:,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fyeau(:,limin:limax,ljminm1:ljmax))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(fxeau,fyeau)

#ifdef key_dyn_wetdry_fct

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                 Part IV : wet-drying purpose                     !
   !  Correction of the sea surface height and the velocities         !
   !  over very small depths (key_dyn_wetdry_fct cpp key)             !
   !                                                                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!$OMP DO SCHEDULE(runtime)
   DO j=ljmin,ljmax
     phimax(:,j)=1.0_rlg
   END DO
!$OMP END DO

   !----------------------------------------------------------------------
   ! For negative water thicknesses, removing of fluxes responsible of the
   ! decreasing of the sea surface height
   !----------------------------------------------------------------------

!$OMP DO SCHEDULE(runtime)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       IF (h0(i,j)+xeuv(i,j) < 0.0_rsh) THEN
         A=fwetp(i,j)*surf(i,j)*(xev(i,j)+h0(i,j))/(MAX(0.0_rlg,SUM(dsigu(1:kmax)*fyeau(1:kmax,i,j  )))*dx_v(i,j  ) &
                                                   -MIN(0.0_rlg,SUM(dsigu(1:kmax)*fyeau(1:kmax,i,j-1)))*dx_v(i,j-1) &
                                                   +MAX(0.0_rlg,SUM(dsigu(1:kmax)*fxeau(1:kmax,i,j  )))*dy_u(i,j  ) &
                                                   -MIN(0.0_rlg,SUM(dsigu(1:kmax)*fxeau(1:kmax,i-1,j)))*dy_u(i-1,j)+epsilon)
         phimax(i,j)=MAX(MIN(1.0_rlg, A),0.0_rlg)
       END IF
     END DO
   END DO
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rlg(0,1,1,limin,limaxp1,ljmin,ljmaxp1,phimax(limin:limaxp1,ljmin:ljmaxp1))
   CALL_MPI ex_j_rlg(0,1,1,limin,limaxp1,ljmin,ljmaxp1,phimax(limin:limaxp1,ljmin:ljmaxp1))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(phimax)

   !----------------------------------------------------------------------
   ! Limitation of the fluxes over small depths (phimax=1 everywhere else)
   !----------------------------------------------------------------------

!$OMP DO SCHEDULE(runtime)
   DO j=ljmin,ljmax
     IF ((j.GE.jmin+2).AND.(j.LE.jmax-1)) THEN
       DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
         IF (SUM(dsigu(1:kmax)*fxeau(1:kmax,i,j)) > 0.0_rlg) THEN
           fxeau(:,i,j)= phimax(i  ,j)*fxeau(:,i,j)
         ELSE
           fxeau(:,i,j)= phimax(i+1,j)*fxeau(:,i,j)
         END IF
       END DO
     END IF
     DO i=liminp2,limaxm1
       IF ((j.GE.jb(i)).AND.(j.LE.jh(i)-1)) THEN
         IF (SUM(dsigu(1:kmax)*fyeau(1:kmax,i,j)) > 0.0_rlg) THEN
           fyeau(:,i,j)= phimax(i,j  )*fyeau(:,i,j)
         ELSE
           fyeau(:,i,j)= phimax(i,j+1)*fyeau(:,i,j)
         END IF
       END IF
     END DO
   END DO
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fxeau(:,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fyeau(:,limin:limax,ljminm1:ljmax))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(fxeau,fyeau)

   !-----------------------------------
   ! Estimate of the sea surface height
   !-----------------------------------

!$OMP DO SCHEDULE(runtime)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       xeuv(i,j)=(xev(i,j)*fwetp(i,j)-(SUM(dsigu(1:kmax)*fxeau(1:kmax,i  ,j))*dy_u(i  ,j) &
                                      -SUM(dsigu(1:kmax)*fxeau(1:kmax,i-1,j))*dy_u(i-1,j) &
                                      +SUM(dsigu(1:kmax)*fyeau(1:kmax,i  ,j))*dx_v(i  ,j) &
                                      -SUM(dsigu(1:kmax)*fyeau(1:kmax,i,j-1))*dx_v(i,j-1) &
                                      -SUM(phieau(1:kmax,i,j))*halfdt )/surf(i,j))/fwet(i,j)
     END DO
   END DO
!$OMP END DO

OMPMPI barrier
OMPMPI master
!utile ?   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(xeuv)

   !---------------------------------------------------------
   ! Apply open boundary conditions to the sea surface height
   !---------------------------------------------------------

!$OMP SINGLE
   IF (l_obc_cycl_y) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc_curv
     CALL_MPI ex_peri_y_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO i=limin,limax ! changed from DO i=imin+2,imax-1
       xeuv(i,jb(i)-1)=xeuv(i,jh(i)-2)
       xeuv(i,jb(i)  )=xeuv(i,jh(i)-1)
       xeuv(i,jh(i)  )=xeuv(i,jb(i)+1)
#ifdef key_shiftobc_curv
       xeuv(i,jb(i)-2)=xeuv(i,jh(i)-3)
       xeuv(i,jh(i)+1)=xeuv(i,jb(i)+2)
       xeuv(i,jh(i)+2)=xeuv(i,jb(i)+3)
#endif
     END DO
#endif
   END IF
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc_curv
     CALL_MPI ex_peri_x_rsh(3,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     CALL_MPI ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif
#else
     DO j=ljmin,ljmax
       xeuv(ig(j)-1,j)=xeuv(id(j)-2,j)
       xeuv(ig(j)  ,j)=xeuv(id(j)-1,j)
       xeuv(id(j)  ,j)=xeuv(ig(j)+1,j)
#ifdef key_shiftobc_curv
       xeuv(ig(j)-2,j)=xeuv(id(j)-3,j)
       xeuv(id(j)+1,j)=xeuv(ig(j)+2,j)
       xeuv(id(j)+2,j)=xeuv(ig(j)+3,j)
#endif
     END DO
#endif
   END IF
!$OMP END SINGLE

   !-------------------------------------------------------
   ! Estimate of the velocities :
   !   - they balance the external pressure gradient
   !   - they do not correspond to the water fluxes (fyeau)
   !-------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=liminp2,limaxm1
       IF (j.GE.jb(i) .AND. j.LE.jh(i)-1) THEN
         IF (ljmin==jmin) av(i,jb(i))=0.0_rlg
         IF (ljmax==jmax) cv(i,jh(i)-1)=0.0_rlg
         uv2(i,j)=(y2dv(i,j)-(cv(i,j)*xeuv(i,j+1)+av(i,j)*xeuv(i,j)))*bv(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

!$OMP SINGLE
   IF (l_obc_cycl_y) THEN

     CALL_MPI ex_peri_y_avcv(liminm1,limax,ljminm1,ljmax,av(liminm1:limax,ljminm1:ljmax),cv(liminm1:limax,ljminm1:ljmax))

     IF_MPI (ljmax == jmax) THEN
     DO i=liminp2,limaxm1
#ifdef key_MPI_2D
       !uv2(i,jmax-1) = uv2(i,jmax-1)-cv(i,jmax-1)*xeuv(i,jmax)*bv(i,jmax-1)
       uv2(i,jh(i)-1) = uv2(i,jh(i)-1)-cv(i,jh(i)-1)*xeuv(i,jh(i))*bv(i,jh(i)-1)
#else
       !uv2(i,jmax-1) = uv2(i,jmax-1)-cv(i,jmin+1)*xeuv(i,jmax)*bv(i,jmax-1)
       uv2(i,jh(i)-1) = uv2(i,jh(i)-1)-cv(i,jb(i))*xeuv(i,jh(i))*bv(i,jh(i)-1)
#endif
     END DO
     ENDIF_MPI

     IF_MPI (ljmin == jmin) THEN
     DO i=liminp2,limaxm1
#ifdef key_MPI_2D
       !uv2(i,jmin+1) = uv2(i,jmin+1)-av(i,jmin+1)*xeuv(i,jmin+1)*bv(i,jmin+1)
       uv2(i,jb(i)) = uv2(i,jb(i))-av(i,jb(i))*xeuv(i,jb(i))*bv(i,jb(i))
#else
       !uv2(i,jmin+1) = uv2(i,jmin+1)-av(i,jmax-1)*xeuv(i,jmin+1)*bv(i,jmin+1)
       uv2(i,jb(i)) = uv2(i,jb(i))-av(i,jh(i)-1)*xeuv(i,jb(i))*bv(i,jb(i))
#endif
     END DO
     ENDIF_MPI

   END IF
!$OMP END SINGLE

   !----------------------------------------------------------
   ! Apply open boundary conditions to the barotropic velocity
   !----------------------------------------------------------

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_char) THEN
     IF (.NOT.l_obc_cycl_y .AND. Agrif_Root_OMP)                   &
       CALL obc_apply_v_mc(u,uv,v,uv2,xev,xeuv,.TRUE.)
     IF_AGRIF (.NOT.l_obc_cycl_y .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_v(uv2)
   ELSE
     IF_AGRIF (.NOT. Agrif_Root_OMP) CALL agrif_obc_v(uv2)
     IF (Agrif_Root_OMP) CALL obc_apply_v(uv2)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uv2)

   ! apply open boundary conditions
   !-------------------------------

#if defined AGRIF_MPI && defined key_agrif_obc_sponge_uv
   IF_AGRIF (.NOT.Agrif_Root_OMP) THEN
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(vz)
   ENDIF_AGRIF
#endif


   IF (l_obc_char) THEN
OMPAGRIF BARRIER
OMPAGRIF MASTER
     IF (.NOT.l_obc_cycl_y .AND. Agrif_Root_OMP)                     &
       CALL obc_apply_vz_mc(uz,vz,boz)
     IF_AGRIF (.NOT.l_obc_cycl_y .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_vz(vz,uv2,xeuv)
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vz)
   ELSE
     IF_AGRIF (.NOT.Agrif_Root_OMP) CALL agrif_obc_vz(vz,uv2,xeuv)
OMPAGRIF BARRIER
OMPAGRIF MASTER
     IF (Agrif_Root_OMP) CALL obc_apply_vz(vz)
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vz)
   END IF


   ! convergence of the 3d velocity towards the barotropic estimate
   ! final estimate of the meridional velocity VZ^n+1/2
   !----------------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF (igdv(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         umoy=SUM(vz(:,i,j)*hzey(:,i,j)*dsigu(:))/hey(i,j)
#else
         umoy=SUM(vz(:,i,j)*dsigu(:))
#endif
         vz(:,i,j)=vz(:,i,j)-umoy+uv2(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_cycl_x .OR. l_obc_cycl_y) THEN
     IF_AGRIF (.NOT.Agrif_Root_OMP) CALL agrif_obc_vz(vz,uv2,xeuv)
     IF (Agrif_Root_OMP) CALL obc_apply_vz(vz)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(vz)

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(vz)

#endif

   ENDIF_AGRIF

!$OMP END PARALLEL

   PRINT_DBG*, 'exit DYN3DXY routine'

  END SUBROUTINE dyn3dxy
