  SUBROUTINE dyn2dyx(v,u,xeu,uv,xeuv,uv2)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dyn2dxy  ***
   !&E
   !&E ** Purpose : Estimate xe^n+1 and u^n+1
   !&E              plus an intermediate variable v_i^n+1
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step2dxy.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables : ustar,vstar,xeuv,uv,uv2,xeint,fxeau,fyeau
   !&E                         htx,hex,igdu,tfg,tfd (dyn2d_height_mask)
   !&E                         hty,hey,igdv,tfb,tfh (dyn2d_height_mask)
   !&E                         ustar,vstar          (dyn2d_uvstar)
   !&E                         y2dxe,y2dv,axe,bxe,cxe,cv,bv,av (dyn2dy_tls_coef dyn2dy_tls_obc)
   !&E                         xeuv,uv              (dyn2dy_*_solve)
   !&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       !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
   !&E       !  2008-08  (T. Odaka) MPI
   !&E       !  2009-01  (V. Garnier) OMP
   !&E       !  2009-02  (A. Mano V. Garnier T. Odaka) MPI2D
   !&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  (N. Bruneau F. Dumas) coupling with swan
   !&E       !  2009-07  (V. Garnier) MPI2d for characteristics method
   !&E       !  2010-04  (T. Odaka) MPI2D (cyclic and transpose_*, tr*_back)
   !&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       !  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
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : jb,jh,ig,id,halfdt,hey,hex,l_obc_cycl_x,l_obc_cycl_y, &
                            h0,dy_u,dx_v,surf,fwet,fwetp
   USE comfluxdo,    ONLY : fxeau,fyeau,phieau
   USE obcapplyu,    ONLY : obc_apply_u,obc_apply_u_mc
   USE obcapplyv,    ONLY : obc_apply_v,obc_apply_v_mc
   USE obccombine,   ONLY : l_obc_char,l_obc_north,l_obc_south,l_obc_east,l_obc_west
   USE comcolig2d,   ONLY : av,bv,cv,y2dv
#ifdef key_waveforcing
   USE waveforcing,  ONLY : div_tens_rad
#endif
   USE_MPI toolmpi,  ONLY : ciminp2,cimaxm1,rjminp2,rjmaxm1
   USE_MPI toolmpi,  ONLY : ex_i_rlg,ex_j_rlg,ex_i_rsh,ex_j_rsh,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_col,trback_col,transpose_row,trback_row
   USE_MPI comcolig_mpi, ONLY : ry2dxe,raxe,rbxe,rcxe,rxeuv
   USE_MPI comcolig_mpi, ONLY : cy2dxe,caxe,cbxe,ccxe,cxeuv

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)    ,INTENT(inout) :: xeu
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)    ,INTENT(inout) :: xeuv
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)    :: v,u
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(inout) :: uv,uv2

   !! * Local declarations
   INTEGER                                        :: i,j,k,i1,i2,j1,j2,test
   LOGICAL                                        :: Agrif_Root_OMP,bidon
   REAL(KIND=rsh), DIMENSION(ljmin:ljmax)         :: tfg,tfd
   REAL(KIND=rsh), DIMENSION(limin:limax)         :: tfb,tfh
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2) :: xeint
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2) :: ustar,vstar
#ifdef key_dyn_wetdry_fct
   REAL(KIND=rlg)                                 :: A
   REAL(KIND=rlg), DIMENSION(limin:limaxp1,ljmin:ljmaxp1)     :: phimax
#endif

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

   PRINT_DBG*, 'enter DYN2DYX routine'

#ifdef key_MPI_2D
   ! To make sur 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()

