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
0539a6785e Jean*0100       _RL    cg2dTolerance_sq
13785d3a37 Jean*0101       _RL    err_sq,  errTile(nSx,nSy)
                0102       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
                0103       _RL    eta_qrNM1, recip_eta_qrNM1
                0104       _RL    cgBeta,  alpha
                0105       _RL    alphaSum,alphaTile(nSx,nSy)
                0106       _RL    sumRHS,  sumRHStile(nSx,nSy)
9f5240b52a Jean*0107       _RL    rhsMax
0539a6785e Jean*0108       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0109       LOGICAL printResidual
d68b36a485 Mart*0110       _RL    cg2d_z(1:sNx,1:sNy,nSx,nSy)
                0111       _RL    cg2d_q(1:sNx,1:sNy,nSx,nSy)
                0112 C     Fields cg2d_r and ct2d_s are exchanged and since there are no AD-exchanges
aecc8b0f47 Mart*0113 C     for fields with smaller overlap, they are defined with full overlap size
                0114       _RL    cg2d_r(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0115       _RL    cg2d_s(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
05b6d6742b Patr*0116 CEOP
                0117
                0118 #ifdef ALLOW_AUTODIFF_TAMC
                0119       IF ( numIters .GT. numItersMax ) THEN
13785d3a37 Jean*0120         WRITE(msgBuf,'(A,I10)')
                0121      &    'CG2D_NSA: numIters > numItersMax =', numItersMax
                0122         CALL PRINT_ERROR( msgBuf, myThid )
                0123         STOP 'ABNORMAL END: S/R CG2D_NSA'
                0124       ENDIF
                0125       IF ( cg2dNormaliseRHS ) THEN
                0126         WRITE(msgBuf,'(A)') 'CG2D_NSA: cg2dNormaliseRHS is disabled'
                0127         CALL PRINT_ERROR( msgBuf, myThid )
                0128         WRITE(msgBuf,'(A)')
                0129      &    'set cg2dTargetResWunit (instead of cg2dTargetResidual)'
                0130         CALL PRINT_ERROR( msgBuf, myThid )
                0131         STOP 'ABNORMAL END: S/R CG2D_NSA'
05b6d6742b Patr*0132       ENDIF
                0133 #endif /* ALLOW_AUTODIFF_TAMC */
                0134
0539a6785e Jean*0135 C--   Initialise auxiliary constant, some output variable and inverter
                0136       cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
                0137       minResidualSq = -1. _d 0
                0138       eta_qrNM1     =  1. _d 0
                0139       recip_eta_qrNM1= 1. _d 0
05b6d6742b Patr*0140
                0141 C--   Normalise RHS
                0142       rhsMax = 0. _d 0
                0143       DO bj=myByLo(myThid),myByHi(myThid)
                0144        DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0145         DO j=1,sNy
                0146          DO i=1,sNx
                0147           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
                0148           rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
05b6d6742b Patr*0149          ENDDO
                0150         ENDDO
                0151        ENDDO
                0152       ENDDO
                0153
13785d3a37 Jean*0154 #ifndef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0155       IF (cg2dNormaliseRHS) THEN
13785d3a37 Jean*0156 C-    Normalise RHS :
                0157        _GLOBAL_MAX_RL( rhsMax, myThid )
05b6d6742b Patr*0158        rhsNorm = 1. _d 0
13785d3a37 Jean*0159        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
                0160        DO bj=myByLo(myThid),myByHi(myThid)
                0161         DO bi=myBxLo(myThid),myBxHi(myThid)
                0162          DO j=1,sNy
                0163           DO i=1,sNx
                0164            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
                0165            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
                0166           ENDDO
05b6d6742b Patr*0167          ENDDO
                0168         ENDDO
                0169        ENDDO
13785d3a37 Jean*0170 C-    end Normalise RHS
05b6d6742b Patr*0171       ENDIF
13785d3a37 Jean*0172 #endif /* ndef ALLOW_AUTODIFF_TAMC */
05b6d6742b Patr*0173
                0174 C--   Update overlaps
08061169a5 Jean*0175       CALL EXCH_XY_RL( cg2d_x, myThid )
05b6d6742b Patr*0176
                0177 C--   Initial residual calculation
                0178 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0179 CADJ STORE cg2d_b = comlev1_cg2d, key = ikey_dynamics, byte = isbyte
                0180 CADJ STORE cg2d_x = comlev1_cg2d, key = ikey_dynamics, byte = isbyte
05b6d6742b Patr*0181 #endif /* ALLOW_AUTODIFF_TAMC */
                0182       DO bj=myByLo(myThid),myByHi(myThid)
                0183        DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0184 C-    Initialise local working arrays:
d68b36a485 Mart*0185         DO j=1,sNy
                0186          DO i=1,sNx
ff51fd05c6 Jean*0187           cg2d_q(i,j,bi,bj) = 0. _d 0
                0188           cg2d_z(i,j,bi,bj) = 0. _d 0
                0189          ENDDO
                0190         ENDDO
                0191         DO j=1-OLy,sNy+OLy
                0192          DO i=1-OLx,sNx+OLx
                0193           cg2d_r(i,j,bi,bj) = 0. _d 0
                0194           cg2d_s(i,j,bi,bj) = 0. _d 0
08061169a5 Jean*0195          ENDDO
                0196         ENDDO
ff51fd05c6 Jean*0197         errTile(bi,bj)    = 0. _d 0
                0198         sumRHStile(bi,bj) = 0. _d 0
0539a6785e Jean*0199         DO j=1,sNy
                0200          DO i=1,sNx
                0201           cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0202      &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0203      &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0204      &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0205      &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0206      &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
05b6d6742b Patr*0207      &    )
13785d3a37 Jean*0208           errTile(bi,bj)    = errTile(bi,bj)
                0209      &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0210           sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
05b6d6742b Patr*0211          ENDDO
                0212         ENDDO
                0213        ENDDO
                0214       ENDDO
                0215
08061169a5 Jean*0216 c     CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
                0217       CALL EXCH_XY_RL ( cg2d_r, myThid )
13785d3a37 Jean*0218       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
                0219       CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
0539a6785e Jean*0220       actualIts = 0
                0221       IF ( err_sq .NE. 0. ) THEN
                0222         firstResidual = SQRT(err_sq)
                0223       ELSE
                0224         firstResidual = 0.
                0225       ENDIF
08061169a5 Jean*0226
0539a6785e Jean*0227       printResidual = .FALSE.
08061169a5 Jean*0228       IF ( debugLevel .GE. debLevZero ) THEN
                0229         _BEGIN_MASTER( myThid )
0539a6785e Jean*0230         printResidual = printResidualFreq.GE.1
08061169a5 Jean*0231         WRITE(standardmessageunit,'(A,1P2E22.14)')
aecc8b0f47 Mart*0232      &  ' cg2d_nsa: Sum(rhs),rhsMax = ', sumRHS,rhsMax
08061169a5 Jean*0233         _END_MASTER( myThid )
                0234       ENDIF
05b6d6742b Patr*0235
                0236 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                0237 Cml   begin main solver loop
                0238 #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
                0239 CADJ LOOP = iteration, cg2d_x = comlev_cg2d_iter
                0240 #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */
                0241       DO it2d=1, numIters
                0242 #ifdef ALLOW_LOOP_DIRECTIVE
                0243 CML      it2d = 0
08061169a5 Jean*0244 CML      DO WHILE ( err_sq .GT. cg2dTolerance_sq .and. it2d .LT. numIters )
05b6d6742b Patr*0245 CML      it2d = it2d+1
                0246 #endif /* ALLOW_LOOP_DIRECTIVE */
                0247
                0248 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0249        ikey = it2d + (ikey_dynamics-1)*numItersMax
                0250 CADJ STORE err_sq = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0251 #endif /* ALLOW_AUTODIFF_TAMC */
aecc8b0f47 Mart*0252        IF ( it2d .LE. cg2dMinItersNSA .OR.
                0253      &      err_sq .GE. cg2dTolerance_sq ) THEN
05b6d6742b Patr*0254
                0255 C--    Solve preconditioning equation and update
                0256 C--    conjugate direction vector "s".
                0257 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0258 CADJ STORE cg2d_r = comlev1_cg2d_iter, key = ikey, byte = isbyte
                0259 CADJ STORE cg2d_s = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0260 #endif /* ALLOW_AUTODIFF_TAMC */
                0261        DO bj=myByLo(myThid),myByHi(myThid)
                0262         DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0263          eta_qrNtile(bi,bj) = 0. _d 0
0539a6785e Jean*0264          DO j=1,sNy
                0265           DO i=1,sNx
                0266            cg2d_z(i,j,bi,bj) =
                0267      &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0268      &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0269      &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0270      &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0271      &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
13785d3a37 Jean*0272            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
0539a6785e Jean*0273      &     +cg2d_z(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
05b6d6742b Patr*0274           ENDDO
                0275          ENDDO
                0276         ENDDO
                0277        ENDDO
                0278
13785d3a37 Jean*0279        CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
05b6d6742b Patr*0280 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0281 CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = ikey, byte = isbyte
                0282 CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key=ikey, byte=isbyte
05b6d6742b Patr*0283 #endif /* ALLOW_AUTODIFF_TAMC */
                0284 CML       cgBeta   = eta_qrN/eta_qrNM1
                0285        cgBeta   = eta_qrN*recip_eta_qrNM1
13785d3a37 Jean*0286 Cml    store normalisation factor for the next interation (in case there is one).
05b6d6742b Patr*0287 CML    store the inverse of the normalization factor for higher precision
                0288 CML       eta_qrNM1 = eta_qrN
aecc8b0f47 Mart*0289        recip_eta_qrNM1 = 0. _d 0
                0290        IF ( eta_qrN .NE. 0. _d 0 ) recip_eta_qrNM1 = 1. _d 0/eta_qrN
05b6d6742b Patr*0291
                0292        DO bj=myByLo(myThid),myByHi(myThid)
                0293         DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0294          DO j=1,sNy
                0295           DO i=1,sNx
                0296            cg2d_s(i,j,bi,bj) = cg2d_z(i,j,bi,bj)
                0297      &                       + cgBeta*cg2d_s(i,j,bi,bj)
05b6d6742b Patr*0298           ENDDO
                0299          ENDDO
                0300         ENDDO
                0301        ENDDO
                0302
                0303 C--    Do exchanges that require messages i.e. between
                0304 C--    processes.
08061169a5 Jean*0305 c      CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
                0306        CALL EXCH_XY_RL ( cg2d_s, myThid )
05b6d6742b Patr*0307
                0308 C==    Evaluate laplace operator on conjugate gradient vector
                0309 C==    q = A.s
                0310 #ifdef ALLOW_AUTODIFF_TAMC
                0311 #ifndef ALLOW_LOOP_DIRECTIVE
edb6656069 Mart*0312 CADJ STORE cg2d_s = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0313 #endif /* not ALLOW_LOOP_DIRECTIVE */
                0314 #endif /* ALLOW_AUTODIFF_TAMC */
                0315        DO bj=myByLo(myThid),myByHi(myThid)
                0316         DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0317          alphaTile(bi,bj) = 0. _d 0
0539a6785e Jean*0318          DO j=1,sNy
                0319           DO i=1,sNx
                0320            cg2d_q(i,j,bi,bj) =
                0321      &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
                0322      &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
                0323      &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
                0324      &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
                0325      &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
13785d3a37 Jean*0326            alphaTile(bi,bj) = alphaTile(bi,bj)
                0327      &                      + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
05b6d6742b Patr*0328           ENDDO
                0329          ENDDO
                0330         ENDDO
                0331        ENDDO
13785d3a37 Jean*0332        CALL GLOBAL_SUM_TILE_RL( alphaTile, alphaSum, myThid )
aecc8b0f47 Mart*0333        alpha = 0. _d 0
                0334        IF ( alphaSum .NE. 0 _d 0 ) alpha = eta_qrN/alphaSum
08061169a5 Jean*0335
0539a6785e Jean*0336 C==    Update simultaneously solution and residual vectors (and Iter number)
05b6d6742b Patr*0337 C      Now compute "interior" points.
                0338 #ifdef ALLOW_AUTODIFF_TAMC
                0339 #ifndef ALLOW_LOOP_DIRECTIVE
edb6656069 Mart*0340 CADJ STORE cg2d_r = comlev1_cg2d_iter, key = ikey, byte = isbyte
05b6d6742b Patr*0341 #endif /* ALLOW_LOOP_DIRECTIVE */
                0342 #endif /* ALLOW_AUTODIFF_TAMC */
                0343        DO bj=myByLo(myThid),myByHi(myThid)
                0344         DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0345          errTile(bi,bj) = 0. _d 0
0539a6785e Jean*0346          DO j=1,sNy
                0347           DO i=1,sNx
                0348            cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
                0349            cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
13785d3a37 Jean*0350            errTile(bi,bj) = errTile(bi,bj)
                0351      &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
05b6d6742b Patr*0352           ENDDO
                0353          ENDDO
                0354         ENDDO
                0355        ENDDO
0539a6785e Jean*0356        actualIts = it2d
05b6d6742b Patr*0357
13785d3a37 Jean*0358        CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq,    myThid )
0539a6785e Jean*0359        IF ( printResidual ) THEN
                0360         IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
                0361          WRITE(msgBuf,'(A,I6,A,1PE21.14)')
aecc8b0f47 Mart*0362      &    ' cg2d_nsa: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
0539a6785e Jean*0363          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0364      &                       SQUEEZE_RIGHT, myThid )
                0365         ENDIF
                0366        ENDIF
