Merge remote-tracking branch 'upstream/master' into genmake2_hybrid

doc/getting_started/getting_started.rst

@@ -2081,6 +2081,9 @@ elliptic solvers are the variables :varlink:`cg2dMaxIters` and
    | :varlink:`cg3dTargetResidual`          | PARM02    | 1.0E-07                                          | 3D conjugate gradient target residual (non-dim. due to RHS normalization );                             |
    |                                        |           |                                                  | requires #define :varlink:`ALLOW_NONHYDROSTATIC`                                                        |
    +----------------------------------------+-----------+--------------------------------------------------+---------------------------------------------------------------------------------------------------------+
+   | :varlink:`cg3dTargetResWunit`          | PARM02    | -1.0E+00                                         | 3D conjugate gradient target residual (:math:`\dot{r}` units);                                          |
+   |                                        |           |                                                  | <0: use RHS normalization, i.e., :varlink:`cg3dTargetResidual` instead                                  |
+   +----------------------------------------+-----------+--------------------------------------------------+---------------------------------------------------------------------------------------------------------+
    | :varlink:`useSRCGSolver`               | PARM02    | FALSE                                            | use conjugate gradient solver with single reduction (single call of mpi_allreduce)                      |
    +----------------------------------------+-----------+--------------------------------------------------+---------------------------------------------------------------------------------------------------------+
    | :varlink:`printResidualFreq`           | PARM02    | 1 unless :varlink:`debugLevel` >4                | frequency (in number of iterations) of printing conjugate gradient residual                             |

doc/tag-index

@@ -1,6 +1,16 @@
     Notes on tags used in MITgcmUV
     ==============================
 
+o model/src:
+  - Add new param "cg3dTargetResWunit" to skip the normalization of cg3d RHS
+    by specifying the CG3D target residual in term of vertical-velocity units ;
+  - Remove header files SOLVE_FOR_PRESSURE{3D}.h ; clean-up "post_cg3d.F" ;
+  - switch test exp. global_ocean.cs32x15.viscA4 & short_surf_wave to use new
+    param "cg3dTargetResWunit".
+o tools/f90mkdepend:
+  - add arguments to f90mkdepend to specify alternative file suffix (e.g.,
+    for MacOS) + minor improvement in error message.
+
 checkpoint69a (2024/09/08)
 o tools/genmake2:
   - list of src files to Differentiate: change from simply "*.list" to commonly

model/inc/CG2D.h

@@ -19,27 +19,27 @@ C     cg2dNormaliseRHS :: flag set to TRUE if normalise RHS in the Solver
       COMMON /CG2D_I_L/ cg2dNormaliseRHS
       LOGICAL cg2dNormaliseRHS
 
-C--   COMMON /CG2D_R/ DEL**2 Laplacian operators
+C--   COMMON /CG2D_I_R/ DEL**2 Laplacian operators
 C     aW2d    :: East-west operator.
 C     aS2d    :: North-south operator.
 C     aC2d    :: 2D operator main diagonal term.
 C     pW      :: East-west off-diagonal term of preconditioner.
 C     pS      :: North-south off-diagonal term of preconditioner.
 C     pC      :: Main diagonal term of preconditioner.
 C     cg2dNorm :: A matrix normalisation factor.
-C     cg2dTolerance :: cg2d solver Tolerance (solver unit = m2/s2 or no unit,
-C                                             depending on cg2dNormaliseRHS)
+C     cg2dTolerance_sq :: square of cg2d solver Tolerance (units depends
+C                 on cg2dNormaliseRHS, solver-unit ^2 = (m2/s2)^2 or no unit)
       COMMON /CG2D_I_RS/
      &      aW2d, aS2d, aC2d,
      &      pW, pS, pC
       COMMON /CG2D_I_RL/
-     &      cg2dNorm, cg2dTolerance
+     &      cg2dNorm, cg2dTolerance_sq
       _RS  aW2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RS  aS2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RS  aC2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RS  pW   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RS  pS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RS  pC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
-      _RL  cg2dNorm, cg2dTolerance
+      _RL  cg2dNorm, cg2dTolerance_sq
 
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

model/inc/CG3D.h

@@ -17,35 +17,42 @@ C     *==========================================================*
 C     \ev
 CEOP
 
+C--   COMMON /CG3D_L/ cg3dNormaliseRHS
+C     cg3dNormaliseRHS :: flag set to TRUE if normalise RHS in the Solver
+      COMMON /CG3D_L/ cg3dNormaliseRHS
+      LOGICAL cg3dNormaliseRHS
+
 C--   COMMON /CG3D_R/ DEL**2 Laplacian operators
 C     aW3d :: East-west operator.
 C     aS3d :: North-south operator.
 C     aV3d :: Vertical operator.
 C     aC3d :: 3D operator main diagonal term.
 C     zMC, zML, zMU :: preconditioner 3D solver
