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

  SUBROUTINE dyn3dyx(v,u,xeu,uv,uv2,xeuv,vz,uz,boz)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn3dyx  ***
!&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 : step3dyx.F90
!&E
!&E ** External calls : turb_smago,dyn3d_height_mask,
!&E                     obc_apply_u,obc_apply_uz,obc_apply_v,obc_apply_vz
!&E                     dyn2dy_tls_obc,dyn2dy_tls_solve,dyn2dy_cyc_solve
!&E                     dyn2dx_tls_obc,dyn2dx_tls_solve,dyn2dx_cyc_solve
!&E                     dyn3dy_tls_coef,dyn3dy_ztls_solve,dyn3dy_epg
!&E                     dyn3dx_tls_coef,dyn3dx_ztls_solve,dyn3dx_epg
!&E
!&E ** Used ij-arrays : igdu,igdv,igdu_old,xeu,u,v,boz
!&E
!&E ** Modified variables : uv,xeuv,uv2,uz,vz,ustar,igdv_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-08  (T. Odaka) 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                                         removing of min,max
!&E       !  2012-02  (F. Vandermeirsch) Complement of previous evolution for key_no_bottom_stress
!&E       !  2016-12  (A. Fiandrino, B. Thouvenin) introduction key_lagoon
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,    ONLY : hex,hey,h0,surf,dx_u,dx_v,dy_u,hty,z0b
   USE comvars2d,    ONLY : halfdt
   USE comvars2d,    ONLY : l_obc_cycl_x,l_obc_cycl_y,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


   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 toolmpi,  ONLY : ex_j_rsh,ex_j_rlg
   USE toolmpi,  ONLY : ex_i_rsh,ex_i_rlg,ex_peri_x_rsh,ex_peri_y_rsh
   USE toolmpi,  ONLY : ex_peri_x_avcv,ex_peri_y_avcv
   USE comcolig_mpi, ONLY : transpose_row,trback_row
   USE comcolig_mpi, ONLY : transpose_col,trback_col
   USE comcolig_mpi, ONLY : cy2dxe,caxe,cbxe,ccxe,cxeuv
   USE comcolig_mpi, ONLY : ry2dxe,raxe,rbxe,rcxe,rxeuv


   IMPLICIT NONE

!! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)     ,INTENT(in   ) :: xeu
   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,ustar
   REAL(KIND=rsh),DIMENSION(limin:limax)                     :: tfh,tfb
   REAL(KIND=rsh),DIMENSION(ljmin:ljmax)                     :: tfg,tfd


   REAL(KIND=rlg)                           :: A
   REAL(KIND=rlg),DIMENSION(limin:limaxp1,ljmin:ljmaxp1) :: phimax

   REAL(kind=rsh) uz_v,vitf,cdf

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

   !PRINT*, 'enter DYN3DYX routine'


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


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

!$OMP PARALLEL default(none) shared(xeu,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,hty     &
!$OMP&,ustar,dx_u,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                                &


!$OMP&,ljmaxp1,limaxp1                                      &

!$OMP&,ry2dxe,raxe,rbxe,rcxe,rxeuv                          &
!$OMP&,cy2dxe,caxe,cbxe,ccxe,cxeuv                          &

!$OMP&,ig,id,jb,jh,av,bv,cv,y2dv                            &
!$OMP&,z0b,cdvit,intz_exp_ustar,intz_exp_vstar              &


!$OMP&,sig                                                  &


!$OMP&,phimax,h0,surf,dx_v,dy_u,phieau,fwet,fwetp           &



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

!$OMP& FIRSTPRIVATE(Agrif_Root_OMP)                         &
!$OMP& PRIVATE(umoy,uz_v,vitf,cdf,                          &

!$OMP& A,                                                   &


!$OMP& i,j,i1,i2,j1,j2)

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

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

   test = 2

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

   !IF ( .NOT. Agrif_GrandMother_OMP) THEN

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

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

     CALL turb_smago(uz,vz,igdv)

   ELSE



   END IF

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




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

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=limin,limax
! strange, the following lines should be sufficient
!DO j=ljminp2,ljmaxm1
!DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       IF (igdu(i,j)==1 .AND. igdu_old(i,j)==1) THEN
!         ustar(i,j)=ustar(i,j)-grav*halfdt*(xeuv(i+1,j)-xeuv(i,j))/dx_u(i,j) 
         ustar(i,j) = u(i,j)+halfdt*(intz_exp_ustar(i,j)-grav*(xeu(i+1,j)-xeu(i,j))/dx_u(i,j) )
       ELSE
         ustar(i,j) = u(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

! At the beginning Ustar contains nothing so we forced it to u
!-------------------------------------------------------------
! comment out this line if you exchange the calling order of stepuv and stepvu
!   IF ( nbouc==1)  ustar(:,:)=u(:,:)

! Apply open boundary conditions to ustar
!----------------------------------------
!!!! BARRIER
!!!! MASTER
   IF (l_obc_char) THEN
     CALL obc_apply_u_mc(u,ustar,v,v,xeu,xeuv,.FALSE.)
   ELSE
     CALL obc_apply_u(ustar)
   END IF
!!!! END MASTER
!!!! BARRIER
!!!! FLUSH(ustar)

! Construction of 2D and 3D matrixes for calculation along y
! 1st part : estimate of right hand side :
! RHS(uz_i^n+1/2,vz^n+1/2) for NL advection, internal pressure gradient,
!                    horizontal viscosity, atmospherical forcing
!                    bottom friction (2d part only)
! RHS(baroclinic uz_i^n+1/2,ustar) 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,xeu,  &

                        ustar)

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


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

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


   IF (l_obc_cycl_y) THEN
     CALL dyn2dy_cyc_solve(cxeuv)
   ELSE
     CALL dyn2dy_tls_solve(cxeuv)
   END IF


!$OMP barrier
!$OMP master
! exchange back to each CPU (limin:limax and ljmi:ljmax)
   CALL trback_col(xeuv)
!$OMP end master
!$OMP barrier
!$OMP flush(xeuv)

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

!$OMP SINGLE
   IF (l_obc_cycl_y) THEN


     CALL ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))


   END IF
   IF (l_obc_cycl_x) THEN


     CALL ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))


   END IF
