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

  SUBROUTINE dyn2d_height_mask(xe,tfg,tfd,tfb,tfh,test)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE dyn2d_height_mask  ***
!&E
!&E ** Purpose : Evaluation of water column height, of the kind
!&E              of cell we deal with and of the type of bound
!&E              for each row and column
!&E
!&E ** Description :
!&E
!&E ** Called by : dyn2dxy.F90,dyn2dyx.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : hx,hy,xe,dhx,goulet,govlet,id,ig,jb,jh
!&E
!&E ** Modified variables : htx,hex,igdu,tfg,tfd,mn,ms,me,mo,
!&E                         hty,hey,igdv,tfb,tfh,whx,why
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2008-07  (T. Duhaut)  Original code
!&E       !  2009-02  (A. Mano, T. Odaka, V. Garnier) MPI2D
!&E       !  2009-03  (V. Garnier) estimate of mn,ms,me,mo
!&E       !  2009-03  (V. Garnier) masks for viscosity and advection terms
!&E       !  2011-03  (F. Dumas)  bottom stress parametrisation + remove case 2 and 3
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,  ONLY  : htx,hty,hx,hy,hey,hex,goulet,govlet,    &
                           h0fond,jb,jh,ig,id,l_obc_cycl_x,l_obc_cycl_y, &
                           dhx,dhy,hm,fwet,fwetmin,dhxy,h0,fwetp
   USE comcolig2d, ONLY : igdu,igdv,mn,ms,me,mo,whx,why

   USE toolmpi,    ONLY : ex_i_rsh,ex_j_rsh,ex_i_int,ex_j_int


   IMPLICIT NONE

!! * Arguments
   INTEGER, INTENT (in)                                        :: test
   REAL(KIND=rsh), DIMENSION(limin:limax), INTENT (out)        :: tfb,tfh
   REAL(KIND=rsh), DIMENSION(ljmin:ljmax), INTENT (out)        :: tfg,tfd
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT (in) :: xe

!! * Local declarations
   INTEGER                                                     :: i,j,i1,i2,j1,j2

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

!-----------------------------------------------------------
! Initialize htx and hty :
! water depth at time n  (1st call from ligne : xev)
! water depth at time nprime (2nd call from ligne : xeint)
! water depth at time n+1/2  (1st call from colonne : xeu)
! water depth at time (n+1/2)prime (2nd call from colonne : xeint)
! htx and hty are used in momentum equations
!-----------------------------------------------------------

   IF (fwetmin/=1.0_rsh) THEN

!$OMP DO schedule(runtime)


   DO j=ljminp2,ljmaxm1
     fwetp(:,j)=fwet(:,j) ! fraction of wet surface of previous time step
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
! test useful ?
       IF (dhxy(i,j) == 0.0_rsh) THEN
         fwet(i,j)=1.0_rsh
       ELSE
         fwet(i,j)=MIN(MAX(fwetmin,(xe(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
       END IF
     END DO
     fwetp(:,j)=fwet(:,j)   ! previous time step inoreoaration for the future
! as this solution is unstable we do not use it but the model is not
! conservtive anymore above tidal flat if fwetmin/=1
   END DO
!$OMP END DO
   END IF


!$OMP DO schedule(runtime)


   DO j=ljmin,ljmaxm1
     DO i=limin,limaxm1
       htx(i,j)=hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j))
!       IF ( htx(i,j)<=0.0 ) With such a test mn and ms are set to 1 more often as htx is smaller when negative
       IF ( htx(i,j)<=0.0_rsh .AND. ( hx(i,j)+xe(i,j) > 0.0_rsh .OR. hx(i,j)+xe(i+1,j) > 0.0_rsh) ) &
         htx(i,j)=hx(i,j)+MAX(xe(i,j),xe(i+1,j))
       hty(i,j)=hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1))
       IF ( hty(i,j)<=0.0_rsh .AND. ( hy(i,j)+xe(i,j) > 0.0_rsh .OR. hy(i,j)+xe(i,j+1) > 0.0_rsh )) &
         hty(i,j)=hy(i,j)+MAX(xe(i,j),xe(i,j+1))
     END DO
   END DO
!$OMP END DO

