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

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

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

revision 1.3 by jmc, Sat Jul 11 21:45:44 2009 UTC revision 1.4 by jmc, Fri May 11 23:29:13 2012 UTC
# Line 13  CBOP Line 13  CBOP
13  C     !ROUTINE: CG2D_NSA  C     !ROUTINE: CG2D_NSA
14  C     !INTERFACE:  C     !INTERFACE:
15        SUBROUTINE CG2D_NSA(        SUBROUTINE CG2D_NSA(
16       I                cg2d_b,       U                cg2d_b, cg2d_x,
17       U                cg2d_x,       O                firstResidual, minResidualSq, lastResidual,
18       O                firstResidual,       U                numIters, nIterMin,
      O                lastResidual,  
      U                numIters,  
19       I                myThid )       I                myThid )
20  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
21  C     *==========================================================*  C     *==========================================================*
# Line 55  C     === Global data === Line 53  C     === Global data ===
53  #include "EEPARAMS.h"  #include "EEPARAMS.h"
54  #include "PARAMS.h"  #include "PARAMS.h"
55  #include "CG2D.h"  #include "CG2D.h"
 c#include "GRID.h"  
 c#include "SURFACE.h"  
56  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
57  # include "tamc.h"  # include "tamc.h"
58  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 64  c#include "SURFACE.h" Line 60  c#include "SURFACE.h"
60    
61  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
62  C     === Routine arguments ===  C     === Routine arguments ===
63  C     cg2d_b    :: The source term or "right hand side"  C     cg2d_b    :: The source term or "right hand side" (Output: normalised RHS)
64  C     cg2d_x    :: The solution  C     cg2d_x    :: The solution (Input: first guess)
65  C     firstResidual :: the initial residual before any iterations  C     firstResidual :: the initial residual before any iterations
66    C     minResidualSq :: the lowest residual reached (squared)
67  C     lastResidual  :: the actual residual reached  C     lastResidual  :: the actual residual reached
68  C     numIters  :: Entry: the maximum number of iterations allowed  C     numIters  :: Inp: the maximum number of iterations allowed
69  C                  Exit:  the actual number of iterations used  C                  Out: the actual number of iterations used
70    C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
71    C                  Out: iteration number corresponding to lowest residual
72  C     myThid    :: Thread on which I am working.  C     myThid    :: Thread on which I am working.
73        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
74        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
75        _RL  firstResidual        _RL  firstResidual
76          _RL  minResidualSq
77        _RL  lastResidual        _RL  lastResidual
78        INTEGER numIters        INTEGER numIters
79          INTEGER nIterMin
80        INTEGER myThid        INTEGER myThid
81    
82  #ifdef ALLOW_CG2D_NSA  #ifdef ALLOW_CG2D_NSA
83  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
84  C     === Local variables ====  C     === Local variables ====
85  C     actualIts      :: Number of iterations taken  C     bi, bj     :: tile index in X and Y.
86  C     actualResidual :: residual  C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
87  C     bi, bj     :: Block index in X and Y.  C     actualIts  :: actual CG iteration number
88  C     eta_qrN    :: Used in computing search directions  C     err_sq     :: Measure of the square of the residual of Ax - b.
89  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrN    :: Used in computing search directions; suffix N and NM1
90  C     cgBeta        previous iterations respectively.  C     eta_qrNM1     denote current and previous iterations respectively.
91  C     recip_eta_qrNM1 :: reciprocal of eta_qrNM1  C     recip_eta_qrNM1 :: reciprocal of eta_qrNM1
92  C     alpha  C     cgBeta     :: coeff used to update conjugate direction vector "s".
93    C     alpha      :: coeff used to update solution & residual
94  C     alpha_aux  :: to avoid the statement: alpha = 1./alpha (for TAMC/TAF)  C     alpha_aux  :: to avoid the statement: alpha = 1./alpha (for TAMC/TAF)
95  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
96  C                   useful debuggin/trouble shooting diagnostic.  C                   debugging/trouble shooting diagnostic. For neumann problems
97  C                   For neumann problems sumRHS needs to be ~0.  C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
98  C                   or they converge at a non-zero residual.  C     cg2d_min   :: used to store solution corresponding to lowest residual.
99  C     err        :: Measure of residual of Ax - b, usually the norm.  C     msgBuf     :: Informational/error message buffer
 C     err_sq     :: square of err (for TAMC/TAF)  
 C     I, J, it2d :: Loop counters ( it2d counts CG iterations )  
       INTEGER actualIts  
       _RL    actualResidual  
