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

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

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

revision 1.2 by heimbach, Wed May 2 02:36:01 2012 UTC revision 1.8 by dgoldberg, Tue May 28 22:32:39 2013 UTC
# Line 7  C---+----1----+----2----+----3----+----4 Line 7  C---+----1----+----2----+----3----+----4
7    
8  CBOP  CBOP
9        SUBROUTINE STREAMICE_CG_SOLVE(        SUBROUTINE STREAMICE_CG_SOLVE(
10       U                               cg_Uin,       U                               cg_Uin,     ! x-velocities
11       U                               cg_Vin,       U                               cg_Vin,     ! y-velocities
12       I                               cg_Bu,       I                               cg_Bu,      ! force in x dir
13       I                               cg_Bv,       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,       I                               tolerance,
19       O                               iters,       O                               iters,
20       I                               myThid )       I                               myThid )
# Line 22  C     | Line 26  C     |
26  C     \============================================================/  C     \============================================================/
27        IMPLICIT NONE        IMPLICIT NONE
28    
 C     === Global variables ===  
29  #include "SIZE.h"  #include "SIZE.h"
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "STREAMICE.h"  #include "STREAMICE.h"
33  #include "STREAMICE_CG.h"  #include "STREAMICE_CG.h"
34    
35    
36    
37    #ifdef ALLOW_PETSC
38    #include "finclude/petsc.h"
39    #include "finclude/petscvec.h"
40    #include "finclude/petscmat.h"
41    #include "finclude/petscksp.h"
42    #include "finclude/petscpc.h"
43    #endif
44    C     === Global variables ===
45    
46                
47  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
48  C     cg_Uin, cg_Vin - input and output velocities  C     cg_Uin, cg_Vin - input and output velocities
# Line 40  C     cg_Bu, cg_Bv - driving stress Line 54  C     cg_Bu, cg_Bv - driving stress
54        _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _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)        _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)        _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57          _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    
63  C     LOCAL VARIABLES  C     LOCAL VARIABLES
64        INTEGER i, j, bi, bj, cg_halo, conv_flag        INTEGER i, j, bi, bj, cg_halo, conv_flag
65        INTEGER iter, is, js, ie, je, colx, coly, k        INTEGER iter, is, js, ie, je, colx, coly, k
66        _RL dot_p1, dot_p2, resid_0, alpha_k, beta_k, resid        _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
67        _RL dot_p1_tile (nSx,nSy)        _RL dot_p1_tile (nSx,nSy)
68        _RL dot_p2_tile (nSx,nSy)        _RL dot_p2_tile (nSx,nSy)
69          CHARACTER*(MAX_LEN_MBUF) msgBuf
70    
71    
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    
       iters = streamice_max_cg_iter  
93    
94  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
95    
96    
97    
98          CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
99    #ifndef STREAMICE_SERIAL_TRISOLVE
100    
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          call MatCreateMPIAIJ (PETSC_COMM_WORLD,
196         &                      local_dofs, local_dofs,
197         &                      global_dofs, global_dofs,
198         &                      18, PETSC_NULL_INTEGER,
199         &                      18, PETSC_NULL_INTEGER,
200         &                      matrix, ierr)
201    
202    
203    ! 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    
226                 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    #else  /* ALLOW_PETSC */
378    
379    
380          iters = streamice_max_cg_iter
381        conv_flag = 0        conv_flag = 0
382    
383    
384        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
385         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
386          DO j=1,sNy          DO j=1,sNy
387           DO i=1,sNx           DO i=1,sNx
388            Zu_SI (i,j,bi,bj) = 0. _d 0            Zu_SI (i,j,bi,bj) = 0. _d 0
389            Zv_SI (i,j,bi,bj) = 0. _d 0            Zv_SI (i,j,bi,bj) = 0. _d 0
           DIAGu_SI (i,j,bi,bj) = 0. _d 0  
           DIAGv_SI (i,j,bi,bj) = 0. _d 0  