!$OMP SINGLE
   IF (ljmax==jmax) THEN
     DO i=limin,limaxm1
       IF (hx(i,jmax) /= -valmanq-h0fond) THEN
         htx(i,jmax)=hx(i,jmax)+0.5_rsh*(xe(i,jmax)+xe(i+1,jmax))
         IF ( htx(i,jmax)<=0.0_rsh .AND. ( hx(i,jmax)+xe(i,jmax) > 0.0_rsh .OR. hx(i,jmax)+xe(i+1,jmax) > 0.0_rsh) ) &
           htx(i,jmax)= hx(i,jmax)+MAX(xe(i,jmax),xe(i+1,jmax))
       END IF
     END DO
   ENDIF

   IF (limax==imax) THEN
     DO j=ljmin,ljmaxm1
       IF (hy(imax,j) /= -valmanq-h0fond) THEN
         hty(imax,j)= hy(imax,j)+0.5_rsh*(xe(imax,j)+xe(imax,j+1))
         IF ( hty(imax,j)<=0.0_rsh .AND. ( hy(imax,j)+xe(imax,j) > 0.0_rsh .OR. hy(imax,j)+xe(imax,j+1) > 0.0_rsh )) &
           hty(imax,j)=hy(imax,j)+MAX(xe(imax,j),xe(imax,j+1))
       END IF
     END DO
   ENDIF
!$OMP END SINGLE

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



!$OMP DO schedule(runtime)

! FVDM pour definir ms au sud j=1 et ameliorer operateur d OBC
! et redescendre le calcul des termes termes NL des schemas quick et quickest
! mais pas compatible avec obc_mars
!DO j=ljminp1,ljmaxm1
   DO j=ljminp2,ljmaxm1

!-----------------------------------------------------------
! Evaluation of whx, variable which determines the location
! of dry points along the x-direction (wall at location hx)
! Evaluation of igdu which determines the nature
! of the cell in the estimate of u
! igdu = 1 water at both sides
! igdu = 5 water at both sides but very shallow depths (2d estimate)
! igdu = 2 water on the left
! igdu = 3 water on the right
! igdu = 4 no flow inside this cell
! Initialize hex :
! water depth at time n  (1st call from ligne : xev)
! water depth at time nprime (2nd call from ligne : xeint)
! water depth at time n+1/2  (1st call from colonne : xeu)
! water depth at time (n+1/2)prime (2nd call from colonne : xeint)
! hex and hey are used in the continuity and for tracers advection
!-----------------------------------------------------------

     i1=ig(j)
     i2=id(j)
     DO i=MAX0(limin,i1),MIN0(limax,i2-1)

       IF (htx(i,j) > 0.0_rsh) THEN
         whx(i,j)=1
         igdu(i,j)=1

         hex(i,j)=(htx(i,j)+dhx(i,j))*goulet(i,j)

       ELSE IF (hx(i,j)+xe(i,j) > 0.0_rsh) THEN
         IF (hx(i,j)+xe(i,j) >= epsilon) igdu(i,j)=1
         IF (hx(i,j)+xe(i,j) <  epsilon) igdu(i,j)=4

         hex(i,j)=(hx(i,j)+xe(i,j)+dhx(i,j))*goulet(i,j)

       ELSE IF (hx(i,j)+xe(i+1,j) > 0.0_rsh) THEN
         IF (hx(i,j)+xe(i+1,j) >= epsilon) igdu(i,j)=1
         IF (hx(i,j)+xe(i+1,j) <  epsilon) igdu(i,j)=4

         hex(i,j)=(hx(i,j)+xe(i+1,j)+dhx(i,j))*goulet(i,j)

       ELSE
         whx(i,j)=0
         igdu(i,j)=4
         hex(i,j)=0.0_rsh
       END IF

!-------------------------------------------
! recherche de la presence de mur au noeud
! phi necessaire au calcul des gradients
! perpendiculaire a la direction de calcul
! look for the presence of a wall at the north and south
! used for estimate of viscosity and advection terms
!-------------------------------------------

       IF (hty(i,j  )+hty(i+1,j  ) > 0.0_rsh .AND. htx(i,j+1) > 0.0_rsh) THEN
         mn(i,j)=1
       ELSE
         mn(i,j)=0
       END IF
       IF (hty(i,j-1)+hty(i+1,j-1) > 0.0_rsh .AND. htx(i,j-1) > 0.0_rsh) THEN
         ms(i,j)=1
       ELSE
         ms(i,j)=0
       END IF


