more Momentum options, mostly related to hFac

model/src/calc_gw.F

@@ -122,7 +122,7 @@ SUBROUTINE CALC_GW(
       PARAMETER( jMin = 1 , jMax = sNy )
 CEOP
 #ifdef ALLOW_DIAGNOSTICS
-      LOGICAL diagDiss, diagAdvec, diag_AB
+      LOGICAL diagDiss, diagAdvec, diagMetric, diag_AB
       LOGICAL  DIAGNOSTICS_IS_ON
       EXTERNAL DIAGNOSTICS_IS_ON
 #endif /* ALLOW_DIAGNOSTICS */
@@ -133,10 +133,12 @@ SUBROUTINE CALC_GW(
       IF ( useDiagnostics ) THEN
         diagDiss  = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid )
         diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid )
+        diagMetric= DIAGNOSTICS_IS_ON( 'Wm_Metr ', myThid )
         diag_AB   = DIAGNOSTICS_IS_ON( 'AB_gW   ', myThid )
       ELSE
         diagDiss  = .FALSE.
         diagAdvec = .FALSE.
+        diagMetric= .FALSE.
         diag_AB   = .FALSE.
       ENDIF
 #endif /* ALLOW_DIAGNOSTICS */
@@ -604,7 +606,7 @@ SUBROUTINE CALC_GW(
 
         IF ( useNHMTerms .AND. k.GT.1 ) THEN
           CALL MOM_W_METRIC_NH(
-     I               bi,bj,k,
+     I               bi, bj, k,
      I               uVel, vVel,
      O               gwAdd,
      I               myThid )
@@ -613,11 +615,20 @@ SUBROUTINE CALC_GW(
              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
            ENDDO
           ENDDO
+#ifdef ALLOW_DIAGNOSTICS
+          IF ( diagMetric ) THEN
+           CALL DIAGNOSTICS_FILL( gwAdd, 'Wm_Metr ',
+     I                            k,1,2, bi,bj, myThid )
+C- note: need to explicitly increment the counter since DIAGNOSTICS_FILL
+C        does it only if k=1 (never the case here)
+           IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Metr ',bi,bj,myThid)
+          ENDIF
+#endif /* ALLOW_DIAGNOSTICS */
         ENDIF
         IF ( select3dCoriScheme.GE.1 .AND. k.GT.1 ) THEN
           CALL MOM_W_CORIOLIS_NH(
-     I               bi,bj,k,
-     I               uVel, vVel,
+     I               bi, bj, k,
+     I               uVel, vVel, recip_rThickC,
      O               gwAdd,
      I               myThid )
           DO j=jMin,jMax
@@ -633,8 +644,6 @@ SUBROUTINE CALC_GW(
         IF ( diagDiss )  THEN
           CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
      &                           k, 1, 2, bi,bj, myThid )
-C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
-C        does it only if k=1 (never the case here)
 c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
         ENDIF
         IF ( diagAdvec ) THEN

model/src/set_parms.F

@@ -63,8 +63,13 @@ SUBROUTINE SET_PARMS( myThid )
      &                     SQUEEZE_RIGHT, myThid )
       ENDIF
       IF ( usingCurvilinearGrid ) THEN
+C-    Horizontal metric terms not implemented for Curvilinear-Grid
         selectMetricTerms = 0
       ENDIF
+      IF ( usingCylindricalGrid ) THEN
+C-    No current alternative metric terms formulation for Cylindrical-Grid
+        selectMetricTerms = MIN( selectMetricTerms, 1 )
+      ENDIF
       IF ( selectCoriMap.EQ.-1 ) THEN
         IF ( usingCartesianGrid.OR.usingCylindricalGrid ) THEN
 C       default is to use Beta-Plane Coriolis
@@ -92,6 +97,8 @@ SUBROUTINE SET_PARMS( myThid )
       momPressureForcing= momStepping .AND. momPressureForcing
       implicitIntGravWave=momPressureForcing .AND. implicitIntGravWave
       momImplVertAdv   = momAdvection .AND. momImplVertAdv
+      useNHMTerms      = momAdvection .AND. useNHMTerms
+      IF ( .NOT.momAdvection ) selectMetricTerms = 0
       implicitViscosity= momViscosity .AND. implicitViscosity
       useSmag3D        = momViscosity .AND. useSmag3D
       use3Dsolver      = nonHydrostatic.OR. implicitIntGravWave

pkg/mom_common/MOM_COMMON_OPTIONS.h

@@ -26,5 +26,10 @@ C as a function of grid cell thickness and roughness length
 C zRoughBot (order 0.01m), assuming a von Karman constant = 0.4.
 #undef ALLOW_BOTTOMDRAG_ROUGHNESS
 
+C-- Deep-Model: use the original vertical advection (with all NH metric terms)
+C   instead of updated version that advects the product deepFac x (u,v)
+C   thus removing the need for NH-metric terms: w.(u,v)/r
+#undef MOM_USE_OLD_DEEP_VERT_ADV
+
 #endif /* ALLOW_MOM_COMMON */
 #endif /* MOM_COMMON_OPTIONS_H */

pkg/mom_common/MOM_VISC.h

@@ -20,6 +20,12 @@ C     useVariableVisc   :: variable (in space or time) viscosity is used
       _RL L4rdt_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL L4rdt_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 
+C--   COMMON /MOM_FLUXFORM_GRID/ just to hold a copy of deepFacC
+C     deepFacAdv :: copy of deepFacC or just Nr-length vector of 1 when
+C                   MOM_USE_OLD_DEEP_VERT_ADV is defined or useNHMTerms=F
+      COMMON /MOM_GRID_COPY/ deepFacAdv
+      _RL deepFacAdv(Nr)
+
 #ifdef ALLOW_SMAG_3D
 C     smag3D_hLsC :: horiz. grid length scale (power 2/3) at grid cell center
 C     smag3D_hLsW :: horiz. grid length scale (power 2/3) at western  edge
@@ -64,7 +70,7 @@ C     viscA4_W  :: Horizontal biharmonic viscosity for vertical momentum
 #ifdef ALLOW_BOTTOMDRAG_ROUGHNESS
 C--   bottom drag coefficents as a function of grid cell thickness
 C     and roughness length
-      COMMON /GRID_DRAGCOEFFS_RS/
+      COMMON /MOM_DRAGCOEFFS_RS/
      &     bottomDragCoeffW, bottomDragCoeffS
       _RS bottomDragCoeffW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RS bottomDragCoeffS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

pkg/mom_common/mom_diagnostics_init.F

@@ -413,13 +413,32 @@ SUBROUTINE MOM_DIAGNOSTICS_INIT( myThid )
 
 #endif /* ALLOW_MOM_VECINV */
 
+C-    diagnostics of metric-terms tendencies:
+      diagName  = 'Um_Metr '
+      diagTitle = 'U momentum tendency from Metric terms'
+      diagCode  = 'UUR     MR      '
+      diagMate  = diagNum + 2
+      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
+     I   diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
+      diagName  = 'Vm_Metr '
+      diagTitle = 'V momentum tendency from Metric terms'
+      diagCode  = 'VVR     MR      '
+      diagMate  = diagNum
+      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
+     I   diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
+
 #ifdef ALLOW_NONHYDROSTATIC
 C-    vertical momentum tendencies
       IF ( usingPCoords ) THEN
         diagUnits = 'Pa/s^2          '
       ELSE
         diagUnits = 'm/s^2           '
       ENDIF
+      diagName  = 'WSidDrag'
+      diagTitle = 'Vertical momentum tendency from Side Drag'
+      diagCode  = 'WMr     LR      '
+      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
+     I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
       diagName  = 'Wm_Diss '
       diagTitle = 'W momentum tendency from Dissipation'
       diagCode  = 'WMr     LR      '
@@ -431,9 +450,8 @@ SUBROUTINE MOM_DIAGNOSTICS_INIT( myThid )
       diagCode  = 'WMr     LR      '
       CALL DIAGNOSTICS_ADDTOLIST( diagNum,
      I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
-
-      diagName  = 'WSidDrag'
-      diagTitle = 'Vertical momentum tendency from Side Drag'
+      diagName  = 'Wm_Metr '
+      diagTitle = 'W momentum tendency from Metric terms'
       diagCode  = 'WMr     LR      '
       CALL DIAGNOSTICS_ADDTOLIST( diagNum,
      I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )

pkg/mom_common/mom_init_fixed.F

@@ -38,11 +38,24 @@ SUBROUTINE MOM_INIT_FIXED( myThid )
       CHARACTER*(MAX_LEN_MBUF) msgBuf
 #endif
 
-      k = 1
       twoThird = 2. _d 0 / 3. _d 0
       recip_dt = 1. _d 0
       IF ( deltaTMom.NE.0. ) recip_dt = 1. _d 0/deltaTMom
 
+      _BEGIN_MASTER(myThid)
+      IF ( useNHMTerms ) THEN
+        DO k=1,Nr
+         deepFacAdv(k) = deepFacC(k)
+        ENDDO
+#ifndef MOM_USE_OLD_DEEP_VERT_ADV
+      ELSE
+#endif
+        DO k=1,Nr
+         deepFacAdv(k) = 1. _d 0
+        ENDDO
+      ENDIF
+      _END_MASTER(myThid)
+
       DO bj=myByLo(myThid), myByHi(myThid)
        DO bi=myBxLo(myThid), myBxHi(myThid)
 

pkg/mom_common/mom_u_coriolis_nh.F

@@ -5,7 +5,7 @@
 
 C !INTERFACE: ==========================================================
       SUBROUTINE MOM_U_CORIOLIS_NH(
-     I               bi,bj,k,wFld,
+     I               bi, bj, k, wFld,
      O               uCoriolisTerm,
      I               myThid )
 
@@ -24,11 +24,11 @@ SUBROUTINE MOM_U_CORIOLIS_NH(
 #include "GRID.h"
 
 C !INPUT PARAMETERS: ===================================================
-C  bi,bj                :: tile indices
+C  bi, bj               :: tile indices
 C  k                    :: vertical level
 C  wFld                 :: vertical flow
 C  myThid               :: thread number
-      INTEGER bi,bj,k
+      INTEGER bi, bj, k
       _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       INTEGER myThid
 
@@ -37,29 +37,48 @@ SUBROUTINE MOM_U_CORIOLIS_NH(
       _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 
 C !LOCAL VARIABLES: ====================================================
-C  i,j                  :: loop indices
-      INTEGER i,j,kp1
+C  i, j                 :: loop indices
+      INTEGER i, j, kp1
       _RL wMsk
 CEOP
 
+      kp1 = MIN(k+1,Nr)
+      wMsk = 1.
+      IF (k.EQ.Nr) wMsk = 0.
 
-      kp1=min(k+1,Nr)
-      wMsk=1.
-      IF (k.EQ.Nr) wMsk=0.
-
-C Energy conserving discretization of 2*Omega*cos(phi)*w
-      DO j=1-Oly,sNy+Oly
-       DO i=2-Olx,sNx+Olx
-        uCoriolisTerm(i,j) =
-     &    0.5*( fCoriCos( i ,j,bi,bj)*angleCosC( i ,j,bi,bj)
-     &         *0.5*( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k )
-     &               +wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
-     &        + fCoriCos(i-1,j,bi,bj)*angleCosC(i-1,j,bi,bj)
-     &         *0.5*( wFld(i-1,j, k ,bi,bj)*rVel2wUnit( k )
-     &               +wFld(i-1,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
-     &        )*gravitySign
+      IF ( select3dCoriScheme.EQ.1 ) THEN
+C-    Original discretization of 2*Omega*cos(phi)*w
+C     documented as Energy conserving
+       DO j=1-OLy,sNy+OLy
+        DO i=2-OLx,sNx+OLx
+         uCoriolisTerm(i,j) = gravitySign*halfRL
+     &       *( fCoriCos( i ,j,bi,bj)*angleCosC( i ,j,bi,bj)*halfRL
+     &           *( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k )
+     &            + wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
+     &        + fCoriCos(i-1,j,bi,bj)*angleCosC(i-1,j,bi,bj)*halfRL
+     &           *( wFld(i-1,j, k ,bi,bj)*rVel2wUnit( k )
+     &            + wFld(i-1,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
+     &        )
+        ENDDO
+       ENDDO
+      ELSE
+C-    Using averaged transport:
+       DO j=1-OLy,sNy+OLy
+        DO i=2-OLx,sNx+OLx
+         uCoriolisTerm(i,j) = gravitySign*halfRL
+     &    *( fCoriCos( i ,j,bi,bj)*angleCosC( i ,j,bi,bj)
+     &      *( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k )
+     &       + wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk
+     &       )*rA( i ,j,bi,bj)*halfRL
+     &     + fCoriCos(i-1,j,bi,bj)*angleCosC(i-1,j,bi,bj)
+     &      *( wFld(i-1,j, k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k )
+     &       + wFld(i-1,j,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk
+     &       )*rA(i-1,j,bi,bj)*halfRL
+     &     )*recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
+     &      *recip_hFacW(i,j,k,bi,bj)
+        ENDDO
        ENDDO
-      ENDDO
+      ENDIF
 
       RETURN
       END

pkg/mom_common/mom_v_coriolis_nh.F

@@ -5,7 +5,7 @@
 
 C !INTERFACE: ==========================================================
       SUBROUTINE MOM_V_CORIOLIS_NH(
-     I               bi,bj,k,wFld,
+     I               bi, bj, k, wFld,
      O               vCoriolisTerm,
      I               myThid )
 
@@ -24,11 +24,11 @@ SUBROUTINE MOM_V_CORIOLIS_NH(
 #include "GRID.h"
 
 C !INPUT PARAMETERS: ===================================================
-C  bi,bj                :: tile indices
+C  bi, bj               :: tile indices
 C  k                    :: vertical level
 C  wFld                 :: vertical flow
 C  myThid               :: thread number
-      INTEGER bi,bj,k
+      INTEGER bi, bj, k
       _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       INTEGER myThid
 
@@ -37,29 +37,48 @@ SUBROUTINE MOM_V_CORIOLIS_NH(
       _RL vCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 
 C !LOCAL VARIABLES: ====================================================
-C  i,j                  :: loop indices
-      INTEGER i,j,kp1
+C  i, j                 :: loop indices
+      INTEGER i, j, kp1
       _RL wMsk
 CEOP
 
+      kp1 = MIN(k+1,Nr)
+      wMsk = 1.
+      IF (k.EQ.Nr) wMsk = 0.
 
-      kp1=min(k+1,Nr)
-      wMsk=1.
-      IF (k.EQ.Nr) wMsk=0.
-
-C Energy conserving discretization of 2*Omega*cos(phi)*w
-      DO j=2-Oly,sNy+Oly
-       DO i=1-Olx,sNx+Olx
-        vCoriolisTerm(i,j) =
-     &   -0.5*( fCoriCos( i ,j,bi,bj)*angleSinC(i, j ,bi,bj)
-     &         *0.5*( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k )
-     &               +wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
-     &        + fCoriCos(i,j-1,bi,bj)*angleSinC(i,j-1,bi,bj)
-     &         *0.5*( wFld(i,j-1, k ,bi,bj)*rVel2wUnit( k )
-     &               +wFld(i,j-1,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
-     &        )*gravitySign
+      IF ( select3dCoriScheme.EQ.1 ) THEN
+C-    Original discretization of 2*Omega*cos(phi)*w
+C     documented as Energy conserving
+       DO j=2-OLy,sNy+OLy
+        DO i=1-OLx,sNx+OLx
+         vCoriolisTerm(i,j) = -gravitySign*halfRL
+     &       *( fCoriCos(i, j ,bi,bj)*angleSinC(i, j ,bi,bj)*halfRL
+     &           *( wFld(i, j , k ,bi,bj)*rVel2wUnit( k )
+     &            + wFld(i, j ,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
+     &        + fCoriCos(i,j-1,bi,bj)*angleSinC(i,j-1,bi,bj)*halfRL
+     &           *( wFld(i,j-1, k ,bi,bj)*rVel2wUnit( k )
+     &            + wFld(i,j-1,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk )
+     &        )
+        ENDDO
+       ENDDO
+      ELSE
+C-    Using averaged transport:
+       DO j=2-OLy,sNy+OLy
+        DO i=1-OLx,sNx+OLx
+         vCoriolisTerm(i,j) = -gravitySign*halfRL
+     &    *( fCoriCos(i, j ,bi,bj)*angleSinC(i, j ,bi,bj)
+     &      *( wFld(i, j , k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k )
+     &       + wFld(i, j ,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk
+     &         )*rA(i, j ,bi,bj)*halfRL
+     &     + fCoriCos(i,j-1,bi,bj)*angleSinC(i,j-1,bi,bj)
+     &      *( wFld(i,j-1, k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k )
+     &       + wFld(i,j-1,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk
+     &         )*rA(i,j-1,bi,bj)*halfRL
+     &     )*recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
+     &      *recip_hFacS(i,j,k,bi,bj)
+        ENDDO
        ENDDO
-      ENDDO
+      ENDIF
 
       RETURN
       END