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

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

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


Revision 1.8 - (hide annotations) (download)
Tue May 28 22:32:39 2013 UTC (12 years, 1 month ago) by dgoldberg
Branch: MAIN
Changes since 1.7: +26 -7 lines
tridiagonal matrix solver for flowline case

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

  ViewVC Help
Powered by ViewVC 1.1.22