1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
4 subroutine buts( ldmx, ldmy, ldmz,
9 > ist, iend, jst, jend,
10 > nx0, ny0, ipt, jpt )
12 c---------------------------------------------------------------------
13 c---------------------------------------------------------------------
15 c---------------------------------------------------------------------
17 c compute the regular-sparse, block upper triangular solution:
21 c---------------------------------------------------------------------
25 c---------------------------------------------------------------------
27 c---------------------------------------------------------------------
28 integer ldmx, ldmy, ldmz
31 double precision omega
32 double precision v( 5, -1:ldmx+2, -1:ldmy+2, *),
34 > d( 5, 5, ldmx, ldmy),
35 > udx( 5, 5, ldmx, ldmy),
36 > udy( 5, 5, ldmx, ldmy),
37 > udz( 5, 5, ldmx, ldmy )
43 c---------------------------------------------------------------------
45 c---------------------------------------------------------------------
48 double precision tmp, tmp1
49 double precision tmat(5,5)
52 c---------------------------------------------------------------------
53 c receive data from south and east
54 c---------------------------------------------------------------------
56 call exchange_1( v,k,iex )
62 > omega * ( udz( m, 1, i, j ) * v( 1, i, j, k+1 )
63 > + udz( m, 2, i, j ) * v( 2, i, j, k+1 )
64 > + udz( m, 3, i, j ) * v( 3, i, j, k+1 )
65 > + udz( m, 4, i, j ) * v( 4, i, j, k+1 )
66 > + udz( m, 5, i, j ) * v( 5, i, j, k+1 ) )
76 tv( m, i, j ) = tv( m, i, j )
77 > + omega * ( udy( m, 1, i, j ) * v( 1, i, j+1, k )
78 > + udx( m, 1, i, j ) * v( 1, i+1, j, k )
79 > + udy( m, 2, i, j ) * v( 2, i, j+1, k )
80 > + udx( m, 2, i, j ) * v( 2, i+1, j, k )
81 > + udy( m, 3, i, j ) * v( 3, i, j+1, k )
82 > + udx( m, 3, i, j ) * v( 3, i+1, j, k )
83 > + udy( m, 4, i, j ) * v( 4, i, j+1, k )
84 > + udx( m, 4, i, j ) * v( 4, i+1, j, k )
85 > + udy( m, 5, i, j ) * v( 5, i, j+1, k )
86 > + udx( m, 5, i, j ) * v( 5, i+1, j, k ) )
89 c---------------------------------------------------------------------
90 c diagonal block inversion
91 c---------------------------------------------------------------------
93 tmat( m, 1 ) = d( m, 1, i, j )
94 tmat( m, 2 ) = d( m, 2, i, j )
95 tmat( m, 3 ) = d( m, 3, i, j )
96 tmat( m, 4 ) = d( m, 4, i, j )
97 tmat( m, 5 ) = d( m, 5, i, j )
100 tmp1 = 1.0d+00 / tmat( 1, 1 )
101 tmp = tmp1 * tmat( 2, 1 )
102 tmat( 2, 2 ) = tmat( 2, 2 )
103 > - tmp * tmat( 1, 2 )
104 tmat( 2, 3 ) = tmat( 2, 3 )
105 > - tmp * tmat( 1, 3 )
106 tmat( 2, 4 ) = tmat( 2, 4 )
107 > - tmp * tmat( 1, 4 )
108 tmat( 2, 5 ) = tmat( 2, 5 )
109 > - tmp * tmat( 1, 5 )
110 tv( 2, i, j ) = tv( 2, i, j )
111 > - tv( 1, i, j ) * tmp
113 tmp = tmp1 * tmat( 3, 1 )
114 tmat( 3, 2 ) = tmat( 3, 2 )
115 > - tmp * tmat( 1, 2 )
116 tmat( 3, 3 ) = tmat( 3, 3 )
117 > - tmp * tmat( 1, 3 )
118 tmat( 3, 4 ) = tmat( 3, 4 )
119 > - tmp * tmat( 1, 4 )
120 tmat( 3, 5 ) = tmat( 3, 5 )
121 > - tmp * tmat( 1, 5 )
122 tv( 3, i, j ) = tv( 3, i, j )
123 > - tv( 1, i, j ) * tmp
125 tmp = tmp1 * tmat( 4, 1 )
126 tmat( 4, 2 ) = tmat( 4, 2 )
127 > - tmp * tmat( 1, 2 )
128 tmat( 4, 3 ) = tmat( 4, 3 )
129 > - tmp * tmat( 1, 3 )
130 tmat( 4, 4 ) = tmat( 4, 4 )
131 > - tmp * tmat( 1, 4 )
132 tmat( 4, 5 ) = tmat( 4, 5 )
133 > - tmp * tmat( 1, 5 )
134 tv( 4, i, j ) = tv( 4, i, j )
135 > - tv( 1, i, j ) * tmp
137 tmp = tmp1 * tmat( 5, 1 )
138 tmat( 5, 2 ) = tmat( 5, 2 )
139 > - tmp * tmat( 1, 2 )
140 tmat( 5, 3 ) = tmat( 5, 3 )
141 > - tmp * tmat( 1, 3 )
142 tmat( 5, 4 ) = tmat( 5, 4 )
143 > - tmp * tmat( 1, 4 )
144 tmat( 5, 5 ) = tmat( 5, 5 )
145 > - tmp * tmat( 1, 5 )
146 tv( 5, i, j ) = tv( 5, i, j )
147 > - tv( 1, i, j ) * tmp
151 tmp1 = 1.0d+00 / tmat( 2, 2 )
152 tmp = tmp1 * tmat( 3, 2 )
153 tmat( 3, 3 ) = tmat( 3, 3 )
154 > - tmp * tmat( 2, 3 )
155 tmat( 3, 4 ) = tmat( 3, 4 )
156 > - tmp * tmat( 2, 4 )
157 tmat( 3, 5 ) = tmat( 3, 5 )
158 > - tmp * tmat( 2, 5 )
159 tv( 3, i, j ) = tv( 3, i, j )
160 > - tv( 2, i, j ) * tmp
162 tmp = tmp1 * tmat( 4, 2 )
163 tmat( 4, 3 ) = tmat( 4, 3 )
164 > - tmp * tmat( 2, 3 )
165 tmat( 4, 4 ) = tmat( 4, 4 )
166 > - tmp * tmat( 2, 4 )
167 tmat( 4, 5 ) = tmat( 4, 5 )
168 > - tmp * tmat( 2, 5 )
169 tv( 4, i, j ) = tv( 4, i, j )
170 > - tv( 2, i, j ) * tmp
172 tmp = tmp1 * tmat( 5, 2 )
173 tmat( 5, 3 ) = tmat( 5, 3 )
174 > - tmp * tmat( 2, 3 )
175 tmat( 5, 4 ) = tmat( 5, 4 )
176 > - tmp * tmat( 2, 4 )
177 tmat( 5, 5 ) = tmat( 5, 5 )
178 > - tmp * tmat( 2, 5 )
179 tv( 5, i, j ) = tv( 5, i, j )
180 > - tv( 2, i, j ) * tmp
184 tmp1 = 1.0d+00 / tmat( 3, 3 )
185 tmp = tmp1 * tmat( 4, 3 )
186 tmat( 4, 4 ) = tmat( 4, 4 )
187 > - tmp * tmat( 3, 4 )
188 tmat( 4, 5 ) = tmat( 4, 5 )
189 > - tmp * tmat( 3, 5 )
190 tv( 4, i, j ) = tv( 4, i, j )
191 > - tv( 3, i, j ) * tmp
193 tmp = tmp1 * tmat( 5, 3 )
194 tmat( 5, 4 ) = tmat( 5, 4 )
195 > - tmp * tmat( 3, 4 )
196 tmat( 5, 5 ) = tmat( 5, 5 )
197 > - tmp * tmat( 3, 5 )
198 tv( 5, i, j ) = tv( 5, i, j )
199 > - tv( 3, i, j ) * tmp
203 tmp1 = 1.0d+00 / tmat( 4, 4 )
204 tmp = tmp1 * tmat( 5, 4 )
205 tmat( 5, 5 ) = tmat( 5, 5 )
206 > - tmp * tmat( 4, 5 )
207 tv( 5, i, j ) = tv( 5, i, j )
208 > - tv( 4, i, j ) * tmp
210 c---------------------------------------------------------------------
212 c---------------------------------------------------------------------
213 tv( 5, i, j ) = tv( 5, i, j )
216 tv( 4, i, j ) = tv( 4, i, j )
217 > - tmat( 4, 5 ) * tv( 5, i, j )
218 tv( 4, i, j ) = tv( 4, i, j )
221 tv( 3, i, j ) = tv( 3, i, j )
222 > - tmat( 3, 4 ) * tv( 4, i, j )
223 > - tmat( 3, 5 ) * tv( 5, i, j )
224 tv( 3, i, j ) = tv( 3, i, j )
227 tv( 2, i, j ) = tv( 2, i, j )
228 > - tmat( 2, 3 ) * tv( 3, i, j )
229 > - tmat( 2, 4 ) * tv( 4, i, j )
230 > - tmat( 2, 5 ) * tv( 5, i, j )
231 tv( 2, i, j ) = tv( 2, i, j )
234 tv( 1, i, j ) = tv( 1, i, j )
235 > - tmat( 1, 2 ) * tv( 2, i, j )
236 > - tmat( 1, 3 ) * tv( 3, i, j )
237 > - tmat( 1, 4 ) * tv( 4, i, j )
238 > - tmat( 1, 5 ) * tv( 5, i, j )
239 tv( 1, i, j ) = tv( 1, i, j )
242 v( 1, i, j, k ) = v( 1, i, j, k ) - tv( 1, i, j )
243 v( 2, i, j, k ) = v( 2, i, j, k ) - tv( 2, i, j )
244 v( 3, i, j, k ) = v( 3, i, j, k ) - tv( 3, i, j )
245 v( 4, i, j, k ) = v( 4, i, j, k ) - tv( 4, i, j )
246 v( 5, i, j, k ) = v( 5, i, j, k ) - tv( 5, i, j )
252 c---------------------------------------------------------------------
253 c send data to north and west
254 c---------------------------------------------------------------------
256 call exchange_1( v,k,iex )