/[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.3 by dgoldberg, Mon May 14 16:52:29 2012 UTC revision 1.9 by dgoldberg, Sat Jun 8 22:15:33 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    ! UNCOMMENT IF V3.0
40    !#include "finclude/petscvec.h"
41    !#include "finclude/petscmat.h"
42    !#include "finclude/petscksp.h"
43    !#include "finclude/petscpc.h"
44    #endif
45    C     === Global variables ===
46    
47                
48  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
49  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 55  C     cg_Bu, cg_Bv - driving stress
55        _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)
56        _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)
57        _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)
58          _RL
59         & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
60         & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
61         & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
62         & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
63    
64  C     LOCAL VARIABLES  C     LOCAL VARIABLES
65        INTEGER i, j, bi, bj, cg_halo, conv_flag        INTEGER i, j, bi, bj, cg_halo, conv_flag
66        INTEGER iter, is, js, ie, je, colx, coly, k        INTEGER iter, is, js, ie, je, colx, coly, k
67        _RL dot_p1, dot_p2, resid_0, alpha_k, beta_k, resid        _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
68        _RL dot_p1_tile (nSx,nSy)        _RL dot_p1_tile (nSx,nSy)
69        _RL dot_p2_tile (nSx,nSy)        _RL dot_p2_tile (nSx,nSy)
70          CHARACTER*(MAX_LEN_MBUF) msgBuf
71    
72    
73    #ifdef ALLOW_PETSC
74          INTEGER indices(2*(snx*nsx*sny*nsy))
75          INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
76          _RL rhs_values(2*(snx*nsx*sny*nsy))
77          _RL solution_values(2*(snx*nsx*sny*nsy))
78    !      _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))
79          _RL mat_values (18,1), mat_val_return(1)
80          INTEGER indices_col(18)
81          INTEGER local_dofs, global_dofs, dof_index, dof_index_col
82          INTEGER local_offset
83          Mat matrix
84          KSP ksp
85          PC  pc
86          Vec rhs
87          Vec solution
88          PetscErrorCode ierr
89    #ifdef ALLOW_USE_MPI
90          integer mpiRC, mpiMyWid
91    #endif
92    #endif
93    
       iters = streamice_max_cg_iter  
