75 WRITE(*,*)
'Riemann problem initialization' 76 WRITE(*,*)
'x_comp(1)',
x_comp(1)
93 EXIT riemann_int_search
97 END DO riemann_int_search
106 IF ( solid_transport_flag )
THEN 110 IF ( temperature_flag ) qp(5,1:i1,:) =
t_w 114 IF ( temperature_flag ) qp(4,1:i1,:) =
t_w 132 WRITE(*,*) j,k,
b_cent(j,k)
149 IF ( solid_transport_flag )
THEN 175 WRITE(*,*) j,k,
b_cent(j,k)
215 REAL*8 :: qp(
n_vars,comp_cells_x,comp_cells_y) , qj(
n_vars)
229 IF ( solid_transport_flag )
THEN 233 IF ( temperature_flag )
THEN 241 IF ( temperature_flag )
THEN 256 qp(1,:,:) = b_cent(:,:) + ( qp(1,:,:)-b_cent(:,:) ) *
released_volume &
257 / (
dx *
dy * sum( qp(1,:,:)-b_cent(:,:) ) )
261 WRITE(*,*)
'Initial volume =',
dx *
dy * sum( qp(1,:,:)-b_cent(:,:) )
264 DO j = 1,comp_cells_x
266 DO k = 1,comp_cells_y
269 CALL qp_to_qc( qp(:,j,k) , b_cent(j,k) , qj )
273 IF (
verbose_level .GE. 1 )
WRITE(*,*) j,k,b_cent(j,k),qp(:,j,k)
308 REAL*8,
ALLOCATABLE :: x_subgrid(:) , y_subgrid(:)
309 INTEGER,
ALLOCATABLE :: check_subgrid(:)
311 INTEGER n_points , n_points2
312 REAL*8 :: source_area
315 n_points2 = n_points**2
317 ALLOCATE( x_subgrid(n_points2) )
318 ALLOCATE( y_subgrid(n_points2) )
319 ALLOCATE( check_subgrid(n_points2) )
326 x_subgrid(h:n_points2:n_points) = h
327 y_subgrid((h-1)*n_points+1:h*n_points) = h
331 x_subgrid = x_subgrid / ( n_points +1 )
332 y_subgrid = y_subgrid / ( n_points +1 )
334 x_subgrid = ( x_subgrid - 0.5d0 ) *
dx 335 y_subgrid = ( y_subgrid - 0.5d0 ) *
dy 345 WHERE ( ( x_comp(j) + x_subgrid - x_source )**2 &
346 + ( y_comp(k) + y_subgrid - y_source )**2 < r_source**2 )
354 source_xy(j,k) =
REAL(sum(check_subgrid))/n_points2
362 WRITE(*,*)
'Source area =',source_area,
' Error =',abs( 1.d0 - &
388 REAL*8,
INTENT(IN) :: x,y
389 REAL*8,
INTENT(IN) :: Bj
391 REAL*8,
PARAMETER :: pig = 4.d0*atan(1.d0)
438 REAL*8,
INTENT(IN) :: x,y
440 REAL*8,
PARAMETER :: pig = 4.0*atan(1.0)
443 INTEGER :: idx1 , idy1
444 REAL*8 :: coeff_x , coeff_y
446 REAL*8 :: dBdx1 , dBdx2 , dBdx
447 REAL*8 :: dBdy1 , dBdy2 , dBdy
449 REAL*8 :: max_slope_angle_rad
450 REAL*8 :: velocity_ang_release_rad
453 r = ( released_volume / pig )**(1.d0/3.d0)
456 idx1 = ceiling( ( x_release - x0 ) / dx )
457 idy1 = ceiling( ( y_release - y0 ) / dy )
460 coeff_x = ( x_release -
x_stag(idx1) ) / dx
461 coeff_y = ( y_release -
y_stag(idy1) ) / dy
465 dbdx1 = (
b_ver(idx1+1,idy1) -
b_ver(idx1,idy1) ) / dx
466 dbdx2 = (
b_ver(idx1+1,idy1+1) -
b_ver(idx1,idy1+1) ) / dx
469 dbdx = coeff_y * dbdx2 + ( 1.d0 - coeff_y ) * dbdx1
472 dbdy1 = (
b_ver(idx1,idy1+1) -
b_ver(idx1,idy1) ) / dy
473 dbdy2 = (
b_ver(idx1+1,idy1+1) -
b_ver(idx1+1,idy1) ) / dx
476 dbdy = coeff_x * dbdy2 + ( 1.d0 - coeff_x ) * dbdy1
480 max_slope_angle_rad = datan2(dbdy,dbdx)
482 IF ( verbose_level .GE. 1 )
THEN 485 WRITE(*,*)
'---',
x_stag(idx1),x_release,
x_stag(idx1+1),coeff_x
486 WRITE(*,*)
'---',
y_stag(idy1),y_release,
y_stag(idy1+1),coeff_y
488 WRITE(*,*)
b_ver(idx1,idy1),
b_ver(idx1+1,idy1)
489 WRITE(*,*)
b_ver(idx1,idy1+1),
b_ver(idx1+1,idy1+1)
490 WRITE(*,*)
'dbdx,dbdy,angle',dbdx,dbdy,max_slope_angle_rad*180.d0/pig
496 IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r )
THEN 510 * cos( max_slope_angle_rad + velocity_ang_release_rad )
542 REAL*8,
INTENT(IN) :: x,y
544 REAL*8,
PARAMETER :: pig = 4.0*atan(1.0)
547 INTEGER :: idx1 , idy1
548 REAL*8 :: coeff_x , coeff_y
550 REAL*8 :: dBdx1 , dBdx2 , dBdx
551 REAL*8 :: dBdy1 , dBdy2 , dBdy
553 REAL*8 :: max_slope_angle_rad
554 REAL*8 :: velocity_ang_release_rad
557 r = ( released_volume / pig )**(1.0/3.0)
560 idx1 = ceiling( ( x_release - x0 ) / dx )
561 idy1 = ceiling( ( y_release - y0 ) / dy )
564 coeff_x = ( x_release -
x_stag(idx1) ) / dx
565 coeff_y = ( y_release -
y_stag(idy1) ) / dy
569 dbdx1 = (
b_ver(idx1+1,idy1) -
b_ver(idx1,idy1) ) / dx
570 dbdx2 = (
b_ver(idx1+1,idy1+1) -
b_ver(idx1,idy1+1) ) / dx
573 dbdx = coeff_y * dbdx2 + ( 1.d0 - coeff_y ) * dbdx1
576 dbdy1 = (
b_ver(idx1,idy1+1) -
b_ver(idx1,idy1) ) / dy
577 dbdy2 = (
b_ver(idx1+1,idy1+1) -
b_ver(idx1+1,idy1) ) / dx
580 dbdy = coeff_x * dbdy2 + ( 1.d0 - coeff_x ) * dbdy1
584 IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r )
THEN 598 * sin( max_slope_angle_rad + velocity_ang_release_rad )
627 REAL*8,
INTENT(IN) :: x,y
629 REAL*8,
PARAMETER :: pig = 4.0*atan(1.0)
632 r = ( released_volume / pig )**(1.0/3.0)
634 IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r )
THEN 665 REAL*8,
INTENT(IN) :: x,y
667 REAL*8,
PARAMETER :: pig = 4.0*atan(1.0)
670 r = ( released_volume / pig )**(1.0/3.0)
672 IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r )
THEN real *8 dy
Control volumes size.
real *8, dimension(:), allocatable x_comp
Location of the centers (x) of the control volume of the domain.
real *8 y0
Bottom of the physical domain.
real *8, dimension(:,:), allocatable b_cent
Topography at the centers of the control volumes.
integer comp_cells_x
Number of control volumes x in the comp. domain.
real *8, dimension(:,:), allocatable thickness_init
real *8, dimension(:), allocatable y_comp
Location of the centers (y) of the control volume of the domain.
logical temperature_flag
Flag to choose if we model temperature transport.
real *8, dimension(:,:), allocatable source_xy
Array defining fraction of cells affected by source term.
real *8 t_init
Initial temperature of the pile of material.
real *8 alphas_e
Right sediment concentration.
real *8 u_e
Right velocity x.
real *8 v_w
Left velocity y.
real *8 dx
Control volumes size.
real *8 y_release
Initial y-coordinate of the pile.
integer comp_cells_y
Number of control volumes y in the comp. domain.
real *8 function sediment_function(x, y)
Sediment function.
real *8 riemann_interface
Riemann problem interface relative position. It is a value between 0 and 1.
real *8 u_w
Left velocity x.
real *8 function thickness_function(x, y, Bj)
Thickness function.
real *8 velocity_ang_release
Initial velocity direction (angle in degree): .
real *8 alphas_init
Initial sediment concentration in the pile of material.
real *8 velocity_mod_release
Initial velocity module of the pile.
real *8 released_volume
Initial volume of the flow.
integer n_vars
Number of conservative variables.
real *8 alphas_ambient
Ambient sediment concentration.
real *8 t_ambient
Ambient temperature.
subroutine initial_conditions
Problem initialization.
logical source_flag
Flag to choose if there is a source of mass within the domain.
real *8 x_release
Initial x-coordiante of the pile.
real *8 x0
Left of the physical domain.
real *8, dimension(:,:), allocatable b_ver
Topography at the vertices of the control volumes.
real *8 t_w
Left temperature.
real *8 v_e
Right velocity y.
real *8 function velocity_u_function(x, y)
Velocity u function.
real *8, dimension(:), allocatable x_stag
Location of the boundaries (x) of the control volumes of the domain.
real *8 function temperature_function(x, y)
Temperature function.
real *8 function velocity_v_function(x, y)
Velocity v function.
real *8 t_e
Right temperature.
subroutine qp_to_qc(qp, B, qc)
Physical to conservative variables.
subroutine riemann_problem
Riemann problem initialization.
real *8, dimension(:), allocatable y_stag
Location of the boundaries (x) of the control volumes of the domain.
real *8, dimension(:,:,:), allocatable q_init
real *8 alphas_w
Left sediment concentration.
logical solid_transport_flag
Flag to choose if we model solid phase transport.
real *8 hb_e
Right height.
subroutine init_source
Source initialization.
real *8, dimension(:,:,:), allocatable q
Conservative variables.