100        INTEGER bi, bj        INTEGER bi, bj
101        INTEGER I, J, it2d        INTEGER i, j, it2d
102          INTEGER actualIts
103        _RL    err        _RL    cg2dTolerance_sq
104        _RL    err_sq        _RL    err_sq
105        _RL    eta_qrN        _RL    eta_qrN
106        _RL    eta_qrNM1        _RL    eta_qrNM1
# Line 113  C     I, J, it2d :: Loop counters ( it2d Line 111  C     I, J, it2d :: Loop counters ( it2d
111        _RL    sumRHS        _RL    sumRHS
112        _RL    rhsMax, rhsMaxGlobal        _RL    rhsMax, rhsMaxGlobal
113        _RL    rhsNorm        _RL    rhsNorm
114        _RL    cg2dTolerance_sq        CHARACTER*(MAX_LEN_MBUF) msgBuf
115          LOGICAL printResidual
116  CEOP  CEOP
117    
118  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 124  CEOP Line 123  CEOP
123        ENDIF        ENDIF
124  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
125    
 CcnhDebugStarts  
 C     CHARACTER*(MAX_LEN_FNAM) suff  
 CcnhDebugEnds  
   
126  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
127            act1 = myThid - 1            act1 = myThid - 1
128            max1 = nTx*nTy            max1 = nTx*nTy
# Line 135  CcnhDebugEnds Line 130  CcnhDebugEnds
130            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
131  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
132    
133  C--   Initialise inverter  C--   Initialise auxiliary constant, some output variable and inverter
134        eta_qrNM1 = 1. _d 0        cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
135        recip_eta_qrNM1 = 1./eta_qrNM1        minResidualSq = -1. _d 0
136          eta_qrNM1     =  1. _d 0
137  CcnhDebugStarts        recip_eta_qrNM1= 1. _d 0
 C     _EXCH_XY_RL( cg2d_b, myThid )  
 C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )  
 C     suff = 'unnormalised'  
 C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)  
 C     STOP  
 CcnhDebugEnds  
138    
139  C--   Normalise RHS  C--   Normalise RHS
140  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 154  CADJ STORE cg2d_b = comlev1_cg2d, key = Line 143  CADJ STORE cg2d_b = comlev1_cg2d, key =
143        rhsMax = 0. _d 0        rhsMax = 0. _d 0
144        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
145         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
146          DO J=1,sNy          DO j=1,sNy
147           DO I=1,sNx           DO i=1,sNx
148            cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*cg2dNorm            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
149            rhsMax = MAX(ABS(cg2d_b(I,J,bi,bj)),rhsMax)            rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
150           ENDDO           ENDDO
151          ENDDO          ENDDO
152         ENDDO         ENDDO
# Line 165  CADJ STORE cg2d_b = comlev1_cg2d, key = Line 154  CADJ STORE cg2d_b = comlev1_cg2d, key =
154    
155        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
156  C     -  Normalise RHS :  C     -  Normalise RHS :
 #ifdef LETS_MAKE_JAM  
 C     _GLOBAL_MAX_RL( rhsMax, myThid )  
       rhsMaxGlobal=1.  
 #else  
157  #ifdef ALLOW_CONST_RHSMAX  #ifdef ALLOW_CONST_RHSMAX
158        rhsMaxGlobal=1.        rhsMaxGlobal=1.
159  #else  #else
160        rhsMaxGlobal=rhsMax        rhsMaxGlobal=rhsMax
161        _GLOBAL_MAX_RL( rhsMaxGlobal, myThid )        _GLOBAL_MAX_RL( rhsMaxGlobal, myThid )
162  #endif /* ALLOW_CONST_RHSMAX */  #endif /* ALLOW_CONST_RHSMAX */
 #endif  