94    
95  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
96    
97    
98    
99          CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
100    #ifndef STREAMICE_SERIAL_TRISOLVE
101    
102    #ifdef ALLOW_PETSC
103    
104    #ifdef ALLOW_USE_MPI
105    
106    
107          CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
108          local_dofs = n_dofs_process (mpiMyWid)
109          global_dofs = 0
110          
111          n_dofs_cum_sum(0) = 0
112          DO i=0,nPx*nPy-1
113           global_dofs = global_dofs + n_dofs_process (i)
114           if (i.ge.1) THEN
115             n_dofs_cum_sum(i) = n_dofs_cum_sum(i-1)+
116         &                     n_dofs_process(i-1)
117           endif
118          ENDDO
119          local_offset = n_dofs_cum_sum(mpimywid)
120    
121    #else
122    
123          local_dofs = n_dofs_process (0)
124          global_dofs = local_dofs
125          local_offset = 0
126    
127    #endif
128    
129    !      call petscInitialize(PETSC_NULL_CHARACTER,ierr)
130    
131    !----------------------
132    
133          call VecCreate(PETSC_COMM_WORLD, rhs, ierr)
134          call VecSetSizes(rhs, local_dofs, global_dofs, ierr)
135          call VecSetType(rhs, VECMPI, ierr)
136    
137          call VecCreate(PETSC_COMM_WORLD, solution, ierr)
138          call VecSetSizes(solution, local_dofs, global_dofs, ierr)
139          call VecSetType(solution, VECMPI, ierr)
140    
141          do i=1,local_dofs
142            indices(i) = i-1 + local_offset
143          end do
144          do i=1,2*nSx*nSy*sNx*sNy
145            rhs_values (i) = 0. _d 0
146            solution_values (i) = 0. _d 0
147          enddo
148    
149    ! gather rhs and initial guess values to populate petsc vectors
150    
151          DO bj = myByLo(myThid), myByHi(myThid)
152           DO bi = myBxLo(myThid), myBxHi(myThid)
153            DO j=1,sNy
154             DO i=1,sNx
155    
156              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
157         &                - local_offset
158    
159              if (dof_index.ge.0) THEN
160              
161               rhs_values(dof_index+1) = cg_Bu(i,j,bi,bj)
162               solution_values(dof_index+1) = cg_Uin(i,j,bi,bj)
163    
164              endif
165    
166    !---------------
167    
168              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
169         &                - local_offset
170    
171              if (dof_index.ge.0) THEN
172    
173               rhs_values(dof_index+1) = cg_Bv(i,j,bi,bj)
174               solution_values(dof_index+1) = cg_Vin(i,j,bi,bj)
175    
176              endif
177    
178             ENDDO
179            ENDDO
180           ENDDO
181          ENDDO
182    
183    
184          call VecSetValues(rhs, local_dofs, indices, rhs_values,
185         &                  INSERT_VALUES, ierr)
186          call VecAssemblyBegin(rhs, ierr)
187          call VecAssemblyEnd(rhs, ierr)
188    
189    
190          call VecSetValues(solution, local_dofs, indices,
191         &                  solution_values, INSERT_VALUES, ierr)
192          call VecAssemblyBegin(solution, ierr)
193          call VecAssemblyEnd(solution, ierr)
194    
195    !     IF USING v3.0 THEN
196    !     call MatCreateMPIAIJ (PETSC_COMM_WORLD,
197          call MatCreateAIJ (PETSC_COMM_WORLD,
198         &                      local_dofs, local_dofs,
199         &                      global_dofs, global_dofs,
200         &                      18, PETSC_NULL_INTEGER,
201         &                      18, PETSC_NULL_INTEGER,
202         &                      matrix, ierr)
203    
204    
205    ! populate petsc matrix
206    
207          DO bj = myByLo(myThid), myByHi(myThid)
208           DO bi = myBxLo(myThid), myBxHi(myThid)
209            DO j=1,sNy
210             DO i=1,sNx
211    
212              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
213    !     &                - local_offset
214    
215              IF (dof_index .ge. 0) THEN
216    
217               DO k=1,18
218                indices_col(k) = 0
219                mat_values(k,1) = 0. _d 0
220               ENDDO
221               k=0
222    
223               DO coly=-1,1
224                DO colx=-1,1
225    
226                 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
227    
228                 if (dof_index_col.ge.0) THEN
229    !               pscal = A_uu(i,j,bi,bj,colx,coly)
230    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
231    !     &              pscal,INSERT_VALUES,ierr)
232                    k=k+1
233                    mat_values (k,1) = A_uu(i,j,bi,bj,colx,coly)
234                    indices_col (k) = dof_index_col
235                 endif
236                  
237                 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
238    
239                 if (dof_index_col.ge.0) THEN
240    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
241    !     &              A_uv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
242                    k=k+1
243                    mat_values (k,1) = A_uv(i,j,bi,bj,colx,coly)
244                    indices_col (k) = dof_index_col
245                 endif
246                
247                ENDDO
248               ENDDO
249    
250               call matSetValues (matrix, 1, dof_index, k, indices_col,
251         &                        mat_values,INSERT_VALUES,ierr)
252    
253    
254              ENDIF
255    
256    ! ----------------------------------------------
257    
258              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
259    !     &                - local_offset
260    
261              IF (dof_index .ge. 0) THEN
262    
263               DO k=1,18
264                indices_col(k) = 0
265                mat_values(k,1) = 0. _d 0
266               ENDDO
267               k=0
268    
269               DO coly=-1,1
270                DO colx=-1,1
271    
272                 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
273    
274                 if (dof_index_col.ge.0) THEN
275    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
276    !     &              A_vu(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
277                    k=k+1
278                    mat_values (k,1) = A_vu(i,j,bi,bj,colx,coly)
279                    indices_col (k) = dof_index_col
280                 endif
281    
282                 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
283    
284                 if (dof_index_col.ge.0) THEN
285    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
286    !     &              A_vv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
287                    k=k+1
288                    mat_values (k,1) = A_vv(i,j,bi,bj,colx,coly)
289                    indices_col (k) = dof_index_col
290                 endif
291    
292                ENDDO
293               ENDDO
294    
295               call matSetValues (matrix, 1, dof_index, k, indices_col,
296         &                        mat_values,INSERT_VALUES,ierr)
297              ENDIF
298    
299             ENDDO
300            ENDDO
301           ENDDO
302          ENDDO
303    
304          call MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY,ierr)
305          call MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY,ierr)
306    
307    
308          call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
309          call KSPSetOperators(ksp, matrix, matrix,
310         &                     DIFFERENT_NONZERO_PATTERN, ierr)
311    
312          SELECT CASE (PETSC_SOLVER_TYPE)
313           CASE ('CG')
314           PRINT *, "PETSC SOLVER: SELECTED CG"
315           call KSPSetType(ksp, KSPCG, ierr)
316           CASE ('GMRES')
317           PRINT *, "PETSC SOLVER: SELECTED GMRES"
318           call KSPSetType(ksp, KSPGMRES, ierr)
319           CASE ('BICG')
320           PRINT *, "PETSC SOLVER: SELECTED BICG"
321           call KSPSetType(ksp, KSPBICG, ierr)
322           CASE DEFAULT
323           PRINT *, "PETSC SOLVER: SELECTED DEFAULT"
324           call KSPSetType(ksp, KSPCG, ierr)
325          END SELECT
326    
327          call KSPGetPC(ksp, pc, ierr)
328          call KSPSetTolerances(ksp,tolerance,
329         &     PETSC_DEFAULT_DOUBLE_PRECISION,    
330         &     PETSC_DEFAULT_DOUBLE_PRECISION,
331         &     streamice_max_cg_iter,ierr)
332    
333          SELECT CASE (PETSC_PRECOND_TYPE)
334           CASE ('BLOCKJACOBI')
335           PRINT *, "PETSC PRECOND: SELECTED BJACOBI"
336           call PCSetType(pc, PCBJACOBI, ierr)
337           CASE ('JACOBI')
338           PRINT *, "PETSC PRECOND: SELECTED JACOBI"
339           call PCSetType(pc, PCJACOBI, ierr)
340           CASE ('ILU')
341           PRINT *, "PETSC PRECOND: SELECTED ILU"
342           call PCSetType(pc, PCILU, ierr)
343           CASE DEFAULT
344           PRINT *, "PETSC PRECOND: SELECTED DEFAULT"
345           call PCSetType(pc, PCBJACOBI, ierr)
346          END SELECT
347    
348    
349          call KSPSolve(ksp, rhs, solution, ierr)
350          call KSPGetIterationNumber(ksp,iters,ierr)
351    
352          call VecGetValues(solution,local_dofs,indices,
353         &      solution_values,ierr)
354    
355          DO bj = myByLo(myThid), myByHi(myThid)
356           DO bi = myBxLo(myThid), myBxHi(myThid)
357            DO j=1,sNy
358             DO i=1,sNx
359    
360              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
361         &                - local_offset
362              if (dof_index.ge.0) THEN
363               cg_Uin(i,j,bi,bj) = solution_values(dof_index+1)
364              endif
365              
366              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
367         &                - local_offset
368              if (dof_index.ge.0) THEN
369               cg_Vin(i,j,bi,bj) = solution_values(dof_index+1)
370              endif
371    
372             ENDDO
373            ENDDO
374           ENDDO
375          ENDDO
376    
377    
378    
379    #else  /* ALLOW_PETSC */
380    
381    
382          iters = streamice_max_cg_iter
383        conv_flag = 0        conv_flag = 0
384    
385    
386        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
387         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
388          DO j=1,sNy          DO j=1,sNy
389           DO i=1,sNx           DO i=1,sNx
390            Zu_SI (i,j,bi,bj) = 0. _d 0            Zu_SI (i,j,bi,bj) = 0. _d 0
391            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  
392            Ru_SI (i,j,bi,bj) = 0. _d 0            Ru_SI (i,j,bi,bj) = 0. _d 0
393            Rv_SI (i,j,bi,bj) = 0. _d 0            Rv_SI (i,j,bi,bj) = 0. _d 0
394            Au_SI (i,j,bi,bj) = 0. _d 0            Au_SI (i,j,bi,bj) = 0. _d 0
395            Av_SI (i,j,bi,bj) = 0. _d 0            Av_SI (i,j,bi,bj) = 0. _d 0
396            Du_SI (i,j,bi,bj) = 0. _d 0            Du_SI (i,j,bi,bj) = 0. _d 0
397            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  
398           ENDDO           ENDDO
399          ENDDO          ENDDO
400         ENDDO         ENDDO
401        ENDDO        ENDDO
402          
403    C     FIND INITIAL RESIDUAL, and initialize r
404    
405  C     PREAMBLE: get bdry contribution, find matrix diagonal,  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
 C     initialize iterates R (residual), Z, and D  
