-
Notifications
You must be signed in to change notification settings - Fork 1
/
initial.f
107 lines (106 loc) · 3.01 KB
/
initial.f
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
subroutine initial(U)
implicit real*8 (a-h,m-z)
parameter (lmx=400,lmy=400)
common /gridI/lx,ly,lsteps
common /gridR/T,dt,dx,dy,xlambda,ylambda
common /conf/config
dimension U(4,-1:lmx+3,-1:lmy+3)
real*8, dimension(4) :: vec1,vec2,vec3,vec4
integer :: config
lx_half=lx/2
ly_half=ly/2
c-----Configuration 3--------------------------------------------------
c
if (config .eq. 3) then
c-----Initial data-----------------------------------------------------
c
c-----First quadrant---------------------------------------------------
rho1 = 1.5d0
p1 = 1.5d0
u1 = 0.0d0
v1 = 0.0d0
e1 = total_energy(rho1,u1,v1,p1)
vec1 = (/ rho1, rho1*u1, rho1*v1, e1/)
c-----Second quadrant--------------------------------------------------
rho2 = 0.5323d0
p2 = 0.3d0
u2 = 1.206d0
v2 = 0.0d0
e2 = total_energy(rho2,u2,v2,p2)
vec2 = (/ rho2, rho2*u2, rho2*v2, e2/)
c-----Third quadrant---------------------------------------------------
rho3 = 0.138d0
p3 = 0.029d0
u3 = 1.206d0
v3 = 1.206d0
e3 = total_energy(rho3,u3,v3,p3)
vec3 = (/ rho3, rho3*u3, rho3*v3, e3/)
c-----Fourth quadrant--------------------------------------------------
rho4 = 0.5323d0
p4 = 0.3d0
u4 = 0.0d0
v4 = 1.206d0
e4 = total_energy(rho4,u4,v4,p4)
vec4 = (/ rho4, rho4*u4, rho4*v4, e4/)
endif
c
if (config .eq. 5) then
c-----Initial data-----------------------------------------------------
c
c-----First quadrant---------------------------------------------------
rho1 = 1.0d0
p1 = 1.0d0
u1 = -0.75d0
v1 = -0.50d0
e1 = total_energy(rho1,u1,v1,p1)
vec1 = (/ rho1, rho1*u1, rho1*v1, e1/)
c-----Second quadrant--------------------------------------------------
rho2 = 2.0d0
p2 = 1.0d0
u2 = -0.75d0
v2 = 0.50d0
e2 = total_energy(rho2,u2,v2,p2)
vec2 = (/ rho2, rho2*u2, rho2*v2, e2/)
c-----Third quadrant---------------------------------------------------
rho3 = 1.0d0
p3 = 1.0d0
u3 = 0.75d0
v3 = 0.50d0
e3 = total_energy(rho3,u3,v3,p3)
vec3 = (/ rho3, rho3*u3, rho3*v3, e3/)
c-----Fourth quadrant--------------------------------------------------
rho4 = 3.0d0
p4 = 1.0d0
u4 = 0.75d0
v4 = -0.50d0
e4 = total_energy(rho4,u4,v4,p4)
vec4 = (/ rho4, rho4*u4, rho4*v4, e4/)
endif
c
c-----Initial solution-------------------------------------------------
c
c-----First quadrant---------------------------------------------------
do j = ly_half+1, ly+3
do i = lx_half+1, lx+3
U(1:4,i,j) = vec1(1:4)
enddo
enddo
c-----Second quadrant--------------------------------------------------
do j = ly_half+1, ly+3
do i = -1, lx_half
U(1:4,i,j) = vec2(1:4)
enddo
enddo
c-----Third quadrant---------------------------------------------------
do j = -1, ly_half
do i = -1, lx_half
U(1:4,i,j) = vec3(1:4)
enddo
enddo
c-----Fourth quadrant--------------------------------------------------
do j = -1, ly_half
do i = lx_half+1, lx+3
U(1:4,i,j) = vec4(1:4)
enddo
enddo
end