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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
| c
real u(nxd,nyd,mn), ux(nxd,nyd,mn), uy(nxd,nyd,mn)
real f(nxd,nyd,mn), fx(nxd,nyd,mn)
real g(nxd,nyd,mn), gy(nxd,nyd,mn)
real v(nxd,nyd), du(nxd,2), df(nxd,2)
xmin(a,b) = 0.5*(sign(1.,a)+sign(1.,b))*min(abs(a),abs(b))
xmic(z,a,b) = xmin( z*xmin(a,b), 0.5*(a+b) )
gm1 = gamma - 1.0
tc = 0.0
istop = 0
do 1000 nt = 1, 10000
do 999 io = 0, 1
c* Periodic boundary condition in both x- & y-direction
do 101 m = 1, mn
do 100 i = 1, md
do 100 j = md + 1, ny + md
u(i, j,m) = u(nx+i,j,m)
u(nx + md+i,j,m) = u(md+i,j,m)
100 continue
do 101 j = 1, md
do 101 i = 1, nx + 2*md
u(i,j, m) = u(i,ny+j,m)
u(i,ny + md+j,m) = u(i,md+j,m)
101 continue
c* Compute f & g and maximum wave speeds "em_x, em_y".
c* See (2.1) & (4.3).
em_x = 1.e-15
em_y = 1.e-15
do 200 j = 1, ny + 2*md
do 200 i = 1, nx + 2*md
den = u(i,j,1)
vex = u(i,j,2) / den
vey = u(i,j,3) / den
eng = u(i,j,4)
pre = gm1 * ( eng - .5*den*( vex*vex + vey*vey ) )
cvel = sqrt( gamma * pre / den )
em_x = max( em_x, abs(vex) + cvel )
em_y = max( em_y, abs(vey) + cvel )
f(i,j,1) = den * vex
f(i,j,2) = den * vex**2 + pre
f(i,j,3) = den * vex * vey
f(i,j,4) = vex * ( pre + eng )
g(i,j,1) = den * vey
g(i,j,2) = den * vex * vey
g(i,j,3) = den * vey**2 + pre
g(i,j,4) = vey * ( pre + eng )
200 continue
c* Compute numerical derivatives "ux", "uy", "fx", "gy".
c*cSee (3.1) & (4.1)
do 330 m = 1, mn
do 310 j = 3, ny + 2*md - 2
do 301 i = 1, nx + 2*md - 1
du(i,1) = u(i+1,j,m) - u(i,j,m)
301 df(i,1) = f(i+1,j,m) - f(i,j,m)
do 302 i = 1, nx + 2*md - 2
du(i,2) = du(i+1,1) - du(i,1)
302 df(i,2) = df(i+1,1) - df(i,1)
if( theta .lt. 2.5 ) then
do 303 i = 3, nx + 2*md - 2
ux(i,j,m) = xmic( theta, du(i-1,1), du(i,1) )
303 fx(i,j,m) = xmic( theta, df(i-1,1), df(i,1) )
else
do 304 i = 3, nx + 2*md - 2
ux(i,j,m)=xmin(du(i-1,1)+.5*xmin(du(i-2,2),du(i-1,2)),
& du(i, 1)-.5*xmin(du(i-1,2),du(i, 2)))
fx(i,j,m)=xmin(df(i-1,1)+.5*xmin(df(i-2,2),df(i-1,2)),
& df(i, 1)-.5*xmin(df(i-1,2),df(i, 2)))
304 continue
endif
310 continue
do 320 i = 3, nx + 2*md - 2
do 311 j = 1, ny + 2*md - 1
du(j,1) = u(i,j+1,m) - u(i,j,m)
311 df(j,1) = g(i,j+1,m) - g(i,j,m)
do 312 j = 1, ny + 2*md - 2
du(j,2) = du(j+1,1) - du(j,1)
312 df(j,2) = df(j+1,1) - df(j,1)
if( theta .lt. 2.5 ) then
do 313 j = 3, ny + 2*md - 2
uy(i,j,m) = xmic( theta, du(j-1,1), du(j,1) )
313 gy(i,j,m) = xmic( theta, df(j-1,1), df(j,1) )
else
do 314 j = 3, ny + 2*md - 2
uy(i,j,m)=xmin(du(j-1,1)+.5*xmin(du(j-2,2),du(j-1,2)),
& du(j, 1)-.5*xmin(du(j-1,2),du(j, 2)))
gy(i,j,m)=xmin(df(j-1,1)+.5*xmin(df(j-2,2),df(j-1,2)),
& df(j, 1)-.5*xmin(df(j-1,2),df(j, 2)))
314 continue
endif
320 continue
330 continue
c* Compute time step size according to the input CFL #.
if(io.eq.0) then
dt = cfl / max( em_x/dx, em_y/dy )
if( ( tc + 2.*dt ) .ge. tf ) then
dt = 0.5 * ( tf - tc )
istop = 1
endif
endif
dtcdx2 = 0.5 * dt / dx
dtcdy2 = 0.5 * dt / dy
c* Compute the flux values of f & g at half time step.
c* See (2.15) & (2.16).
do 400 j = 3, ny + 2*md - 2
do 400 i = 3, nx + 2*md - 2
den = u(i,j,1) - dtcdx2*fx(i,j,1) - dtcdy2*gy(i,j,1)
xmt = u(i,j,2) - dtcdx2*fx(i,j,2) - dtcdy2*gy(i,j,2)
ymt = u(i,j,3) - dtcdx2*fx(i,j,3) - dtcdy2*gy(i,j,3)
eng = u(i,j,4) - dtcdx2*fx(i,j,4) - dtcdy2*gy(i,j,4)
pre = gm1 * ( eng - .5 * ( xmt*xmt + ymt*ymt ) / den )
f(i,j,1) = xmt
f(i,j,2) = xmt * xmt / den + pre
f(i,j,3) = xmt * ymt / den
f(i,j,4) = xmt / den * ( pre + eng )
g(i,j,1) = ymt
g(i,j,2) = xmt * ymt / den
g(i,j,3) = ymt * ymt / den + pre
g(i,j,4) = ymt / den * ( pre + eng )
400 continue
c* Compute the values of "u" at the next time level. See (2.16).
do 510 m = 1, mn
do 501 j = md + 1 - io, ny + md - io
do 501 i = md + 1 - io, nx + md - io
v(i,j) = 0.25 * ( u(i,j, m) + u(i+1,j, m)
& + u(i,j+1,m) + u(i+1,j+1,m) )
& + 0.0625 * ( ux(i,j, m) - ux(i+1,j, m)
& + ux(i,j+1,m) - ux(i+1,j+1,m)
& + uy(i,j, m) + uy(i+1,j, m)
& - uy(i,j+1,m) - uy(i+1,j+1,m) )
& + dtcdx2 * ( f(i,j, m) - f(i+1,j, m)
& + f(i,j+1,m) - f(i+1,j+1,m) )
& + dtcdy2 * ( g(i,j, m) + g(i+1,j, m)
& - g(i,j+1,m) - g(i+1,j+1,m) )
501 continue
do 502 j = md + 1, ny + md
do 502 i = md + 1, nx + md
502 u(i,j,m) = v(i-io,j-io)
510 continue
tc = tc + dt
999 continue
if(istop.eq.1 ) goto 1001
1000 continue
1001 return
end |
