model/src/cg2d_nsa.F

1eadaea85b Jean*0001 #include "PACKAGES_CONFIG.h"
05b6d6742b Patr*0002 #include "CPP_OPTIONS.h"
517dbdc414 Jean*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
aecc8b0f47 Mart*0006 #ifdef ALLOW_CTRL
                0007 # include "CTRL_OPTIONS.h"
                0008 #endif
05b6d6742b Patr*0009
                0010 CML THIS DOES NOT WORK +++++
                0011 #undef ALLOW_LOOP_DIRECTIVE
                0012 CBOP
                0013 C     !ROUTINE: CG2D_NSA
                0014 C     !INTERFACE:
08061169a5 Jean*0015       SUBROUTINE CG2D_NSA(
0539a6785e Jean*0016      U                cg2d_b, cg2d_x,
                0017      O                firstResidual, minResidualSq, lastResidual,
                0018      U                numIters, nIterMin,
05b6d6742b Patr*0019      I                myThid )
                0020 C     !DESCRIPTION: \bv
aecc8b0f47 Mart*0021 C     *================================================================*
08061169a5 Jean*0022 C     | SUBROUTINE CG2D_NSA
aecc8b0f47 Mart*0023 C     | o Two-dimensional grid problem conjugate-gradient inverter
                0024 C     |   (with preconditioner).
05b6d6742b Patr*0025 C     | o This version is used only in the case when the matrix
aecc8b0f47 Mart*0026 C     |   operator is Not "Self-Adjoint" (NSA). Any remaining residuals
                0027 C     |   will be immediately reported to the National Security Agency
                0028 C     *================================================================*
08061169a5 Jean*0029 C     | Con. grad is an iterative procedure for solving Ax = b.
                0030 C     | It requires the A be symmetric.
aecc8b0f47 Mart*0031 C     | This implementation assumes A is a five-diagonal matrix
                0032 C     | of the form that arises in the discrete representation of
                0033 C     | the del^2 operator in a two-dimensional space.
                0034 C     *================================================================*
05b6d6742b Patr*0035 C     \ev
                0036
                0037 C     !USES:
                0038       IMPLICIT NONE
                0039 C     === Global data ===
                0040 #include "SIZE.h"
                0041 #include "EEPARAMS.h"
                0042 #include "PARAMS.h"
                0043 #include "CG2D.h"
7c50f07931 Mart*0044 #ifdef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0045 # include "tamc.h"
                0046 #endif
                0047
                0048 C     !INPUT/OUTPUT PARAMETERS:
0539a6785e Jean*0049 C     cg2d_b    :: The source term or "right hand side" (Output: normalised RHS)
                0050 C     cg2d_x    :: The solution (Input: first guess)
08061169a5 Jean*0051 C     firstResidual :: the initial residual before any iterations
0539a6785e Jean*0052 C     minResidualSq :: the lowest residual reached (squared)
08061169a5 Jean*0053 C     lastResidual  :: the actual residual reached
0539a6785e Jean*0054 C     numIters  :: Inp: the maximum number of iterations allowed
                0055 C                  Out: the actual number of iterations used
                0056 C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
                0057 C                  Out: iteration number corresponding to lowest residual
9f5240b52a Jean*0058 C     myThid    :: my Thread Id number
05b6d6742b Patr*0059       _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0060       _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0061       _RL  firstResidual
0539a6785e Jean*0062       _RL  minResidualSq
05b6d6742b Patr*0063       _RL  lastResidual
                0064       INTEGER numIters
0539a6785e Jean*0065       INTEGER nIterMin
05b6d6742b Patr*0066       INTEGER myThid
                0067
                0068 #ifdef ALLOW_CG2D_NSA
                0069 C     !LOCAL VARIABLES:
0539a6785e Jean*0070 C     bi, bj     :: tile index in X and Y.
                0071 C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
                0072 C     actualIts  :: actual CG iteration number
                0073 C     err_sq     :: Measure of the square of the residual of Ax - b.
                0074 C     eta_qrN    :: Used in computing search directions; suffix N and NM1
                0075 C     eta_qrNM1     denote current and previous iterations respectively.
