geometry_2d.f90

Go to the documentation of this file.

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

  USE variables, ONLY : wp , sp
  USE parameters_2d, ONLY : verbose_level

  IMPLICIT NONE

  REAL(wp), ALLOCATABLE :: x_comp(:)

  REAL(wp), ALLOCATABLE :: x_stag(:)

  REAL(wp), ALLOCATABLE :: y_comp(:)

  REAL(wp), ALLOCATABLE :: y_stag(:)

  REAL(wp), ALLOCATABLE :: grid_output(:,:)

  INTEGER, ALLOCATABLE :: source_cell(:,:)
  LOGICAL, ALLOCATABLE :: sourcee(:,:)
  LOGICAL, ALLOCATABLE :: sourcew(:,:)
  LOGICAL, ALLOCATABLE :: sources(:,:)
  LOGICAL, ALLOCATABLE :: sourcen(:,:)

  REAL(wp), ALLOCATABLE :: dist_sourcee(:,:)
  REAL(wp), ALLOCATABLE :: dist_sourcew(:,:)
  REAL(wp), ALLOCATABLE :: dist_sources(:,:)
  REAL(wp), ALLOCATABLE :: dist_sourcen(:,:)

  REAL(wp), ALLOCATABLE :: sourcee_vect_x(:,:)
  REAL(wp), ALLOCATABLE :: sourcee_vect_y(:,:)

  REAL(wp), ALLOCATABLE :: sourcew_vect_x(:,:)
  REAL(wp), ALLOCATABLE :: sourcew_vect_y(:,:)

  REAL(wp), ALLOCATABLE :: sources_vect_x(:,:)
  REAL(wp), ALLOCATABLE :: sources_vect_y(:,:)

  REAL(wp), ALLOCATABLE :: sourcen_vect_x(:,:)
  REAL(wp), ALLOCATABLE :: sourcen_vect_y(:,:)

  REAL(wp) :: pi_g

  INTEGER :: n_topography_profile_x, n_topography_profile_y

  REAL(wp) :: dx
  REAL(wp) :: x0
  REAL(wp) :: xn
  REAL(wp) :: dy
  REAL(wp) :: y0
  REAL(wp) :: yn
  REAL(wp) :: dx2
  REAL(wp) :: dy2
  INTEGER :: comp_cells_x
  INTEGER :: comp_interfaces_x
  INTEGER :: comp_cells_y
  INTEGER :: comp_interfaces_y
  REAL(wp) :: cell_size
  INTEGER :: comp_cells_xy

  INTEGER :: j_source , k_source
  INTEGER :: j_source_w , k_source_s

  REAL(wp), ALLOCATABLE :: cell_fract(:,:)
  REAL(wp), ALLOCATABLE :: dx_rel(:,:)
  REAL(wp), ALLOCATABLE :: dy_rel(:,:)

  REAL(wp), ALLOCATABLE :: upwind_dist(:,:)
  REAL(wp), ALLOCATABLE :: crosswind_dist(:,:)

