/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.7 - (show annotations) (download)
Sat Apr 6 17:43:41 2013 UTC (12 years, 3 months ago) by dgoldberg
Branch: MAIN
Changes since 1.6: +315 -10 lines
use PETSc solver for matrix

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F,v 1.6 2012/07/26 16:22:58 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 myThid )
21 C /============================================================\
22 C | SUBROUTINE |
23 C | o |
24 C |============================================================|
25 C | |
26 C \============================================================/
27 IMPLICIT NONE
28
29 #include "SIZE.h"
30 #include "EEPARAMS.h"
31 #include "PARAMS.h"
32 #include "STREAMICE.h"
33 #include "STREAMICE_CG.h"
34 #ifdef ALLOW_PETSC
35 #include "finclude/petsc.h"
36 !#include "finclude/petscvec.h"
37 !#include "finclude/petscmat.h"
38 !#include "finclude/petscksp.h"
39 !#include "finclude/petscpc.h"
40 #endif
41 C === Global variables ===
42
43
44 C !INPUT/OUTPUT ARGUMENTS
45 C cg_Uin, cg_Vin - input and output velocities
46 C cg_Bu, cg_Bv - driving stress
47 INTEGER myThid
48 INTEGER iters
49 _RL tolerance
50 _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51 _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52 _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
53 _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54 _RL
55 & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
56 & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
57 & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
58 & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
59
60 C LOCAL VARIABLES
61 INTEGER i, j, bi, bj, cg_halo, conv_flag
62 INTEGER iter, is, js, ie, je, colx, coly, k
63 _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
64 _RL dot_p1_tile (nSx,nSy)
65 _RL dot_p2_tile (nSx,nSy)
66 CHARACTER*(MAX_LEN_MBUF) msgBuf
67
68
69 #ifdef ALLOW_PETSC
70 INTEGER indices(2*(snx*nsx*sny*nsy))
71 INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
72 _RL rhs_values(2*(snx*nsx*sny*nsy))
73 _RL solution_values(2*(snx*nsx*sny*nsy))
74 ! _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))
75 _RL mat_values (18,1), mat_val_return(1)
76 INTEGER indices_col(18)
77 INTEGER local_dofs, global_dofs, dof_index, dof_index_col
78 INTEGER local_offset
79 Mat matrix
80 KSP ksp
81 PC pc
82 Vec rhs
83 Vec solution
84 PetscErrorCode ierr
85 #ifdef ALLOW_USE_MPI
86 integer mpiRC, mpiMyWid
87 #endif
88 #endif
89
90
91 #ifdef ALLOW_STREAMICE
92
93 CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
94
95 #ifdef ALLOW_PETSC
96
97 #ifdef ALLOW_USE_MPI
98
99
100 CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
101 local_dofs = n_dofs_process (mpiMyWid)
102 global_dofs = 0
103
104 n_dofs_cum_sum(0) = 0
105 DO i=0,nPx*nPy-1
106 global_dofs = global_dofs + n_dofs_process (i)
107 if (i.ge.1) THEN
108 n_dofs_cum_sum(i) = n_dofs_cum_sum(i-1)+
109 & n_dofs_process(i-1)
110 endif
111 ENDDO
112 local_offset = n_dofs_cum_sum(mpimywid)
113
114 #else
115
116 local_dofs = n_dofs_process (0)
117 global_dofs = local_dofs
118 local_offset = 0
119
120 #endif
121
122 ! call petscInitialize(PETSC_NULL_CHARACTER,ierr)
123
124 !----------------------
125
126 call VecCreate(PETSC_COMM_WORLD, rhs, ierr)
127 call VecSetSizes(rhs, local_dofs, global_dofs, ierr)
128 call VecSetType(rhs, VECMPI, ierr)
129
130 call VecCreate(PETSC_COMM_WORLD, solution, ierr)
131 call VecSetSizes(solution, local_dofs, global_dofs, ierr)
132 call VecSetType(solution, VECMPI, ierr)
133
134 do i=1,local_dofs
135 indices(i) = i-1 + local_offset
136 end do
137 do i=1,2*nSx*nSy*sNx*sNy
138 rhs_values (i) = 0. _d 0
139 solution_values (i) = 0. _d 0
140 enddo
141
142 ! gather rhs and initial guess values to populate petsc vectors
143
144 DO bj = myByLo(myThid), myByHi(myThid)
145 DO bi = myBxLo(myThid), myBxHi(myThid)
146 DO j=1,sNy
147 DO i=1,sNx
148
149 dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
150 & - local_offset
151
152 if (dof_index.ge.0) THEN
153
154 rhs_values(dof_index+1) = cg_Bu(i,j,bi,bj)
155 solution_values(dof_index+1) = cg_Uin(i,j,bi,bj)
156
157 endif
158
159 !---------------
160
161 dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
162 & - local_offset
163
164 if (dof_index.ge.0) THEN
165
166 rhs_values(dof_index+1) = cg_Bv(i,j,bi,bj)
167 solution_values(dof_index+1) = cg_Vin(i,j,bi,bj)
168
169 endif
170
171 ENDDO
172 ENDDO
173 ENDDO
174 ENDDO
175
176
177 call VecSetValues(rhs, local_dofs, indices, rhs_values,
178 & INSERT_VALUES, ierr)
179 call VecAssemblyBegin(rhs, ierr)
180 call VecAssemblyEnd(rhs, ierr)
181
182
183 call VecSetValues(solution, local_dofs, indices,
184 & solution_values, INSERT_VALUES, ierr)
185 call VecAssemblyBegin(solution, ierr)
186 call VecAssemblyEnd(solution, ierr)
187
188
189 call MatCreateAIJ (PETSC_COMM_WORLD,
190 & local_dofs, local_dofs,
191 & global_dofs, global_dofs,
192 & 18, PETSC_NULL_INTEGER,
193 & 18, PETSC_NULL_INTEGER,
194 & matrix, ierr)
195
196 ! populate petsc matrix
197
198 DO bj = myByLo(myThid), myByHi(myThid)
199 DO bi = myBxLo(myThid), myBxHi(myThid)
200 DO j=1,sNy
201 DO i=1,sNx
202
203 dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
204 ! & - local_offset
205
206 IF (dof_index .ge. 0) THEN
207
208 DO k=1,18
209 indices_col(k) = 0
210 mat_values(k,1) = 0. _d 0
211 ENDDO
212 k=0
213
214 DO coly=-1,1
215 DO colx=-1,1
216
217 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
218
219 if (dof_index_col.ge.0) THEN
220 ! pscal = A_uu(i,j,bi,bj,colx,coly)
221 ! CALL MatSetValue (matrix,dof_index, dof_index_col,
222 ! & pscal,INSERT_VALUES,ierr)
223 k=k+1
224 mat_values (k,1) = A_uu(i,j,bi,bj,colx,coly)
225 indices_col (k) = dof_index_col
226 endif
227
228 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
229
230 if (dof_index_col.ge.0) THEN
231 ! CALL MatSetValue (matrix,dof_index, dof_index_col,
232 ! & A_uv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
233 k=k+1
234 mat_values (k,1) = A_uv(i,j,bi,bj,colx,coly)
235 indices_col (k) = dof_index_col
236 endif
237
238 ENDDO
239 ENDDO
240
241 call matSetValues (matrix, 1, dof_index, k, indices_col,
242 & mat_values,INSERT_VALUES,ierr)
243
244
245 ENDIF
246
247 ! ----------------------------------------------
248
249 dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
250 ! & - local_offset
251
252 IF (dof_index .ge. 0) THEN
253
254 DO k=1,18
255 indices_col(k) = 0
256 mat_values(k,1) = 0. _d 0
257 ENDDO
258 k=0
259
260 DO coly=-1,1
261 DO colx=-1,1
262
263 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
264
265 if (dof_index_col.ge.0) THEN
266 ! CALL MatSetValue (matrix,dof_index, dof_index_col,
267 ! & A_vu(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
268 k=k+1
269 mat_values (k,1) = A_vu(i,j,bi,bj,colx,coly)
270 indices_col (k) = dof_index_col
271 endif
272
273 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
274
275 if (dof_index_col.ge.0) THEN
276 ! CALL MatSetValue (matrix,dof_index, dof_index_col,
277 ! & A_vv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
278 k=k+1
279 mat_values (k,1) = A_vv(i,j,bi,bj,colx,coly)
280 indices_col (k) = dof_index_col
281 endif
282
283 ENDDO
284 ENDDO
285
286 call matSetValues (matrix, 1, dof_index, k, indices_col,
287 & mat_values,INSERT_VALUES,ierr)
288 ENDIF
289
290 ENDDO
291 ENDDO
292 ENDDO
293 ENDDO
294
295 call MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY,ierr)
296 call MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY,ierr)
297
298
299 call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
300 call KSPSetOperators(ksp, matrix, matrix,
301 & DIFFERENT_NONZERO_PATTERN, ierr)
302
303 SELECT CASE (PETSC_SOLVER_TYPE)
304 CASE ('CG')
305 PRINT *, "PETSC SOLVER: SELECTED CG"
306 call KSPSetType(ksp, KSPCG, ierr)
307 CASE ('GMRES')
308 PRINT *, "PETSC SOLVER: SELECTED GMRES"
309 call KSPSetType(ksp, KSPGMRES, ierr)
310 CASE ('BICG')
311 PRINT *, "PETSC SOLVER: SELECTED BICG"
312 call KSPSetType(ksp, KSPBICG, ierr)
313 CASE DEFAULT
314 PRINT *, "PETSC SOLVER: SELECTED DEFAULT"
315 call KSPSetType(ksp, KSPCG, ierr)
316 END SELECT
317
318 call KSPGetPC(ksp, pc, ierr)
319 call KSPSetTolerances(ksp,tolerance,
320 & PETSC_DEFAULT_DOUBLE_PRECISION,
321 & PETSC_DEFAULT_DOUBLE_PRECISION,
322 & streamice_max_cg_iter,ierr)
323
324 SELECT CASE (PETSC_PRECOND_TYPE)
325 CASE ('BLOCKJACOBI')
326 PRINT *, "PETSC PRECOND: SELECTED BJACOBI"
327 call PCSetType(pc, PCBJACOBI, ierr)
328 CASE ('JACOBI')
329 PRINT *, "PETSC PRECOND: SELECTED JACOBI"
330 call PCSetType(pc, PCJACOBI, ierr)
331 CASE ('ILU')
332 PRINT *, "PETSC PRECOND: SELECTED ILU"
333 call PCSetType(pc, PCILU, ierr)
334 CASE DEFAULT
335 PRINT *, "PETSC PRECOND: SELECTED DEFAULT"
336 call PCSetType(pc, PCBJACOBI, ierr)
337 END SELECT
338
339
340 call KSPSolve(ksp, rhs, solution, ierr)
341 call KSPGetIterationNumber(ksp,iters,ierr)
342
343 call VecGetValues(solution,local_dofs,indices,
344 & solution_values,ierr)
345
346 DO bj = myByLo(myThid), myByHi(myThid)
347 DO bi = myBxLo(myThid), myBxHi(myThid)
348 DO j=1,sNy
349 DO i=1,sNx
350
351 dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
352 & - local_offset
353 if (dof_index.ge.0) THEN
354 cg_Uin(i,j,bi,bj) = solution_values(dof_index+1)
355 endif
356
357 dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
358 & - local_offset
359 if (dof_index.ge.0) THEN
360 cg_Vin(i,j,bi,bj) = solution_values(dof_index+1)
361 endif
362
363 ENDDO
364 ENDDO
365 ENDDO
366 ENDDO
367
368
369
370 #else
371
372
373 iters = streamice_max_cg_iter
374 conv_flag = 0
375
376
377 DO bj = myByLo(myThid), myByHi(myThid)
378 DO bi = myBxLo(myThid), myBxHi(myThid)
379 DO j=1,sNy
380 DO i=1,sNx
381 Zu_SI (i,j,bi,bj) = 0. _d 0
382 Zv_SI (i,j,bi,bj) = 0. _d 0
383 Ru_SI (i,j,bi,bj) = 0. _d 0
384 Rv_SI (i,j,bi,bj) = 0. _d 0
385 Au_SI (i,j,bi,bj) = 0. _d 0
386 Av_SI (i,j,bi,bj) = 0. _d 0
387 Du_SI (i,j,bi,bj) = 0. _d 0
388 Dv_SI (i,j,bi,bj) = 0. _d 0
389 ENDDO
390 ENDDO
391 ENDDO
392 ENDDO
393
394 C FIND INITIAL RESIDUAL, and initialize r
395
396 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
397
398
399
400 DO bj = myByLo(myThid), myByHi(myThid)
401 DO bi = myBxLo(myThid), myBxHi(myThid)
402 DO j=1,sNy
403 DO i=1,sNx
404 DO colx=-1,1
405 DO coly=-1,1
406 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
407 & A_uu(i,j,bi,bj,colx,coly)*
408 & cg_Uin(i+colx,j+coly,bi,bj)+
409 & A_uv(i,j,bi,bj,colx,coly)*
410 & cg_Vin(i+colx,j+coly,bi,bj)
411
412
413 Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
414 & A_vu(i,j,bi,bj,colx,coly)*
415 & cg_Uin(i+colx,j+coly,bi,bj)+
416 & A_vv(i,j,bi,bj,colx,coly)*
417 & cg_Vin(i+colx,j+coly,bi,bj)
418 ENDDO
419 ENDDO
420 ENDDO
421 ENDDO
422 ENDDO
423 ENDDO
424
425
426 _EXCH_XY_RL( Au_SI, myThid )
427 _EXCH_XY_RL( Av_SI, myThid )
428
429
430 DO bj = myByLo(myThid), myByHi(myThid)
431 DO bi = myBxLo(myThid), myBxHi(myThid)
432 DO j=1-OLy,sNy+OLy
433 DO i=1-OLx,sNx+OLx
434 Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
435 & Au_SI(i,j,bi,bj)
436 Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
437 & Av_SI(i,j,bi,bj)
438 ENDDO
439 ENDDO
440 dot_p1_tile(bi,bj) = 0. _d 0
441 dot_p2_tile(bi,bj) = 0. _d 0
442 ENDDO
443 ENDDO
444
445
446
447 DO bj = myByLo(myThid), myByHi(myThid)
448 DO bi = myBxLo(myThid), myBxHi(myThid)
449 DO j=1,sNy
450 DO i=1,sNx
451 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
452 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
453 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
454 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
455 ENDDO
456 ENDDO
457 ENDDO
458 ENDDO
459
460 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
461 resid_0 = sqrt(dot_p1)
462
463 WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
464 & resid_0
465 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
466 & SQUEEZE_RIGHT , 1)
467
468 C CCCCCCCCCCCCCCCCCCCC
469
470 DO bj = myByLo(myThid), myByHi(myThid)
471 DO bi = myBxLo(myThid), myBxHi(myThid)
472 DO j=1-OLy,sNy+OLy
473 DO i=1-OLx,sNx+OLx
474 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
475 & Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)
476 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
477 & Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)
478 ENDDO
479 ENDDO
480 ENDDO
481 ENDDO
482
483 cg_halo = min(OLx-1,OLy-1)
484 conv_flag = 0
485
486 DO bj = myByLo(myThid), myByHi(myThid)
487 DO bi = myBxLo(myThid), myBxHi(myThid)
488 DO j=1-OLy,sNy+OLy
489 DO i=1-OLx,sNx+OLx
490 Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
491 Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
492 ENDDO
493 ENDDO
494 ENDDO
495 ENDDO
496
497 resid = resid_0
498 iters = 0
499
500 c !!!!!!!!!!!!!!!!!!
501 c !! !!
502 c !! MAIN CG LOOP !!
503 c !! !!
504 c !!!!!!!!!!!!!!!!!!
505
506
507
508
509 c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
510
511 WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
512 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
513 & SQUEEZE_RIGHT , 1)
514
515 ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
516
517
518 do iter = 1, streamice_max_cg_iter
519 if (resid .gt. tolerance*resid_0) then
520
521 c to avoid using "exit"
522 iters = iters + 1
523
524 is = 1 - cg_halo
525 ie = sNx + cg_halo
526 js = 1 - cg_halo
527 je = sNy + cg_halo
528
529 DO bj = myByLo(myThid), myByHi(myThid)
530 DO bi = myBxLo(myThid), myBxHi(myThid)
531 DO j=1-OLy,sNy+OLy
532 DO i=1-OLx,sNx+OLx
533 Au_SI(i,j,bi,bj) = 0. _d 0
534 Av_SI(i,j,bi,bj) = 0. _d 0
535 ENDDO
536 ENDDO
537 ENDDO
538 ENDDO
539
540 ! IF (STREAMICE_construct_matrix) THEN
541
542 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
543
544 DO bj = myByLo(myThid), myByHi(myThid)
545 DO bi = myBxLo(myThid), myBxHi(myThid)
546 DO j=js,je
547 DO i=is,ie
548 DO colx=-1,1
549 DO coly=-1,1
550 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
551 & A_uu(i,j,bi,bj,colx,coly)*
552 & Du_SI(i+colx,j+coly,bi,bj)+
553 & A_uv(i,j,bi,bj,colx,coly)*
554 & Dv_SI(i+colx,j+coly,bi,bj)
555 Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
556 & A_vu(i,j,bi,bj,colx,coly)*
557 & Du_SI(i+colx,j+coly,bi,bj)+
558 & A_vv(i,j,bi,bj,colx,coly)*
559 & Dv_SI(i+colx,j+coly,bi,bj)
560 ENDDO
561 ENDDO
562 ENDDO
563 ENDDO
564 ENDDO
565 ENDDO
566
567 ! else
568 ! #else
569 !
570 ! CALL STREAMICE_CG_ACTION( myThid,
571 ! O Au_SI,
572 ! O Av_SI,
573 ! I Du_SI,
574 ! I Dv_SI,
575 ! I is,ie,js,je)
576 !
577 ! ! ENDIF
578 !
579 ! #endif
580
581
582 DO bj = myByLo(myThid), myByHi(myThid)
583 DO bi = myBxLo(myThid), myBxHi(myThid)
584 dot_p1_tile(bi,bj) = 0. _d 0
585 dot_p2_tile(bi,bj) = 0. _d 0
586 ENDDO
587 ENDDO
588
589 DO bj = myByLo(myThid), myByHi(myThid)
590 DO bi = myBxLo(myThid), myBxHi(myThid)
591 DO j=1,sNy
592 DO i=1,sNx
593 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
594 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
595 & Ru_SI(i,j,bi,bj)
596 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
597 & Au_SI(i,j,bi,bj)
598 ENDIF
599 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
600 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
601 & Rv_SI(i,j,bi,bj)
602 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
603 & Av_SI(i,j,bi,bj)
604 ENDIF
605 ENDDO
606 ENDDO
607 ENDDO
608 ENDDO
609
610 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
611 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
612 alpha_k = dot_p1/dot_p2
613
614 DO bj = myByLo(myThid), myByHi(myThid)
615 DO bi = myBxLo(myThid), myBxHi(myThid)
616 DO j=1-OLy,sNy+OLy
617 DO i=1-OLx,sNx+OLx
618
619 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
620 cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
621 & alpha_k*Du_SI(i,j,bi,bj)
622 Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
623 Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
624 Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
625 & alpha_k*Au_SI(i,j,bi,bj)
626 Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
627 & DIAGu_SI(i,j,bi,bj)
628 ENDIF
629
630 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
631 cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
632 & alpha_k*Dv_SI(i,j,bi,bj)
633 Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
634 Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
635 Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
636 & alpha_k*Av_SI(i,j,bi,bj)
637 Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
638 & DIAGv_SI(i,j,bi,bj)
639
640 ENDIF
641 ENDDO
642 ENDDO
643 ENDDO
644 ENDDO
645
646 DO bj = myByLo(myThid), myByHi(myThid)
647 DO bi = myBxLo(myThid), myBxHi(myThid)
648 dot_p1_tile(bi,bj) = 0. _d 0
649 dot_p2_tile(bi,bj) = 0. _d 0
650 ENDDO
651 ENDDO
652
653 DO bj = myByLo(myThid), myByHi(myThid)
654 DO bi = myBxLo(myThid), myBxHi(myThid)
655 DO j=1,sNy
656 DO i=1,sNx
657
658 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
659 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
660 & Ru_SI(i,j,bi,bj)
661 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
662 & Ru_old_SI(i,j,bi,bj)
663 ENDIF
664
665 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
666 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
667 & Rv_SI(i,j,bi,bj)
668 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
669 & Rv_old_SI(i,j,bi,bj)
670 ENDIF
671
672 ENDDO
673 ENDDO
674 ENDDO
675 ENDDO
676
677 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
678 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
679
680 beta_k = dot_p1/dot_p2
681
682 DO bj = myByLo(myThid), myByHi(myThid)
683 DO bi = myBxLo(myThid), myBxHi(myThid)
684 DO j=1-OLy,sNy+OLy
685 DO i=1-OLx,sNx+OLx
686 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
687 & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
688 & Zu_SI(i,j,bi,bj)
689 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
690 & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
691 & Zv_SI(i,j,bi,bj)
692 ENDDO
693 ENDDO
694 ENDDO
695 ENDDO
696
697 DO bj = myByLo(myThid), myByHi(myThid)
698 DO bi = myBxLo(myThid), myBxHi(myThid)
699 dot_p1_tile(bi,bj) = 0. _d 0
700 ENDDO
701 ENDDO
702
703 DO bj = myByLo(myThid), myByHi(myThid)
704 DO bi = myBxLo(myThid), myBxHi(myThid)
705 DO j=1,sNy
706 DO i=1,sNx
707 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
708 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
709 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
710 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
711 ENDDO
712 ENDDO
713 ENDDO
714 ENDDO
715
716 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
717 resid = sqrt(dot_p1)
718
719 ! IF (iter .eq. 1) then
720 ! print *, alpha_k, beta_k, resid
721 ! ENDIF
722
723 cg_halo = cg_halo - 1
724
725 if (cg_halo .eq. 0) then
726 cg_halo = min(OLx-1,OLy-1)
727 _EXCH_XY_RL( Du_SI, myThid )
728 _EXCH_XY_RL( Dv_SI, myThid )
729 _EXCH_XY_RL( Ru_SI, myThid )
730 _EXCH_XY_RL( Rv_SI, myThid )
731 _EXCH_XY_RL( cg_Uin, myThid )
732 _EXCH_XY_RL( cg_Vin, myThid )
733 endif
734
735
736 endif
737 enddo ! end of CG loop
738
739 c to avoid using "exit"
740 c if iters has reached max_iters there is no convergence
741
742 IF (iters .lt. streamice_max_cg_iter) THEN
743 conv_flag = 1
744 ENDIF
745
746 ! DO bj = myByLo(myThid), myByHi(myThid)
747 ! DO bi = myBxLo(myThid), myBxHi(myThid)
748 ! DO j=1-OLy,sNy+OLy
749 ! DO i=1-OLy,sNx+OLy
750 ! IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
751 ! & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
752 ! IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
753 ! & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
754 ! ENDDO
755 ! ENDDO
756 ! ENDDO
757 ! ENDDO
758 !
759 ! _EXCH_XY_RL( cg_Uin, myThid )
760 ! _EXCH_XY_RL( cg_Vin, myThid )
761
762 #endif
763
764 CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
765
766
767 #endif
768 RETURN
769 END

  ViewVC Help
Powered by ViewVC 1.1.22