!-------------------------------------------
! type of beginning and end of the line for boundary conditions
! tfg = 1 water at both sides = flow at open boundary (westwards or eastwards
! tfg = 0 no flow from the boundary (open or shore)
! tfg = 2 flow towards the east only
! tfd = 1 water at both sides = flow at open boundary
! tfd = 0 no flow from the boundary (open or shore)
! tfd = 2 flow towards the west only
!-------------------------------------------

       IF (test == 1) THEN
       IF (.NOT.l_obc_cycl_x) THEN
         IF (i.EQ. i1) THEN
           IF (igdu(i1,j) == 1)                      tfg(j) = 1.0_rsh
           IF (igdu(i1,j) == 4)                      tfg(j) = 0.0_rsh
         END IF

         IF (i.EQ.i2-1) THEN
           IF (igdu(i2-1,j) == 1)                    tfd(j) = 1.0_rsh
           IF (igdu(i2-1,j) == 4)                    tfd(j) = 0.0_rsh
         END IF
       END IF
       END IF

!-----------------------------------------------------------
! Evaluation of whx, variable which determines the location
! of dry points along the x-direction (wall at location hx)
! hypothesis : straight walls at obc
!-----------------------------------------------------------
       IF (i.EQ. i1) THEN
         whx(i1-1,j)=whx(i1,j)
       END IF
       IF (i.EQ.i2-1) THEN
         whx(i2,j)=whx(i2-1,j)
       END IF

     END DO  ! ends loop on i
   END DO
!$OMP END DO


!$OMP DO SCHEDULE(runtime)

   DO j=ljmin,ljmax

     DO i=liminp2,limaxm1
       j1=jb(i)
       j2=jh(i)
! FVDM pour definir mo a ouest i=1 et ameliorer operateur d OBC
! et redescendre le calcul des termes termes NL des schemas quick et quickest
!IF ( (i.GE.liminp1.AND.i.LE.limaxm1) .AND. (j>=j1 .AND. j<=j2-1) ) THEN
! mais pas compatible avec obc_mars

       IF ( j>=j1 .AND. j<=j2-1 ) THEN

!-----------------------------------------------------------
! Evaluation of why variable which determines the location
! of dry points along the y-direction (wall at location hy)
! Evaluation of igdv which determines the nature
! of the cell in the estimate of v
! igdv = 1 water at both sides
! igdv = 5 water at both sides but very shallow depths (2d estimate)
! igdv = 2 water at the south
! igdv = 3 water at the north
! igdv = 4 no flow inside this cell
! Initialize hey :
! water depth at time n  (1st call from ligne : xev)
! water depth at time nprime (2nd call from ligne : xeint)
! water depth at time n+1/2  (1st call from colonne : xeu)
! water depth at time (n+1/2)prime (2nd call from colonne : xeint)
! hex and hey are used in the continuity and for tracers advection
!-----------------------------------------------------------

         IF (hty(i,j) > 0.0_rsh) THEN
           why(i,j)=1
           igdv(i,j)=1

           hey(i,j)=(hty(i,j)+dhy(i,j))*govlet(i,j)

         ELSE IF (hy(i,j)+xe(i,j) > 0.0_rsh) THEN
           IF (hy(i,j)+xe(i,j) >= epsilon) igdv(i,j)=1
           IF (hy(i,j)+xe(i,j) <  epsilon) igdv(i,j)=4

           hey(i,j)=(hy(i,j)+xe(i,j)+dhy(i,j))*govlet(i,j)

         ELSE IF (hy(i,j)+xe(i,j+1) > 0.0_rsh) THEN
           IF (hy(i,j)+xe(i,j+1) >= epsilon) igdv(i,j)=1
           IF (hy(i,j)+xe(i,j+1) <  epsilon) igdv(i,j)=4

           hey(i,j)=(hy(i,j)+xe(i,j+1)+dhy(i,j))*govlet(i,j)

         ELSE
           why(i,j)=0
           igdv(i,j)=4
           hey(i,j)=0.0_rsh
         END IF

!---------------------------------------------------------
!  detection des murs aux noeuds phi pour l evaluation des
!  gradients perpendiculaires a la direction de
!  calcul (ici les gradients selon Ox)
! look for the presence of a wall at the west and east
! used for estimate of viscosity and advection terms
!-------------------------------------------

         IF (htx(i  ,j)+htx(i  ,j+1) > 0.0_rsh .AND. hty(i+1,j) > 0.0_rsh) THEN
           me(i,j)=1
         ELSE
           me(i,j)=0
         END IF
         IF (htx(i-1,j)+htx(i-1,j+1) > 0.0_rsh .AND. hty(i-1,j) > 0.0_rsh) THEN
           mo(i,j)=1
         ELSE
           mo(i,j)=0
         END IF

!-------------------------------------------
! type of beginning and end of the line for boundary conditions
! tfb = 1 water at both sides = flow at open boundary (northwards or southwards)
! tfb = 0 no flow from the boundary (open or shore)
! tfb = 2 flow towards the north only
! tfh = 1 water at both sides = flow at open boundary
! tfh = 0 no flow from the boundary (open or shore)
! tfh = 2 flow towards the south only
!-------------------------------------------

         IF (test == 2) THEN
           IF (.NOT.l_obc_cycl_y) THEN
             IF (j.EQ. j1) THEN
               IF (igdv(i,j1) == 1)                      tfb(i) = 1.0_rsh
               IF (igdv(i,j1) == 4)                      tfb(i) = 0.0_rsh
             END IF

             IF (j.EQ. j2-1) THEN
               IF (igdv(i,j2-1) == 1)                    tfh(i) = 1.0_rsh
               IF (igdv(i,j2-1) == 4)                    tfh(i) = 0.0_rsh
             END IF
           END IF
         END IF


!-----------------------------------------------------------
! Evaluation of why, variable which determines the location
! of dry points along the x-direction (wall at location hy)
! hypothesis : straight walls at obc
!-----------------------------------------------------------
         IF (j.EQ. j1) THEN
           why(i,j1-1)=why(i,j1)
         END IF
         IF (j.EQ. j2-1) THEN
           why(i,j2)=why(i,j2-1)
         END IF

       END IF

     END DO
   END DO
!$OMP END DO

!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljmin,ljmax,hex(liminm2:limaxp2,ljmin:ljmax))
   CALL ex_j_rsh(-2,2,1,limin,limax,ljminm2,ljmaxp2,hey(limin:limax,ljminm2:ljmaxp2))
   CALL ex_i_int(-1,1,1,liminm1,limaxp1,ljmin,ljmax,whx(liminm1:limaxp1,ljmin:ljmax))
   CALL ex_j_int(-1,1,1,limin,limax,ljminm1,ljmaxp1,why(limin:limax,ljminm1:ljmaxp1))
   CALL ex_i_int(-1,0,1,liminm1,limax  ,ljmin,ljmax,mo(liminm1:limax  ,ljmin:ljmax))
   CALL ex_i_int(-1,1,1,liminm1,limaxp1,ljminm1,ljmax,me(liminm1:limaxp1,ljminm1:ljmax))
   CALL ex_j_int(-1,0,1,liminm1,limaxp1,ljminm1,ljmax,me(liminm1:limaxp1,ljminm1:ljmax))
   CALL ex_j_int(-1,0,1,limin,limax,ljminm1,ljmax  ,ms(limin:limax,ljminm1:ljmax  ))
   CALL ex_i_int(-1,0,1,liminm1,limax,ljminm1,ljmaxp1,mn(liminm1:limax,ljminm1:ljmaxp1))
   CALL ex_j_int(-1,1,1,liminm1,limax,ljminm1,ljmaxp1,mn(liminm1:limax,ljminm1:ljmaxp1))
!$OMP end master
!$OMP barrier
!$OMP flush(hex,hey,whx,why,mn,ms,me,mo)
 
  END SUBROUTINE dyn2d_height_mask
