1 |
C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_cg_solve.F,v 1.8 2015/02/16 16:46:44 dgoldberg Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "STREAMICE_OPTIONS.h" |
5 |
|
6 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
7 |
|
8 |
CBOP |
9 |
SUBROUTINE STREAMICE_CG_SOLVE( |
10 |
U cg_Uin, ! x-velocities |
11 |
U cg_Vin, ! y-velocities |
12 |
I cg_Bu, ! force in x dir |
13 |
I cg_Bv, ! force in y dir |
14 |
I A_uu, ! section of matrix that multiplies u and projects on u |
15 |
I A_uv, ! section of matrix that multiplies v and projects on u |
16 |
I A_vu, ! section of matrix that multiplies u and projects on v |
17 |
I A_vv, ! section of matrix that multiplies v and projects on v |
18 |
I tolerance, |
19 |
O iters, |
20 |
I maxIter, |
21 |
I myThid ) |
22 |
C /============================================================\ |
23 |
C | SUBROUTINE | |
24 |
C | o | |
25 |
C |============================================================| |
26 |
C | | |
27 |
C \============================================================/ |
28 |
IMPLICIT NONE |
29 |
|
30 |
#include "SIZE.h" |
31 |
#include "EEPARAMS.h" |
32 |
#include "PARAMS.h" |
33 |
#include "STREAMICE.h" |
34 |
#include "STREAMICE_CG.h" |
35 |
|
36 |
|
37 |
|
38 |
!#ifdef ALLOW_PETSC |
39 |
!#include "finclude/petsc.h" |
40 |
! UNCOMMENT IF V3.0 |
41 |
!#include "finclude/petscvec.h" |
42 |
!#include "finclude/petscmat.h" |
43 |
!#include "finclude/petscksp.h" |
44 |
!#include "finclude/petscpc.h" |
45 |
!#endif |
46 |
C === Global variables === |
47 |
|
48 |
|
49 |
C !INPUT/OUTPUT ARGUMENTS |
50 |
C cg_Uin, cg_Vin - input and output velocities |
51 |
C cg_Bu, cg_Bv - driving stress |
52 |
INTEGER myThid |
53 |
INTEGER iters |
54 |
INTEGER maxIter |
55 |
_RL tolerance |
56 |
_RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
57 |
_RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
58 |
_RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
59 |
_RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
60 |
_RL |
61 |
& A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), |
62 |
& A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), |
63 |
& A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), |
64 |
& A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1) |
65 |
|
66 |
C LOCAL VARIABLES |
67 |
INTEGER i, j, bi, bj, cg_halo, conv_flag |
68 |
INTEGER iter, is, js, ie, je, colx, coly, k |
69 |
_RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0 |
70 |
_RL dot_p1_tile (nSx,nSy) |
71 |
_RL dot_p2_tile (nSx,nSy) |
72 |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
73 |
|
74 |
|
75 |
!#ifdef ALLOW_PETSC |
76 |
! INTEGER indices(2*(snx*nsx*sny*nsy)) |
77 |
! INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1) |
78 |
! _RL rhs_values(2*(snx*nsx*sny*nsy)) |
79 |
! _RL solution_values(2*(snx*nsx*sny*nsy)) |
80 |
! _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy)) |
81 |
! _RL mat_values (18,1), mat_val_return(1) |
82 |
! INTEGER indices_col(18) |
83 |
! INTEGER local_dofs, global_dofs, dof_index, dof_index_col |
84 |
! INTEGER local_offset |
85 |
! Mat matrix |
86 |
! KSP ksp |
87 |
! PC pc |
88 |
! Vec rhs |
89 |
! Vec solution |
90 |
! PetscErrorCode ierr |
91 |
!#ifdef ALLOW_USE_MPI |
92 |
! integer mpiRC, mpiMyWid |
93 |
!#endif |
94 |
!#endif |
95 |
|
96 |
|
97 |
#ifdef ALLOW_STREAMICE |
98 |
|
99 |
|
100 |
|
101 |
CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid) |
102 |
#ifndef STREAMICE_SERIAL_TRISOLVE |
103 |
|
104 |
#ifdef ALLOW_PETSC |
105 |
|
106 |
if (streamice_use_petsc) then |
107 |
|
108 |
CALL STREAMICE_CG_SOLVE_PETSC( |
109 |
U cg_Uin, ! x-velocities |
110 |
U cg_Vin, ! y-velocities |
111 |
I cg_Bu, ! force in x dir |
112 |
I cg_Bv, ! force in y dir |
113 |
I A_uu, ! section of matrix that multiplies u and projects on u |
114 |
I A_uv, ! section of matrix that multiplies v and projects on u |
115 |
I A_vu, ! section of matrix that multiplies u and projects on v |
116 |
I A_vv, ! section of matrix that multiplies v and projects on v |
117 |
I tolerance, |
118 |
I iters, |
119 |
O maxiter, |
120 |
I myThid ) |
121 |
|
122 |
|
123 |
else |
124 |
|
125 |
#endif /* ALLOW_PETSC */ |
126 |
|
127 |
|
128 |
iters = maxIter |
129 |
conv_flag = 0 |
130 |
|
131 |
|
132 |
DO bj = myByLo(myThid), myByHi(myThid) |
133 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
134 |
DO j=1,sNy |
135 |
DO i=1,sNx |
136 |
Zu_SI (i,j,bi,bj) = 0. _d 0 |
137 |
Zv_SI (i,j,bi,bj) = 0. _d 0 |
138 |
Ru_SI (i,j,bi,bj) = 0. _d 0 |
139 |
Rv_SI (i,j,bi,bj) = 0. _d 0 |
140 |
Au_SI (i,j,bi,bj) = 0. _d 0 |
141 |
Av_SI (i,j,bi,bj) = 0. _d 0 |
142 |
Du_SI (i,j,bi,bj) = 0. _d 0 |
143 |
Dv_SI (i,j,bi,bj) = 0. _d 0 |
144 |
ENDDO |
145 |
ENDDO |
146 |
ENDDO |
147 |
ENDDO |
148 |
|
149 |
C FIND INITIAL RESIDUAL, and initialize r |
150 |
|
151 |
! #ifdef STREAMICE_CONSTRUCT_MATRIX |
152 |
|
153 |
|
154 |
|
155 |
DO bj = myByLo(myThid), myByHi(myThid) |
156 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
157 |
DO j=1,sNy |
158 |
DO i=1,sNx |
159 |
DO colx=-1,1 |
160 |
DO coly=-1,1 |
161 |
Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) + |
162 |
& A_uu(i,j,bi,bj,colx,coly)* |
163 |
& cg_Uin(i+colx,j+coly,bi,bj)+ |
164 |
& A_uv(i,j,bi,bj,colx,coly)* |
165 |
& cg_Vin(i+colx,j+coly,bi,bj) |
166 |
|
167 |
|
168 |
Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) + |
169 |
& A_vu(i,j,bi,bj,colx,coly)* |
170 |
& cg_Uin(i+colx,j+coly,bi,bj)+ |
171 |
& A_vv(i,j,bi,bj,colx,coly)* |
172 |
& cg_Vin(i+colx,j+coly,bi,bj) |
173 |
ENDDO |
174 |
ENDDO |
175 |
ENDDO |
176 |
ENDDO |
177 |
ENDDO |
178 |
ENDDO |
179 |
|
180 |
|
181 |
_EXCH_XY_RL( Au_SI, myThid ) |
182 |
_EXCH_XY_RL( Av_SI, myThid ) |
183 |
|
184 |
|
185 |
DO bj = myByLo(myThid), myByHi(myThid) |
186 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
187 |
DO j=1-OLy,sNy+OLy |
188 |
DO i=1-OLx,sNx+OLx |
189 |
Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)- |
190 |
& Au_SI(i,j,bi,bj) |
191 |
Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)- |
192 |
& Av_SI(i,j,bi,bj) |
193 |
ENDDO |
194 |
ENDDO |
195 |
dot_p1_tile(bi,bj) = 0. _d 0 |
196 |
dot_p2_tile(bi,bj) = 0. _d 0 |
197 |
ENDDO |
198 |
ENDDO |
199 |
|
200 |
DO bj = myByLo(myThid), myByHi(myThid) |
201 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
202 |
DO j=1,sNy |
203 |
DO i=1,sNx |
204 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) |
205 |
& dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2 |
206 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) |
207 |
& dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2 |
208 |
ENDDO |
209 |
ENDDO |
210 |
ENDDO |
211 |
ENDDO |
212 |
|
213 |
|
214 |
CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) |
215 |
resid_0 = sqrt(dot_p1) |
216 |
|
217 |
DO bj = myByLo(myThid), myByHi(myThid) |
218 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
219 |
|
220 |
WRITE(msgBuf,'(A,I1,I1,E14.7)') 'CONJ GRAD INIT RESID LOCAL, ', |
221 |
& bi,bj, dot_p1_tile(bi,bj) |
222 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
223 |
& SQUEEZE_RIGHT , 1) |
224 |
|
225 |
enddo |
226 |
enddo |
227 |
|
228 |
WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ', |
229 |
& resid_0 |
230 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
231 |
& SQUEEZE_RIGHT , 1) |
232 |
|
233 |
C CCCCCCCCCCCCCCCCCCCC |
234 |
|
235 |
DO bj = myByLo(myThid), myByHi(myThid) |
236 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
237 |
DO j=1-OLy,sNy+OLy |
238 |
DO i=1-OLx,sNx+OLx |
239 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) |
240 |
& Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj) |
241 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) |
242 |
& Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj) |
243 |
ENDDO |
244 |
ENDDO |
245 |
ENDDO |
246 |
ENDDO |
247 |
|
248 |
cg_halo = min(OLx-1,OLy-1) |
249 |
conv_flag = 0 |
250 |
|
251 |
DO bj = myByLo(myThid), myByHi(myThid) |
252 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
253 |
DO j=1-OLy,sNy+OLy |
254 |
DO i=1-OLx,sNx+OLx |
255 |
Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj) |
256 |
Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj) |
257 |
ENDDO |
258 |
ENDDO |
259 |
ENDDO |
260 |
ENDDO |
261 |
|
262 |
resid = resid_0 |
263 |
iters = 0 |
264 |
|
265 |
c !!!!!!!!!!!!!!!!!! |
266 |
c !! !! |
267 |
c !! MAIN CG LOOP !! |
268 |
c !! !! |
269 |
c !!!!!!!!!!!!!!!!!! |
270 |
|
271 |
|
272 |
|
273 |
|
274 |
c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!) |
275 |
|
276 |
WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP' |
277 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
278 |
& SQUEEZE_RIGHT , 1) |
279 |
|
280 |
! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid) |
281 |
|
282 |
|
283 |
do iter = 1, maxIter |
284 |
if (resid .gt. tolerance*resid_0) then |
285 |
|
286 |
c to avoid using "exit" |
287 |
iters = iters + 1 |
288 |
|
289 |
is = 1 - cg_halo |
290 |
ie = sNx + cg_halo |
291 |
js = 1 - cg_halo |
292 |
je = sNy + cg_halo |
293 |
|
294 |
DO bj = myByLo(myThid), myByHi(myThid) |
295 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
296 |
DO j=1-OLy,sNy+OLy |
297 |
DO i=1-OLx,sNx+OLx |
298 |
Au_SI(i,j,bi,bj) = 0. _d 0 |
299 |
Av_SI(i,j,bi,bj) = 0. _d 0 |
300 |
ENDDO |
301 |
ENDDO |
302 |
ENDDO |
303 |
ENDDO |
304 |
|
305 |
! IF (STREAMICE_construct_matrix) THEN |
306 |
|
307 |
! #ifdef STREAMICE_CONSTRUCT_MATRIX |
308 |
|
309 |
DO bj = myByLo(myThid), myByHi(myThid) |
310 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
311 |
DO j=js,je |
312 |
DO i=is,ie |
313 |
DO colx=-1,1 |
314 |
DO coly=-1,1 |
315 |
Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) + |
316 |
& A_uu(i,j,bi,bj,colx,coly)* |
317 |
& Du_SI(i+colx,j+coly,bi,bj)+ |
318 |
& A_uv(i,j,bi,bj,colx,coly)* |
319 |
& Dv_SI(i+colx,j+coly,bi,bj) |
320 |
Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) + |
321 |
& A_vu(i,j,bi,bj,colx,coly)* |
322 |
& Du_SI(i+colx,j+coly,bi,bj)+ |
323 |
& A_vv(i,j,bi,bj,colx,coly)* |
324 |
& Dv_SI(i+colx,j+coly,bi,bj) |
325 |
ENDDO |
326 |
ENDDO |
327 |
ENDDO |
328 |
ENDDO |
329 |
ENDDO |
330 |
ENDDO |
331 |
|
332 |
! else |
333 |
! #else |
334 |
! |
335 |
! CALL STREAMICE_CG_ACTION( myThid, |
336 |
! O Au_SI, |
337 |
! O Av_SI, |
338 |
! I Du_SI, |
339 |
! I Dv_SI, |
340 |
! I is,ie,js,je) |
341 |
! |
342 |
! ! ENDIF |
343 |
! |
344 |
! #endif |
345 |
|
346 |
|
347 |
DO bj = myByLo(myThid), myByHi(myThid) |
348 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
349 |
dot_p1_tile(bi,bj) = 0. _d 0 |
350 |
dot_p2_tile(bi,bj) = 0. _d 0 |
351 |
ENDDO |
352 |
ENDDO |
353 |
|
354 |
DO bj = myByLo(myThid), myByHi(myThid) |
355 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
356 |
DO j=1,sNy |
357 |
DO i=1,sNx |
358 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN |
359 |
dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)* |
360 |
& Ru_SI(i,j,bi,bj) |
361 |
dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)* |
362 |
& Au_SI(i,j,bi,bj) |
363 |
ENDIF |
364 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN |
365 |
dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)* |
366 |
& Rv_SI(i,j,bi,bj) |
367 |
dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)* |
368 |
& Av_SI(i,j,bi,bj) |
369 |
ENDIF |
370 |
ENDDO |
371 |
ENDDO |
372 |
ENDDO |
373 |
ENDDO |
374 |
|
375 |
CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) |
376 |
CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid ) |
377 |
alpha_k = dot_p1/dot_p2 |
378 |
|
379 |
DO bj = myByLo(myThid), myByHi(myThid) |
380 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
381 |
DO j=1-OLy,sNy+OLy |
382 |
DO i=1-OLx,sNx+OLx |
383 |
|
384 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN |
385 |
cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+ |
386 |
& alpha_k*Du_SI(i,j,bi,bj) |
387 |
Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) |
388 |
Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj) |
389 |
Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)- |
390 |
& alpha_k*Au_SI(i,j,bi,bj) |
391 |
Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) / |
392 |
& DIAGu_SI(i,j,bi,bj) |
393 |
ENDIF |
394 |
|
395 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN |
396 |
cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+ |
397 |
& alpha_k*Dv_SI(i,j,bi,bj) |
398 |
Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) |
399 |
Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj) |
400 |
Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)- |
401 |
& alpha_k*Av_SI(i,j,bi,bj) |
402 |
Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) / |
403 |
& DIAGv_SI(i,j,bi,bj) |
404 |
|
405 |
ENDIF |
406 |
ENDDO |
407 |
ENDDO |
408 |
ENDDO |
409 |
ENDDO |
410 |
|
411 |
DO bj = myByLo(myThid), myByHi(myThid) |
412 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
413 |
dot_p1_tile(bi,bj) = 0. _d 0 |
414 |
dot_p2_tile(bi,bj) = 0. _d 0 |
415 |
ENDDO |
416 |
ENDDO |
417 |
|
418 |
DO bj = myByLo(myThid), myByHi(myThid) |
419 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
420 |
DO j=1,sNy |
421 |
DO i=1,sNx |
422 |
|
423 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN |
424 |
dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)* |
425 |
& Ru_SI(i,j,bi,bj) |
426 |
dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)* |
427 |
& Ru_old_SI(i,j,bi,bj) |
428 |
ENDIF |
429 |
|
430 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN |
431 |
dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)* |
432 |
& Rv_SI(i,j,bi,bj) |
433 |
dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)* |
434 |
& Rv_old_SI(i,j,bi,bj) |
435 |
ENDIF |
436 |
|
437 |
ENDDO |
438 |
ENDDO |
439 |
ENDDO |
440 |
ENDDO |
441 |
|
442 |
CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) |
443 |
CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid ) |
444 |
|
445 |
beta_k = dot_p1/dot_p2 |
446 |
|
447 |
DO bj = myByLo(myThid), myByHi(myThid) |
448 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
449 |
DO j=1-OLy,sNy+OLy |
450 |
DO i=1-OLx,sNx+OLx |
451 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) |
452 |
& Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+ |
453 |
& Zu_SI(i,j,bi,bj) |
454 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) |
455 |
& Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+ |
456 |
& Zv_SI(i,j,bi,bj) |
457 |
ENDDO |
458 |
ENDDO |
459 |
ENDDO |
460 |
ENDDO |
461 |
|
462 |
DO bj = myByLo(myThid), myByHi(myThid) |
463 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
464 |
dot_p1_tile(bi,bj) = 0. _d 0 |
465 |
ENDDO |
466 |
ENDDO |
467 |
|
468 |
DO bj = myByLo(myThid), myByHi(myThid) |
469 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
470 |
DO j=1,sNy |
471 |
DO i=1,sNx |
472 |
IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) |
473 |
& dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2 |
474 |
IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) |
475 |
& dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2 |
476 |
ENDDO |
477 |
ENDDO |
478 |
ENDDO |
479 |
ENDDO |
480 |
|
481 |
CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) |
482 |
resid = sqrt(dot_p1) |
483 |
|
484 |
! IF (iter .eq. 1) then |
485 |
! print *, alpha_k, beta_k, resid |
486 |
! ENDIF |
487 |
|
488 |
cg_halo = cg_halo - 1 |
489 |
|
490 |
if (cg_halo .eq. 0) then |
491 |
cg_halo = min(OLx-1,OLy-1) |
492 |
_EXCH_XY_RL( Du_SI, myThid ) |
493 |
_EXCH_XY_RL( Dv_SI, myThid ) |
494 |
_EXCH_XY_RL( Ru_SI, myThid ) |
495 |
_EXCH_XY_RL( Rv_SI, myThid ) |
496 |
_EXCH_XY_RL( cg_Uin, myThid ) |
497 |
_EXCH_XY_RL( cg_Vin, myThid ) |
498 |
endif |
499 |
|
500 |
|
501 |
endif |
502 |
enddo ! end of CG loop |
503 |
|
504 |
c to avoid using "exit" |
505 |
c if iters has reached max_iters there is no convergence |
506 |
|
507 |
IF (iters .lt. maxIter) THEN |
508 |
conv_flag = 1 |
509 |
ENDIF |
510 |
PRINT *, "GOT HERE CG ITERATIONS", iters |
511 |
|
512 |
! DO bj = myByLo(myThid), myByHi(myThid) |
513 |
! DO bi = myBxLo(myThid), myBxHi(myThid) |
514 |
! DO j=1-OLy,sNy+OLy |
515 |
! DO i=1-OLy,sNx+OLy |
516 |
! IF (STREAMICE_umask(i,j,bi,bj).eq.3.0) |
517 |
! & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj) |
518 |
! IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0) |
519 |
! & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj) |
520 |
! ENDDO |
521 |
! ENDDO |
522 |
! ENDDO |
523 |
! ENDDO |
524 |
! |
525 |
! _EXCH_XY_RL( cg_Uin, myThid ) |
526 |
! _EXCH_XY_RL( cg_Vin, myThid ) |
527 |
|
528 |
|
529 |
#ifdef ALLOW_PETSC |
530 |
endif !if (streamice_use_petsc) |
531 |
#endif |
532 |
|
533 |
#else /* STREAMICE_SERIAL_TRISOLVE */ |
534 |
|
535 |
iters = 0 |
536 |
|
537 |
CALL STREAMICE_TRIDIAG_SOLVE( |
538 |
U cg_Uin, ! x-velocities |
539 |
U cg_Vin, |
540 |
U cg_Bu, ! force in x dir |
541 |
I A_uu, ! section of matrix that multiplies u and projects on u |
542 |
I STREAMICE_umask, |
543 |
I myThid ) |
544 |
|
545 |
#endif |
546 |
|
547 |
CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid) |
548 |
|
549 |
|
550 |
#endif |
551 |
RETURN |
552 |
END |