05b6d6742b Patr*0367
08061169a5 Jean*0368 c      CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
                0369        CALL EXCH_XY_RL ( cg2d_r, myThid )
05b6d6742b Patr*0370
0539a6785e Jean*0371 Cml   end of if "err >= cg2dTolerance" block ; end main solver loop
05b6d6742b Patr*0372       ENDIF
                0373       ENDDO
                0374
13785d3a37 Jean*0375 #ifndef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0376       IF (cg2dNormaliseRHS) THEN
                0377 C--   Un-normalise the answer
                0378         DO bj=myByLo(myThid),myByHi(myThid)
                0379          DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0380           DO j=1,sNy
                0381            DO i=1,sNx
                0382             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
05b6d6742b Patr*0383            ENDDO
                0384           ENDDO
                0385          ENDDO
                0386         ENDDO
                0387       ENDIF
13785d3a37 Jean*0388 #endif /* ndef ALLOW_AUTODIFF_TAMC */
05b6d6742b Patr*0389
                0390 C--   Return parameters to caller
0539a6785e Jean*0391       IF ( err_sq .NE. 0. ) THEN
                0392         lastResidual = SQRT(err_sq)
                0393       ELSE
                0394         lastResidual = 0.
                0395       ENDIF
                0396       numIters = actualIts