163  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
164  CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte
165  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 191  CADJ STORE cg2d_x = comlev1_cg2d, key = Line 175  CADJ STORE cg2d_x = comlev1_cg2d, key =
175  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
176        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
177         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
178          DO J=1,sNy          DO j=1,sNy
179           DO I=1,sNx           DO i=1,sNx
180            cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*rhsNorm            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
181            cg2d_x(I,J,bi,bj) = cg2d_x(I,J,bi,bj)*rhsNorm            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
182           ENDDO           ENDDO
183          ENDDO          ENDDO
184         ENDDO         ENDDO
# Line 203  C- end Normalise RHS Line 187  C- end Normalise RHS
187        ENDIF        ENDIF
188    
189  C--   Update overlaps  C--   Update overlaps
 c     CALL EXCH_XY_RL( cg2d_b, myThid )  
190        CALL EXCH_XY_RL( cg2d_x, myThid )        CALL EXCH_XY_RL( cg2d_x, myThid )
 CcnhDebugStarts  
 C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )  
 C     suff = 'normalised'  
 C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)  
 CcnhDebugEnds  
191    
192  C--   Initial residual calculation  C--   Initial residual calculation
       err = 0. _d 0  
193        err_sq = 0. _d 0        err_sq = 0. _d 0
194        sumRHS = 0. _d 0        sumRHS = 0. _d 0
195  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 221  CADJ STORE cg2d_x = comlev1_cg2d, key = Line 198  CADJ STORE cg2d_x = comlev1_cg2d, key =
198  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
199        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
200         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
201          DO J=1-1,sNy+1          DO j=0,sNy+1
202           DO I=1-1,sNx+1           DO i=0,sNx+1
203            cg2d_s(I,J,bi,bj) = 0.            cg2d_s(i,j,bi,bj) = 0.
204           ENDDO           ENDDO
205          ENDDO          ENDDO
206          DO J=1,sNy          DO j=1,sNy
207           DO I=1,sNx           DO i=1,sNx
208            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -            cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
209       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)       &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
210       &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)       &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
211       &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)       &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
212       &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)       &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
213       &    +aC2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
214       &    )       &    )
 c    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 c    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 c    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 c    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 c    &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  
 c    &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  
 cML     &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm  
 c    &    )  
215            err_sq         = err_sq         +            err_sq         = err_sq         +
216       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &     cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
217            sumRHS            = sumRHS            +            sumRHS            = sumRHS            +
218       &     cg2d_b(I,J,bi,bj)       &     cg2d_b(i,j,bi,bj)
219           ENDDO           ENDDO
220          ENDDO          ENDDO
221         ENDDO         ENDDO
# Line 256  c     CALL EXCH_S3D_RL( cg2d_r, 1, myThi Line 225  c     CALL EXCH_S3D_RL( cg2d_r, 1, myThi
225        CALL EXCH_XY_RL ( cg2d_r, myThid )        CALL EXCH_XY_RL ( cg2d_r, myThid )
226        _GLOBAL_SUM_RL( sumRHS, myThid )        _GLOBAL_SUM_RL( sumRHS, myThid )
227        _GLOBAL_SUM_RL( err_sq, myThid )        _GLOBAL_SUM_RL( err_sq, myThid )
228         IF ( err_sq .NE. 0. ) THEN        actualIts = 0
229            err = SQRT(err_sq)        IF ( err_sq .NE. 0. ) THEN
230         ELSE          firstResidual = SQRT(err_sq)
231            err = 0.        ELSE
232         ENDIF          firstResidual = 0.
233         actualIts      = 0        ENDIF
        actualResidual = err  
234    
235          printResidual = .FALSE.
236        IF ( debugLevel .GE. debLevZero ) THEN        IF ( debugLevel .GE. debLevZero ) THEN
237          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
238            printResidual = printResidualFreq.GE.1
239          WRITE(standardmessageunit,'(A,1P2E22.14)')          WRITE(standardmessageunit,'(A,1P2E22.14)')
240       &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMaxGlobal       &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMaxGlobal
241          _END_MASTER( myThid )          _END_MASTER( myThid )
242        ENDIF        ENDIF
 C     _BARRIER  
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',  
 c    & actualIts, actualResidual  
 c     _END_MASTER( myThid )  
       firstResidual=actualResidual  
       cg2dTolerance_sq = cg2dTolerance**2  
243    
244  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
245  Cml   begin main solver loop  Cml   begin main solver loop
# Line 292  CML      it2d = it2d+1 Line 255  CML      it2d = it2d+1
255    
256  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
257         icg2dkey = (ikey-1)*numItersMax + it2d         icg2dkey = (ikey-1)*numItersMax + it2d
 CMLCADJ STORE err = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  
258  CADJ STORE err_sq = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE err_sq = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
259  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
260  CML      IF ( err .LT. cg2dTolerance ) THEN        IF ( err_sq .GE. cg2dTolerance_sq ) THEN
       IF ( err_sq .LT. cg2dTolerance_sq ) THEN  
 Cml      DO NOTHING  
       ELSE  
261    
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' residual = ',  
 C    &  actualResidual  
 CcnhDebugEnds  
262  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
263  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
264         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
# Line 313  CADJ STORE cg2d_s = comlev1_cg2d_iter, k Line 268  CADJ STORE cg2d_s = comlev1_cg2d_iter, k
268  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
269         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
270          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
271           DO J=1,sNy           DO j=1,sNy
272            DO I=1,sNx            DO i=1,sNx
273             cg2d_z(I,J,bi,bj) =             cg2d_z(i,j,bi,bj) =
274       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)       &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
275       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)       &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
276       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)       &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
277       &     +pS(I  ,J  ,bi,bj)*cg2d_r(I  ,J-1,bi,bj)       &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
278       &     +pS(I  ,J+1,bi,bj)*cg2d_r(I  ,J+1,bi,bj)       &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
 CcnhDebugStarts  
 C          cg2d_z(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)  
 CcnhDebugEnds  
279             eta_qrN = eta_qrN             eta_qrN = eta_qrN
280       &     +cg2d_z(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &     +cg2d_z(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
281            ENDDO            ENDDO
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284         ENDDO         ENDDO
285    
286         _GLOBAL_SUM_RL(eta_qrN, myThid)         _GLOBAL_SUM_RL(eta_qrN, myThid)
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' eta_qrN = ',eta_qrN  
 CcnhDebugEnds  
287  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
288  CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
289  CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
290  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
291  CML       cgBeta   = eta_qrN/eta_qrNM1  CML       cgBeta   = eta_qrN/eta_qrNM1
292         cgBeta   = eta_qrN*recip_eta_qrNM1         cgBeta   = eta_qrN*recip_eta_qrNM1
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' beta = ',cgBeta  
 CcnhDebugEnds  
293  Cml    store normalisation factor for the next interation  Cml    store normalisation factor for the next interation
294  Cml    (in case there is one).  Cml    (in case there is one).
295  CML    store the inverse of the normalization factor for higher precision  CML    store the inverse of the normalization factor for higher precision
# Line 352  CML       eta_qrNM1 = eta_qrN Line 298  CML       eta_qrNM1 = eta_qrN
298    
299         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
300          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
301           DO J=1,sNy           DO j=1,sNy
302            DO I=1,sNx            DO i=1,sNx
303             cg2d_s(I,J,bi,bj) = cg2d_z(I,J,bi,bj)             cg2d_s(i,j,bi,bj) = cg2d_z(i,j,bi,bj)
304       &                       + cgBeta*cg2d_s(I,J,bi,bj)       &                       + cgBeta*cg2d_s(i,j,bi,bj)
305            ENDDO            ENDDO
306           ENDDO           ENDDO
307          ENDDO          ENDDO
# Line 377  CADJ STORE cg2d_s = comlev1_cg2d_iter, k Line 323  CADJ STORE cg2d_s = comlev1_cg2d_iter, k
323  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
324         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
325          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
326           DO J=1,sNy           DO j=1,sNy
327            DO I=1,sNx            DO i=1,sNx
328             cg2d_q(I,J,bi,bj) =             cg2d_q(i,j,bi,bj) =
329       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)       &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
330       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)       &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
331       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)       &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
332       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)       &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
333       &    +aC2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
334  c    &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)             alpha_aux = alpha_aux+cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
 c    &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  
 c    &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  
 c    &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  
 c    &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  
 c    &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  
 cML     &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm  
            alpha_aux = alpha_aux+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)  
335            ENDDO            ENDDO
336           ENDDO           ENDDO
337          ENDDO          ENDDO
338         ENDDO         ENDDO
339         _GLOBAL_SUM_RL(alpha_aux,myThid)         _GLOBAL_SUM_RL(alpha_aux,myThid)
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' SUM(s*q)= ',alpha_aux  
 CcnhDebugEnds  
