geometry_2d.f90

Go to the documentation of this file.

!*********************************************************************
!
!*********************************************************************
MODULE geometry_2d

  USE parameters_2d, ONLY : verbose_level

  IMPLICIT NONE

  REAL*8, ALLOCATABLE :: x_comp(:)

  REAL*8, ALLOCATABLE :: x_stag(:)

  REAL*8, ALLOCATABLE :: y_comp(:)

  REAL*8, ALLOCATABLE :: y_stag(:)

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

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

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

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

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

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

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

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

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

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

  INTEGER :: n_topography_profile_x, n_topography_profile_y

  REAL*8 :: dx
  REAL*8 :: x0
  REAL*8 :: xn
  REAL*8 :: dy
  REAL*8 :: y0
  REAL*8 :: yn
  REAL*8 :: dx2
  REAL*8 :: dy2
  INTEGER :: comp_cells_x
  INTEGER :: comp_interfaces_x
  INTEGER :: comp_cells_y
  INTEGER :: comp_interfaces_y
  REAL*8 :: cell_size

CONTAINS

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

  SUBROUTINE init_grid

    USE parameters_2d, ONLY: eps_sing, topography_demfile

    IMPLICIT none

    INTEGER j,k

    comp_interfaces_x = comp_cells_x+1
    comp_interfaces_y = comp_cells_y+1

    ALLOCATE( x_comp(comp_cells_x) )
    ALLOCATE( x_stag(comp_interfaces_x) )
    ALLOCATE( y_comp(comp_cells_y) )
    ALLOCATE( y_stag(comp_interfaces_y) )

    ALLOCATE( b_cent(comp_cells_x,comp_cells_y) )
    ALLOCATE( b_prime_x(comp_cells_x,comp_cells_y) )
    ALLOCATE( b_prime_y(comp_cells_x,comp_cells_y) )
    ALLOCATE( b_ver(comp_interfaces_x,comp_interfaces_y) )
    ALLOCATE( b_stag_x(comp_interfaces_x,comp_cells_y) )
    ALLOCATE( b_stag_y(comp_cells_x,comp_interfaces_y) )

    ALLOCATE( grid_output(comp_cells_x,comp_cells_y) )

    ALLOCATE( grav_surf(comp_cells_x,comp_cells_y) )

    IF ( comp_cells_x .GT. 1 ) THEN

       dx = cell_size

    ELSE
       
       dx = 1.d0

    END IF

    
    IF ( comp_cells_y .GT. 1 ) THEN

       dy = cell_size
    
    ELSE

       dy = 1.d0

    END IF

    xn = x0 + comp_cells_x * dx
    yn = y0 + comp_cells_y * dy

    dx2 = dx / 2.d0
    dy2 = dy / 2.d0

    ! eps_sing = MIN( dx ** 4.D0,dy ** 4.D0 )
    eps_sing=min(min( dx ** 4.d0,dy ** 4.d0 ),1.d-20)

    IF ( verbose_level .GE. 1 ) WRITE(*,*) 'eps_sing = ',eps_sing

    x_comp(1) = x0 + 0.5d0 * dx
    x_stag(1) = x0
    y_comp(1) = y0 + 0.5d0 * dy
    y_stag(1) = y0

    ! if topography is defined in file .inp we do a rescaling
    IF ( .NOT.topography_demfile ) THEN

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

       topography_profile(2,:,:) = y0 + ( yn - y0 ) * topography_profile(2,:,:)

       b_ver(1,1) = topography_profile(3,1,1)

       ! bottom row
       DO j=1,comp_cells_x

          x_stag(j+1) = x_stag(j) + dx

          IF(k.EQ.comp_cells_y)THEN

             ! right-bottom vertex
             b_ver(j+1,1)=topography_profile(3,n_topography_profile_x,1)

          ELSE

             CALL interp_1d_scalar( topography_profile(1,:,1) ,                 &
                  topography_profile(3,:,1) , x_stag(j+1) , b_ver(j+1,1) ) 

          ENDIF

          b_stag_y(j,1)=0.5*(b_ver(j+1,1)+b_ver(j,1))

       ENDDO

       ! left column
       DO k=1,comp_cells_y

          y_stag(k+1) = y_stag(k) + dy

          IF(k.EQ.comp_cells_y)THEN

             ! left-top vertex
             b_ver(1,k+1) = topography_profile(3,1,n_topography_profile_y)

          ELSE

             CALL interp_1d_scalar( topography_profile(2,1,:) ,                 &
                  topography_profile(3,1,:) , y_stag(k+1) , b_ver(1,k+1) ) 

          ENDIF

          b_stag_x(1,k)=0.5*(b_ver(1,k+1)+b_ver(1,k))

       ENDDO

       ! all the other cells
       DO j = 1,comp_cells_x

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

          DO k = 1,comp_cells_y

             y_comp(k) = 0.5 * ( y_stag(k) + y_stag(k+1) )

             ! right column
             IF ( j.EQ.comp_cells_x .AND. k.NE.comp_cells_y ) THEN

                CALL interp_1d_scalar(                                          &
                     topography_profile(2,n_topography_profile_x,:) ,           &
                     topography_profile(3,n_topography_profile_x,:) ,           &
                     y_stag(k+1) , b_ver(j+1,k+1) )
                
                ! top row
             ELSEIF ( j.NE.comp_cells_x .AND. k.EQ.comp_cells_y ) THEN

                CALL interp_1d_scalar(                                          &
                     topography_profile(1,:,n_topography_profile_y) ,           &
                     topography_profile(3,:,n_topography_profile_y) ,           &
                     x_stag(j+1) , b_ver(j+1,k+1) ) 

                ! right-top vertex
             ELSEIF ( j.EQ.comp_cells_x .AND. k.EQ.comp_cells_y ) THEN

                b_ver(j+1,k+1) = topography_profile( 3, n_topography_profile_x ,&
                     n_topography_profile_y)

                ! internal cells
             ELSE

                CALL interp_2d_scalar( topography_profile(1,:,:) ,              &
                     & topography_profile(2,:,:), topography_profile(3,:,:) ,   &
                     & x_stag(j+1), y_stag(k+1) , b_ver(j+1,k+1) )

             ENDIF

             ! Eq. 3.12 K&P
             b_cent(j,k) = 0.25 * ( b_ver(j,k) + b_ver(j+1,k) + b_ver(j,k+1)    &
                  + b_ver(j+1,k+1) )

             ! Eq. 3.13 K&P
             b_stag_x(j+1,k) = 0.5d0 * (b_ver(j+1,k+1)+b_ver(j+1,k))

             ! Eq. 3.14 K&P
             b_stag_y(j,k+1) = 0.5d0 * (b_ver(j+1,k+1)+b_ver(j,k+1))

             ! Second factor in RHS 1st Eq. 3.16 K&P
             b_prime_x(j,k) = ( b_stag_x(j+1,k) - b_stag_x(j,k) ) /             &
                  (  x_stag(j+1) - x_stag(j) )

             ! Second factor in RHS 2nd Eq. 3.16 K&P
             b_prime_y(j,k) = ( b_stag_y(j,k+1) - b_stag_y(j,k) ) /             &
                  (  y_stag(k+1) - y_stag(k) )

             IF ( verbose_level .GE. 2 ) THEN

                WRITE(*,*) topography_profile(1,:,:) 
                WRITE(*,*) topography_profile(2,:,:)
                WRITE(*,*) topography_profile(3,:,:)
                WRITE(*,*) x_stag(j+1) , y_stag(k+1) , b_stag_x(j+1,k) ,        &
                     b_stag_y(j,k+1) , x_comp(j) , y_comp(k) , b_cent(j,k) ,    &
                     b_ver(j+1,k+1) ,  b_prime_x(j,k) , b_prime_y(j,k) 
                READ(*,*)

             END IF

          END DO

       ENDDO

    ! topography from larger dem
    ELSE

       DO j=1,comp_interfaces_x

          x_stag(j) = x0 + (j-1) * dx

       END DO
          
       DO k=1,comp_interfaces_y

          y_stag(k) = y0 + (k-1) * dy
          
       END DO

       DO j=1,comp_cells_x

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

       END DO

       DO k=1,comp_cells_y

          y_comp(k) = 0.5d0 * ( y_stag(k) + y_stag(k+1) )

       END DO

       DO j=1,comp_interfaces_x
          
          DO k=1,comp_interfaces_y

             CALL interp_2d_scalar( topography_profile(1,:,:) ,                 &
                  topography_profile(2,:,:), topography_profile(3,:,:) ,        &
                  x_stag(j), y_stag(k) , b_ver(j,k) )

          END DO

       END DO
       
       DO j=1,comp_cells_x
          
          DO k=1,comp_interfaces_y
             
             ! Eq. 3.14 K&P
             b_stag_y(j,k) = 0.5d0 * ( b_ver(j+1,k) + b_ver(j,k) )
             
          END DO
          
       END DO
       
       DO j=1,comp_interfaces_x
          
          DO k=1,comp_cells_y
             
             ! Eq. 3.13 K&P
             b_stag_x(j,k) = 0.5d0 * ( b_ver(j,k+1) + b_ver(j,k) )
             
          END DO
          
       END DO

       DO j=1,comp_cells_x
          
          DO k=1,comp_cells_y

             ! Eq. 3.12 K&P
             b_cent(j,k) = 0.25d0 * ( b_ver(j,k) + b_ver(j+1,k) + b_ver(j,k+1)  &
                  + b_ver(j+1,k+1) )
             
             ! Second factor in RHS 1st Eq. 3.16 K&P
             b_prime_x(j,k) = ( b_stag_x(j+1,k) - b_stag_x(j,k) ) /             &
                  (  x_stag(j+1) - x_stag(j) )

             ! Second factor in RHS 2nd Eq. 3.16 K&P
             b_prime_y(j,k) = ( b_stag_y(j,k+1) - b_stag_y(j,k) ) /             &
                  (  y_stag(k+1) - y_stag(k) )

          END DO

       ENDDO

    ENDIF

    ! this coefficient is used when the the scalar dot between the normal to the 
    ! topography and gravity is computed
    DO j = 1,comp_cells_x

       DO k=1,comp_cells_y

          grav_surf(j,k) = - ( 1.d0/ dsqrt( 1.d0 + b_prime_x(j,k)**2            & 
               + b_prime_y(j,k)**2) )

       ENDDO

    ENDDO

  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

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

  SUBROUTINE interp_2d_scalar(x1, y1, f1, x2, y2, f2)
    IMPLICIT NONE

    REAL*8, INTENT(IN), DIMENSION(:,:) :: x1, y1, f1
    REAL*8, INTENT(IN) :: x2, y2
    REAL*8, INTENT(OUT) :: f2

    INTEGER :: ix , iy
    REAL*8 :: alfa_x , alfa_y

    IF ( size(x1,1) .GT. 1 ) THEN

       ix = floor( ( x2 - x1(1,1) ) / ( x1(2,1) - x1(1,1) ) ) + 1
       ix = min( ix , SIZE(x1,1)-1 )
       alfa_x = ( x1(ix+1,1) - x2 ) / (  x1(ix+1,1) - x1(ix,1) )

    ELSE

       ix = 1
       alfa_x = 0.d0
       
    END IF
    
    IF ( size(x1,2) .GT. 1 ) THEN

       iy = floor( ( y2 - y1(1,1) ) / ( y1(1,2) - y1(1,1) ) ) + 1
       iy = min( iy , SIZE(x1,2)-1 )
       alfa_y = ( y1(1,iy+1) - y2 ) / (  y1(1,iy+1) - y1(1,iy) )
    
    ELSE

       iy = 1
       alfa_y = 0.d0
       
    END IF
       
    IF ( size(x1,1) .EQ. 1 ) THEN
       
       f2 = alfa_y * f1(ix,iy) + ( 1.d0 - alfa_y ) * f1(ix,iy+1)
       
    ELSEIF ( size(x1,2) .EQ. 1 ) THEN
       
       f2 = alfa_x * f1(ix,iy)  + ( 1.d0 - alfa_x ) * f1(ix+1,iy)
       
    ELSE
       
       f2 = alfa_x * ( alfa_y * f1(ix,iy) + ( 1.d0 - alfa_y ) * f1(ix,iy+1) )   &
            + ( 1.d0 - alfa_x ) *  ( alfa_y * f1(ix+1,iy) + ( 1.d0 - alfa_y )   &
            * f1(ix+1,iy+1) )
       
    END IF
    
  END SUBROUTINE interp_2d_scalar


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

  SUBROUTINE interp_2d_scalarb(x1, y1, f1, x2, y2, f2)
    IMPLICIT NONE

    REAL*8, INTENT(IN), DIMENSION(:) :: x1, y1
    REAL*8, INTENT(IN), DIMENSION(:,:) :: f1
    REAL*8, INTENT(IN) :: x2, y2
    REAL*8, INTENT(OUT) :: f2

    INTEGER :: ix , iy
    REAL*8 :: alfa_x , alfa_y

    IF ( size(x1) .GT. 1 ) THEN

       ix = floor( ( x2 - x1(1) ) / ( x1(2) - x1(1) ) ) + 1
       ix = min( ix , SIZE(x1)-1 )
       alfa_x = ( x1(ix+1) - x2 ) / (  x1(ix+1) - x1(ix) )

    ELSE

       ix = 1
       alfa_x = 0.d0
       
    END IF
    
    IF ( size(y1) .GT. 1 ) THEN

       iy = floor( ( y2 - y1(1) ) / ( y1(2) - y1(1) ) ) + 1
       iy = min( iy , SIZE(x1)-1 )
       alfa_y = ( y1(iy+1) - y2 ) / (  y1(iy+1) - y1(iy) )
    
    ELSE

       iy = 1
       alfa_y = 0.d0
       
    END IF
       
    IF ( size(x1) .EQ. 1 ) THEN
       
       f2 = alfa_y * f1(ix,iy) + ( 1.d0 - alfa_y ) * f1(ix,iy+1)
       
    ELSEIF ( size(y1) .EQ. 1 ) THEN
       
       f2 = alfa_x * f1(ix,iy)  + ( 1.d0 - alfa_x ) * f1(ix+1,iy)
       
    ELSE
       
       f2 = alfa_x * ( alfa_y * f1(ix,iy) + ( 1.d0 - alfa_y ) * f1(ix,iy+1) )   &
            + ( 1.d0 - alfa_x ) *  ( alfa_y * f1(ix+1,iy) + ( 1.d0 - alfa_y )   &
            * f1(ix+1,iy+1) )
       
    END IF
    
  END SUBROUTINE interp_2d_scalarb


  !------------------------------------------------------------------------------
  !
  !------------------------------------------------------------------------------
  REAL*8 FUNCTION topography_function(x,y)
    IMPLICIT NONE

    REAL*8, INTENT(IN) :: x,y

    REAL*8, PARAMETER :: pig = 4.0*atan(1.0)
    REAL*8, PARAMETER :: eps_dis = 1.d-8
    REAL*8 :: a

    ! example 1D from Kurganov and Petrova 2007    
    !IF(y.LT.0.0)THEN
    !
    !  topography_function = 1.d0
    !
    !ELSEIF(y.GE.0.0.AND.y.LE.0.4)THEN
    !
    !  topography_function = COS(pig*y)**2
    !
    !ELSEIF(y.GT.0.4.AND.y.LE.0.5)THEN 
    !
    !  topography_function = COS(pig*y)**2+0.25*(COS(10.0*pig*(y-0.5))+1)
    !
    !ELSEIF(y.GT.0.5.AND.y.LE.0.6)THEN
    !
    !  topography_function = 0.5*COS(pig*y)**4+0.25*(COS(10.0*pig*(y-0.5))+1)
    !
    !ELSEIF(y.GT.0.6.AND.y.LT.1.0-eps_dis)THEN
    !
    !  topography_function = 0.5*COS(pig*y)**4
    !
    !ELSEIF(y.GE.1.0-eps_dis.AND.y.LE.1.0+eps_dis)THEN
    !
    !  topography_function = 0.25
    !
    !ELSEIF(y.GT.1.0+eps_dis.AND.y.LE.1.5)THEN
    !
    !  topography_function = 0.25*SIN(2*pig*(y-1))
    !
    !ELSE
    !
    !  topography_function = 0.d0
    !
    !ENDIF


    ! example 2D from Kurganov and Petrova 2007    
    IF(abs(y).LE.0.5.AND.x.LE.(y-1.0)/2.0)THEN

       a=y**2

    ELSEIF(abs(y).GT.0.5.AND.x.LE.(y-1.0)/2.0)THEN

       a=y**2+0.1*sin(pig*x)

    ELSE

       a=max(0.125,y**2+0.1*sin(pig*x))

    ENDIF

    topography_function=7.0/32.0*exp(-8.0*(x-0.3)**2-60.0*(y-0.1)**2)- &
         & 1.0/8.0*exp(-30.0*(x+0.1)**2-90.0*(y+0.2)**2) + a

  END FUNCTION topography_function


END MODULE geometry_2d