!$OMP END SINGLE

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP 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 ex_peri_y_avcv(liminm1,limax,ljminm1,ljmax,av(liminm1:limax,ljminm1:ljmax),cv(liminm1:limax,ljminm1:ljmax))

     IF (ljmax == jmax) THEN
     DO i=liminp2,limaxm1

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

     END DO
     ENDIF

     IF (ljmin == jmin) THEN
     DO i=liminp2,limaxm1

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

     END DO
     ENDIF

   END IF
!$OMP END SINGLE

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

!!!! BARRIER
!!!! 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 (.NOT.l_obc_cycl_y .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_v(uv)
   ELSE
     !IF (.NOT. Agrif_Root_OMP) CALL agrif_obc_v(uv)
     IF (Agrif_Root_OMP) CALL obc_apply_v(uv)
   END IF
!!!! END MASTER
!!!! BARRIER
!!!! FLUSH(uv)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Intermediate estimate of the 3d meridional velocity              !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

! 28- gradient de pression externe a l instant n
   CALL dyn3dy_epg(xeu,xeuv)

! complete estimate of the 3d meridional velocity : VZ_i^n+1 taking
! all the explicit terms (yz_exp)
! and the implicit vertical diffusivity or mixing into account
!----------------------------------------------------------

   CALL dyn3dy_ztls_solve(vz,uv)

! Apply open boundary conditions to intermediate estimate of VZ_i^n+1
!--------------------------------------------------------------------



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

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

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF (igdv(i,j)/=4) THEN

         umoy=SUM(vz(:,i,j)*dsigu(:))

         vz(:,i,j)=vz(:,i,j)-umoy+uv(i,j)
       ELSE
         vz(:,i,j)=uv(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

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

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
!part really usefull CALL ex_j_rsh(-2,2,kmax,limin,limaxp1,ljminm2,ljmaxp2,vz(:,limin:limaxp1,ljminm2:ljmaxp2))
!larger exchanges to make transfer in step2d* and step3d* easier
   CALL ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(vz)

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

! Start building Vstar from the integrated right hand side
! At this step, Vstar does not contains the external
! pressure gradient term :
! v*^n+1 = v_i^n+1 + RHS(uz_i^n+1/2,vz^n+1/2)
!---------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmaxm1
     DO i = liminp2, limaxm1
       IF (igdv(i,j)==1) THEN

         uz_v=0.25_rsh*(uz(1,i,j)+uz(1,i-1,j)+uz(1,i,j+1)+uz(1,i-1,j+1))
         vitf=SQRT(uz_v*uz_v+vz(1,i,j)*vz(1,i,j))

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



         cdvit(i,j)=cdf*vitf

         intz_exp(i,j)=intz_exp(i,j)-cdvit(i,j)*vz(1,i,j)/hty(i,j)*cfroi
         intz_exp_vstar(i,j) = intz_exp(i,j)
!         vstar(i,j)=uv(i,j)+halfdt*intz_exp(i,j)
       END IF
       igdv_old(i,j)=igdv(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 row (tfg,tfd) used in dyn2dxtlsobc
! the mask used for the viscosity term (mn and ms)
!-----------------------------------------------------------------------

   test = 1

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminm2,ljmaxp2
     xeint(liminm2:limaxp2,j) = alphap*xeu(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 x
! 1st part : estimate of right hand side :
! RHS(uz_i^n+1/2,vz_i^n+1) 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 dyn3dx_tls_coef(vz,boz,uz,u,v,xeu,  &

                    uv)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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 ( imin==limin .OR. imax==limax ) THEN

!$OMP DO SCHEDULE(RUNTIME)

     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
   END IF

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

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

   CALL dyn2dx_tls_presolve

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

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


   IF (l_obc_cycl_x) THEN
     CALL dyn2dx_cyc_solve(rxeuv)
   ELSE
     CALL dyn2dx_tls_solve(rxeuv)
   END IF


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

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

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN


     CALL ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))


   END IF
   IF (l_obc_cycl_y) THEN


     CALL ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))


   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 ex_peri_x_avcv(liminm1,limax,ljminm1,ljmax,av(liminm1:limax,ljminm1:ljmax),cv(liminm1:limax,ljminm1:ljmax))

     IF (limax == imax) THEN
     DO j=ljminp2,ljmaxm1

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

     END DO
     ENDIF

     IF (limin == imin) THEN
     DO j=ljminp2,ljmaxm1

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

     END DO
     ENDIF

   END IF