-C     cg3dNorm - A matrix normalisation factor.
-      COMMON /CG3D_R/
-     &      aW3d, aS3d, aV3d,
-     &      aC3d,
-     &      zMC, zML, zMU,
-     &      cg3dNorm
+C     cg3dNorm :: A matrix normalisation factor.
+C     cg3dTolerance_sq :: square of cg3d solver Tolerance (units depends
+C                 on cg3dNormaliseRHS, solver-unit ^2 = (m2/s2)^2 or no unit)
+      COMMON /CG3D_RS/
+     &      aW3d, aS3d, aV3d, aC3d,
+     &      zMC, zML, zMU
+      COMMON /CG3D_RL/
+     &      cg3dNorm, cg3dTolerance_sq
       _RS  aW3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS  aS3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS  aV3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS  aC3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS  zMC  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS  zML  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS  zMU  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
-      _RL  cg3dNorm
+      _RL  cg3dNorm, cg3dTolerance_sq
 
 C--   COMMON /CG3D_WK_R/  Work array common block
-C     cg3d_q - Intermediate matrix-vector product term
-C     cg3d_r -   idem
-C     cg3d_s -   idem
+C     cg3d_q :: Intermediate matrix-vector product term
+C     cg3d_r ::   idem
+C     cg3d_s ::   idem
       COMMON /CG3D_WK_R/
      & cg3d_q, cg3d_r, cg3d_s
-      _RL  cg3d_q(1-1  :sNx+1  ,1-1  :sNy+1  ,Nr,nSx,nSy)
-      _RL  cg3d_r(1-1  :sNx+1  ,1-1  :sNy+1  ,Nr,nSx,nSy)
-      _RL  cg3d_s(1-1  :sNx+1  ,1-1  :sNy+1  ,Nr,nSx,nSy)
+      _RL  cg3d_q( 0:sNx+1, 0:sNy+1, Nr,nSx,nSy)
+      _RL  cg3d_r( 0:sNx+1, 0:sNy+1, Nr,nSx,nSy)
+      _RL  cg3d_s( 0:sNx+1, 0:sNy+1, Nr,nSx,nSy)
 
 #endif /* ALLOW_NONHYDROSTATIC */

model/inc/PARAMS.h

@@ -562,11 +562,13 @@ C                        & Last iteration, in addition multiple of dumpFreq iter
 
 C--   COMMON /PARM_R/ "Real" valued parameters used by the model.
 C     cg2dTargetResidual
-C          :: Target residual for cg2d solver; no unit (RHS normalisation)
+C          :: Target residual for cg2d solver ; no unit (RHS normalisation)
 C     cg2dTargetResWunit
-C          :: Target residual for cg2d solver; W unit (No RHS normalisation)
+C          :: Target residual for cg2d solver ; W unit (No RHS normalisation)
 C     cg3dTargetResidual
-C               :: Target residual for cg3d solver.
+C          :: Target residual for cg3d solver ; no unit (RHS normalisation)
+C     cg3dTargetResWunit
+C          :: Target residual for cg3d solver ; W unit (No RHS normalisation)
 C     cg2dpcOffDFac :: Averaging weight for preconditioner off-diagonal.
 C     Note. 20th May 1998
 C           I made a weird discovery! In the model paper we argue
@@ -815,7 +817,7 @@ C     phiEuler      :: Euler angle, rotation about original z-axis
 C     thetaEuler    :: Euler angle, rotation about new x-axis
 C     psiEuler      :: Euler angle, rotation about new z-axis
       COMMON /PARM_R/ cg2dTargetResidual, cg2dTargetResWunit,
-     & cg2dpcOffDFac, cg3dTargetResidual,
+     & cg2dpcOffDFac, cg3dTargetResidual, cg3dTargetResWunit,
      & delR, delRc, xgOrigin, ygOrigin, rSphere, recip_rSphere,
      & radius_fromHorizGrid, seaLev_Z, top_Pres, rSigmaBnd,
      & deltaT, deltaTMom, dTtracerLev, deltaTFreeSurf, deltaTClock,
@@ -864,6 +866,7 @@ C     psiEuler      :: Euler angle, rotation about new z-axis
       _RL cg2dTargetResidual
       _RL cg2dTargetResWunit
       _RL cg3dTargetResidual
+      _RL cg3dTargetResWunit
       _RL cg2dpcOffDFac
       _RL delR(Nr)
       _RL delRc(Nr+1)

model/src/cg2d.F