05b6d6742b Patr*0397
                0398 #endif /* ALLOW_CG2D_NSA */
                0399       RETURN
                0400       END
                0401
0539a6785e Jean*0402 #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
13785d3a37 Jean*0403
05b6d6742b Patr*0404 C     These routines are routinely part of the TAMC/TAF library that is
                0405 C     not included in the MITcgm, therefore they are mimicked here.
13785d3a37 Jean*0406
aecc8b0f47 Mart*0407       SUBROUTINE ADSTORE(chardum,int1,idow,int2,int3,icount)
05b6d6742b Patr*0408
aecc8b0f47 Mart*0409       IMPLICIT NONE
05b6d6742b Patr*0410
                0411 #include "SIZE.h"
                0412 #include "tamc.h"
                0413
aecc8b0f47 Mart*0414       CHARACTER*(*) chardum
                0415       INTEGER int1, int2, int3, idow, icount
05b6d6742b Patr*0416
08061169a5 Jean*0417 C     the length of this vector must be greater or equal
05b6d6742b Patr*0418 C     twice the number of timesteps
aecc8b0f47 Mart*0419       INTEGER nidow
05b6d6742b Patr*0420 #ifdef ALLOW_TAMC_CHECKPOINTING
aecc8b0f47 Mart*0421       PARAMETER ( nidow = 2*nchklev_1*nchklev_2*nchklev_3 )
05b6d6742b Patr*0422 #else
aecc8b0f47 Mart*0423       PARAMETER ( nidow = 1000000 )
05b6d6742b Patr*0424 #endif /* ALLOW_TAMC_CHECKPOINTING */
aecc8b0f47 Mart*0425       INTEGER istoreidow(nidow)
                0426       COMMON /istorecommon/ istoreidow