390            Ru_SI (i,j,bi,bj) = 0. _d 0            Ru_SI (i,j,bi,bj) = 0. _d 0
391            Rv_SI (i,j,bi,bj) = 0. _d 0            Rv_SI (i,j,bi,bj) = 0. _d 0
392            Au_SI (i,j,bi,bj) = 0. _d 0            Au_SI (i,j,bi,bj) = 0. _d 0
393            Av_SI (i,j,bi,bj) = 0. _d 0            Av_SI (i,j,bi,bj) = 0. _d 0
394            Du_SI (i,j,bi,bj) = 0. _d 0            Du_SI (i,j,bi,bj) = 0. _d 0
395            Dv_SI (i,j,bi,bj) = 0. _d 0            Dv_SI (i,j,bi,bj) = 0. _d 0
           ubd_SI (i,j,bi,bj) = 0. _d 0  
           vbd_SI (i,j,bi,bj) = 0. _d 0  
396           ENDDO           ENDDO
397          ENDDO          ENDDO
398         ENDDO         ENDDO
399        ENDDO        ENDDO
400          
401    C     FIND INITIAL RESIDUAL, and initialize r
402    
403  C     PREAMBLE: get bdry contribution, find matrix diagonal,  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
 C     initialize iterates R (residual), Z, and D  
404    
405        CALL STREAMICE_CG_BOUND_VALS( myThid,              
      O    ubd_SI,  
      O    vbd_SI)  
406    
407        DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
408         DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
409          DO j=1-OLy,sNy+OLy            DO j=1,sNy
410           DO i=1-OLx,sNx+OLx             DO i=1,sNx
411            RHSu_SI (i,j,bi,bj) = cg_Bu(i,j,bi,bj)              DO colx=-1,1
412       &     - ubd_SI(i,j,bi,bj)               DO coly=-1,1
413            RHSv_SI (i,j,bi,bj) = cg_Bv(i,j,bi,bj)                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
414       &     - vbd_SI(i,j,bi,bj)       &         A_uu(i,j,bi,bj,colx,coly)*
415           ENDDO       &         cg_Uin(i+colx,j+coly,bi,bj)+
416          ENDDO       &         A_uv(i,j,bi,bj,colx,coly)*    
417         ENDDO       &         cg_Vin(i+colx,j+coly,bi,bj)
       ENDDO  
         
       _EXCH_XY_RL( RHSu_SI, myThid )  
       _EXCH_XY_RL( RHSv_SI, myThid )  
418    
       CALL STREAMICE_CG_ADIAG( myThid,  
      O    DIAGu_SI,  
      O    DIAGv_SI)  
419    
420        _EXCH_XY_RL( DIAGu_SI, myThid )                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
421        _EXCH_XY_RL( DIAGv_SI, myThid )       &         A_vu(i,j,bi,bj,colx,coly)*
422               &         cg_Uin(i+colx,j+coly,bi,bj)+
423  C     FIX PROBLEM WITH PRECOND LATER       &         A_vv(i,j,bi,bj,colx,coly)*    
424         &         cg_Vin(i+colx,j+coly,bi,bj)
425  !       DO bj = myByLo(myThid), myByHi(myThid)               ENDDO
426  !        DO bi = myBxLo(myThid), myBxHi(myThid)              ENDDO
427  !         DO j=1-OLy,sNy+OLy             ENDDO
428  !          DO i=1-OLx,sNx+OLx            ENDDO
429  !           DIAGu_SI(i,j,bi,bj)=1.0           ENDDO
430  !           DIAGv_SI(i,j,bi,bj)=1.0          ENDDO
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
         
       CALL STREAMICE_CG_ACTION( myThid,  
      O    Au_SI,  
      O    Av_SI,  
      I    cg_Uin,  
      I    cg_Vin,  
      I    0, sNx+1, 0, sNy+1 )  