!$OMP PARALLEL DEFAULT(NONE) SHARED(u,v,xeu,ustar,vstar,uv,uv2,xeuv,     &
!$OMP& xeint,hex,hey,fxeau,fyeau,halfdt,av,bv,cv,y2dv,                   &
#ifdef key_MPI_2D
!$OMP& ry2dxe,raxe,rbxe,rcxe,rxeuv,cy2dxe,caxe,cbxe,ccxe,cxeuv,          &
!$OMP& ljminm1,liminm1,                                                  &
#ifdef key_dyn_wetdry_fct
!$OMP& limaxp1,                                                          &
#endif
#endif
#ifdef key_dyn_wetdry_fct
!$OMP& phimax,h0,dy_u,dx_v,surf,                                         &
!$OMP& phieau,fwet,fwetp,                                                &
#endif
!$OMP& l_obc_char,l_obc_cycl_x,l_obc_cycl_y,ig,id,jb,jh,tfg,tfd,tfb,tfh, &
!$OMP& l_obc_north,l_obc_south,l_obc_east,l_obc_west,                    &
!$OMP& ljmin,ljmax,ljminp2,ljmaxp1,ljmaxm1,ljminm2,ljmaxp2,              &
!$OMP& limin,limax,liminp2,limaxm1,liminm2,limaxp2)                      &
#ifdef key_agrif
!$OMP&  SHARED(imin,imax,jmin,jmax)                                      &
#endif
!$OMP& FIRSTPRIVATE(Agrif_Root_OMP)                                      &
!$OMP& PRIVATE(i,j,j1,j2,i1,i2,                                          &
#ifdef key_dyn_wetdry_fct
!$OMP& A,                                                                &
#endif
!$OMP& test)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                             Part I                               !
   ! Estimate of intermediate variables :                             !
   !                         XE^n+1* and v_i^n+1                      !
   !                                                                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

   test = 2

   CALL dyn2d_height_mask(xeu,tfg,tfd,tfb,tfh,test)

#ifdef key_waveforcing
   ! compute the divergence of the radiation stress components
   ! ---------------------------------------------------------
   CALL div_tens_rad
#endif

   !-----------------------------------------------------------------------
   ! Estimate of ustar and vstar for the first guess
   ! of the barotropic velocity
   !-----------------------------------------------------------------------

! why is it useful ? for nbouc=1 only ?
!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmaxp1
     ustar(liminm2:limaxp2,j)=u(liminm2:limaxp2,j)
     vstar(liminm2:limaxp2,j)=v(liminm2:limaxp2,j)
   END DO
!$OMP END DO
#ifdef key_MPI_2D
   ustar(liminm2:limaxp2,ljmaxp2)=u(liminm2:limaxp2,ljmaxp2)
   vstar(liminm2:limaxp2,ljmaxp2)=v(liminm2:limaxp2,ljmaxp2)
#endif

   CALL dyn2d_uvstar(u,v,xeu,ustar,vstar)

   !-----------------------------------------------------------------------
   ! Apply open boundary conditions to ustar,vstar
   !-----------------------------------------------------------------------

OMPAGRIF BARRIER
OMPAGRIF MASTER
   IF (l_obc_char) THEN
     CALL obc_apply_u_mc(u,ustar,v,vstar,xeu,xeuv,.FALSE.)
     CALL obc_apply_v_mc(u,ustar,v,vstar,xeu,xeuv,.FALSE.)
   ELSE
     CALL obc_apply_u(ustar)
     CALL obc_apply_v(vstar)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(ustar,vstar)

   !-----------------------------------------------------------------------
   ! Construction of 2D matrix for calculation along x
   ! Estimate of tridiagonal linear system coefficients :
   ! av,bv,cv,y2dv,axe,bxe,cxe,y2dxe
   !-----------------------------------------------------------------------

   CALL dyn2dy_tls_coef(v,u,xeu,vstar,ustar)

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

   !-----------------------------------------------------------------------
   ! 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(RUNTIME)
     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,xeu,xeuv)
     END DO
!$OMP END DO
     ENDIF_MPI
   END IF

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

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

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 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
         uv(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
       !uv(i,jmax-1) = uv(i,jmax-1)-cv(i,jmax-1)*xeuv(i,jmax)*bv(i,jmax-1)
       uv(i,jh(i)-1) = uv(i,jh(i)-1)-cv(i,jh(i)-1)*xeuv(i,jh(i))*bv(i,jh(i)-1)
#else
       !uv(i,jmax-1) = uv(i,jmax-1)-cv(i,jmin+1)*xeuv(i,jmax)*bv(i,jmax-1)
       uv(i,jh(i)-1) = uv(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
       !uv(i,jmin+1) = uv(i,jmin+1)-av(i,jmin+1)*xeuv(i,jmin+1)*bv(i,jmin+1)
       uv(i,jb(i)) = uv(i,jb(i))-av(i,jb(i))*xeuv(i,jb(i))*bv(i,jb(i))
#else
       !uv(i,jmin+1) = uv(i,jmin+1)-av(i,jmax-1)*xeuv(i,jmin+1)*bv(i,jmin+1)
       uv(i,jb(i)) = uv(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 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,ustar,v,uv,xeu,xeuv,.TRUE.)
     IF_AGRIF (.NOT.l_obc_cycl_y .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_v(uv)
   ELSE
     IF_AGRIF (.NOT. Agrif_Root_OMP) CALL agrif_obc_v(uv)
     IF (Agrif_Root_OMP) CALL obc_apply_v(uv)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uv)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                             Part II                              !
   ! Estimate of final variables :                                    !
   !                         XE^n+1 and u^n+1                         !
   !                                                                  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   !-----------------------------------------------------------------------
   ! Update of heights hex and hey
   ! used in continuity equation (and for tracers advection)
   ! Update of heights htx and hty in momentum equation
   ! Definition of the type of cell (igdu and igdv),
   ! the type of extremity of the column (tfb,tfh) used in dyn2dxtlsobc
   !-----------------------------------------------------------------------

   test = 1

#ifdef key_MPI_2D
!$OMP DO SCHEDULE(RUNTIME)
#else
!$OMP DO SCHEDULE(dynamic,5)
#endif
   DO j=ljminm2,ljmaxp2
     xeint(liminm2:limaxp2,j) = 0.0_rsh
     xeint(liminm2:limaxp2,j) = alphap*xeu(liminm2:limaxp2,j) + alphan*xeuv(liminm2:limaxp2,j)
   END DO
!$OMP END DO

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

   !-----------------------------------------------------------------------
   ! Construction of 2D matrix for calculation along x
   ! Estimate of tridiagonal linear system coefficients :
   ! av,bv,cv,y2dv,axe,bxe,cxe,y2dxe
   !-----------------------------------------------------------------------

   ! update coefficients of the tridiagonal linear system
   CALL dyn2dxctymtm2(u,xeu,xeuv,v,uv,ustar,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,5)
#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,xeu,xeuv)
     END DO
!$OMP END DO
     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 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)
   CALL_MPI ex_i_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_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

   ! --------------------------------------------------
   ! 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
       uv2(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
       !uv2(imax-1,j) = uv2(imax-1,j)-cv(imax-1,j)*xeuv(imax,j)*bv(imax-1,j)
       uv2(id(j)-1,j) = uv2(id(j)-1,j)-cv(id(j)-1,j)*xeuv(id(j),j)*bv(id(j)-1,j)
#else
       !uv2(imax-1,j) = uv2(imax-1,j)-cv(imin+1,j)*xeuv(imax,j)*bv(imax-1,j)
       uv2(id(j)-1,j) = uv2(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
       !uv2(imin+1,j) = uv2(imin+1,j)-av(imin+1,j)*xeuv(imin+1,j)*bv(imin+1,j)
       uv2(ig(j),j) = uv2(ig(j),j)-av(ig(j),j)*xeuv(ig(j),j)*bv(ig(j),j)
#else
       !uv2(imin+1,j) = uv2(imin+1,j)-av(imax-1,j)*xeuv(imin+1,j)*bv(imin+1,j)
       uv2(ig(j),j) = uv2(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 the 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,uv2,v,uv,xeu,xeuv,.TRUE.)
     IF_AGRIF (.NOT.l_obc_cycl_x .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_u(uv2)
   ELSE
     IF_AGRIF (.NOT. Agrif_Root_OMP) CALL agrif_obc_u(uv2)
     IF (Agrif_Root_OMP) CALL obc_apply_u(uv2)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uv2)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !                                                                  !
   !                             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)
         fxeau(1,i,j)=hex(i,j)*(alphap*u(i,j)+alphan*uv2(i,j))*halfdt
       END DO
     END IF
     DO i=liminp2,limaxm1
       IF ( j>=jb(i) .AND. j<=jh(i)-1 ) fyeau(1,i,j)=hey(i,j)*(alphap*v(i,j)+alphan*uv(i,j))*halfdt
     END DO
   END DO
!$OMP END DO


#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

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

   !----------------------------------------------------------------------
   ! 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)*(xeu(i,j)+h0(i,j))/(MAX(0.0_rlg,fyeau(1,i,j))*dx_v(i,j  )-MIN(0.0_rlg,fyeau(1,i,j-1))*dx_v(i,j-1) &
                                           +MAX(0.0_rlg,fxeau(1,i,j))*dy_u(i,j  )-MIN(0.0_rlg,fxeau(1,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 (fxeau(1,i,j) > 0.0_rlg) THEN
           fxeau(1,i,j)= phimax(i  ,j)*fxeau(1,i,j)
         ELSE
           fxeau(1,i,j)= phimax(i+1,j)*fxeau(1,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 (fyeau(1,i,j) > 0.0_rlg) THEN
           fyeau(1,i,j)= phimax(i,j  )*fyeau(1,i,j)
         ELSE
           fyeau(1,i,j)= phimax(i,j+1)*fyeau(1,i,j)
         END IF
       END IF
     END DO
   END DO
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rlg(-1,0,1,liminm1,limax,ljmin,ljmax,fxeau(1,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_j_rlg(-1,0,1,limin,limax,ljminm1,ljmax,fyeau(1,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)=(fwetp(i,j)*xeu(i,j)-(fxeau(1,i,j)*dy_u(i  ,j)-fxeau(1,i-1,j)*dy_u(i-1,j)+ &
                                       fyeau(1,i,j)*dx_v(i  ,j)-fyeau(1,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
   CALL_MPI ex_i_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_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

   !-------------------------------------------------------
   ! 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=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
       uv2(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
       !uv2(imax-1,j) = uv2(imax-1,j)-cv(imax-1,j)*xeuv(imax,j)*bv(imax-1,j)
       uv2(id(j)-1,j) = uv2(id(j)-1,j)-cv(id(j)-1,j)*xeuv(id(j),j)*bv(id(j)-1,j)
#else
       !uv2(imax-1,j) = uv2(imax-1,j)-cv(imin+1,j)*xeuv(imax,j)*bv(imax-1,j)
       uv2(id(j)-1,j) = uv2(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
       !uv2(imin+1,j) = uv2(imin+1,j)-av(imin+1,j)*xeuv(imin+1,j)*bv(imin+1,j)
       uv2(ig(j),j) = uv2(ig(j),j)-av(ig(j),j)*xeuv(ig(j),j)*bv(ig(j),j)
#else
       !uv2(imin+1,j) = uv2(imin+1,j)-av(imax-1,j)*xeuv(imin+1,j)*bv(imin+1,j)
       uv2(ig(j),j) = uv2(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 the 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,uv2,v,uv,xeu,xeuv,.TRUE.)
     IF_AGRIF (.NOT.l_obc_cycl_x .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_u(uv2)
   ELSE
     IF_AGRIF (.NOT. Agrif_Root_OMP) CALL agrif_obc_u(uv2)
     IF (Agrif_Root_OMP) CALL obc_apply_u(uv2)
   END IF
OMPAGRIF END MASTER
OMPAGRIF BARRIER
OMPAGRIF FLUSH(uv2)

#endif

!$OMP END PARALLEL

   PRINT_DBG*, 'exit DYN2DYX routine'

  END SUBROUTINE dyn2dyx