05b6d6742b Patr*0427
aecc8b0f47 Mart*0428       PRINT *, 'ADSTORE: ', chardum, int1, idow, int2, int3, icount
05b6d6742b Patr*0429
aecc8b0f47 Mart*0430       IF ( icount .GT. nidow ) THEN
                0431          PRINT *, 'adstore: error: icount > nidow = ', nidow
                0432          STOP 'ABNORMAL STOP in ADSTORE'
                0433       ENDIF
05b6d6742b Patr*0434
                0435       istoreidow(icount) = idow
                0436
aecc8b0f47 Mart*0437       RETURN
                0438       END
05b6d6742b Patr*0439
aecc8b0f47 Mart*0440       SUBROUTINE ADRESTO(chardum,int1,idow,int2,int3,icount)
05b6d6742b Patr*0441
aecc8b0f47 Mart*0442       IMPLICIT NONE
05b6d6742b Patr*0443
                0444 #include "SIZE.h"
                0445 #include "tamc.h"
                0446
aecc8b0f47 Mart*0447       CHARACTER*(*) chardum
                0448       INTEGER int1, int2, int3, idow, icount
05b6d6742b Patr*0449
08061169a5 Jean*0450 C     the length of this vector must be greater or equal
05b6d6742b Patr*0451 C     twice the number of timesteps
aecc8b0f47 Mart*0452       INTEGER nidow
05b6d6742b Patr*0453 #ifdef ALLOW_TAMC_CHECKPOINTING
aecc8b0f47 Mart*0454       PARAMETER ( nidow = 2*nchklev_1*nchklev_2*nchklev_3 )
05b6d6742b Patr*0455 #else
aecc8b0f47 Mart*0456       PARAMETER ( nidow = 1000000 )
05b6d6742b Patr*0457 #endif /* ALLOW_TAMC_CHECKPOINTING */
aecc8b0f47 Mart*0458       INTEGER istoreidow(nidow)
                0459       COMMON /istorecommon/ istoreidow
05b6d6742b Patr*0460
aecc8b0f47 Mart*0461       PRINT *, 'ADRESTO: ', chardum, int1, idow, int2, int3, icount
05b6d6742b Patr*0462
aecc8b0f47 Mart*0463       IF ( icount .GT. nidow ) THEN
                0464          PRINT *, 'ADRESTO: error: icount > nidow = ', nidow
                0465          STOP 'ABNORMAL STOP in ADRESTO'
                0466       ENDIF
05b6d6742b Patr*0467
                0468       idow = istoreidow(icount)
                0469
aecc8b0f47 Mart*0470       RETURN
                0471       END
0539a6785e Jean*0472 #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */