/[MITgcm]/MITgcm/pkg/diagnostics/diag_cg2d.F
ViewVC logotype

Diff of /MITgcm/pkg/diagnostics/diag_cg2d.F

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

revision 1.2 by jmc, Fri Jul 22 19:53:39 2011 UTC revision 1.3 by jmc, Mon Aug 1 20:40:43 2011 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "DIAG_OPTIONS.h"  #include "DIAG_OPTIONS.h"
5    #undef DEBUG_DIAG_CG2D
6    
7  CBOP  CBOP
8  C     !ROUTINE: DIAG_CG2D  C     !ROUTINE: DIAG_CG2D
# Line 34  C     === Global data === Line 35  C     === Global data ===
35  #include "SIZE.h"  #include "SIZE.h"
36  #include "EEPARAMS.h"  #include "EEPARAMS.h"
37  #include "PARAMS.h"  #include "PARAMS.h"
 c#include "CG2D.h"  
 c#include "GRID.h"  
 c#include "SURFACE.h"  
38    
39  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
40  C     === Routine arguments ===  C     === Routine arguments ===
# Line 85  C     i, j, it2d :: Loop counters ( it2d Line 83  C     i, j, it2d :: Loop counters ( it2d
83        _RL    alpha,  alphaTile(nSx,nSy)        _RL    alpha,  alphaTile(nSx,nSy)
84        _RL    sumRHS, sumRHStile(nSx,nSy)        _RL    sumRHS, sumRHStile(nSx,nSy)
85        _RL    pW_tmp, pS_tmp        _RL    pW_tmp, pS_tmp
86          _RL    diagCG_pcOffDFac
87        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
88        LOGICAL printResidual        LOGICAL printResidual
89  CEOP  CEOP
# Line 93  CEOP Line 92  CEOP
92        _RS  pS  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  pS  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
93        _RS  pC  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  pC  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
94        _RL  q2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)        _RL  q2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
95    #ifdef DEBUG_DIAG_CG2D
96          CHARACTER*(10) sufx
97          _RL  r2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
98    #else
99        _RL  r2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)        _RL  r2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
100    #endif
101        _RL  s2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)        _RL  s2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
102        _RL  x2dm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  x2dm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
103    
# Line 112  C--   Set matrice main diagnonal: Line 116  C--   Set matrice main diagnonal:
116        CALL EXCH_XY_RS(aC2d, myThid)        CALL EXCH_XY_RS(aC2d, myThid)
117    
118  C--   Initialise preconditioner  C--   Initialise preconditioner
119          diagCG_pcOffDFac = cg2dpcOffDFac
120        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
121         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
122          DO j=1,sNy+1          DO j=1,sNy+1
# Line 126  C--   Initialise preconditioner Line 131  C--   Initialise preconditioner
131             pW(i,j,bi,bj) = 0.             pW(i,j,bi,bj) = 0.
132            ELSE            ELSE
133             pW(i,j,bi,bj) = -aW2d(i,j,bi,bj)             pW(i,j,bi,bj) = -aW2d(i,j,bi,bj)
134       &                   /( (cg2dpcOffDFac*pW_tmp)**2 )       &                   /( (diagCG_pcOffDFac*pW_tmp)**2 )
135            ENDIF            ENDIF
136            pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)            pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
137            IF ( pS_tmp .EQ. 0. ) THEN            IF ( pS_tmp .EQ. 0. ) THEN
138             pS(i,j,bi,bj) = 0.             pS(i,j,bi,bj) = 0.
139            ELSE            ELSE
140             pS(i,j,bi,bj) = -aS2d(i,j,bi,bj)             pS(i,j,bi,bj) = -aS2d(i,j,bi,bj)
141       &                   /( (cg2dpcOffDFac*pS_tmp)**2 )       &                   /( (diagCG_pcOffDFac*pS_tmp)**2 )
142            ENDIF            ENDIF
143    c         pC(i,j,bi,bj) = 1.
144    c         pW(i,j,bi,bj) = 0.
145    c         pS(i,j,bi,bj) = 0.
146           ENDDO           ENDDO
147          ENDDO          ENDDO
148         ENDDO         ENDDO
# Line 172  C--   Initial residual calculation Line 180  C--   Initial residual calculation
180          ENDDO          ENDDO
181         ENDDO         ENDDO
182        ENDDO        ENDDO
183    #ifdef DEBUG_DIAG_CG2D
184          CALL EXCH_XY_RL ( r2d, myThid )
185    #else
186        CALL EXCH_S3D_RL( r2d, 1, myThid )        CALL EXCH_S3D_RL( r2d, 1, myThid )
187    #endif
188        CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )        CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
189        CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )        IF ( printResidFrq.GE.1 )
190         &  CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
191        err = SQRT(err)        err = SQRT(err)
192        it2d = 0        it2d = 0
193        firstResidual = err        firstResidual = err
# Line 185  C--   Initial residual calculation Line 198  C--   Initial residual calculation
198        IF ( debugLevel .GE. debLevZero ) THEN        IF ( debugLevel .GE. debLevZero ) THEN
199          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
200          printResidual = printResidFrq.GE.1          printResidual = printResidFrq.GE.1
201  c       WRITE(standardmessageunit,'(A,1P2E22.14)')          IF ( printResidual ) THEN
202  c    &  ' diag_cg2d: Sum(rhs),rhsMax = ', sumRHS, rhsMax           WRITE(msgBuf,'(2A,I6,A,1PE17.9,A,1PE14.6)')' diag_cg2d:',
203         &    ' iter=', it2d, ' ; resid.=', err, ' ; sumRHS=', sumRHS
204             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
205         &                       SQUEEZE_RIGHT, myThid )
206            ENDIF
207          _END_MASTER( myThid )          _END_MASTER( myThid )
208        ENDIF        ENDIF
209    #ifdef DEBUG_DIAG_CG2D
210          IF ( printResidFrq.GE.1 ) THEN
211            WRITE(sufx,'(I10.10)') 0
212            CALL WRITE_FLD_XY_RL( 'r2d.',sufx, r2d, 1, myThid )
213          ENDIF
214    #endif
215    
216        IF ( err .LT. residCriter ) GOTO 11        IF ( err .LT. residCriter ) GOTO 11
217    
# Line 275  C      Now compute "interior" points. Line 298  C      Now compute "interior" points.
298         err = SQRT(err)         err = SQRT(err)
299         IF ( printResidual ) THEN         IF ( printResidual ) THEN
300          IF ( MOD( it2d-1, printResidFrq ).EQ.0 ) THEN          IF ( MOD( it2d-1, printResidFrq ).EQ.0 ) THEN
301           WRITE(msgBuf,'(A,I6,A,1PE21.14)')           WRITE(msgBuf,'(A,I6,A,1PE17.9)')
302       &    ' diag_cg2d: iter=', it2d, ' ; resid.= ', err       &    ' diag_cg2d: iter=', it2d, ' ; resid.=', err
303           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
304       &                       SQUEEZE_RIGHT, myThid )       &                       SQUEEZE_RIGHT, myThid )
305          ENDIF          ENDIF
# Line 297  C-     Store lowest residual solution Line 320  C-     Store lowest residual solution
320           ENDDO           ENDDO
321         ENDIF         ENDIF
322    
323    #ifdef DEBUG_DIAG_CG2D
324           CALL EXCH_XY_RL( r2d, myThid )
325           IF ( printResidFrq.GE.1 ) THEN
326            WRITE(sufx,'(I10.10)') it2d
327            CALL WRITE_FLD_XY_RL( 'r2d.',sufx, r2d, 1, myThid )
328            CALL WRITE_FLD_XY_RL( 'x2d.',sufx, x2d, 1, myThid )
329           ENDIF
330    #else
331         CALL EXCH_S3D_RL( r2d, 1, myThid )         CALL EXCH_S3D_RL( r2d, 1, myThid )
332    #endif
333    
334     10 CONTINUE     10 CONTINUE
335        it2d = numIters        it2d = numIters
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22