!$OMP END SINGLE


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

!!!! BARRIER
!!!! 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 (.NOT.l_obc_cycl_x .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_u(uv2)
   ELSE
     !IF (.NOT. Agrif_Root_OMP) CALL agrif_obc_u(uv2)
     IF (Agrif_Root_OMP) CALL obc_apply_u(uv2)
   END IF
!!!! END MASTER
!!!! BARRIER
!!!! FLUSH(uv2)

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

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

    CALL dyn3dx_epg(xeu,xeuv)

! Estimate of the zonal velocity UZ^n+1 taking
! all the explicit terms (yz_exp)
! and the implicit vertical diffusivity or mixing into account
!----------------------------------------------------------

   CALL dyn3dx_ztls_solve(uz,uv2)

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




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

! convergence of the 3d velocity towards the barotropic estimate
! final estimate of the zonal velocity UZ^n+1
!----------------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=liminp1,limaxm1
       IF (igdu(i,j)/=4) THEN

         umoy=SUM(uz(:,i,j)*dsigu(:))

         uz(:,i,j)=uz(:,i,j)-umoy+uv2(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

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

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(uz)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                  !
!                             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(:,i,j)=hex(i,j)*uz(:,i,j)*halfdt


       END DO
     END IF
     DO i=liminp2,limaxm1


       IF ( j>=jb(i) .AND. j<=jh(i)-1 ) fyeau(:,i,j)=hey(i,j)*vz(:,i,j)*halfdt


     END DO
   END DO
!$OMP END DO

!$OMP barrier
!$OMP master
   CALL ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fxeau(:,liminm1:limax,ljmin:ljmax))
   CALL ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fyeau(:,limin:limax,ljminm1:ljmax))
!$OMP end master
!$OMP barrier
!$OMP flush(fxeau,fyeau)



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                  !
!                 Part IV : wet-drying purpose                     !
!  Correction of the sea surface height and the velocities         !
!  over very small depths (1 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)*(xeu(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

!$OMP barrier
!$OMP master
   CALL ex_i_rlg(0,1,1,limin,limaxp1,ljmin,ljmaxp1,phimax(limin:limaxp1,ljmin:ljmaxp1))
   CALL ex_j_rlg(0,1,1,limin,limaxp1,ljmin,ljmaxp1,phimax(limin:limaxp1,ljmin:ljmaxp1))
!$OMP end master
!$OMP barrier
!$OMP 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

!$OMP barrier
!$OMP master
   CALL ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fxeau(:,liminm1:limax,ljmin:ljmax))
   CALL ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fyeau(:,limin:limax,ljminm1:ljmax))
!$OMP end master
!$OMP barrier
!$OMP 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)=(xeu(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

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
!utile   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(xeuv)

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

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN


     CALL ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))


   END IF
   IF (l_obc_cycl_y) THEN


     CALL ex_peri_y_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))


   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 ex_peri_x_avcv(liminm1,limax,ljminm1,ljmax,av(liminm1:limax,ljminm1:ljmax),cv(liminm1:limax,ljminm1:ljmax))

     IF (limax == imax) THEN
     DO j=ljminp2,ljmaxm1

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

     END DO
     ENDIF

     IF (limin == imin) THEN
     DO j=ljminp2,ljmaxm1

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

     END DO
     ENDIF

   END IF
!$OMP END SINGLE

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

!!!! BARRIER
!!!! 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 (.NOT.l_obc_cycl_x .AND. .NOT.Agrif_Root_OMP) CALL agrif_obc_u(uv2)
   ELSE
     !IF (.NOT. Agrif_Root_OMP) CALL agrif_obc_u(uv2)
     IF (Agrif_Root_OMP) CALL obc_apply_u(uv2)
   END IF
!!!! END MASTER
!!!! BARRIER
!!!! FLUSH(uv2)

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




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


! convergence of the 3d velocity towards the barotropic estimate
! final estimate of the zonal velocity UZ^n+1
!----------------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=liminp1,limaxm1
       IF (igdu(i,j)/=4) THEN

         umoy=SUM(uz(:,i,j)*dsigu(:))

         uz(:,i,j)=uz(:,i,j)-umoy+uv2(i,j)
       END IF
     END DO
   END DO
!$OMP END DO

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

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(uz)



   !ENDIF

!$OMP END PARALLEL

   !PRINT*, 'exit DYN3DYX routine'

  END SUBROUTINE dyn3dyx
