1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
program tsunami
! This version solves the linearized 1-d advection equation:
!
! du/dt + c du/dx = 0
implicit none
integer :: i, n
integer, parameter :: grid_size = 100 ! grid size in x
integer, parameter :: num_time_steps = 100 ! number of time steps
real, parameter :: dt = 1 ! time step [s]
real, parameter :: dx = 1 ! grid spacing [m]
real, parameter :: c = 1 ! background flow speed [m/s]
real :: h(grid_size), dh(grid_size)
integer, parameter :: icenter = 25
real, parameter :: decay = 0.02
character(*), parameter :: fmt = '(i0,*(1x,es15.8e2))'
! check input parameter values
if (grid_size <= 0) stop 'grid_size must be > 0'
if (dt <= 0) stop 'time step dt must be > 0'
if (dx <= 0) stop 'grid spacing dx must be > 0'
if (c <= 0) stop 'background flow speed c must be > 0'
! initialize water height to a Gaussian shape
do concurrent(i = 1:grid_size)
h(i) = exp(-decay * (i - icenter)**2)
end do
! write initial state to screen
print fmt, 0, h
time_loop: do n = 1, num_time_steps
! apply the periodic boundary condition
dh(1) = h(1) - h(grid_size)
! calculate the difference of u in space
do concurrent (i = 2:grid_size)
dh(i) = h(i) - h(i-1)
end do
! compute u at next time step
do concurrent (i = 1:grid_size)
h(i) = h(i) - c * dh(i) / dx * dt
end do
! write current state to screen
print fmt, n, h
end do time_loop
print *, "sum(h): ", sum(h)
if (abs(sum(h) - 12.5331345) > 1e-7) error stop
end program tsunami
|