@@ -78,7 +78,6 @@ SUBROUTINE CG2D(
       INTEGER bi, bj
       INTEGER i, j, it2d
       INTEGER actualIts
-      _RL    cg2dTolerance_sq
       _RL    err_sq,  errTile(nSx,nSy)
       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
       _RL    eta_qrNM1
@@ -99,7 +98,6 @@ SUBROUTINE CG2D(
 CEOP
 
 C--   Initialise auxiliary constant, some output variable and inverter
-      cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
       minResidualSq = -1. _d 0
       eta_qrNM1     =  1. _d 0
 
@@ -118,19 +116,19 @@ SUBROUTINE CG2D(
 
       IF (cg2dNormaliseRHS) THEN
 C-  Normalise RHS :
-      _GLOBAL_MAX_RL( rhsMax, myThid )
-      rhsNorm = 1. _d 0
-      IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
-      DO bj=myByLo(myThid),myByHi(myThid)
-       DO bi=myBxLo(myThid),myBxHi(myThid)
-        DO j=1,sNy
-         DO i=1,sNx
-          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
-          cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
+       _GLOBAL_MAX_RL( rhsMax, myThid )
+       rhsNorm = 1. _d 0
+       IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
+       DO bj=myByLo(myThid),myByHi(myThid)
+        DO bi=myBxLo(myThid),myBxHi(myThid)
+         DO j=1,sNy
+          DO i=1,sNx
+           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
+           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
+          ENDDO
          ENDDO
         ENDDO
        ENDDO
-      ENDDO
 C- end Normalise RHS
       ENDIF
 
@@ -198,7 +196,7 @@ SUBROUTINE CG2D(
       IF ( debugLevel .GE. debLevZero ) THEN
         _BEGIN_MASTER( myThid )
         printResidual = printResidualFreq.GE.1
-        WRITE(standardmessageunit,'(A,1P2E22.14)')
+        WRITE(standardMessageUnit,'(A,1P2E22.14)')
      &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
         _END_MASTER( myThid )
       ENDIF
@@ -372,11 +370,13 @@ SUBROUTINE CG2D(
 
       IF (cg2dNormaliseRHS) THEN
 C--   Un-normalise the answer
+c       rhsMax = 1. _d 0 / rhsNorm
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1,sNy
            DO i=1,sNx
             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
+c           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsMax
            ENDDO
           ENDDO
          ENDDO

model/src/cg2d_ex0.F

@@ -82,7 +82,6 @@ SUBROUTINE CG2D_EX0(
       INTEGER i, j, it2d
       INTEGER actualIts(nSx,nSy)
       INTEGER minResIter(nSx,nSy)
-      _RL    cg2dTolerance_sq
       _RL    err_sq,  errTile(nSx,nSy)
       _RL    eta_qrNtile(nSx,nSy)
       _RL    eta_qrNM1(nSx,nSy)
@@ -101,7 +100,6 @@ SUBROUTINE CG2D_EX0(
 CEOP
 
 C--   Initialise auxiliary constant, some output variable
-      cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
 
 C--   Initialise inverter and Normalise RHS
       rhsMax = 0. _d 0
@@ -194,7 +192,7 @@ SUBROUTINE CG2D_EX0(
       IF ( debugLevel .GE. debLevZero ) THEN
         _BEGIN_MASTER( myThid )
         printResidual = printResidualFreq.GE.1
-        WRITE(standardmessageunit,'(A,1P2E22.14)')
+        WRITE(standardMessageUnit,'(A,1P2E22.14)')
      &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
         _END_MASTER( myThid )
       ENDIF

model/src/cg2d_nsa.F

@@ -97,7 +97,6 @@ SUBROUTINE CG2D_NSA(
 #else
       _RL    rhsNorm
 #endif
-      _RL    cg2dTolerance_sq
       _RL    err_sq,  errTile(nSx,nSy)
       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
       _RL    eta_qrNM1, recip_eta_qrNM1
@@ -133,7 +132,6 @@ SUBROUTINE CG2D_NSA(
 #endif /* ALLOW_AUTODIFF_TAMC */
 
 C--   Initialise auxiliary constant, some output variable and inverter
-      cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
       minResidualSq = -1. _d 0
       eta_qrNM1     =  1. _d 0
       recip_eta_qrNM1= 1. _d 0
@@ -228,7 +226,7 @@ SUBROUTINE CG2D_NSA(
       IF ( debugLevel .GE. debLevZero ) THEN
         _BEGIN_MASTER( myThid )
         printResidual = printResidualFreq.GE.1
-        WRITE(standardmessageunit,'(A,1P2E22.14)')
+        WRITE(standardMessageUnit,'(A,1P2E22.14)')
      &  ' cg2d_nsa: Sum(rhs),rhsMax = ', sumRHS,rhsMax
         _END_MASTER( myThid )
       ENDIF

model/src/cg2d_sr.F

@@ -87,7 +87,6 @@ SUBROUTINE CG2D_SR(
       INTEGER bi, bj
       INTEGER i, j, it2d
 c     INTEGER actualIts
-      _RL    cg2dTolerance_sq
       _RL    err_sq,  errTile(nSx,nSy)
       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
       _RL    eta_qrNM1
@@ -110,7 +109,6 @@ SUBROUTINE CG2D_SR(
 CEOP
 
 C--   Initialise auxiliary constant, some output variable and inverter
-      cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
       minResidualSq = -1. _d 0
       eta_qrNM1     =  1. _d 0