431    
432                
433        _EXCH_XY_RL( Au_SI, myThid )        _EXCH_XY_RL( Au_SI, myThid )
434        _EXCH_XY_RL( Av_SI, myThid )        _EXCH_XY_RL( Av_SI, myThid )
435    
436    
437        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
438         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
439          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
440           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
441            Ru_SI(i,j,bi,bj)=RHSu_SI(i,j,bi,bj)-            Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
442       &     Au_SI(i,j,bi,bj)       &     Au_SI(i,j,bi,bj)
443            Rv_SI(i,j,bi,bj)=RHSv_SI(i,j,bi,bj)-            Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
444       &     Av_SI(i,j,bi,bj)       &     Av_SI(i,j,bi,bj)
445           ENDDO           ENDDO
446          ENDDO          ENDDO
# Line 143  C     FIX PROBLEM WITH PRECOND LATER Line 449  C     FIX PROBLEM WITH PRECOND LATER
449         ENDDO         ENDDO
450        ENDDO        ENDDO
451    
452          
453    
454        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
455         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
456          DO j=1,sNy          DO j=1,sNy
# Line 159  C     FIX PROBLEM WITH PRECOND LATER Line 467  C     FIX PROBLEM WITH PRECOND LATER
467        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
468        resid_0 = sqrt(dot_p1)        resid_0 = sqrt(dot_p1)
469    
470          WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
471         &         resid_0
472           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
473         &                     SQUEEZE_RIGHT , 1)
474    
475    C    CCCCCCCCCCCCCCCCCCCC
476    
477        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
478         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
479          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 172  C     FIX PROBLEM WITH PRECOND LATER Line 487  C     FIX PROBLEM WITH PRECOND LATER
487         ENDDO         ENDDO
488        ENDDO        ENDDO
489    
   
490        cg_halo = min(OLx-1,OLy-1)        cg_halo = min(OLx-1,OLy-1)
491        conv_flag = 0        conv_flag = 0
492    
# Line 195  c  !!              !! Line 509  c  !!              !!
509  c  !! MAIN CG LOOP !!  c  !! MAIN CG LOOP !!
510  c  !!              !!  c  !!              !!
511  c  !!!!!!!!!!!!!!!!!!  c  !!!!!!!!!!!!!!!!!!
512    
513    
514        
515        
516  c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)  c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
517    
518        print *, "BEGINNING MAIN CG LOOP"        WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
519           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
520         &                     SQUEEZE_RIGHT , 1)
521    
522    !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
523    
       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)  
524    
525        do iter = 1, streamice_max_cg_iter        do iter = 1, streamice_max_cg_iter
526         if (resid .gt. tolerance*resid_0) then         if (resid .gt. tolerance*resid_0) then
# Line 219  c      to avoid using "exit" Line 539  c      to avoid using "exit"
539            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
540             Au_SI(i,j,bi,bj) = 0. _d 0             Au_SI(i,j,bi,bj) = 0. _d 0
541             Av_SI(i,j,bi,bj) = 0. _d 0             Av_SI(i,j,bi,bj) = 0. _d 0
 !            Du_SI(i,j,bi,bj) = Real(i)  
 !            Dv_SI(i,j,bi,bj) = 0.0  