CONTAINS

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

  SUBROUTINE init_grid

    USE parameters_2d, ONLY: eps_sing

    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( source_cell(comp_cells_x,comp_cells_y) )

    ! cell where are equations are solved
    source_cell(1:comp_cells_x,1:comp_cells_y) = 0


    ALLOCATE( sourcee(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcew(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcen(comp_cells_x,comp_cells_y) )
    ALLOCATE( sources(comp_cells_x,comp_cells_y) )

    ALLOCATE( dist_sourcee(comp_cells_x,comp_cells_y) )
    ALLOCATE( dist_sourcew(comp_cells_x,comp_cells_y) )
    ALLOCATE( dist_sourcen(comp_cells_x,comp_cells_y) )
    ALLOCATE( dist_sources(comp_cells_x,comp_cells_y) )

    ALLOCATE( sourcee_vect_x(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcee_vect_y(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcew_vect_x(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcew_vect_y(comp_cells_x,comp_cells_y) )
    ALLOCATE( sources_vect_x(comp_cells_x,comp_cells_y) )
    ALLOCATE( sources_vect_y(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcen_vect_x(comp_cells_x,comp_cells_y) )
    ALLOCATE( sourcen_vect_y(comp_cells_x,comp_cells_y) )

    ALLOCATE( grid_output(comp_cells_x,comp_cells_y) )

    WRITE(*,*) 'init_grid: dx = ',dx
    WRITE(*,*) 'init_grid: dy = ',dy
    
    IF ( comp_cells_x .GT. 1 ) THEN

       dx = cell_size

    ELSE
       
       dx = 1.0_wp

    END IF

    
    IF ( comp_cells_y .GT. 1 ) THEN

       dy = cell_size
    
    ELSE

       dy = 1.0_wp

    END IF

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

    dx2 = dx / 2.0_wp
    dy2 = dy / 2.0_wp


    IF ( wp .EQ. sp ) THEN

       eps_sing=min(min( dx ** 4.0_wp,dy ** 4.0_wp ),1.0e-6_wp)

    ELSE

       eps_sing=min(min( dx ** 4.0_wp,dy ** 4.0_wp ),1.0e-10_wp)

    END IF

    IF ( verbose_level .GE. 1 ) WRITE(*,*) 'eps_sing = ',eps_sing
    
    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.5_wp * ( x_stag(j) + x_stag(j+1) )

    END DO

    DO k=1,comp_cells_y

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

    END DO

    ALLOCATE( cell_fract(comp_cells_x,comp_cells_y) )
    ALLOCATE( dx_rel(comp_cells_x,comp_cells_y) )
    ALLOCATE( dy_rel(comp_cells_x,comp_cells_y) )

    ALLOCATE( upwind_dist(comp_cells_x,comp_cells_y) )
    ALLOCATE( crosswind_dist(comp_cells_x,comp_cells_y) )
    
    RETURN

  END SUBROUTINE init_grid

  SUBROUTINE deallocate_grid

    DEALLOCATE( x_comp )
    DEALLOCATE( x_stag )
    DEALLOCATE( y_comp )
    DEALLOCATE( y_stag )

    DEALLOCATE( source_cell )

    DEALLOCATE( sourcee )
    DEALLOCATE( sourcew )
    DEALLOCATE( sourcen )
    DEALLOCATE( sources )

    DEALLOCATE( dist_sourcee )
    DEALLOCATE( dist_sourcew )
    DEALLOCATE( dist_sourcen )
    DEALLOCATE( dist_sources )

    DEALLOCATE( sourcee_vect_x )
    DEALLOCATE( sourcee_vect_y )
    DEALLOCATE( sourcew_vect_x )
    DEALLOCATE( sourcew_vect_y )
    DEALLOCATE( sources_vect_x )
    DEALLOCATE( sources_vect_y )
    DEALLOCATE( sourcen_vect_x )
    DEALLOCATE( sourcen_vect_y )

    DEALLOCATE( grid_output )

    DEALLOCATE( cell_fract )
    DEALLOCATE( dx_rel )
    DEALLOCATE( dy_rel )

    DEALLOCATE( upwind_dist )
    DEALLOCATE( crosswind_dist )

    
  END SUBROUTINE deallocate_grid
  

  SUBROUTINE init_source

    USE parameters_2d, ONLY : x_source , y_source , r_source

    IMPLICIT NONE
    
    INTEGER :: j,k

    REAL(wp) :: total_source

    j_source = floor( ( x_source - x0 ) / cell_size )
    k_source = floor( ( y_source - y0 ) / cell_size )

    j_source_w = floor( ( x_source - r_source - x0 ) / cell_size ) 
    k_source_s = floor( ( y_source - r_source - y0 ) / cell_size ) 

    ! cell where are equations are solved
    source_cell(1:comp_cells_x,1:comp_cells_y) = 0

    sourcee(1:comp_cells_x,1:comp_cells_y) = .false.
    sourcew(1:comp_cells_x,1:comp_cells_y) = .false.
    sourcen(1:comp_cells_x,1:comp_cells_y) = .false.
    sources(1:comp_cells_x,1:comp_cells_y) = .false.

    dist_sourcee(1:comp_cells_x,1:comp_cells_y) = 0.0_wp
    dist_sourcew(1:comp_cells_x,1:comp_cells_y) = 0.0_wp
    dist_sourcen(1:comp_cells_x,1:comp_cells_y) = 0.0_wp
    dist_sources(1:comp_cells_x,1:comp_cells_y) = 0.0_wp

    total_source = 0.0_wp

    DO k = 2,comp_cells_y-1

       DO j = 2,comp_cells_x-1

          IF ( ( x_comp(j) - x_source )**2 + ( y_comp(k) - y_source )**2 .LE.   &
               r_source**2 ) THEN

             ! cells where equations are not solved
             source_cell(j,k) = 1 

             ! check on west cell
             IF ( ( x_comp(j-1) - x_source )**2 + ( y_comp(k) - y_source )**2   &
                  .GE. r_source**2 ) THEN

                ! cells where radial source boundary condition are applied
                source_cell(j-1,k) = 2
                sourcee(j-1,k) = .true.
                dist_sourcee(j-1,k) = sqrt( ( x_stag(j) - x_source )**2         &
                     + ( y_comp(k) - y_source )**2 )

                sourcee_vect_x(j-1,k) = ( x_stag(j) - x_source ) * r_source     &
                     / dist_sourcee(j-1,k)**2

                sourcee_vect_y(j-1,k) = ( y_comp(k) - y_source ) * r_source     &
                     / dist_sourcee(j-1,k)**2

                total_source = total_source + dx * abs( sourcee_vect_x(j-1,k) )
         
             ELSEIF ( ( x_comp(j+1) - x_source )**2 + ( y_comp(k)-y_source )**2 &
                  .GE. r_source**2 ) THEN
                ! check on east cell
       
                ! cells where radial source boundary condition are applied
                source_cell(j+1,k) = 2
                sourcew(j+1,k) = .true.
                dist_sourcew(j+1,k) = sqrt( ( x_stag(j+1) - x_source )**2       &
                     + ( y_comp(k) - y_source )**2 )

                sourcew_vect_x(j+1,k) = ( x_stag(j+1) - x_source ) * r_source   &
                     / dist_sourcew(j+1,k)**2

                sourcew_vect_y(j+1,k) = ( y_comp(k) - y_source ) * r_source     &
                     / dist_sourcew(j+1,k)**2

                total_source = total_source + dx * abs( sourcew_vect_x(j+1,k) )

             END IF
    
             ! check on south cell
             IF ( ( x_comp(j) - x_source )**2 + ( y_comp(k-1) - y_source )**2   &
                  .GE. r_source**2 ) THEN

                ! cells where radial source boundary condition are applied
                source_cell(j,k-1) = 2
                sourcen(j,k-1) = .true.
                dist_sourcen(j,k-1) = sqrt( ( x_comp(j) - x_source )**2         &
                     + ( y_stag(k) - y_source )**2 )

                sourcen_vect_x(j,k-1) = ( x_comp(j) - x_source ) * r_source     &
                     / dist_sourcen(j,k-1)**2

                sourcen_vect_y(j,k-1) = ( y_stag(k) - y_source ) * r_source     &
                     / dist_sourcen(j,k-1)**2

                total_source = total_source + dy * abs( sourcen_vect_y(j,k-1) )

             ELSEIF ( ( x_comp(j)-x_source )**2 + ( y_comp(k+1) - y_source )**2 &
                  .GE. r_source**2 ) THEN

                ! cells where radial source boundary condition are applied
                source_cell(j,k+1) = 2
                sources(j,k+1) = .true.
                dist_sources(j,k+1) = sqrt( ( x_comp(j) - x_source )**2         &
                     + ( y_stag(k+1) - y_source )**2 )
                
                sources_vect_x(j,k+1) = ( x_comp(j) - x_source ) * r_source     &
                     / dist_sources(j,k+1)**2

                sources_vect_y(j,k+1) = ( y_stag(k+1) - y_source ) * r_source   &
                     / dist_sources(j,k+1)**2

                total_source = total_source + dy * abs( sources_vect_y(j,k+1) )

             END IF

          END IF

       END DO

    END DO
    
    RETURN

  END SUBROUTINE init_source

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

  SUBROUTINE interp_1d_scalar(x1, f1, x2, f2)
    IMPLICIT NONE

    REAL(wp), INTENT(IN), DIMENSION(:) :: x1, f1
    REAL(wp), INTENT(IN) :: x2
    REAL(wp), INTENT(OUT) :: f2
    INTEGER :: n, n1x, t
    REAL(wp) :: 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.0_wp ) .AND. ( rel_pos .LE. 1.0_wp ) ) 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.0_wp ) 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(wp), INTENT(IN), DIMENSION(:,:) :: x1, y1, f1
    REAL(wp), INTENT(IN) :: x2, y2
    REAL(wp), INTENT(OUT) :: f2

    INTEGER :: ix , iy
    REAL(wp) :: 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.0_wp
       
    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.0_wp
       
    END IF
       
    IF ( size(x1,1) .EQ. 1 ) THEN
       
       f2 = alfa_y * f1(ix,iy) + ( 1.0_wp - alfa_y ) * f1(ix,iy+1)
       
    ELSEIF ( size(x1,2) .EQ. 1 ) THEN
       
       f2 = alfa_x * f1(ix,iy)  + ( 1.0_wp - alfa_x ) * f1(ix+1,iy)
       
    ELSE
       
       f2 = alfa_x * ( alfa_y * f1(ix,iy) + ( 1.0_wp - alfa_y ) * f1(ix,iy+1) ) &
            + ( 1.0_wp - alfa_x ) * ( alfa_y * f1(ix+1,iy) + ( 1.0_wp - 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(wp), INTENT(IN), DIMENSION(:) :: x1, y1
    REAL(wp), INTENT(IN), DIMENSION(:,:) :: f1
    REAL(wp), INTENT(IN) :: x2, y2
    REAL(wp), INTENT(OUT) :: f2

    INTEGER :: ix , iy
    REAL(wp) :: alfa_x , alfa_y

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

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

    ELSE

       ix = 1
       alfa_x = 0.0_wp
       
    END IF
    
    IF ( size(y1) .GT. 1 ) THEN

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

       iy = 1
       alfa_y = 0.0_wp
       
    END IF

    IF ( ( alfa_x .LT. 0.0_wp ) .OR. ( alfa_x .GT. 1.0_wp )                     &
         .OR. ( alfa_y .LT. 0.0_wp ) .OR. ( alfa_y .GT. 1.0_wp ) ) THEN

       f2 = 0.0_wp
       RETURN

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

    RETURN
    
  END SUBROUTINE interp_2d_scalarb

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

  SUBROUTINE interp_2d_slope(x1, y1, f1, x2, y2, f_x, f_y)
    IMPLICIT NONE

    REAL(wp), INTENT(IN), DIMENSION(:,:) :: x1, y1, f1
    REAL(wp), INTENT(IN) :: x2, y2
    REAL(wp), INTENT(OUT) :: f_x , f_y

    INTEGER :: ix , iy
    REAL(wp) :: alfa_x , alfa_y
    REAL(wp) :: f_x1 , f_x2 , f_y1 , f_y2


    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 = 1.0_wp
       
    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 = 1.0_wp
  
    END IF

    f_x1 = 0.0_wp
    f_x2 = 0.0_wp

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

       f_x1 = ( f1(ix+1,iy) -  f1(ix,iy) ) / (  x1(2,1) - x1(1,1) )

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

          f_x2 = ( f1(ix+1,iy+1) -  f1(ix,iy+1) ) / (  x1(2,1) - x1(1,1) )

       END IF

    END IF

    f_x = alfa_y * f_x1 + ( 1.0_wp - alfa_y ) * f_x2

    f_y1 = 0.0_wp
    f_y2 = 0.0_wp

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

       f_y1 = ( f1(ix,iy+1) - f1(ix,iy) ) / (  y1(1,2) - y1(1,1) )

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

          f_y2 = ( f1(ix+1,iy+1) - f1(ix+1,iy) ) / (  y1(1,2) - y1(1,1) )

       END IF

    END IF

    f_y = alfa_x * f_y1 + ( 1.0_wp - alfa_x ) * f_y2


    RETURN
    
  END SUBROUTINE interp_2d_slope

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

  SUBROUTINE regrid_scalar(xin, yin, fin, xl, xr , yl, yr, fout)
    IMPLICIT NONE

    REAL(wp), INTENT(IN), DIMENSION(:) :: xin, yin
    REAL(wp), INTENT(IN), DIMENSION(:,:) :: fin
    REAL(wp), INTENT(IN) :: xl, xr , yl , yr
    REAL(wp), INTENT(OUT) :: fout

    INTEGER :: ix , iy
    INTEGER :: ix1 , ix2 , iy1 , iy2
    REAL(wp) :: alfa_x , alfa_y
    REAL(wp) :: dXin , dYin
    
    INTEGER nXin,nYin
    
    nxin = size(xin)-1
    nyin = size(yin)-1

    dxin = xin(2) - xin(1)
    dyin = yin(2) - yin(1)
    
    ix1 = max(1,ceiling( ( xl - xin(1) ) / dxin ))
    ix2 = min(nxin,ceiling( ( xr -xin(1) ) / dxin )+1)
        
    iy1 = max(1,ceiling( ( yl - yin(1) ) / dyin ))
    iy2 = min(nyin,ceiling( ( yr - yin(1) ) / dyin ) + 1)

    fout = 0.0_wp
    
    DO ix=ix1,ix2-1
        
       alfa_x = ( min(xr,xin(ix+1)) - max(xl,xin(ix)) ) / ( xr - xl )

       DO iy=iy1,iy2-1
                
          alfa_y = ( min(yr,yin(iy+1)) - max(yl,yin(iy)) ) / ( yr - yl )

          fout = fout + alfa_x * alfa_y * fin(ix,iy)
                
       END DO
            
    END DO
        
  END SUBROUTINE regrid_scalar
    
  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE limit( v , z , limiter , slope_lim )

    USE parameters_2d, ONLY : theta

    IMPLICIT none

    REAL(wp), INTENT(IN) :: v(3)
    REAL(wp), INTENT(IN) :: z(3)
    INTEGER, INTENT(IN) :: limiter

    REAL(wp), INTENT(OUT) :: slope_lim

    REAL(wp) :: a , b , c

    REAL(wp) :: sigma1 , sigma2

    a = ( v(3) - v(2) ) / ( z(3) - z(2) )
    b = ( v(2) - v(1) ) / ( z(2) - z(1) )
    c = ( v(3) - v(1) ) / ( z(3) - z(1) )

    SELECT CASE (limiter)

    CASE ( 0 )

       slope_lim = 0.0_wp

    CASE ( 1 )

       ! minmod
       slope_lim = minmod(a,b)

    CASE ( 2 )

       ! superbee
       sigma1 = minmod( a , 2.0_wp*b )
       sigma2 = minmod( 2.0_wp*a , b )
       slope_lim = maxmod( sigma1 , sigma2 )

    CASE ( 3 )

       ! generalized minmod
       slope_lim = minmod( c , theta * minmod( a , b ) )

    CASE ( 4 )

       ! monotonized central-difference (MC, LeVeque p.112)
       slope_lim = minmod( c , 2.0 * minmod( a , b ) )

    CASE ( 5 )

       ! centered
       slope_lim = c

    CASE (6) 

       ! backward
       slope_lim = a

    CASE (7)
       
       !forward
       slope_lim = b

    END SELECT

  END SUBROUTINE limit


  REAL(wp) FUNCTION minmod(a,b)

    IMPLICIT none

    REAL(wp) :: a , b , sa , sb 

    IF ( a*b .EQ. 0.0_wp ) THEN

       minmod = 0.0_wp

    ELSE

       sa = a / abs(a)
       sb = b / abs(b)

       minmod = 0.5_wp * ( sa+sb ) * min( abs(a) , abs(b) )

    END IF

  END FUNCTION minmod

  REAL(wp) function maxmod(a,b)

    IMPLICIT none

    REAL(wp) :: a , b , sa , sb 

    IF ( a*b .EQ. 0.0_wp ) THEN

       maxmod = 0.0_wp

    ELSE

       sa = a / abs(a)
       sb = b / abs(b)

       maxmod = 0.5_wp * ( sa+sb ) * max( abs(a) , abs(b) )

    END IF

  END function maxmod

  SUBROUTINE compute_cell_fract

    USE parameters_2d, ONLY: x_source, y_source, r_source , dr_dz
  
    IMPLICIT NONE

    REAL(wp), ALLOCATABLE :: x_subgrid(:) , y_subgrid(:)

    REAL(wp) :: r_rel
    
    INTEGER, ALLOCATABLE :: check_subgrid(:)

    INTEGER n_points , n_points2

    REAL(wp) :: source_area
    REAL(wp) :: delta_x , delta_y
    
    INTEGER :: h ,j, k

    
    j_source = floor( ( x_source - x0 ) / cell_size )
    k_source = floor( ( y_source - y0 ) / cell_size )
    
    j_source_w = floor( ( x_source - r_source - x0 ) / cell_size ) 
    k_source_s = floor( ( y_source - r_source - y0 ) / cell_size )

    n_points = 200
    n_points2 = n_points**2

    ALLOCATE( x_subgrid(n_points2) )
    ALLOCATE( y_subgrid(n_points2) )
    ALLOCATE( check_subgrid(n_points2) )

    x_subgrid = 0.0_wp
    y_subgrid = 0.0_wp
    
    DO h = 1,n_points

       x_subgrid(h:n_points2:n_points) = dble(h)
       y_subgrid((h-1)*n_points+1:h*n_points) = dble(h)

    END DO

    x_subgrid = ( 2.0_wp * x_subgrid - 1.0_wp ) / ( 2.0_wp * dble(n_points) )
    y_subgrid = ( 2.0_wp * y_subgrid - 1.0_wp ) / ( 2.0_wp * dble(n_points) )
    
    x_subgrid = ( x_subgrid - 0.5_wp ) * dx
    y_subgrid = ( y_subgrid - 0.5_wp ) * dy
    
    DO j=1,comp_cells_x

      DO k=1,comp_cells_y

         IF ( ( x_stag(j+1) .LT. ( x_source - r_source ) ) .OR.                             &
              ( x_stag(j) .GT. ( x_source + r_source ) ) .OR.                               &
              ( y_stag(k+1) .LT. ( y_source - r_source ) ) .OR.                             &
              ( y_stag(k) .GT. ( y_source + r_source ) ) ) THEN

            cell_fract(j,k) = 0.0_wp 

         ELSE

            check_subgrid = 0
            
            WHERE ( ( x_comp(j) + x_subgrid - x_source )**2                           &
                 + ( y_comp(k) + y_subgrid - y_source )**2 < r_source**2 )
               
               check_subgrid = 1
               
            END WHERE
            
            cell_fract(j,k) = REAL(sum(check_subgrid))/n_points2

            delta_x = x_comp(j) - x_source
            delta_y = y_comp(k) - y_source
            r_rel = min( 1.0_wp , sqrt( delta_x**2 + delta_y**2 ) / r_source )

            dx_rel(j,k) = delta_x / sqrt( delta_x**2 + delta_y**2 ) * r_rel * dr_dz
            dy_rel(j,k) = delta_y / sqrt( delta_x**2 + delta_y**2 ) * r_rel * dr_dz
            
         END IF
         
      ENDDO

    ENDDO

    source_area = dx*dy*sum(cell_fract)

    IF ( verbose_level .GE. 0 ) THEN
       
       WRITE(*,*) 'Source area =',source_area,' Error =',abs( 1.0_wp -          &
            dx*dy*sum(cell_fract) / ( 4.0_wp*atan(1.0_wp)*r_source**2 ) )

    END IF
       
    DEALLOCATE( x_subgrid )
    DEALLOCATE( y_subgrid )
    DEALLOCATE( check_subgrid )

    RETURN

  END SUBROUTINE compute_cell_fract


END MODULE geometry_2d