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