geometry.f90

Go to the documentation of this file.

!*********************************************************************
!
!*********************************************************************
MODULE geometry

  USE parameters, ONLY : verbose_level

  IMPLICIT NONE

  REAL*8, ALLOCATABLE :: x_comp(:)

  REAL*8, ALLOCATABLE :: x_stag(:)

  REAL*8, ALLOCATABLE :: B_stag(:)

  REAL*8, ALLOCATABLE :: B_cent(:)

  REAL*8, ALLOCATABLE :: B_prime(:)

  REAL*8, ALLOCATABLE :: batimetry_profile(:,:)

  INTEGER :: n_batimetry_profile

  REAL*8 :: dx
  REAL*8 :: x0
  REAL*8 :: xN
  INTEGER :: comp_cells
  INTEGER :: comp_interfaces

CONTAINS

  !*********************************************************************
  !
  !*********************************************************************

  SUBROUTINE init_grid

    USE parameters, ONLY: eps_sing, batimetry_function_flag

    IMPLICIT none

    INTEGER j      !> loop counter

    comp_interfaces = comp_cells+1

    ALLOCATE( x_comp(comp_cells) )
    ALLOCATE( x_stag(comp_interfaces) )

    ALLOCATE( b_cent(comp_cells) )
    ALLOCATE( b_prime(comp_cells) )
    ALLOCATE( b_stag(comp_interfaces) )
    
    dx = ( xn - x0 ) / comp_cells
    
    eps_sing = dx ** 4

    WRITE(*,*) 'eps_sing = ',eps_sing

    x_comp(1) = x0 + 0.5d0 * dx
    x_stag(1) = x0

    ! rescaling in the case of batimetry defined in input file
    IF(batimetry_function_flag.EQV..false.)THEN

       batimetry_profile(1,:) = x0 + ( xn - x0 ) * batimetry_profile(1,:) 

    ENDIF

    ! definition of batimetry by input file
    IF(batimetry_function_flag.EQV..false.)THEN

      b_stag(1) = batimetry_profile(2,1)

      DO j=1,comp_cells

         x_stag(j+1) = x_stag(j) + dx
       
         x_comp(j) = 0.5 * ( x_stag(j) + x_stag(j+1) )

         CALL interp_1d_scalar( batimetry_profile(1,:) , batimetry_profile(2,:) , &
              x_stag(j+1) , b_stag(j+1) )


         b_cent(j) = 0.5 * ( b_stag(j) + b_stag(j+1) )
       
         b_prime(j) = ( b_stag(j+1) - b_stag(j) ) / (  x_stag(j+1) - x_stag(j) )

         IF ( verbose_level .GE. 2 ) THEN

            WRITE(*,*) batimetry_profile(1,:) 
            WRITE(*,*) batimetry_profile(2,:)
            WRITE(*,*) x_stag(j+1) , b_stag(j+1) ,x_comp(j) , b_cent(j) , b_prime(j) 
            READ(*,*)

         END IF

      END DO

    ! definition of batimetry by a function
    ELSE

      b_stag(1) = batimetry_function(x_stag(1))

      DO j=1,comp_cells

         x_stag(j+1) = x_stag(j) + dx
       
         x_comp(j) = 0.5 * ( x_stag(j) + x_stag(j+1) )

         b_stag(j+1) = batimetry_function(x_stag(j+1))

         b_cent(j) = 0.5 * ( b_stag(j) + b_stag(j+1) )
       
         b_prime(j) = ( b_stag(j+1) - b_stag(j) ) / (  x_stag(j+1) - x_stag(j) )

         IF ( verbose_level .GE. 2 ) THEN

            WRITE(*,*) batimetry_profile(1,:) 
            WRITE(*,*) batimetry_profile(2,:)
            WRITE(*,*) x_stag(j+1) , b_stag(j+1) ,x_comp(j) , b_cent(j) , b_prime(j) 
            READ(*,*)

         END IF

      END DO

    ENDIF

  END SUBROUTINE init_grid

!---------------------------------------------------------------------------
!
!---------------------------------------------------------------------------

  SUBROUTINE interp_1d_scalar(x1, f1, x2, f2)
    IMPLICIT NONE
    
    REAL*8, INTENT(IN), DIMENSION(:) :: x1, f1
    REAL*8, INTENT(IN) :: x2
    REAL*8, INTENT(OUT) :: f2
    INTEGER :: n, n1x, t
    REAL*8 :: grad , rel_pos
    
    n1x = SIZE(x1)
  
    !
    ! ... locate the grid points near the topographic points
    ! ... and interpolate linearly the profile
    !
    t = 1

    search:DO n = 1, n1x-1

       rel_pos = ( x2 - x1(n) ) / ( x1(n+1) - x1(n) )

       IF ( ( rel_pos .GE. 0.d0 ) .AND. ( rel_pos .LE. 1.d0 ) ) THEN

          grad = ( f1(n+1)-f1(n) ) / ( x1(n+1)-x1(n) )
          f2 = f1(n) + ( x2-x1(n) ) * grad

          EXIT search
          
       ELSEIF  ( rel_pos .LT. 0.d0 ) THEN

          f2 = f1(n)
          
       ELSE

          f2 = f1(n+1)

       END IF

    END DO search

    RETURN
    
  END SUBROUTINE interp_1d_scalar


!---------------------------------------------------------------------------
!
!---------------------------------------------------------------------------
  REAL*8 FUNCTION batimetry_function(x)
    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: x

    REAL*8, PARAMETER :: pig = 4.0*atan(1.0)
    REAL*8, PARAMETER :: eps_dis = 10.0**(-8)

    ! example from Kurganov and Petrova 2007    
    IF(x.LT.0.0)THEN

      batimetry_function = 1.d0

    ELSEIF(x.GE.0.0.AND.x.LE.0.4)THEN

      batimetry_function = cos(pig*x)**2

    ELSEIF(x.GT.0.4.AND.x.LE.0.5)THEN

      batimetry_function = cos(pig*x)**2+0.25*(cos(10.0*pig*(x-0.5))+1)

    ELSEIF(x.GT.0.5.AND.x.LE.0.6)THEN

      batimetry_function = 0.5*cos(pig*x)**4+0.25*(cos(10.0*pig*(x-0.5))+1)

    ELSEIF(x.GT.0.6.AND.x.LT.1.0-eps_dis)THEN

      batimetry_function = 0.5*cos(pig*x)**4

    ELSEIF(x.GE.1.0-eps_dis.AND.x.LE.1.0+eps_dis)THEN

      batimetry_function = 0.25

    ELSEIF(x.GT.1.0+eps_dis.AND.x.LE.1.5)THEN

      batimetry_function = 0.25*sin(2*pig*(x-1))

    ELSE

      batimetry_function = 0.d0

    ENDIF

  END FUNCTION batimetry_function
  
  
END MODULE geometry