08061169a5 Jean*0076 C     recip_eta_qrNM1 :: reciprocal of eta_qrNM1
0539a6785e Jean*0077 C     cgBeta     :: coeff used to update conjugate direction vector "s".
                0078 C     alpha      :: coeff used to update solution & residual
13785d3a37 Jean*0079 C     alphaSum   :: to avoid the statement: alpha = 1./alpha (for TAMC/TAF)
0539a6785e Jean*0080 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
                0081 C                   debugging/trouble shooting diagnostic. For neumann problems
                0082 C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
                0083 C     msgBuf     :: Informational/error message buffer
aecc8b0f47 Mart*0084 C--   local working array (used to be in CG2D.h common block:
                0085 C     cg2d_q     :: Intermediate matrix-vector product term
                0086 C     cg2d_r     ::   *same*
                0087 C     cg2d_s     ::   *same*
                0088 C     cg2d_z     :: Intermediate matrix-vector product term
                0089 C                :: reduces the number of recomputation in adjoint mode
                0090 C                :: this field is superfluous if your cg2d is self-adjoint.
08061169a5 Jean*0091       INTEGER bi, bj
0539a6785e Jean*0092       INTEGER i, j, it2d
                0093       INTEGER actualIts
7c50f07931 Mart*0094 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0095 C     ikey       :: tape key (depends on iteration it2d)
                0096       INTEGER ikey
9f5240b52a Jean*0097 #else
                0098       _RL    rhsNorm
7c50f07931 Mart*0099 #endif
13785d3a37 Jean*0100       _RL    err_sq,  errTile(nSx,nSy)
                0101       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
                0102       _RL    eta_qrNM1, recip_eta_qrNM1
                0103       _RL    cgBeta,  alpha
                0104       _RL    alphaSum,alphaTile(nSx,nSy)
                0105       _RL    sumRHS,  sumRHStile(nSx,nSy)
9f5240b52a Jean*0106       _RL    rhsMax
0539a6785e Jean*0107       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0108       LOGICAL printResidual
d68b36a485 Mart*0109       _RL    cg2d_z(1:sNx,1:sNy,nSx,nSy)
                0110       _RL    cg2d_q(1:sNx,1:sNy,nSx,nSy)
                0111 C     Fields cg2d_r and ct2d_s are exchanged and since there are no AD-exchanges
aecc8b0f47 Mart*0112 C     for fields with smaller overlap, they are defined with full overlap size
                0113       _RL    cg2d_r(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0114       _RL    cg2d_s(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
05b6d6742b Patr*0115 CEOP
                0116
                0117 #ifdef ALLOW_AUTODIFF_TAMC
                0118       IF ( numIters .GT. numItersMax ) THEN
13785d3a37 Jean*0119         WRITE(msgBuf,'(A,I10)')
                0120      &    'CG2D_NSA: numIters > numItersMax =', numItersMax
                0121         CALL PRINT_ERROR( msgBuf, myThid )
                0122         STOP 'ABNORMAL END: S/R CG2D_NSA'
                0123       ENDIF
                0124       IF ( cg2dNormaliseRHS ) THEN
                0125         WRITE(msgBuf,'(A)') 'CG2D_NSA: cg2dNormaliseRHS is disabled'
                0126         CALL PRINT_ERROR( msgBuf, myThid )
                0127         WRITE(msgBuf,'(A)')
                0128      &    'set cg2dTargetResWunit (instead of cg2dTargetResidual)'
                0129         CALL PRINT_ERROR( msgBuf, myThid )
                0130         STOP 'ABNORMAL END: S/R CG2D_NSA'
05b6d6742b Patr*0131       ENDIF
                0132 #endif /* ALLOW_AUTODIFF_TAMC */
                0133
0539a6785e Jean*0134 C--   Initialise auxiliary constant, some output variable and inverter
                0135       minResidualSq = -1. _d 0
                0136       eta_qrNM1     =  1. _d 0
                0137       recip_eta_qrNM1= 1. _d 0
05b6d6742b Patr*0138
                0139 C--   Normalise RHS
                0140       rhsMax = 0. _d 0
                0141       DO bj=myByLo(myThid),myByHi(myThid)
                0142        DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0143         DO j=1,sNy
                0144          DO i=1,sNx
                0145           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
                0146           rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
05b6d6742b Patr*0147          ENDDO
                0148         ENDDO
                0149        ENDDO
                0150       ENDDO
                0151
13785d3a37 Jean*0152 #ifndef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0153       IF (cg2dNormaliseRHS) THEN
13785d3a37 Jean*0154 C-    Normalise RHS :
                0155        _GLOBAL_MAX_RL( rhsMax, myThid )
05b6d6742b Patr*0156        rhsNorm = 1. _d 0
13785d3a37 Jean*0157        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
                0158        DO bj=myByLo(myThid),myByHi(myThid)
                0159         DO bi=myBxLo(myThid),myBxHi(myThid)
                0160          DO j=1,sNy
                0161           DO i=1,sNx
                0162            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
                0163            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
                0164           ENDDO
05b6d6742b Patr*0165          ENDDO
                0166         ENDDO
                0167        ENDDO
13785d3a37 Jean*0168 C-    end Normalise RHS
05b6d6742b Patr*0169       ENDIF
13785d3a37 Jean*0170 #endif /* ndef ALLOW_AUTODIFF_TAMC */
05b6d6742b Patr*0171
                0172 C--   Update overlaps
08061169a5 Jean*0173       CALL EXCH_XY_RL( cg2d_x, myThid )
05b6d6742b Patr*0174
                0175 C--   Initial residual calculation
                0176 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0177 CADJ STORE cg2d_b = comlev1_cg2d, key = ikey_dynamics, byte = isbyte
                0178 CADJ STORE cg2d_x = comlev1_cg2d, key = ikey_dynamics, byte = isbyte
05b6d6742b Patr*0179 #endif /* ALLOW_AUTODIFF_TAMC */
                0180       DO bj=myByLo(myThid),myByHi(myThid)
                0181        DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0182 C-    Initialise local working arrays:
d68b36a485 Mart*0183         DO j=1,sNy
                0184          DO i=1,sNx
ff51fd05c6 Jean*0185           cg2d_q(i,j,bi,bj) = 0. _d 0
                0186           cg2d_z(i,j,bi,bj) = 0. _d 0
                0187          ENDDO
                0188         ENDDO
                0189         DO j=1-OLy,sNy+OLy
                0190          DO i=1-OLx,sNx+OLx
                0191           cg2d_r(i,j,bi,bj) = 0. _d 0
                0192           cg2d_s(i,j,bi,bj) = 0. _d 0
08061169a5 Jean*0193          ENDDO
                0194         ENDDO
ff51fd05c6 Jean*0195         errTile(bi,bj)    = 0. _d 0
                0196         sumRHStile(bi,bj) = 0. _d 0
0539a6785e Jean*0197         DO j=1,sNy
                0198          DO i=1,sNx
                0199           cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0200      &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0201      &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0202      &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0203      &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0204      &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
05b6d6742b Patr*0205      &    )
13785d3a37 Jean*0206           errTile(bi,bj)    = errTile(bi,bj)
                0207      &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0208           sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
05b6d6742b Patr*0209          ENDDO
                0210         ENDDO
                0211        ENDDO
                0212       ENDDO
                0213
08061169a5 Jean*0214 c     CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
                0215       CALL EXCH_XY_RL ( cg2d_r, myThid )
13785d3a37 Jean*0216       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
                0217       CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
0539a6785e Jean*0218       actualIts = 0
                0219       IF ( err_sq .NE. 0. ) THEN
                0220         firstResidual = SQRT(err_sq)
                0221       ELSE
                0222         firstResidual = 0.
                0223       ENDIF
08061169a5 Jean*0224
0539a6785e Jean*0225       printResidual = .FALSE.
08061169a5 Jean*0226       IF ( debugLevel .GE. debLevZero ) THEN
                0227         _BEGIN_MASTER( myThid )
0539a6785e Jean*0228         printResidual = printResidualFreq.GE.1
e6e223b277 Jean*0229         WRITE(standardMessageUnit,'(A,1P2E22.14)')
aecc8b0f47 Mart*0230      &  ' cg2d_nsa: Sum(rhs),rhsMax = ', sumRHS,rhsMax
08061169a5 Jean*0231         _END_MASTER( myThid )
                0232       ENDIF
05b6d6742b Patr*0233
                0234 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                0235 Cml   begin main solver loop
                0236 #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
                0237 CADJ LOOP = iteration, cg2d_x = comlev_cg2d_iter
                0238 #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */
                0239       DO it2d=1, numIters
                0240 #ifdef ALLOW_LOOP_DIRECTIVE
                0241 CML      it2d = 0
08061169a5 Jean*0242 CML      DO WHILE ( err_sq .GT. cg2dTolerance_sq .and. it2d .LT. numIters )
05b6d6742b Patr*0243 CML      it2d = it2d+1
                0244 #endif /* ALLOW_LOOP_DIRECTIVE */
                0245
                0246 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0247        ikey = it2d + (ikey_dynamics-1)*numItersMax
                0248 CADJ STORE err_sq = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0249 #endif /* ALLOW_AUTODIFF_TAMC */
aecc8b0f47 Mart*0250        IF ( it2d .LE. cg2dMinItersNSA .OR.
                0251      &      err_sq .GE. cg2dTolerance_sq ) THEN
05b6d6742b Patr*0252
                0253 C--    Solve preconditioning equation and update
                0254 C--    conjugate direction vector "s".
                0255 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0256 CADJ STORE cg2d_r = comlev1_cg2d_iter, key = ikey, byte = isbyte
                0257 CADJ STORE cg2d_s = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0258 #endif /* ALLOW_AUTODIFF_TAMC */
                0259        DO bj=myByLo(myThid),myByHi(myThid)
                0260         DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0261          eta_qrNtile(bi,bj) = 0. _d 0
0539a6785e Jean*0262          DO j=1,sNy
                0263           DO i=1,sNx
                0264            cg2d_z(i,j,bi,bj) =
                0265      &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0266      &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0267      &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0268      &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0269      &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
13785d3a37 Jean*0270            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
0539a6785e Jean*0271      &     +cg2d_z(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
05b6d6742b Patr*0272           ENDDO
                0273          ENDDO
                0274         ENDDO
                0275        ENDDO
                0276
13785d3a37 Jean*0277        CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
05b6d6742b Patr*0278 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0279 CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = ikey, byte = isbyte
                0280 CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key=ikey, byte=isbyte
05b6d6742b Patr*0281 #endif /* ALLOW_AUTODIFF_TAMC */
                0282 CML       cgBeta   = eta_qrN/eta_qrNM1
                0283        cgBeta   = eta_qrN*recip_eta_qrNM1
13785d3a37 Jean*0284 Cml    store normalisation factor for the next interation (in case there is one).
05b6d6742b Patr*0285 CML    store the inverse of the normalization factor for higher precision
                0286 CML       eta_qrNM1 = eta_qrN
aecc8b0f47 Mart*0287        recip_eta_qrNM1 = 0. _d 0
                0288        IF ( eta_qrN .NE. 0. _d 0 ) recip_eta_qrNM1 = 1. _d 0/eta_qrN
05b6d6742b Patr*0289
                0290        DO bj=myByLo(myThid),myByHi(myThid)
                0291         DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0292          DO j=1,sNy
                0293           DO i=1,sNx
                0294            cg2d_s(i,j,bi,bj) = cg2d_z(i,j,bi,bj)
                0295      &                       + cgBeta*cg2d_s(i,j,bi,bj)
05b6d6742b Patr*0296           ENDDO
                0297          ENDDO
                0298         ENDDO
                0299        ENDDO
                0300
                0301 C--    Do exchanges that require messages i.e. between
                0302 C--    processes.
08061169a5 Jean*0303 c      CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
                0304        CALL EXCH_XY_RL ( cg2d_s, myThid )
05b6d6742b Patr*0305
                0306 C==    Evaluate laplace operator on conjugate gradient vector
                0307 C==    q = A.s
                0308 #ifdef ALLOW_AUTODIFF_TAMC
                0309 #ifndef ALLOW_LOOP_DIRECTIVE
edb6656069 Mart*0310 CADJ STORE cg2d_s = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0311 #endif /* not ALLOW_LOOP_DIRECTIVE */
                0312 #endif /* ALLOW_AUTODIFF_TAMC */
                0313        DO bj=myByLo(myThid),myByHi(myThid)
                0314         DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0315          alphaTile(bi,bj) = 0. _d 0
0539a6785e Jean*0316          DO j=1,sNy
                0317           DO i=1,sNx
                0318            cg2d_q(i,j,bi,bj) =
                0319      &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
                0320      &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
                0321      &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
                0322      &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
                0323      &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
13785d3a37 Jean*0324            alphaTile(bi,bj) = alphaTile(bi,bj)
                0325      &                      + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
05b6d6742b Patr*0326           ENDDO
                0327          ENDDO
                0328         ENDDO
                0329        ENDDO
13785d3a37 Jean*0330        CALL GLOBAL_SUM_TILE_RL( alphaTile, alphaSum, myThid )
aecc8b0f47 Mart*0331        alpha = 0. _d 0
                0332        IF ( alphaSum .NE. 0 _d 0 ) alpha = eta_qrN/alphaSum
08061169a5 Jean*0333
0539a6785e Jean*0334 C==    Update simultaneously solution and residual vectors (and Iter number)
05b6d6742b Patr*0335 C      Now compute "interior" points.
                0336 #ifdef ALLOW_AUTODIFF_TAMC
                0337 #ifndef ALLOW_LOOP_DIRECTIVE
edb6656069 Mart*0338 CADJ STORE cg2d_r = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0339 #endif /* ALLOW_LOOP_DIRECTIVE */
                0340 #endif /* ALLOW_AUTODIFF_TAMC */
                0341        DO bj=myByLo(myThid),myByHi(myThid)
                0342         DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0343          errTile(bi,bj) = 0. _d 0
0539a6785e Jean*0344          DO j=1,sNy
                0345           DO i=1,sNx
                0346            cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
                0347            cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
13785d3a37 Jean*0348            errTile(bi,bj) = errTile(bi,bj)
                0349      &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
05b6d6742b Patr*0350           ENDDO
                0351          ENDDO
                0352         ENDDO
                0353        ENDDO
0539a6785e Jean*0354        actualIts = it2d
05b6d6742b Patr*0355
13785d3a37 Jean*0356        CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq,    myThid )
0539a6785e Jean*0357        IF ( printResidual ) THEN
                0358         IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
                0359          WRITE(msgBuf,'(A,I6,A,1PE21.14)')
aecc8b0f47 Mart*0360      &    ' cg2d_nsa: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
0539a6785e Jean*0361          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0362      &                       SQUEEZE_RIGHT, myThid )
                0363         ENDIF
                0364        ENDIF
05b6d6742b Patr*0365
08061169a5 Jean*0366 c      CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
                0367        CALL EXCH_XY_RL ( cg2d_r, myThid )
05b6d6742b Patr*0368
0539a6785e Jean*0369 Cml   end of if "err >= cg2dTolerance" block ; end main solver loop
05b6d6742b Patr*0370       ENDIF
                0371       ENDDO
                0372
13785d3a37 Jean*0373 #ifndef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0374       IF (cg2dNormaliseRHS) THEN
                0375 C--   Un-normalise the answer
                0376         DO bj=myByLo(myThid),myByHi(myThid)
                0377          DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0378           DO j=1,sNy
                0379            DO i=1,sNx
                0380             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
05b6d6742b Patr*0381            ENDDO
                0382           ENDDO
                0383          ENDDO
                0384         ENDDO
                0385       ENDIF
13785d3a37 Jean*0386 #endif /* ndef ALLOW_AUTODIFF_TAMC */
05b6d6742b Patr*0387
                0388 C--   Return parameters to caller
0539a6785e Jean*0389       IF ( err_sq .NE. 0. ) THEN
                0390         lastResidual = SQRT(err_sq)
                0391       ELSE
                0392         lastResidual = 0.
                0393       ENDIF
                0394       numIters = actualIts
05b6d6742b Patr*0395
                0396 #endif /* ALLOW_CG2D_NSA */
                0397       RETURN
                0398       END
                0399
0539a6785e Jean*0400 #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
13785d3a37 Jean*0401
05b6d6742b Patr*0402 C     These routines are routinely part of the TAMC/TAF library that is
                0403 C     not included in the MITcgm, therefore they are mimicked here.
13785d3a37 Jean*0404
aecc8b0f47 Mart*0405       SUBROUTINE ADSTORE(chardum,int1,idow,int2,int3,icount)
05b6d6742b Patr*0406
aecc8b0f47 Mart*0407       IMPLICIT NONE
05b6d6742b Patr*0408
                0409 #include "SIZE.h"
                0410 #include "tamc.h"
                0411
aecc8b0f47 Mart*0412       CHARACTER*(*) chardum
                0413       INTEGER int1, int2, int3, idow, icount
05b6d6742b Patr*0414
08061169a5 Jean*0415 C     the length of this vector must be greater or equal
05b6d6742b Patr*0416 C     twice the number of timesteps
aecc8b0f47 Mart*0417       INTEGER nidow
05b6d6742b Patr*0418 #ifdef ALLOW_TAMC_CHECKPOINTING
aecc8b0f47 Mart*0419       PARAMETER ( nidow = 2*nchklev_1*nchklev_2*nchklev_3 )
05b6d6742b Patr*0420 #else
aecc8b0f47 Mart*0421       PARAMETER ( nidow = 1000000 )
05b6d6742b Patr*0422 #endif /* ALLOW_TAMC_CHECKPOINTING */
aecc8b0f47 Mart*0423       INTEGER istoreidow(nidow)
                0424       COMMON /istorecommon/ istoreidow
05b6d6742b Patr*0425
aecc8b0f47 Mart*0426       PRINT *, 'ADSTORE: ', chardum, int1, idow, int2, int3, icount
05b6d6742b Patr*0427
aecc8b0f47 Mart*0428       IF ( icount .GT. nidow ) THEN
                0429          PRINT *, 'adstore: error: icount > nidow = ', nidow
                0430          STOP 'ABNORMAL STOP in ADSTORE'
                0431       ENDIF
05b6d6742b Patr*0432
                0433       istoreidow(icount) = idow
                0434
aecc8b0f47 Mart*0435       RETURN
                0436       END
05b6d6742b Patr*0437
aecc8b0f47 Mart*0438       SUBROUTINE ADRESTO(chardum,int1,idow,int2,int3,icount)
05b6d6742b Patr*0439
aecc8b0f47 Mart*0440       IMPLICIT NONE
05b6d6742b Patr*0441
                0442 #include "SIZE.h"
                0443 #include "tamc.h"
                0444
aecc8b0f47 Mart*0445       CHARACTER*(*) chardum
                0446       INTEGER int1, int2, int3, idow, icount
05b6d6742b Patr*0447
08061169a5 Jean*0448 C     the length of this vector must be greater or equal
05b6d6742b Patr*0449 C     twice the number of timesteps
aecc8b0f47 Mart*0450       INTEGER nidow
05b6d6742b Patr*0451 #ifdef ALLOW_TAMC_CHECKPOINTING
aecc8b0f47 Mart*0452       PARAMETER ( nidow = 2*nchklev_1*nchklev_2*nchklev_3 )
05b6d6742b Patr*0453 #else
aecc8b0f47 Mart*0454       PARAMETER ( nidow = 1000000 )
05b6d6742b Patr*0455 #endif /* ALLOW_TAMC_CHECKPOINTING */
aecc8b0f47 Mart*0456       INTEGER istoreidow(nidow)
                0457       COMMON /istorecommon/ istoreidow
05b6d6742b Patr*0458
aecc8b0f47 Mart*0459       PRINT *, 'ADRESTO: ', chardum, int1, idow, int2, int3, icount
05b6d6742b Patr*0460
aecc8b0f47 Mart*0461       IF ( icount .GT. nidow ) THEN
                0462          PRINT *, 'ADRESTO: error: icount > nidow = ', nidow
                0463          STOP 'ABNORMAL STOP in ADRESTO'
                0464       ENDIF
05b6d6742b Patr*0465
                0466       idow = istoreidow(icount)
                0467
aecc8b0f47 Mart*0468       RETURN
                0469       END
0539a6785e Jean*0470 #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */