 SUBROUTINE  fct_bathy(x,y,z,zb,zbx,zby,cf)
!=========================== DOCUMENTATION =============================
!DOCDEBh*    PHY/USER/fct_bathy  *
!
!    NOM
! fct_bathy(x,y,z,zb,zbx,zby,cf)
!
!    DESCRIPTION
!*
! Initialisation du coefficient de friction, de la bathymetrie et eventuellement son gradient par une fonction utilisateur.
! Initialise the friction coefficient, the bathymetry and possibly its gradient using a user function.
!*
! on connait egalement le temps courant tc / the current time tc is also known
!
!    ENTREES / INPUT
!    x,y,z      : coordonnees du point ou effectuer le calcul 
!                / coordinates of the point where to perform the calculation
!    num_bathy  : numero de la fonction definissant la bathymetrie (variable globale)
!               / number of the function defining the bathymetry (global variable)
!    
!    SORTIES / OUTPUT
!    zb      : bathymetrie au point considere /bathymetry at the point under consideration
!    zbx,zby : gradient de la bathymetrie au point considere / gradient of bathymetry at the point under consideration
!    cf      : coefficient de friction / coefficient of friction
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ====================
!.1-----  Implicit, Use
   use phy_typ
   IMPLICIT NONE
!.2-----  Declaration
   real(kind=kind(0.d0)), intent(in)                :: x,y,z !! Coordonnees / Coordinates
   real(kind=kind(0.d0)), intent(out)               :: zb      
!! bathymetrie au point considere /bathymetry at the point under consideration
   real(kind=kind(0.d0)), intent(out)               :: zbx,zby
!! gradient de la bathymetrie au point considere / gradient of bathymetry at the point under consideration
   real(kind=kind(0.d0)), intent(out)               :: cf      
!! coefficient de friction / coefficient of friction

   real(kind=kind(0.d0))  :: h,hp,xi
   integer ::i,n

!=========================== DEBUT DU CODE EXECUTABLE ==================

  zb  = 0.d0*x*y*z
  zbx = 0.d0
  zby = 0.d0
  cf  = 0.d0
  
   SELECT CASE (num_bathy)

!------ Cas test "Bump" de Swashes
    CASE (1)
      IF(x.gt.8.D0 .and. x.lt.12.d0)then
          zb  = 0.2d0-0.05d0*(x-10.d0)**2
          zbx = -0.1d0*(x-10.d0)
      endif

!------ Cas test "McDonald" de Swashes
    CASE (2)
      n=int((1000.d0-x)*10.d0)
      zb=0.d0
      do i=1,n
          xi=1000.d0-(i*1.d0-0.5d0)/10.d0
          h=(4.d0/9.81d0)**(1.d0/3.d0)*(1.d0+0.5d0*exp(-16.d0*(xi*1.d-3-0.5d0)**2))
          hp=-(4.d0/9.81d0)**(1.d0/3.d0)*2.d0/125.d0*(xi*1.d-3-0.5d0)*exp(-16.d0*(xi*1.d-3-0.5d0)**2)
          cf=0.033d0**2
          zbx=(1.d0-4.d0/(9.81d0*h**3))*hp+4*cf/h**(10.d0/3.d0)
          zb=zb+zbx/10.d0
      enddo
      
!------ Cas "Montfaucon"
    CASE (3)
      if(x.le.2.d0)then
        zb = 6.2d0
        zbx= 0.d0
        cf = 0.d0
      endif
      if(x.gt.2.d0.and.x.le.16.d0)then
        zb =6.2d0-0.25d0*(x-2.d0)
        zbx=-0.25d0
        cf = 1.8d-4
      endif
      if(x.gt.16.d0.and.x.le.22.d0)then
        zb = 2.7d0
        zbx= 0.d0
        cf = 3.3d-4
      endif
      if(x.gt.22.d0)then
        zb =2.7d0-0.36666667d0*(x-22.d0)
        zbx=-0.25d0
        cf = 4.d-4
      endif
      

    CASE DEFAULT
      call print_err('fct_bathy', 'No user function defined')

   END SELECT
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE fct_bathy
!DOCFIN