406    
407        CALL STREAMICE_CG_BOUND_VALS( myThid,              
      O    ubd_SI,  
      O    vbd_SI)  
408    
409        DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
410         DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
411          DO j=1-OLy,sNy+OLy            DO j=1,sNy
412           DO i=1-OLx,sNx+OLx             DO i=1,sNx
413            RHSu_SI (i,j,bi,bj) = cg_Bu(i,j,bi,bj)              DO colx=-1,1
414       &     - ubd_SI(i,j,bi,bj)               DO coly=-1,1
415            RHSv_SI (i,j,bi,bj) = cg_Bv(i,j,bi,bj)                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
416       &     - vbd_SI(i,j,bi,bj)       &         A_uu(i,j,bi,bj,colx,coly)*
417           ENDDO       &         cg_Uin(i+colx,j+coly,bi,bj)+
418          ENDDO       &         A_uv(i,j,bi,bj,colx,coly)*    
419         ENDDO       &         cg_Vin(i+colx,j+coly,bi,bj)
       ENDDO  
         
       _EXCH_XY_RL( RHSu_SI, myThid )  
       _EXCH_XY_RL( RHSv_SI, myThid )  
420    
       CALL STREAMICE_CG_ADIAG( myThid,  
      O    DIAGu_SI,  
      O    DIAGv_SI)  
421    
422        _EXCH_XY_RL( DIAGu_SI, myThid )                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
423        _EXCH_XY_RL( DIAGv_SI, myThid )       &         A_vu(i,j,bi,bj,colx,coly)*
424               &         cg_Uin(i+colx,j+coly,bi,bj)+
425  C     FIX PROBLEM WITH PRECOND LATER       &         A_vv(i,j,bi,bj,colx,coly)*    
426         &         cg_Vin(i+colx,j+coly,bi,bj)
427  !       DO bj = myByLo(myThid), myByHi(myThid)               ENDDO
428  !        DO bi = myBxLo(myThid), myBxHi(myThid)              ENDDO
429  !         DO j=1-OLy,sNy+OLy             ENDDO
430  !          DO i=1-OLx,sNx+OLx            ENDDO
431  !           DIAGu_SI(i,j,bi,bj)=1.0           ENDDO
432  !           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 )  
433    
434                
435        _EXCH_XY_RL( Au_SI, myThid )        _EXCH_XY_RL( Au_SI, myThid )
436        _EXCH_XY_RL( Av_SI, myThid )        _EXCH_XY_RL( Av_SI, myThid )
437    
438    
439        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
440         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
441          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
442           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
443            Ru_SI(i,j,bi,bj)=RHSu_SI(i,j,bi,bj)-            Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
444       &     Au_SI(i,j,bi,bj)       &     Au_SI(i,j,bi,bj)
445            Rv_SI(i,j,bi,bj)=RHSv_SI(i,j,bi,bj)-            Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
446       &     Av_SI(i,j,bi,bj)       &     Av_SI(i,j,bi,bj)
447           ENDDO           ENDDO
448          ENDDO          ENDDO
# Line 143  C     FIX PROBLEM WITH PRECOND LATER Line 451  C     FIX PROBLEM WITH PRECOND LATER
451         ENDDO         ENDDO
452        ENDDO        ENDDO
453    
454          
455    
456        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
457         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
458          DO j=1,sNy          DO j=1,sNy
# Line 159  C     FIX PROBLEM WITH PRECOND LATER Line 469  C     FIX PROBLEM WITH PRECOND LATER
469        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
470        resid_0 = sqrt(dot_p1)        resid_0 = sqrt(dot_p1)
471    
472          WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
473         &         resid_0
474           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
475         &                     SQUEEZE_RIGHT , 1)
476    
477    C    CCCCCCCCCCCCCCCCCCCC
478    
479        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
480         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
481          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 172  C     FIX PROBLEM WITH PRECOND LATER Line 489  C     FIX PROBLEM WITH PRECOND LATER
489         ENDDO         ENDDO
490        ENDDO        ENDDO
491    
   
492        cg_halo = min(OLx-1,OLy-1)        cg_halo = min(OLx-1,OLy-1)
493        conv_flag = 0        conv_flag = 0
494    
# Line 195  c  !!              !! Line 511  c  !!              !!
511  c  !! MAIN CG LOOP !!  c  !! MAIN CG LOOP !!
512  c  !!              !!  c  !!              !!
513  c  !!!!!!!!!!!!!!!!!!  c  !!!!!!!!!!!!!!!!!!
514    
515    
516        
517        
518  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!!)
519    
520        print *, "BEGINNING MAIN CG LOOP"        WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
521           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
522         &                     SQUEEZE_RIGHT , 1)
523    
524  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
 #ifdef STREAMICE_CONSTRUCT_MATRIX  
   
       CALL STREAMICE_CG_MAKE_A(myThid)  
525    
 #endif  
526    
527        do iter = 1, streamice_max_cg_iter        do iter = 1, streamice_max_cg_iter
528         if (resid .gt. tolerance*resid_0) then         if (resid .gt. tolerance*resid_0) then
# Line 224  c      to avoid using "exit" Line 541  c      to avoid using "exit"
541            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
542             Au_SI(i,j,bi,bj) = 0. _d 0             Au_SI(i,j,bi,bj) = 0. _d 0
543             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  
544            ENDDO            ENDDO
545           ENDDO           ENDDO
546          ENDDO          ENDDO
# Line 233  c      to avoid using "exit" Line 548  c      to avoid using "exit"
548    
549  !        IF (STREAMICE_construct_matrix) THEN  !        IF (STREAMICE_construct_matrix) THEN
550    
551  #ifdef STREAMICE_CONSTRUCT_MATRIX  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
552    
553          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
554           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 242  c      to avoid using "exit" Line 557  c      to avoid using "exit"
557              DO colx=-1,1              DO colx=-1,1
558               DO coly=-1,1               DO coly=-1,1
559                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
560       &         streamice_cg_A1(i,j,bi,bj,colx,coly)*       &         A_uu(i,j,bi,bj,colx,coly)*
561       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
562       &         streamice_cg_A2(i,j,bi,bj,colx,coly)*           &         A_uv(i,j,bi,bj,colx,coly)*    
563       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
564                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
565       &         streamice_cg_A3(i,j,bi,bj,colx,coly)*       &         A_vu(i,j,bi,bj,colx,coly)*
566       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
567       &         streamice_cg_A4(i,j,bi,bj,colx,coly)*           &         A_vv(i,j,bi,bj,colx,coly)*    
568       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
569               ENDDO               ENDDO
570              ENDDO              ENDDO
# Line 259  c      to avoid using "exit" Line 574  c      to avoid using "exit"
574          ENDDO          ENDDO
575    
576  !        else  !        else
577  #else  ! #else
578    !
579          CALL STREAMICE_CG_ACTION( myThid,  !         CALL STREAMICE_CG_ACTION( myThid,
580       O     Au_SI,  !      O     Au_SI,
581       O     Av_SI,  !      O     Av_SI,
582       I     Du_SI,  !      I     Du_SI,
583       I     Dv_SI,  !      I     Dv_SI,
584       I     is,ie,js,je)  !      I     is,ie,js,je)
585    !
586  !        ENDIF  ! !        ENDIF
587    !
588  #endif  ! #endif
589                
 ! !        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  
   
   
590    
591         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
592          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 454  c      to avoid using "exit" Line 741  c      to avoid using "exit"
741          _EXCH_XY_RL( cg_Vin, myThid )          _EXCH_XY_RL( cg_Vin, myThid )
742         endif         endif
743    
744    
745         endif         endif
746        enddo ! end of CG loop        enddo ! end of CG loop
747                
# Line 464  c     if iters has reached max_iters the Line 752  c     if iters has reached max_iters the
752         conv_flag = 1         conv_flag = 1
753        ENDIF        ENDIF
754    
755        DO bj = myByLo(myThid), myByHi(myThid)  !       DO bj = myByLo(myThid), myByHi(myThid)
756         DO bi = myBxLo(myThid), myBxHi(myThid)  !        DO bi = myBxLo(myThid), myBxHi(myThid)
757          DO j=1-OLy,sNy+OLy  !         DO j=1-OLy,sNy+OLy
758           DO i=1-OLy,sNx+OLy  !          DO i=1-OLy,sNx+OLy
759            IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
760       &     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)
761            IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
762       &     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)
763           ENDDO  !          ENDDO
764          ENDDO  !         ENDDO
765         ENDDO  !        ENDDO
766        ENDDO        !       ENDDO      
767    !
768        _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
769        _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
770    
771    #endif /* ifndef ALLOW_PETSC */
772    
773    #else /* STREAMICE_SERIAL_TRISOLVE */
774    
775          CALL STREAMICE_TRIDIAG_SOLVE(
776         U                               cg_Uin,     ! x-velocities
777         U                               cg_Vin,
778         U                               cg_Bu,      ! force in x dir
779         I                               A_uu,       ! section of matrix that multiplies u and projects on u
780         I                               STREAMICE_umask,
781         I                               myThid )
782    
783  #endif  #endif
       RETURN  
       END  
   
784    
785                CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
   
   
   
786    
787    
788    #endif
789          RETURN
790          END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22