1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
4 subroutine blts ( ldmx, ldmy, ldmz,
9 > ist, iend, jst, jend,
12 c---------------------------------------------------------------------
13 c---------------------------------------------------------------------
15 c---------------------------------------------------------------------
17 c compute the regular-sparse, block lower 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, *),
33 > ldz( 5, 5, ldmx, ldmy),
34 > ldy( 5, 5, ldmx, ldmy),
35 > ldx( 5, 5, ldmx, ldmy),
36 > d( 5, 5, ldmx, ldmy)
42 c---------------------------------------------------------------------
44 c---------------------------------------------------------------------
47 double precision tmp, tmp1
48 double precision tmat(5,5)
51 c---------------------------------------------------------------------
52 c receive data from north and west
53 c---------------------------------------------------------------------
55 call exchange_1( v,k,iex )
62 v( m, i, j, k ) = v( m, i, j, k )
63 > - omega * ( ldz( m, 1, i, j ) * v( 1, i, j, k-1 )
64 > + ldz( m, 2, i, j ) * v( 2, i, j, k-1 )
65 > + ldz( m, 3, i, j ) * v( 3, i, j, k-1 )
66 > + ldz( m, 4, i, j ) * v( 4, i, j, k-1 )
67 > + ldz( m, 5, i, j ) * v( 5, i, j, k-1 ) )
79 v( m, i, j, k ) = v( m, i, j, k )
80 > - omega * ( ldy( m, 1, i, j ) * v( 1, i, j-1, k )
81 > + ldx( m, 1, i, j ) * v( 1, i-1, j, k )
82 > + ldy( m, 2, i, j ) * v( 2, i, j-1, k )
83 > + ldx( m, 2, i, j ) * v( 2, i-1, j, k )
84 > + ldy( m, 3, i, j ) * v( 3, i, j-1, k )
85 > + ldx( m, 3, i, j ) * v( 3, i-1, j, k )
86 > + ldy( m, 4, i, j ) * v( 4, i, j-1, k )
87 > + ldx( m, 4, i, j ) * v( 4, i-1, j, k )
88 > + ldy( m, 5, i, j ) * v( 5, i, j-1, k )
89 > + ldx( m, 5, i, j ) * v( 5, i-1, j, k ) )
93 c---------------------------------------------------------------------
94 c diagonal block inversion
97 c---------------------------------------------------------------------
99 tmat( m, 1 ) = d( m, 1, i, j )
100 tmat( m, 2 ) = d( m, 2, i, j )
101 tmat( m, 3 ) = d( m, 3, i, j )
102 tmat( m, 4 ) = d( m, 4, i, j )
103 tmat( m, 5 ) = d( m, 5, i, j )
106 tmp1 = 1.0d+00 / tmat( 1, 1 )
107 tmp = tmp1 * tmat( 2, 1 )
108 tmat( 2, 2 ) = tmat( 2, 2 )
109 > - tmp * tmat( 1, 2 )
110 tmat( 2, 3 ) = tmat( 2, 3 )
111 > - tmp * tmat( 1, 3 )
112 tmat( 2, 4 ) = tmat( 2, 4 )
113 > - tmp * tmat( 1, 4 )
114 tmat( 2, 5 ) = tmat( 2, 5 )
115 > - tmp * tmat( 1, 5 )
116 v( 2, i, j, k ) = v( 2, i, j, k )
117 > - v( 1, i, j, k ) * tmp
119 tmp = tmp1 * tmat( 3, 1 )
120 tmat( 3, 2 ) = tmat( 3, 2 )
121 > - tmp * tmat( 1, 2 )
122 tmat( 3, 3 ) = tmat( 3, 3 )
123 > - tmp * tmat( 1, 3 )
124 tmat( 3, 4 ) = tmat( 3, 4 )
125 > - tmp * tmat( 1, 4 )
126 tmat( 3, 5 ) = tmat( 3, 5 )
127 > - tmp * tmat( 1, 5 )
128 v( 3, i, j, k ) = v( 3, i, j, k )
129 > - v( 1, i, j, k ) * tmp
131 tmp = tmp1 * tmat( 4, 1 )
132 tmat( 4, 2 ) = tmat( 4, 2 )
133 > - tmp * tmat( 1, 2 )
134 tmat( 4, 3 ) = tmat( 4, 3 )
135 > - tmp * tmat( 1, 3 )
136 tmat( 4, 4 ) = tmat( 4, 4 )
137 > - tmp * tmat( 1, 4 )
138 tmat( 4, 5 ) = tmat( 4, 5 )
139 > - tmp * tmat( 1, 5 )
140 v( 4, i, j, k ) = v( 4, i, j, k )
141 > - v( 1, i, j, k ) * tmp
143 tmp = tmp1 * tmat( 5, 1 )
144 tmat( 5, 2 ) = tmat( 5, 2 )
145 > - tmp * tmat( 1, 2 )
146 tmat( 5, 3 ) = tmat( 5, 3 )
147 > - tmp * tmat( 1, 3 )
148 tmat( 5, 4 ) = tmat( 5, 4 )
149 > - tmp * tmat( 1, 4 )
150 tmat( 5, 5 ) = tmat( 5, 5 )
151 > - tmp * tmat( 1, 5 )
152 v( 5, i, j, k ) = v( 5, i, j, k )
153 > - v( 1, i, j, k ) * tmp
157 tmp1 = 1.0d+00 / tmat( 2, 2 )
158 tmp = tmp1 * tmat( 3, 2 )
159 tmat( 3, 3 ) = tmat( 3, 3 )
160 > - tmp * tmat( 2, 3 )
161 tmat( 3, 4 ) = tmat( 3, 4 )
162 > - tmp * tmat( 2, 4 )
163 tmat( 3, 5 ) = tmat( 3, 5 )
164 > - tmp * tmat( 2, 5 )
165 v( 3, i, j, k ) = v( 3, i, j, k )
166 > - v( 2, i, j, k ) * tmp
168 tmp = tmp1 * tmat( 4, 2 )
169 tmat( 4, 3 ) = tmat( 4, 3 )
170 > - tmp * tmat( 2, 3 )
171 tmat( 4, 4 ) = tmat( 4, 4 )
172 > - tmp * tmat( 2, 4 )
173 tmat( 4, 5 ) = tmat( 4, 5 )
174 > - tmp * tmat( 2, 5 )
175 v( 4, i, j, k ) = v( 4, i, j, k )
176 > - v( 2, i, j, k ) * tmp
178 tmp = tmp1 * tmat( 5, 2 )
179 tmat( 5, 3 ) = tmat( 5, 3 )
180 > - tmp * tmat( 2, 3 )
181 tmat( 5, 4 ) = tmat( 5, 4 )
182 > - tmp * tmat( 2, 4 )
183 tmat( 5, 5 ) = tmat( 5, 5 )
184 > - tmp * tmat( 2, 5 )
185 v( 5, i, j, k ) = v( 5, i, j, k )
186 > - v( 2, i, j, k ) * tmp
190 tmp1 = 1.0d+00 / tmat( 3, 3 )
191 tmp = tmp1 * tmat( 4, 3 )
192 tmat( 4, 4 ) = tmat( 4, 4 )
193 > - tmp * tmat( 3, 4 )
194 tmat( 4, 5 ) = tmat( 4, 5 )
195 > - tmp * tmat( 3, 5 )
196 v( 4, i, j, k ) = v( 4, i, j, k )
197 > - v( 3, i, j, k ) * tmp
199 tmp = tmp1 * tmat( 5, 3 )
200 tmat( 5, 4 ) = tmat( 5, 4 )
201 > - tmp * tmat( 3, 4 )
202 tmat( 5, 5 ) = tmat( 5, 5 )
203 > - tmp * tmat( 3, 5 )
204 v( 5, i, j, k ) = v( 5, i, j, k )
205 > - v( 3, i, j, k ) * tmp
209 tmp1 = 1.0d+00 / tmat( 4, 4 )
210 tmp = tmp1 * tmat( 5, 4 )
211 tmat( 5, 5 ) = tmat( 5, 5 )
212 > - tmp * tmat( 4, 5 )
213 v( 5, i, j, k ) = v( 5, i, j, k )
214 > - v( 4, i, j, k ) * tmp
216 c---------------------------------------------------------------------
218 c---------------------------------------------------------------------
219 v( 5, i, j, k ) = v( 5, i, j, k )
222 v( 4, i, j, k ) = v( 4, i, j, k )
223 > - tmat( 4, 5 ) * v( 5, i, j, k )
224 v( 4, i, j, k ) = v( 4, i, j, k )
227 v( 3, i, j, k ) = v( 3, i, j, k )
228 > - tmat( 3, 4 ) * v( 4, i, j, k )
229 > - tmat( 3, 5 ) * v( 5, i, j, k )
230 v( 3, i, j, k ) = v( 3, i, j, k )
233 v( 2, i, j, k ) = v( 2, i, j, k )
234 > - tmat( 2, 3 ) * v( 3, i, j, k )
235 > - tmat( 2, 4 ) * v( 4, i, j, k )
236 > - tmat( 2, 5 ) * v( 5, i, j, k )
237 v( 2, i, j, k ) = v( 2, i, j, k )
240 v( 1, i, j, k ) = v( 1, i, j, k )
241 > - tmat( 1, 2 ) * v( 2, i, j, k )
242 > - tmat( 1, 3 ) * v( 3, i, j, k )
243 > - tmat( 1, 4 ) * v( 4, i, j, k )
244 > - tmat( 1, 5 ) * v( 5, i, j, k )
245 v( 1, i, j, k ) = v( 1, i, j, k )
252 c---------------------------------------------------------------------
253 c send data to east and south
254 c---------------------------------------------------------------------
256 call exchange_1( v,k,iex )