/[MITgcm]/MITgcm/model/src/cg3d.F
ViewVC logotype

Diff of /MITgcm/model/src/cg3d.F

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

revision 1.16 by jmc, Tue Nov 8 06:18:10 2005 UTC revision 1.17 by jmc, Tue Dec 20 20:31:28 2005 UTC
# Line 50  C     === Global data === Line 50  C     === Global data ===
50  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
51  C     === Routine arguments ===  C     === Routine arguments ===
52  C     myThid    - Thread on which I am working.  C     myThid    - Thread on which I am working.
53  C     cg2d_b    - The source term or "right hand side"  C     cg3d_b    - The source term or "right hand side"
54  C     cg2d_x    - The solution  C     cg3d_x    - The solution
55  C     firstResidual - the initial residual before any iterations  C     firstResidual - the initial residual before any iterations
56  C     lastResidual  - the actual residual reached  C     lastResidual  - the actual residual reached
57  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  - Entry: the maximum number of iterations allowed
# Line 81  C                   useful debuggin/trou Line 81  C                   useful debuggin/trou
81  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
82  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
83  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
84  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, K, N  - Loop counters ( N counts CG iterations )
85        INTEGER actualIts        INTEGER actualIts
86        _RL    actualResidual        _RL    actualResidual
87        INTEGER bi, bj                      INTEGER bi, bj              
88        INTEGER I, J, K, it3d        INTEGER I, J, K, it3d
89        INTEGER KM1, KP1        INTEGER Km1, Kp1
90          _RL    maskM1, maskP1
91        _RL    err, errTile        _RL    err, errTile
92        _RL    eta_qrN, eta_qrNtile        _RL    eta_qrN, eta_qrNtile
93        _RL    eta_qrNM1        _RL    eta_qrNM1
# Line 95  C     I, J, N     - Loop counters ( N co Line 96  C     I, J, N     - Loop counters ( N co
96        _RL    sumRHS, sumRHStile        _RL    sumRHS, sumRHStile
97        _RL    rhsMax        _RL    rhsMax
98        _RL    rhsNorm        _RL    rhsNorm
       _RL    topLevTerm  
99  CEOP  CEOP
100    
101    
# Line 110  C--   Normalise RHS Line 110  C--   Normalise RHS
110           DO J=1,sNy           DO J=1,sNy
111            DO I=1,sNx            DO I=1,sNx
112             cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*cg3dNorm             cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*cg3dNorm
113         &                          * maskC(I,J,K,bi,bj)
114             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)
115            ENDDO            ENDDO
116           ENDDO           ENDDO
# Line 144  C--   Initial residual calculation (with Line 145  C--   Initial residual calculation (with
145          errTile    = 0. _d 0          errTile    = 0. _d 0
146          sumRHStile = 0. _d 0          sumRHStile = 0. _d 0
147          DO K=1,Nr          DO K=1,Nr
148           KM1 = K-1           Km1 = MAX(K-1, 1 )
149           IF ( K .EQ. 1 ) KM1 = 1           Kp1 = MIN(K+1, Nr)
150           KP1 = K+1           maskM1 = 1. _d 0
151           IF ( K .EQ. Nr ) KP1 = 1           maskP1 = 1. _d 0
152           topLevTerm = 0.           IF ( K .EQ. 1 ) maskM1 = 0. _d 0
153           IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*           IF ( K .EQ. Nr) maskP1 = 0. _d 0
154       &      (horiVertRatio/gravity)/deltaTMom/deltaTMom            
155           DO J=1,sNy           DO J=1,sNy
156            DO I=1,sNx            DO I=1,sNx
157             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.
# Line 158  C--   Initial residual calculation (with Line 159  C--   Initial residual calculation (with
159       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)
160       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)
161       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)
162       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,KM1,bi,bj)       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,Km1,bi,bj)*maskM1
163       &     +aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,KP1,bi,bj)       &     +aV3d(I  ,J  ,Kp1,bi,bj)*cg3d_x(I  ,J  ,Kp1,bi,bj)*maskP1
164       &     -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)
      &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)  
165       &     )       &     )
166             errTile = errTile             errTile = errTile
167       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
# Line 184  C--   Initial residual calculation (with Line 179  C--   Initial residual calculation (with
179          sumRHS = sumRHS + sumRHStile          sumRHS = sumRHS + sumRHStile
180         ENDDO         ENDDO
181        ENDDO        ENDDO
 C     _EXCH_XYZ_R8( cg3d_r, myThid )  
182         CALL EXCH_S3D_RL( cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, myThid )
 C     _EXCH_XYZ_R8( cg3d_s, myThid )  
183  c      CALL EXCH_S3D_RL( cg3d_s, myThid )  c      CALL EXCH_S3D_RL( cg3d_s, myThid )
184        _GLOBAL_SUM_R8( sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHS, myThid )
185        _GLOBAL_SUM_R8( err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
# Line 208  c     _END_MASTER( myThid ) Line 201  c     _END_MASTER( myThid )
201        firstResidual=actualResidual        firstResidual=actualResidual
202    
203  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
204        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, numIters
205    
206  CcnhDebugStarts  CcnhDebugStarts
207  c      IF ( mod(it3d-1,10).EQ.0)  c      IF ( mod(it3d-1,10).EQ.0)
# Line 320  CcnhDebugEnds Line 313  CcnhDebugEnds
313  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
314  C==    q = A.s  C==    q = A.s
315         alpha = 0. _d 0         alpha = 0. _d 0
        topLevTerm = freeSurfFac*cg3dNorm*  
      &      (horiVertRatio/gravity)/deltaTMom/deltaTMom  
316         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
317          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
318           alphaTile = 0. _d 0           alphaTile = 0. _d 0
# Line 335  C==    q = A.s Line 326  C==    q = A.s
326       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
327       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
328       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)
329       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
      &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)  
330               alphaTile = alphaTile               alphaTile = alphaTile
331       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
332              ENDDO              ENDDO
# Line 355  C==    q = A.s Line 341  C==    q = A.s
341       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)
342       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
343       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
344       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
      &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)  
345               alphaTile = alphaTile               alphaTile = alphaTile
346       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
347              ENDDO              ENDDO
# Line 376  C==    q = A.s Line 358  C==    q = A.s
358       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
359       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)
360       &     +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)       &     +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)
361       &     -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
      &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
362              alphaTile = alphaTile              alphaTile = alphaTile
363       &                +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &                +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
364             ENDDO             ENDDO
# Line 397  C==    q = A.s Line 374  C==    q = A.s
374       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
375       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
376       &      +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &      +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)
377       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
      &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
378               alphaTile = alphaTile               alphaTile = alphaTile
379       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
380              ENDDO              ENDDO
# Line 433  C      Now compute "interior" points. Line 406  C      Now compute "interior" points.
406       &            +alpha*cg3d_s(I,J,K,bi,bj)       &            +alpha*cg3d_s(I,J,K,bi,bj)
407              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)
408       &            -alpha*cg3d_q(I,J,K,bi,bj)       &            -alpha*cg3d_q(I,J,K,bi,bj)
409             errTile = errTile              errTile = errTile
410       &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
411             ENDDO             ENDDO
412            ENDDO            ENDDO
# Line 447  C      Now compute "interior" points. Line 420  C      Now compute "interior" points.
420         actualIts      = it3d         actualIts      = it3d
421         actualResidual = err         actualResidual = err
422         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
 C      _EXCH_XYZ_R8(cg3d_r, myThid )  
423         CALL EXCH_S3D_RL( cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, myThid )
424    
425     10 CONTINUE     10 CONTINUE
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22