1 |
C $Header: /u/gcmpack/MITgcm/model/src/solve_tridiagonal.F,v 1.10 2012/03/15 15:25:18 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "CPP_OPTIONS.h" |
5 |
|
6 |
CBOP |
7 |
C !ROUTINE: SOLVE_TRIDIAGONAL |
8 |
C !INTERFACE: |
9 |
SUBROUTINE SOLVE_TRIDIAGONAL( |
10 |
I iMin,iMax, jMin,jMax, |
11 |
I a3d, b3d, c3d, |
12 |
U y3d, |
13 |
O errCode, |
14 |
I bi, bj, myThid ) |
15 |
C !DESCRIPTION: \bv |
16 |
C *==========================================================* |
17 |
C | S/R SOLVE_TRIDIAGONAL |
18 |
C | o Solve a tri-diagonal system A*X=Y (dimension Nr) |
19 |
C *==========================================================* |
20 |
C | o Used to solve implicitly vertical advection & diffusion |
21 |
C *==========================================================* |
22 |
C \ev |
23 |
|
24 |
C !USES: |
25 |
IMPLICIT NONE |
26 |
C == Global data == |
27 |
#include "SIZE.h" |
28 |
#include "EEPARAMS.h" |
29 |
|
30 |
C !INPUT/OUTPUT PARAMETERS: |
31 |
C == Routine Arguments == |
32 |
C a3d :: matrix lower diagnonal |
33 |
C b3d :: matrix main diagnonal |
34 |
C c3d :: matrix upper diagnonal |
35 |
C y3d :: Input = Y vector ; Output = X = solution of A*X=Y |
36 |
C errCode :: > 0 if singular matrix |
37 |
INTEGER iMin,iMax,jMin,jMax |
38 |
_RL a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
39 |
_RL b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
40 |
_RL c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
41 |
_RL y3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
42 |
INTEGER errCode |
43 |
INTEGER bi, bj, myThid |
44 |
|
45 |
C !LOCAL VARIABLES: |
46 |
C == Local variables == |
47 |
INTEGER i,j,k |
48 |
_RL tmpVar |
49 |
#ifndef SOLVE_DIAGONAL_LOWMEMORY |
50 |
_RL recVar |
51 |
_RL c3d_m1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
52 |
_RL y3d_m1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
53 |
# ifdef SOLVE_DIAGONAL_KINNER |
54 |
_RL c3d_prime(Nr) |
55 |
_RL y3d_prime(Nr) |
56 |
_RL y3d_update(Nr) |
57 |
# else |
58 |
_RL c3d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
59 |
_RL y3d_prime(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
60 |
_RL y3d_update(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
61 |
# endif |
62 |
#endif /* SOLVE_DIAGONAL_LOWMEMORY */ |
63 |
CEOP |
64 |
|
65 |
errCode = 0 |
66 |
|
67 |
#ifdef SOLVE_DIAGONAL_LOWMEMORY |
68 |
|
69 |
C-- Beginning of forward sweep (top level) |
70 |
DO j=jMin,jMax |
71 |
DO i=iMin,iMax |
72 |
IF ( b3d(i,j,1).NE.0. _d 0 ) THEN |
73 |
b3d(i,j,1) = 1. _d 0 / b3d(i,j,1) |
74 |
ELSE |
75 |
b3d(i,j,1) = 0. _d 0 |
76 |
errCode = 1 |
77 |
ENDIF |
78 |
ENDDO |
79 |
ENDDO |
80 |
|
81 |
C-- Middle of forward sweep |
82 |
DO k=2,Nr |
83 |
DO j=jMin,jMax |
84 |
DO i=iMin,iMax |
85 |
tmpVar = b3d(i,j,k) - a3d(i,j,k)*c3d(i,j,k-1)*b3d(i,j,k-1) |
86 |
IF ( tmpVar .NE. 0. _d 0 ) THEN |
87 |
b3d(i,j,k) = 1. _d 0 / tmpVar |
88 |
ELSE |
89 |
b3d(i,j,k) = 0. _d 0 |
90 |
errCode = 1 |
91 |
ENDIF |
92 |
ENDDO |
93 |
ENDDO |
94 |
ENDDO |
95 |
|
96 |
DO j=jMin,jMax |
97 |
DO i=iMin,iMax |
98 |
y3d(i,j,1) = y3d(i,j,1)*b3d(i,j,1) |
99 |
ENDDO |
100 |
ENDDO |
101 |
DO k=2,Nr |
102 |
DO j=jMin,jMax |
103 |
DO i=iMin,iMax |
104 |
y3d(i,j,k) = ( y3d(i,j,k) |
105 |
& - a3d(i,j,k)*y3d(i,j,k-1) |
106 |
& )*b3d(i,j,k) |
107 |
ENDDO |
108 |
ENDDO |
109 |
ENDDO |
110 |
|
111 |
C-- Backward sweep |
112 |
DO k=Nr-1,1,-1 |
113 |
DO j=jMin,jMax |
114 |
DO i=iMin,iMax |
115 |
y3d(i,j,k) = y3d(i,j,k) |
116 |
& - c3d(i,j,k)*b3d(i,j,k)*y3d(i,j,k+1) |
117 |
ENDDO |
118 |
ENDDO |
119 |
ENDDO |
120 |
|
121 |
#else /* ndef SOLVE_DIAGONAL_LOWMEMORY */ |
122 |
|
123 |
C-- Temporary array |
124 |
DO k=1,Nr |
125 |
DO j=1-OLy,sNy+OLy |
126 |
DO i=1-OLx,sNx+OLx |
127 |
c DO j=jMin,jMax |
128 |
c DO i=iMin,iMax |
129 |
c3d_m1(i,j,k) = c3d(i,j,k) |
130 |
y3d_m1(i,j,k) = y3d(i,j,k) |
131 |
ENDDO |
132 |
ENDDO |
133 |
ENDDO |
134 |
|
135 |
#ifdef SOLVE_DIAGONAL_KINNER |
136 |
|
137 |
C-- Main loop |
138 |
DO j=1-OLy,sNy+OLy |
139 |
DO i=1-OLx,sNx+OLx |
140 |
|
141 |
DO k=1,Nr |
142 |
c3d_prime(k) = 0. _d 0 |
143 |
y3d_prime(k) = 0. _d 0 |
144 |
y3d_update(k) = 0. _d 0 |
145 |
ENDDO |
146 |
|
147 |
C-- Forward sweep |
148 |
DO k=1,Nr |
149 |
IF ( k.EQ.1 ) THEN |
150 |
IF ( b3d(i,j,1).NE.0. _d 0 ) THEN |
151 |
c3d_prime(1) = c3d_m1(i,j,1) / b3d(i,j,1) |
152 |
y3d_prime(1) = y3d_m1(i,j,1) / b3d(i,j,1) |
153 |
ELSE |
154 |
c3d_prime(1) = 0. _d 0 |
155 |
y3d_prime(1) = 0. _d 0 |
156 |
errCode = 1 |
157 |
ENDIF |
158 |
ELSE |
159 |
tmpVar = b3d(i,j,k) - a3d(i,j,k)*c3d_prime(k-1) |
160 |
IF ( tmpVar .NE. 0. _d 0 ) THEN |
161 |
recVar = 1. _d 0 / tmpVar |
162 |
c3d_prime(k) = c3d_m1(i,j,k)*recVar |
163 |
y3d_prime(k) = (y3d_m1(i,j,k) - y3d_prime(k-1)*a3d(i,j,k)) |
164 |
& *recVar |
165 |
ELSE |
166 |
c3d_prime(k) = 0. _d 0 |
167 |
y3d_prime(k) = 0. _d 0 |
168 |
errCode = 1 |
169 |
ENDIF |
170 |
ENDIF |
171 |
ENDDO |
172 |
|
173 |
C-- Backward sweep |
174 |
DO k=Nr,1,-1 |
175 |
IF ( k.EQ.Nr ) THEN |
176 |
y3d_update(k)=y3d_prime(k) |
177 |
ELSE |
178 |
y3d_update(k)=y3d_prime(k)-c3d_prime(k)*y3d_update(k+1) |
179 |
ENDIF |
180 |
ENDDO |
181 |
|
182 |
C-- Update array |
183 |
DO k=1,Nr |
184 |
y3d(i,j,k) = y3d_update(k) |
185 |
ENDDO |
186 |
|
187 |
ENDDO |
188 |
ENDDO |
189 |
|
190 |
#else /* ndef SOLVE_DIAGONAL_KINNER */ |
191 |
|
192 |
C-- Init. |
193 |
DO k=1,Nr |
194 |
DO j=1-OLy,sNy+OLy |
195 |
DO i=1-OLx,sNx+OLx |
196 |
c DO j=jMin,jMax |
197 |
c DO i=iMin,iMax |
198 |
c3d_prime(i,j,k) = 0. _d 0 |
199 |
y3d_prime(i,j,k) = 0. _d 0 |
200 |
y3d_update(i,j,k) = 0. _d 0 |
201 |
ENDDO |
202 |
ENDDO |
203 |
ENDDO |
204 |
|
205 |
CADJ loop = sequential |
206 |
C-- Forward sweep |
207 |
DO k=1,Nr |
208 |
|
209 |
DO j=1-OLy,sNy+OLy |
210 |
DO i=1-OLx,sNx+OLx |
211 |
c DO j=jMin,jMax |
212 |
c DO i=iMin,iMax |
213 |
IF ( k.EQ.1 ) THEN |
214 |
IF ( b3d(i,j,1).NE.0. _d 0 ) THEN |
215 |
recVar = 1. _d 0 / b3d(i,j,1) |
216 |
c3d_prime(i,j,1) = c3d_m1(i,j,1)*recVar |
217 |
y3d_prime(i,j,1) = y3d_m1(i,j,1)*recVar |
218 |
ELSE |
219 |
c3d_prime(i,j,1) = 0. _d 0 |
220 |
y3d_prime(i,j,1) = 0. _d 0 |
221 |
errCode = 1 |
222 |
ENDIF |
223 |
ELSE |
224 |
tmpVar = b3d(i,j,k) - a3d(i,j,k)*c3d_prime(i,j,k-1) |
225 |
IF ( tmpVar .NE. 0. _d 0 ) THEN |
226 |
recVar = 1. _d 0 / tmpVar |
227 |
c3d_prime(i,j,k) = c3d_m1(i,j,k)*recVar |
228 |
y3d_prime(i,j,k) = ( y3d_m1(i,j,k) |
229 |
& - a3d(i,j,k)*y3d_prime(i,j,k-1) |
230 |
& )*recVar |
231 |
ELSE |
232 |
c3d_prime(i,j,k) = 0. _d 0 |
233 |
y3d_prime(i,j,k) = 0. _d 0 |
234 |
errCode = 1 |
235 |
ENDIF |
236 |
ENDIF |
237 |
ENDDO |
238 |
ENDDO |
239 |
C-- k-loop |
240 |
ENDDO |
241 |
|
242 |
CADJ loop = sequential |
243 |
C-- Backward sweep |
244 |
DO k=Nr,1,-1 |
245 |
DO j=1-OLy,sNy+OLy |
246 |
DO i=1-OLx,sNx+OLx |
247 |
c DO j=jMin,jMax |
248 |
c DO i=iMin,iMax |
249 |
IF ( k.EQ.Nr ) THEN |
250 |
y3d_update(i,j,k) = y3d_prime(i,j,k) |
251 |
ELSE |
252 |
y3d_update(i,j,k) = y3d_prime(i,j,k) |
253 |
& - c3d_prime(i,j,k)*y3d_update(i,j,k+1) |
254 |
ENDIF |
255 |
ENDDO |
256 |
ENDDO |
257 |
C-- k-loop |
258 |
ENDDO |
259 |
|
260 |
C-- Update array |
261 |
DO k=1,Nr |
262 |
DO j=1-OLy,sNy+OLy |
263 |
DO i=1-OLx,sNx+OLx |
264 |
c DO j=jMin,jMax |
265 |
c DO i=iMin,iMax |
266 |
y3d(i,j,k) = y3d_update(i,j,k) |
267 |
ENDDO |
268 |
ENDDO |
269 |
C-- k-loop |
270 |
ENDDO |
271 |
|
272 |
#endif /* SOLVE_DIAGONAL_KINNER */ |
273 |
|
274 |
#endif /* SOLVE_DIAGONAL_LOWMEMORY */ |
275 |
|
276 |
RETURN |
277 |
END |