340         alpha = eta_qrN/alpha_aux         alpha = eta_qrN/alpha_aux
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' alpha= ',alpha  
 CcnhDebugEnds  
341    
342  C==    Update solution and residual vectors  C==    Update simultaneously solution and residual vectors (and Iter number)
343  C      Now compute "interior" points.  C      Now compute "interior" points.
        err    = 0. _d 0  
344         err_sq = 0. _d 0         err_sq = 0. _d 0
345  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
346  #ifndef ALLOW_LOOP_DIRECTIVE  #ifndef ALLOW_LOOP_DIRECTIVE
# Line 417  CADJ STORE cg2d_r = comlev1_cg2d_iter, k Line 349  CADJ STORE cg2d_r = comlev1_cg2d_iter, k
349  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
350         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
351          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
352           DO J=1,sNy           DO j=1,sNy
353            DO I=1,sNx            DO i=1,sNx
354             cg2d_x(I,J,bi,bj)=cg2d_x(I,J,bi,bj)+alpha*cg2d_s(I,J,bi,bj)             cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
355             cg2d_r(I,J,bi,bj)=cg2d_r(I,J,bi,bj)-alpha*cg2d_q(I,J,bi,bj)             cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
356             err_sq = err_sq+cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)             err_sq = err_sq+cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
357            ENDDO            ENDDO
358           ENDDO           ENDDO
359          ENDDO          ENDDO
360         ENDDO         ENDDO
361           actualIts = it2d
362    
363         _GLOBAL_SUM_RL( err_sq   , myThid )         _GLOBAL_SUM_RL( err_sq   , myThid )
364         if ( err_sq .ne. 0. ) then         IF ( printResidual ) THEN
365            err = SQRT(err_sq)          IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
366         else           WRITE(msgBuf,'(A,I6,A,1PE21.14)')
367            err = 0.       &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
368         end if           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
369         actualIts      = it2d       &                       SQUEEZE_RIGHT, myThid )
370         actualResidual = err          ENDIF
371           ENDIF
372    
373  c      CALL EXCH_S3D_RL( cg2d_r, 1, myThid )  c      CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
374         CALL EXCH_XY_RL ( cg2d_r, myThid )         CALL EXCH_XY_RL ( cg2d_r, myThid )
375    
376  Cml   end of IF ( err .LT. cg2dTolerance ) THEN; ELSE  Cml   end of if "err >= cg2dTolerance" block ; end main solver loop
377        ENDIF        ENDIF
 Cml     end main solver loop  
378        ENDDO        ENDDO
379    
380        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
# Line 452  CADJ STORE cg2d_x = comlev1_cg2d, key = Line 385  CADJ STORE cg2d_x = comlev1_cg2d, key =
385  C--   Un-normalise the answer  C--   Un-normalise the answer
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              cg2d_x(I  ,J  ,bi,bj) = cg2d_x(I  ,J  ,bi,bj)/rhsNorm              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
391             ENDDO             ENDDO
392            ENDDO            ENDDO
393           ENDDO           ENDDO
394          ENDDO          ENDDO
395        ENDIF        ENDIF
396    
 C     The following exchange was moved up to solve_for_pressure  
 C     for compatibility with TAMC.  
 C     _EXCH_XY_RL(cg2d_x, myThid )  
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',  
 c    & actualIts, actualResidual  
 c     _END_MASTER( )  
   
397  C--   Return parameters to caller  C--   Return parameters to caller
398        lastResidual=actualResidual        IF ( err_sq .NE. 0. ) THEN
399        numIters=actualIts          lastResidual = SQRT(err_sq)
400          ELSE
401            lastResidual = 0.
402          ENDIF
403          numIters = actualIts
404    
405  #endif /* ALLOW_CG2D_NSA */  #endif /* ALLOW_CG2D_NSA */
406        RETURN        RETURN
407        END        END
408    
409  # if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))  #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
410  C  C
411  C     These routines are routinely part of the TAMC/TAF library that is  C     These routines are routinely part of the TAMC/TAF library that is
412  C     not included in the MITcgm, therefore they are mimicked here.  C     not included in the MITcgm, therefore they are mimicked here.
# Line 548  C     twice the number of timesteps Line 477  C     twice the number of timesteps
477    
478        return        return
479        end        end
480  # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */  #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */
   
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22