/[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.7 - (hide 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 dgoldberg 1.7 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F,v 1.6 2012/07/26 16:22:58 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.7 #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 heimbach 1.1
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 dgoldberg 1.4 _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 heimbach 1.1
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 dgoldberg 1.4 _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
64 heimbach 1.1 _RL dot_p1_tile (nSx,nSy)
65     _RL dot_p2_tile (nSx,nSy)
66 dgoldberg 1.5 CHARACTER*(MAX_LEN_MBUF) msgBuf
67 heimbach 1.1
68 dgoldberg 1.7
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 heimbach 1.1
91     #ifdef ALLOW_STREAMICE
92    
93 dgoldberg 1.7 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 dgoldberg 1.4
219 dgoldberg 1.7 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 heimbach 1.1 conv_flag = 0
375    
376 dgoldberg 1.4
377 heimbach 1.1 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 dgoldberg 1.4
394     C FIND INITIAL RESIDUAL, and initialize r
395 heimbach 1.1
396 dgoldberg 1.4 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
397 heimbach 1.1
398 dgoldberg 1.4
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 heimbach 1.1
412 dgoldberg 1.4
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 heimbach 1.1 ENDDO
423     ENDDO
424    
425 dgoldberg 1.4
426 heimbach 1.1 _EXCH_XY_RL( Au_SI, myThid )
427     _EXCH_XY_RL( Av_SI, myThid )
428    
429 dgoldberg 1.4
430 heimbach 1.1 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 dgoldberg 1.4 Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
435 heimbach 1.1 & Au_SI(i,j,bi,bj)
436 dgoldberg 1.4 Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
437 heimbach 1.1 & 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 dgoldberg 1.4
446    
447 heimbach 1.1 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 dgoldberg 1.5 WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
464     & resid_0
465     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
466     & SQUEEZE_RIGHT , 1)
467    
468 dgoldberg 1.4 C CCCCCCCCCCCCCCCCCCCC
469    
470 heimbach 1.1 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 dgoldberg 1.4
506    
507    
508 heimbach 1.1
509     c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
510    
511 dgoldberg 1.6 WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
512     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
513     & SQUEEZE_RIGHT , 1)
514 heimbach 1.1
515 dgoldberg 1.3 ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
516    
517 heimbach 1.1
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 dgoldberg 1.3 ! IF (STREAMICE_construct_matrix) THEN
541    
542 dgoldberg 1.4 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
543 dgoldberg 1.3
544 heimbach 1.1 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 dgoldberg 1.4 & A_uu(i,j,bi,bj,colx,coly)*
552 heimbach 1.1 & Du_SI(i+colx,j+coly,bi,bj)+
553 dgoldberg 1.4 & A_uv(i,j,bi,bj,colx,coly)*
554 heimbach 1.1 & Dv_SI(i+colx,j+coly,bi,bj)
555     Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
556 dgoldberg 1.4 & A_vu(i,j,bi,bj,colx,coly)*
557 heimbach 1.1 & Du_SI(i+colx,j+coly,bi,bj)+
558 dgoldberg 1.4 & A_vv(i,j,bi,bj,colx,coly)*
559 heimbach 1.1 & Dv_SI(i+colx,j+coly,bi,bj)
560     ENDDO
561     ENDDO
562     ENDDO
563     ENDDO
564     ENDDO
565     ENDDO
566    
567 dgoldberg 1.3 ! else
568 dgoldberg 1.4 ! #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 heimbach 1.1
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 dgoldberg 1.4
736 heimbach 1.1 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 dgoldberg 1.4 ! 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 heimbach 1.1
762     #endif
763    
764 dgoldberg 1.7 CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
765 heimbach 1.1
766    
767 dgoldberg 1.7 #endif
768     RETURN
769     END

  ViewVC Help
Powered by ViewVC 1.1.22