542            ENDDO            ENDDO
543           ENDDO           ENDDO
544          ENDDO          ENDDO
545         ENDDO         ENDDO
546    
547         IF (STREAMICE_construct_matrix) THEN  !        IF (STREAMICE_construct_matrix) THEN
548    
549    ! #ifdef STREAMICE_CONSTRUCT_MATRIX
550    
551          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
552           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
553            DO j=js,je            DO j=js,je
# Line 234  c      to avoid using "exit" Line 555  c      to avoid using "exit"
555              DO colx=-1,1              DO colx=-1,1
556               DO coly=-1,1               DO coly=-1,1
557                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
558       &         streamice_cg_A1(i,j,bi,bj,colx,coly)*       &         A_uu(i,j,bi,bj,colx,coly)*
559       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
560       &         streamice_cg_A2(i,j,bi,bj,colx,coly)*           &         A_uv(i,j,bi,bj,colx,coly)*    
561       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
562                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
563       &         streamice_cg_A3(i,j,bi,bj,colx,coly)*       &         A_vu(i,j,bi,bj,colx,coly)*
564       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
565       &         streamice_cg_A4(i,j,bi,bj,colx,coly)*           &         A_vv(i,j,bi,bj,colx,coly)*    
566       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
567               ENDDO               ENDDO
568              ENDDO              ENDDO
# Line 250  c      to avoid using "exit" Line 571  c      to avoid using "exit"
571           ENDDO           ENDDO
572          ENDDO          ENDDO
573    
574         else  !        else
575    ! #else
576          CALL STREAMICE_CG_ACTION( myThid,  !
577       O     Au_SI,  !         CALL STREAMICE_CG_ACTION( myThid,
578       O     Av_SI,  !      O     Au_SI,
579       I     Du_SI,  !      O     Av_SI,
580       I     Dv_SI,  !      I     Du_SI,
581       I     is,ie,js,je)  !      I     Dv_SI,
582    !      I     is,ie,js,je)
583         ENDIF  !
584    ! !        ENDIF
585    !
586    ! #endif
587                
 ! !        if (iter.eq.1) then  
 ! !         CALL WRITE_FLD_XY_RL ("Au2","",Au_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL ("Av2","",Av_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL ("Du","",Du_SI,0,myThid)  
 ! !          
 ! !         CALL WRITE_FLD_XY_RL("Dv","",Au_SI,0,myThid)  
 ! !          
 ! !         CALL WRITE_FLD_XY_RL("DiagU1","",Diagu_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagV1","",Diagv_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagU2","",  
 ! !      &    streamice_cg_A1(i,j,bi,bj,0,0),0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagV2","",  
 ! !      &    streamice_cg_A4(i,j,bi,bj,0,0),0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagV2","",Diagv_SI,0,myThid)  
 ! !  
 ! !        endif  
   
   
   
   
   
 ! !        if (iter.eq.1) then  
 ! !         CALL WRITE_FLD_XY_RL ("Au1","",Au_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL ("Av1","",Av_SI,0,myThid)  
 ! !          
 ! !        endif  
   
   
588    
589         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
590          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 443  c      to avoid using "exit" Line 739  c      to avoid using "exit"
739          _EXCH_XY_RL( cg_Vin, myThid )          _EXCH_XY_RL( cg_Vin, myThid )
740         endif         endif
741    
742    
743         endif         endif
744        enddo ! end of CG loop        enddo ! end of CG loop
745                
# Line 453  c     if iters has reached max_iters the Line 750  c     if iters has reached max_iters the
750         conv_flag = 1         conv_flag = 1
751        ENDIF        ENDIF
752    
753        DO bj = myByLo(myThid), myByHi(myThid)  !       DO bj = myByLo(myThid), myByHi(myThid)
754         DO bi = myBxLo(myThid), myBxHi(myThid)  !        DO bi = myBxLo(myThid), myBxHi(myThid)
755          DO j=1-OLy,sNy+OLy  !         DO j=1-OLy,sNy+OLy
756           DO i=1-OLy,sNx+OLy  !          DO i=1-OLy,sNx+OLy
757            IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)  !           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)  !      &     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)  !           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)  !      &     cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
761           ENDDO  !          ENDDO
762          ENDDO  !         ENDDO
763         ENDDO  !        ENDDO
764        ENDDO        !       ENDDO      
765    !
766        _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
767        _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
768    
769    #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  #endif  #endif
       RETURN  
       END  
   
782    
783                CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
   
   
   
784    
785    
786    #endif
787          RETURN
788          END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22