set_constants.f
4.28 KB
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
c---------------------------------------------------------------------
c---------------------------------------------------------------------
subroutine set_constants
c---------------------------------------------------------------------
c---------------------------------------------------------------------
include 'header.h'
ce(1,1) = 2.0d0
ce(1,2) = 0.0d0
ce(1,3) = 0.0d0
ce(1,4) = 4.0d0
ce(1,5) = 5.0d0
ce(1,6) = 3.0d0
ce(1,7) = 0.5d0
ce(1,8) = 0.02d0
ce(1,9) = 0.01d0
ce(1,10) = 0.03d0
ce(1,11) = 0.5d0
ce(1,12) = 0.4d0
ce(1,13) = 0.3d0
ce(2,1) = 1.0d0
ce(2,2) = 0.0d0
ce(2,3) = 0.0d0
ce(2,4) = 0.0d0
ce(2,5) = 1.0d0
ce(2,6) = 2.0d0
ce(2,7) = 3.0d0
ce(2,8) = 0.01d0
ce(2,9) = 0.03d0
ce(2,10) = 0.02d0
ce(2,11) = 0.4d0
ce(2,12) = 0.3d0
ce(2,13) = 0.5d0
ce(3,1) = 2.0d0
ce(3,2) = 2.0d0
ce(3,3) = 0.0d0
ce(3,4) = 0.0d0
ce(3,5) = 0.0d0
ce(3,6) = 2.0d0
ce(3,7) = 3.0d0
ce(3,8) = 0.04d0
ce(3,9) = 0.03d0
ce(3,10) = 0.05d0
ce(3,11) = 0.3d0
ce(3,12) = 0.5d0
ce(3,13) = 0.4d0
ce(4,1) = 2.0d0
ce(4,2) = 2.0d0
ce(4,3) = 0.0d0
ce(4,4) = 0.0d0
ce(4,5) = 0.0d0
ce(4,6) = 2.0d0
ce(4,7) = 3.0d0
ce(4,8) = 0.03d0
ce(4,9) = 0.05d0
ce(4,10) = 0.04d0
ce(4,11) = 0.2d0
ce(4,12) = 0.1d0
ce(4,13) = 0.3d0
ce(5,1) = 5.0d0
ce(5,2) = 4.0d0
ce(5,3) = 3.0d0
ce(5,4) = 2.0d0
ce(5,5) = 0.1d0
ce(5,6) = 0.4d0
ce(5,7) = 0.3d0
ce(5,8) = 0.05d0
ce(5,9) = 0.04d0
ce(5,10) = 0.03d0
ce(5,11) = 0.1d0
ce(5,12) = 0.3d0
ce(5,13) = 0.2d0
c1 = 1.4d0
c2 = 0.4d0
c3 = 0.1d0
c4 = 1.0d0
c5 = 1.4d0
bt = dsqrt(0.5d0)
dnxm1 = 1.0d0 / dble(grid_points(1)-1)
dnym1 = 1.0d0 / dble(grid_points(2)-1)
dnzm1 = 1.0d0 / dble(grid_points(3)-1)
c1c2 = c1 * c2
c1c5 = c1 * c5
c3c4 = c3 * c4
c1345 = c1c5 * c3c4
conz1 = (1.0d0-c1c5)
tx1 = 1.0d0 / (dnxm1 * dnxm1)
tx2 = 1.0d0 / (2.0d0 * dnxm1)
tx3 = 1.0d0 / dnxm1
ty1 = 1.0d0 / (dnym1 * dnym1)
ty2 = 1.0d0 / (2.0d0 * dnym1)
ty3 = 1.0d0 / dnym1
tz1 = 1.0d0 / (dnzm1 * dnzm1)
tz2 = 1.0d0 / (2.0d0 * dnzm1)
tz3 = 1.0d0 / dnzm1
dx1 = 0.75d0
dx2 = 0.75d0
dx3 = 0.75d0
dx4 = 0.75d0
dx5 = 0.75d0
dy1 = 0.75d0
dy2 = 0.75d0
dy3 = 0.75d0
dy4 = 0.75d0
dy5 = 0.75d0
dz1 = 1.0d0
dz2 = 1.0d0
dz3 = 1.0d0
dz4 = 1.0d0
dz5 = 1.0d0
dxmax = dmax1(dx3, dx4)
dymax = dmax1(dy2, dy4)
dzmax = dmax1(dz2, dz3)
dssp = 0.25d0 * dmax1(dx1, dmax1(dy1, dz1) )
c4dssp = 4.0d0 * dssp
c5dssp = 5.0d0 * dssp
dttx1 = dt*tx1
dttx2 = dt*tx2
dtty1 = dt*ty1
dtty2 = dt*ty2
dttz1 = dt*tz1
dttz2 = dt*tz2
c2dttx1 = 2.0d0*dttx1
c2dtty1 = 2.0d0*dtty1
c2dttz1 = 2.0d0*dttz1
dtdssp = dt*dssp
comz1 = dtdssp
comz4 = 4.0d0*dtdssp
comz5 = 5.0d0*dtdssp
comz6 = 6.0d0*dtdssp
c3c4tx3 = c3c4*tx3
c3c4ty3 = c3c4*ty3
c3c4tz3 = c3c4*tz3
dx1tx1 = dx1*tx1
dx2tx1 = dx2*tx1
dx3tx1 = dx3*tx1
dx4tx1 = dx4*tx1
dx5tx1 = dx5*tx1
dy1ty1 = dy1*ty1
dy2ty1 = dy2*ty1
dy3ty1 = dy3*ty1
dy4ty1 = dy4*ty1
dy5ty1 = dy5*ty1
dz1tz1 = dz1*tz1
dz2tz1 = dz2*tz1
dz3tz1 = dz3*tz1
dz4tz1 = dz4*tz1
dz5tz1 = dz5*tz1
c2iv = 2.5d0
con43 = 4.0d0/3.0d0
con16 = 1.0d0/6.0d0
xxcon1 = c3c4tx3*con43*tx3
xxcon2 = c3c4tx3*tx3
xxcon3 = c3c4tx3*conz1*tx3
xxcon4 = c3c4tx3*con16*tx3
xxcon5 = c3c4tx3*c1c5*tx3
yycon1 = c3c4ty3*con43*ty3
yycon2 = c3c4ty3*ty3
yycon3 = c3c4ty3*conz1*ty3
yycon4 = c3c4ty3*con16*ty3
yycon5 = c3c4ty3*c1c5*ty3
zzcon1 = c3c4tz3*con43*tz3
zzcon2 = c3c4tz3*tz3
zzcon3 = c3c4tz3*conz1*tz3
zzcon4 = c3c4tz3*con16*tz3
zzcon5 = c3c4tz3*c1c5*tz3
return
end