Back to home page

MITgcm
 
 

    

File indexing completed on 2025-05-05 05:07:59 UTC

view on githubraw file Latest commit 31fb0e0e on 2025-05-05 02:15:14 UTC
a25a2def2e Jean*0001 #include "PACKAGES_CONFIG.h"
88830be691 Alis*0002 #include "CPP_OPTIONS.h"
0bd9d32119 Jean*0003 #ifdef ALLOW_MOM_COMMON
                0004 # include "MOM_COMMON_OPTIONS.h"
                0005 #endif
53169f2ab0 Jean*0006 #define CALC_GW_NEW_THICK
88830be691 Alis*0007 
9366854e02 Chri*0008 CBOP
                0009 C     !ROUTINE: CALC_GW
                0010 C     !INTERFACE:
52c561f327 Jean*0011       SUBROUTINE CALC_GW(
efd2c352d2 Jean*0012      I               bi, bj, kappaRU, kappaRV,
efc81681f0 Jean*0013      I               str13, str23, str33,
                0014      I               viscAh3d_00, viscAh3d_13, viscAh3d_23,
d455cbf76e Jean*0015      I               myTime, myIter, myThid )
9366854e02 Chri*0016 C     !DESCRIPTION: \bv
                0017 C     *==========================================================*
52c561f327 Jean*0018 C     | S/R CALC_GW
597720f5c5 Jean*0019 C     | o Calculate vertical velocity tendency terms
                0020 C     |   ( Non-Hydrostatic only )
9366854e02 Chri*0021 C     *==========================================================*
597720f5c5 Jean*0022 C     | In NH, the vertical momentum tendency must be
52c561f327 Jean*0023 C     | calculated explicitly and included as a source term
9366854e02 Chri*0024 C     | for a 3d pressure eqn. Calculate that term here.
                0025 C     | This routine is not used in HYD calculations.
                0026 C     *==========================================================*
52c561f327 Jean*0027 C     \ev
88830be691 Alis*0028 
9366854e02 Chri*0029 C     !USES:
                0030       IMPLICIT NONE
88830be691 Alis*0031 C     == Global variables ==
                0032 #include "SIZE.h"
                0033 #include "EEPARAMS.h"
                0034 #include "PARAMS.h"
                0035 #include "GRID.h"
117ee419f5 Jean*0036 #include "RESTART.h"
1487584115 Jean*0037 #include "SURFACE.h"
d1b81ea0bc Jean*0038 #include "DYNVARS.h"
                0039 #include "NH_VARS.h"
229ac9feb6 Jean*0040 #ifdef ALLOW_ADDFLUID
                0041 # include "FFIELDS.h"
                0042 #endif
0bd9d32119 Jean*0043 #ifdef ALLOW_MOM_COMMON
                0044 # include "MOM_VISC.h"
                0045 #endif
88830be691 Alis*0046 
9366854e02 Chri*0047 C     !INPUT/OUTPUT PARAMETERS:
88830be691 Alis*0048 C     == Routine arguments ==
d455cbf76e Jean*0049 C     bi,bj   :: current tile indices
efd2c352d2 Jean*0050 C     kappaRU :: vertical viscosity at U points
                0051 C     kappaRV :: vertical viscosity at V points
d455cbf76e Jean*0052 C     myTime  :: Current time in simulation
                0053 C     myIter  :: Current iteration number in simulation
                0054 C     myThid  :: Thread number for this instance of the routine.
                0055       INTEGER bi,bj
efd2c352d2 Jean*0056       _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0057       _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
efc81681f0 Jean*0058 #ifdef ALLOW_SMAG_3D
                0059       _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0060       _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0061       _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0062       _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0063       _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0064       _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0065 #else /* ALLOW_SMAG_3D */
                0066       _RL str13(1), str23(1), str33(1)
                0067       _RL viscAh3d_00(1), viscAh3d_13(1), viscAh3d_23(1)
                0068 #endif /* ALLOW_SMAG_3D */
c1b1e715b8 Jean*0069       _RL     myTime
                0070       INTEGER myIter
88830be691 Alis*0071       INTEGER myThid
                0072 
87f515096d Alis*0073 #ifdef ALLOW_NONHYDROSTATIC
0bd9d32119 Jean*0074 #ifdef ALLOW_MOM_COMMON
87f515096d Alis*0075 
9366854e02 Chri*0076 C     !LOCAL VARIABLES:
88830be691 Alis*0077 C     == Local variables ==
0bd9d32119 Jean*0078 C     biharmonicVisc:: use horizontal biharmonic viscosity for vertical momentum
efc81681f0 Jean*0079 C     iMin, iMax    :: Ranges and sub-block indices on which calculations
                0080 C     jMin, jMax       are applied.
597720f5c5 Jean*0081 C     xA            :: W-Cell face area normal to X
                0082 C     yA            :: W-Cell face area normal to Y
                0083 C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
                0084 C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
1487584115 Jean*0085 C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
427e24e121 Jean*0086 C     recip_rThickC :: reciprocal thickness of W-Cell (centered on W-point)
597720f5c5 Jean*0087 C     flx_NS        :: vertical momentum flux, meridional direction
                0088 C     flx_EW        :: vertical momentum flux, zonal direction
                0089 C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
                0090 C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
                0091 C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
                0092 C     gwDiss        :: vertical momentum dissipation tendency
50603588ee Jean*0093 C     gwAdd         :: other tendencies (Coriolis, Metric-terms)
b88f9ec5cb Jean*0094 C     gw_AB         :: tendency increment from Adams-Bashforth
50603588ee Jean*0095 C     del2w         :: laplacian of wVel
                0096 C     wFld          :: local copy of wVel
597720f5c5 Jean*0097 C     i,j,k         :: Loop counters
0bd9d32119 Jean*0098       LOGICAL biharmonicVisc
d455cbf76e Jean*0099       INTEGER iMin,iMax,jMin,jMax
597720f5c5 Jean*0100       _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101       _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0102       _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0103       _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1487584115 Jean*0104       _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
597720f5c5 Jean*0105       _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d455cbf76e Jean*0106       _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
597720f5c5 Jean*0109       _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110       _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111       _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112       _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b88f9ec5cb Jean*0113       _RL    gw_AB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d455cbf76e Jean*0114       _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50603588ee Jean*0115       _RL    wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
982e105a17 Jean*0116       INTEGER i,j,k, km1, kp1
                0117       _RL  mskM1, mskP1
597720f5c5 Jean*0118       _RL  tmp_WbarZ
                0119       _RL  uTrans, vTrans, rTrans
d455cbf76e Jean*0120       _RL  viscLoc
53169f2ab0 Jean*0121       PARAMETER( iMin = 1 , iMax = sNx )
                0122       PARAMETER( jMin = 1 , jMax = sNy )
9366854e02 Chri*0123 CEOP
41d43b2d0d Jean*0124 #ifdef ALLOW_DIAGNOSTICS
31fb0e0e6d Jean*0125       LOGICAL diagDiss, diagAdvec, diagMetric, diag_AB
41d43b2d0d Jean*0126       LOGICAL  DIAGNOSTICS_IS_ON
                0127       EXTERNAL DIAGNOSTICS_IS_ON
                0128 #endif /* ALLOW_DIAGNOSTICS */
88830be691 Alis*0129 
982e105a17 Jean*0130 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
52c561f327 Jean*0131 
41d43b2d0d Jean*0132 #ifdef ALLOW_DIAGNOSTICS
                0133       IF ( useDiagnostics ) THEN
                0134         diagDiss  = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid )
                0135         diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid )
31fb0e0e6d Jean*0136         diagMetric= DIAGNOSTICS_IS_ON( 'Wm_Metr ', myThid )
67662de2f2 Jean*0137         diag_AB   = DIAGNOSTICS_IS_ON( 'AB_gW   ', myThid )
41d43b2d0d Jean*0138       ELSE
                0139         diagDiss  = .FALSE.
                0140         diagAdvec = .FALSE.
31fb0e0e6d Jean*0141         diagMetric= .FALSE.
b88f9ec5cb Jean*0142         diag_AB   = .FALSE.
41d43b2d0d Jean*0143       ENDIF
                0144 #endif /* ALLOW_DIAGNOSTICS */
                0145 
0bd9d32119 Jean*0146       biharmonicVisc = viscA4W.NE.zeroRL
                0147      &           .OR. ( useVariableVisc .AND. useBiharmonicVisc )
                0148 
d455cbf76e Jean*0149 C--   Initialise gW to zero
                0150       DO k=1,Nr
                0151         DO j=1-OLy,sNy+OLy
                0152          DO i=1-OLx,sNx+OLx
88830be691 Alis*0153            gW(i,j,k,bi,bj) = 0.
                0154          ENDDO
                0155         ENDDO
d455cbf76e Jean*0156       ENDDO
597720f5c5 Jean*0157 C-    Initialise gwDiss to zero
                0158       DO j=1-OLy,sNy+OLy
                0159        DO i=1-OLx,sNx+OLx
                0160          gwDiss(i,j) = 0.
                0161        ENDDO
                0162       ENDDO
50603588ee Jean*0163       IF (momViscosity) THEN
                0164 C-    Initialize del2w to zero:
935bd32a21 Jean*0165         DO j=1-OLy,sNy+OLy
                0166          DO i=1-OLx,sNx+OLx
50603588ee Jean*0167            del2w(i,j) = 0. _d 0
                0168          ENDDO
                0169         ENDDO
                0170       ENDIF
88830be691 Alis*0171 
9091e25e07 Jean*0172 C--   Boundaries condition at top (vertical advection of vertical momentum):
982e105a17 Jean*0173       DO j=1-OLy,sNy+OLy
                0174        DO i=1-OLx,sNx+OLx
                0175          flxAdvUp(i,j) = 0.
                0176 c        flxDisUp(i,j) = 0.
                0177        ENDDO
                0178       ENDDO
                0179 
d455cbf76e Jean*0180 C---  Sweep down column
982e105a17 Jean*0181       DO k=1,Nr
                0182         km1 = MAX( k-1, 1 )
                0183         kp1 = MIN( k+1,Nr )
                0184         mskM1 = 1.
                0185         mskP1 = 1.
                0186         IF ( k.EQ. 1 ) mskM1 = 0.
                0187         IF ( k.EQ.Nr ) mskP1 = 0.
                0188         IF ( k.GT.1 ) THEN
597720f5c5 Jean*0189 C--   Compute grid factor arround a W-point:
53169f2ab0 Jean*0190 #ifdef CALC_GW_NEW_THICK
935bd32a21 Jean*0191          DO j=1-OLy,sNy+OLy
                0192           DO i=1-OLx,sNx+OLx
53169f2ab0 Jean*0193            IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
                0194      &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
                0195              recip_rThickC(i,j) = 0.
                0196            ELSE
                0197 C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
                0198              recip_rThickC(i,j) = 1. _d 0 /
                0199      &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
                0200      &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
                0201      &        )
                0202            ENDIF
982e105a17 Jean*0203           ENDDO
53169f2ab0 Jean*0204          ENDDO
982e105a17 Jean*0205          IF (momViscosity) THEN
935bd32a21 Jean*0206          DO j=1-OLy,sNy+OLy
                0207           DO i=1-OLx,sNx+OLx
53169f2ab0 Jean*0208            rThickC_C(i,j) = MAX( zeroRS,
                0209      &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
                0210      &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
                0211      &                         )
                0212           ENDDO
                0213          ENDDO
935bd32a21 Jean*0214          DO j=1-OLy,sNy+OLy
                0215           DO i=1-OLx+1,sNx+OLx
53169f2ab0 Jean*0216            rThickC_W(i,j) = MAX( zeroRS,
50603588ee Jean*0217      &                           MIN( rSurfW(i,j,bi,bj), rC(k-1) )
                0218      &                          -MAX(  rLowW(i,j,bi,bj), rC(k)   )
53169f2ab0 Jean*0219      &                         )
                0220 C     W-Cell Western face area:
                0221            xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
                0222 c    &                              *deepFacF(k)
                0223           ENDDO
                0224          ENDDO
935bd32a21 Jean*0225          DO j=1-OLy+1,sNy+OLy
                0226           DO i=1-OLx,sNx+OLx
53169f2ab0 Jean*0227            rThickC_S(i,j) = MAX( zeroRS,
50603588ee Jean*0228      &                           MIN( rSurfS(i,j,bi,bj), rC(k-1) )
                0229      &                          -MAX(  rLowS(i,j,bi,bj), rC(k)   )
53169f2ab0 Jean*0230      &                         )
                0231 C     W-Cell Southern face area:
                0232            yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
                0233 c    &                              *deepFacF(k)
                0234 C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
                0235 C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
                0236           ENDDO
                0237          ENDDO
982e105a17 Jean*0238          ENDIF
53169f2ab0 Jean*0239 #else /* CALC_GW_NEW_THICK */
935bd32a21 Jean*0240          DO j=1-OLy,sNy+OLy
                0241           DO i=1-OLx,sNx+OLx
597720f5c5 Jean*0242 C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
                0243            IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
                0244              recip_rThickC(i,j) = 0.
                0245            ELSE
                0246              recip_rThickC(i,j) = 1. _d 0 /
85162e9502 Jean*0247      &        ( drF(k-1)*halfRS
597720f5c5 Jean*0248      &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
                0249      &        )
                0250            ENDIF
1487584115 Jean*0251 c          IF (momViscosity) THEN
                0252 #ifdef NONLIN_FRSURF
4606c28752 Jean*0253            rThickC_C(i,j) =
1487584115 Jean*0254      &          drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0255      &        + drF( k )*MIN( h0FacC(i,j,k  ,bi,bj), halfRS )
                0256 #else
4606c28752 Jean*0257            rThickC_C(i,j) =
1487584115 Jean*0258      &          drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0259      &        + drF( k )*MIN( _hFacC(i,j,k  ,bi,bj), halfRS )
                0260 #endif
                0261            rThickC_W(i,j) =
                0262      &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0263      &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
                0264            rThickC_S(i,j) =
                0265      &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0266      &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
597720f5c5 Jean*0267 C     W-Cell Western face area:
                0268            xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
4606c28752 Jean*0269 c    &                              *deepFacF(k)
597720f5c5 Jean*0270 C     W-Cell Southern face area:
                0271            yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
4606c28752 Jean*0272 c    &                              *deepFacF(k)
                0273 C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
                0274 C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
1487584115 Jean*0275 c          ENDIF
982e105a17 Jean*0276           ENDDO
597720f5c5 Jean*0277          ENDDO
53169f2ab0 Jean*0278 #endif /* CALC_GW_NEW_THICK */
982e105a17 Jean*0279         ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
935bd32a21 Jean*0280          DO j=1-OLy,sNy+OLy
                0281           DO i=1-OLx,sNx+OLx
982e105a17 Jean*0282            recip_rThickC(i,j) = recip_drC(k)
                0283 c          rThickC_C(i,j) = drC(k)
                0284 c          rThickC_W(i,j) = drC(k)
                0285 c          rThickC_S(i,j) = drC(k)
                0286 c          xA(i,j) = _dyG(i,j,bi,bj)*drC(k)
                0287 c          yA(i,j) = _dxG(i,j,bi,bj)*drC(k)
                0288           ENDDO
                0289          ENDDO
                0290         ENDIF
                0291 
                0292 C--   local copy of wVel:
935bd32a21 Jean*0293         DO j=1-OLy,sNy+OLy
                0294           DO i=1-OLx,sNx+OLx
982e105a17 Jean*0295             wFld(i,j) = wVel(i,j,k,bi,bj)
                0296           ENDDO
                0297         ENDDO
597720f5c5 Jean*0298 
d455cbf76e Jean*0299 C--   horizontal bi-harmonic dissipation
0bd9d32119 Jean*0300         IF ( momViscosity .AND. k.GT.1 .AND. biharmonicVisc ) THEN
50603588ee Jean*0301 
d455cbf76e Jean*0302 C-    calculate the horizontal Laplacian of vertical flow
0b6cbae535 Mart*0303 C     Zonal flux d/dx W
50603588ee Jean*0304           IF ( useCubedSphereExchange ) THEN
                0305 C     to compute d/dx(W), fill corners with appropriate values:
                0306             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
                0307      &                                 wFld, bi,bj, myThid )
                0308           ENDIF
935bd32a21 Jean*0309           DO j=1-OLy,sNy+OLy
                0310            flx_EW(1-OLx,j)=0.
                0311            DO i=1-OLx+1,sNx+OLx
597720f5c5 Jean*0312             flx_EW(i,j) =
50603588ee Jean*0313      &               ( wFld(i,j) - wFld(i-1,j) )
597720f5c5 Jean*0314      &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
0b6cbae535 Mart*0315 #ifdef COSINEMETH_III
4606c28752 Jean*0316      &              *sqCosFacU(j,bi,bj)
0b6cbae535 Mart*0317 #endif
b09ee35a04 Jean*0318 #ifdef ALLOW_OBCS
                0319      &              *maskInW(i,j,bi,bj)
                0320 #endif
0b6cbae535 Mart*0321            ENDDO
d455cbf76e Jean*0322           ENDDO
50603588ee Jean*0323 
0b6cbae535 Mart*0324 C     Meridional flux d/dy W
50603588ee Jean*0325           IF ( useCubedSphereExchange ) THEN
                0326 C     to compute d/dy(W), fill corners with appropriate values:
                0327             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
                0328      &                                 wFld, bi,bj, myThid )
                0329           ENDIF
935bd32a21 Jean*0330           DO i=1-OLx,sNx+OLx
                0331            flx_NS(i,1-OLy)=0.
0b6cbae535 Mart*0332           ENDDO
935bd32a21 Jean*0333           DO j=1-OLy+1,sNy+OLy
                0334            DO i=1-OLx,sNx+OLx
597720f5c5 Jean*0335             flx_NS(i,j) =
50603588ee Jean*0336      &               ( wFld(i,j) - wFld(i,j-1) )
597720f5c5 Jean*0337      &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
0b6cbae535 Mart*0338 #ifdef ISOTROPIC_COS_SCALING
                0339 #ifdef COSINEMETH_III
597720f5c5 Jean*0340      &              *sqCosFacV(j,bi,bj)
0b6cbae535 Mart*0341 #endif
                0342 #endif
b09ee35a04 Jean*0343 #ifdef ALLOW_OBCS
                0344      &              *maskInS(i,j,bi,bj)
                0345 #endif
0b6cbae535 Mart*0346            ENDDO
                0347           ENDDO
52c561f327 Jean*0348 
0b6cbae535 Mart*0349 C     del^2 W
50603588ee Jean*0350 C     Divergence of horizontal fluxes
935bd32a21 Jean*0351           DO j=1-OLy,sNy+OLy-1
                0352            DO i=1-OLx,sNx+OLx-1
50603588ee Jean*0353             del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
                0354      &                    +( flx_NS(i,j+1)-flx_NS(i,j) )
597720f5c5 Jean*0355      &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
4606c28752 Jean*0356      &                    *recip_deepFac2F(k)
0b6cbae535 Mart*0357            ENDDO
                0358           ENDDO
50603588ee Jean*0359 C end if biharmonic viscosity
d455cbf76e Jean*0360         ENDIF
0b6cbae535 Mart*0361 
982e105a17 Jean*0362         IF ( momViscosity .AND. k.GT.1 ) THEN
597720f5c5 Jean*0363 C Viscous Flux on Western face
                0364           DO j=jMin,jMax
                0365            DO i=iMin,iMax+1
                0366              flx_EW(i,j)=
                0367      &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
                0368      &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
41a8c7589b Jean*0369      &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
1487584115 Jean*0370      &              *cosFacU(j,bi,bj)
597720f5c5 Jean*0371      &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
                0372      &              *(del2w(i,j)-del2w(i-1,j))
41a8c7589b Jean*0373      &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
0b6cbae535 Mart*0374 #ifdef COSINEMETH_III
597720f5c5 Jean*0375      &              *sqCosFacU(j,bi,bj)
0b6cbae535 Mart*0376 #else
1487584115 Jean*0377      &              *cosFacU(j,bi,bj)
0b6cbae535 Mart*0378 #endif
597720f5c5 Jean*0379            ENDDO
                0380           ENDDO
                0381 C Viscous Flux on Southern face
                0382           DO j=jMin,jMax+1
                0383            DO i=iMin,iMax
                0384              flx_NS(i,j)=
                0385      &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
                0386      &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
41a8c7589b Jean*0387      &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
1487584115 Jean*0388 #ifdef ISOTROPIC_COS_SCALING
                0389      &              *cosFacV(j,bi,bj)
                0390 #endif
597720f5c5 Jean*0391      &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
                0392      &              *(del2w(i,j)-del2w(i,j-1))
41a8c7589b Jean*0393      &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
597720f5c5 Jean*0394 #ifdef ISOTROPIC_COS_SCALING
0b6cbae535 Mart*0395 #ifdef COSINEMETH_III
597720f5c5 Jean*0396      &              *sqCosFacV(j,bi,bj)
52c561f327 Jean*0397 #else
1487584115 Jean*0398      &              *cosFacV(j,bi,bj)
0b6cbae535 Mart*0399 #endif
597720f5c5 Jean*0400 #endif
                0401            ENDDO
                0402           ENDDO
                0403 C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
                0404           DO j=jMin,jMax
                0405            DO i=iMin,iMax
                0406 C     Interpolate vert viscosity to center of tracer-cell (level k):
efd2c352d2 Jean*0407              viscLoc = ( kappaRU(i,j,k)  +kappaRU(i+1,j,k)
                0408      &                  +kappaRU(i,j,k+1)+kappaRU(i+1,j,k+1)
                0409      &                  +kappaRV(i,j,k)  +kappaRV(i,j+1,k)
                0410      &                  +kappaRV(i,j,k+1)+kappaRV(i,j+1,k+1)
597720f5c5 Jean*0411      &                 )*0.125 _d 0
                0412              flx_Dn(i,j) =
982e105a17 Jean*0413      &          - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
597720f5c5 Jean*0414      &                     -wVel(i,j, k ,bi,bj) )*rkSign
41a8c7589b Jean*0415      &                   *recip_drF(k)*rA(i,j,bi,bj)
4606c28752 Jean*0416      &                   *deepFac2C(k)*rhoFacC(k)
597720f5c5 Jean*0417            ENDDO
                0418           ENDDO
9091e25e07 Jean*0419           IF ( k.EQ.2 ) THEN
                0420 C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)
                0421            DO j=jMin,jMax
                0422             DO i=iMin,iMax
                0423 C     Interpolate horizontally (but not vertically) vert viscosity to center:
                0424 C     Although background visc. might be defined at k=1, this is not
                0425 C     generally true when using variable visc. (from vertical mixing scheme).
                0426 C     Therefore, no vert. interp. and only horizontal interpolation.
efd2c352d2 Jean*0427              viscLoc = ( kappaRU(i,j,k) + kappaRU(i+1,j,k)
                0428      &                  +kappaRV(i,j,k) + kappaRV(i,j+1,k)
9091e25e07 Jean*0429      &                 )*0.25 _d 0
                0430              flxDisUp(i,j) =
                0431      &          - viscLoc*( wVel(i,j, k ,bi,bj)
                0432      &                     -wVel(i,j,k-1,bi,bj) )*rkSign
                0433      &                   *recip_drF(k-1)*rA(i,j,bi,bj)
                0434      &                   *deepFac2C(k-1)*rhoFacC(k-1)
                0435 C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero)
                0436 c            flxDisUp(i,j) = 0.
                0437             ENDDO
                0438            ENDDO
                0439           ENDIF
597720f5c5 Jean*0440 C     Tendency is minus divergence of viscous fluxes:
4606c28752 Jean*0441 C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
597720f5c5 Jean*0442           DO j=jMin,jMax
                0443            DO i=iMin,iMax
                0444              gwDiss(i,j) =
41a8c7589b Jean*0445      &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
                0446      &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
                0447      &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
4606c28752 Jean*0448      &                                          *recip_rhoFacF(k)
41a8c7589b Jean*0449      &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
4606c28752 Jean*0450      &          *recip_deepFac2F(k)
d455cbf76e Jean*0451 C--        prepare for next level (k+1)
597720f5c5 Jean*0452              flxDisUp(i,j)=flx_Dn(i,j)
                0453            ENDDO
                0454           ENDDO
                0455         ENDIF
                0456 
982e105a17 Jean*0457         IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
597720f5c5 Jean*0458 C-     No-slip BCs impose a drag at walls...
1487584115 Jean*0459           CALL MOM_W_SIDEDRAG(
                0460      I               bi,bj,k,
                0461      I               wVel, del2w,
                0462      I               rThickC_C, recip_rThickC,
                0463      I               viscAh_W, viscA4_W,
                0464      O               gwAdd,
                0465      I               myThid )
                0466           DO j=jMin,jMax
                0467            DO i=iMin,iMax
                0468             gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
                0469            ENDDO
                0470           ENDDO
597720f5c5 Jean*0471         ENDIF
                0472 
efc81681f0 Jean*0473 #ifdef ALLOW_SMAG_3D
                0474         IF ( useSmag3D .AND. k.GT.1 ) THEN
                0475              CALL MOM_W_SMAG_3D(
                0476      I         str13, str23, str33,
                0477      I         viscAh3d_00, viscAh3d_13, viscAh3d_23,
                0478      I         rThickC_W, rThickC_S, rThickC_C, recip_rThickC,
                0479      O         gwAdd,
                0480      I         iMin,iMax,jMin,jMax, k, bi, bj, myThid )
                0481           DO j = jMin,jMax
                0482            DO i = iMin,iMax
                0483             gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
                0484            ENDDO
                0485           ENDDO
                0486         ENDIF
                0487 #endif /* ALLOW_SMAG_3D */
                0488 
597720f5c5 Jean*0489 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0490 
                0491         IF ( momAdvection ) THEN
982e105a17 Jean*0492 
                0493          IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
597720f5c5 Jean*0494 C Advective Flux on Western face
                0495           DO j=jMin,jMax
                0496            DO i=iMin,iMax+1
                0497 C     transport through Western face area:
                0498              uTrans = (
982e105a17 Jean*0499      &          drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
                0500      &                  *rhoFacC(km1)*mskM1
597720f5c5 Jean*0501      &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
4606c28752 Jean*0502      &                  *rhoFacC(k)
                0503      &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
982e105a17 Jean*0504              flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
                0505 c            flx_EW(i,j)=
                0506 c    &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
597720f5c5 Jean*0507            ENDDO
                0508           ENDDO
                0509 C Advective Flux on Southern face
                0510           DO j=jMin,jMax+1
                0511            DO i=iMin,iMax
                0512 C     transport through Southern face area:
                0513              vTrans = (
982e105a17 Jean*0514      &          drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
                0515      &                  *rhoFacC(km1)*mskM1
597720f5c5 Jean*0516      &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
4606c28752 Jean*0517      &                  *rhoFacC(k)
                0518      &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
982e105a17 Jean*0519              flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
                0520 c            flx_NS(i,j)=
                0521 c    &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
597720f5c5 Jean*0522            ENDDO
                0523           ENDDO
982e105a17 Jean*0524          ENDIF
597720f5c5 Jean*0525 C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
982e105a17 Jean*0526 c        IF (.TRUE.) THEN
597720f5c5 Jean*0527           DO j=jMin,jMax
                0528            DO i=iMin,iMax
53169f2ab0 Jean*0529 C     NH in p-coord.: advect wSpeed [m/s] with rTrans
                0530              tmp_WbarZ = halfRL*
982e105a17 Jean*0531      &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k )
                0532      &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 )
597720f5c5 Jean*0533 C     transport through Lower face area:
4606c28752 Jean*0534              rTrans = halfRL*
                0535      &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
                0536      &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
982e105a17 Jean*0537      &                                   *mskP1
4606c28752 Jean*0538      &              )*rA(i,j,bi,bj)
41a8c7589b Jean*0539              flx_Dn(i,j) = rTrans*tmp_WbarZ
597720f5c5 Jean*0540            ENDDO
                0541           ENDDO
982e105a17 Jean*0542 c        ENDIF
                0543          IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
                0544 C Advective Flux on Upper face of W-Cell (= at surface)
9091e25e07 Jean*0545            DO j=jMin,jMax
                0546             DO i=iMin,iMax
982e105a17 Jean*0547              tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
                0548              rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
                0549      &               *rA(i,j,bi,bj)
9091e25e07 Jean*0550              flxAdvUp(i,j) = rTrans*tmp_WbarZ
                0551 c            flxAdvUp(i,j) = 0.
                0552             ENDDO
                0553            ENDDO
982e105a17 Jean*0554          ENDIF
935bd32a21 Jean*0555 
982e105a17 Jean*0556          IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
597720f5c5 Jean*0557 C     Tendency is minus divergence of advective fluxes:
4606c28752 Jean*0558 C     anelastic: all transports & advect. fluxes are scaled by rhoFac
597720f5c5 Jean*0559           DO j=jMin,jMax
                0560            DO i=iMin,iMax
89a3c0c31d Jean*0561 C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)
                0562 c            IF (k.EQ.2) flxAdvUp(i,j) = 0.
597720f5c5 Jean*0563              gW(i,j,k,bi,bj) =
41a8c7589b Jean*0564      &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
                0565      &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
53169f2ab0 Jean*0566      &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
41a8c7589b Jean*0567      &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
4606c28752 Jean*0568      &          *recip_deepFac2F(k)*recip_rhoFacF(k)
982e105a17 Jean*0569            ENDDO
                0570           ENDDO
935bd32a21 Jean*0571 #ifdef ALLOW_ADDFLUID
                0572           IF ( selectAddFluid.GE.1 ) THEN
                0573            DO j=jMin,jMax
                0574             DO i=iMin,iMax
                0575              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
                0576      &        + wVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
                0577      &          *( addMass(i,j,k,bi,bj)
                0578      &            +addMass(i,j,km1,bi,bj)*mskM1 )
                0579      &          *recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
                0580      &          *recip_deepFac2F(k)*recip_rhoFacF(k)
                0581             ENDDO
                0582            ENDDO
                0583           ENDIF
                0584 #endif /* ALLOW_ADDFLUID */
982e105a17 Jean*0585          ENDIF
                0586 
                0587          DO j=jMin,jMax
                0588            DO i=iMin,iMax
597720f5c5 Jean*0589 C--          prepare for next level (k+1)
                0590              flxAdvUp(i,j)=flx_Dn(i,j)
                0591            ENDDO
982e105a17 Jean*0592          ENDDO
                0593 
                0594 c       ELSE
                0595 C-    if momAdvection / else
                0596 c         DO j=1-OLy,sNy+OLy
                0597 c          DO i=1-OLx,sNx+OLx
                0598 c            gW(i,j,k,bi,bj) = 0. _d 0
                0599 c          ENDDO
                0600 c         ENDDO
                0601 
                0602 C-    endif momAdvection.
597720f5c5 Jean*0603         ENDIF
                0604 
982e105a17 Jean*0605 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0606 
                0607         IF ( useNHMTerms .AND. k.GT.1 ) THEN
597720f5c5 Jean*0608           CALL MOM_W_METRIC_NH(
31fb0e0e6d Jean*0609      I               bi, bj, k,
597720f5c5 Jean*0610      I               uVel, vVel,
                0611      O               gwAdd,
                0612      I               myThid )
                0613           DO j=jMin,jMax
                0614            DO i=iMin,iMax
                0615              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
                0616            ENDDO
                0617           ENDDO
31fb0e0e6d Jean*0618 #ifdef ALLOW_DIAGNOSTICS
                0619           IF ( diagMetric ) THEN
                0620            CALL DIAGNOSTICS_FILL( gwAdd, 'Wm_Metr ',
                0621      I                            k,1,2, bi,bj, myThid )
                0622 C- note: need to explicitly increment the counter since DIAGNOSTICS_FILL
                0623 C        does it only if k=1 (never the case here)
                0624            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Metr ',bi,bj,myThid)
                0625           ENDIF
                0626 #endif /* ALLOW_DIAGNOSTICS */
597720f5c5 Jean*0627         ENDIF
427e24e121 Jean*0628         IF ( select3dCoriScheme.GE.1 .AND. k.GT.1 ) THEN
597720f5c5 Jean*0629           CALL MOM_W_CORIOLIS_NH(
31fb0e0e6d Jean*0630      I               bi, bj, k,
                0631      I               uVel, vVel, recip_rThickC,
597720f5c5 Jean*0632      O               gwAdd,
                0633      I               myThid )
                0634           DO j=jMin,jMax
                0635            DO i=iMin,iMax
                0636              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
                0637            ENDDO
                0638           ENDDO
                0639         ENDIF
88830be691 Alis*0640 
52c561f327 Jean*0641 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ac07cf9a5a Jean*0642 
41d43b2d0d Jean*0643 #ifdef ALLOW_DIAGNOSTICS
50603588ee Jean*0644         IF ( diagDiss )  THEN
                0645           CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
                0646      &                           k, 1, 2, bi,bj, myThid )
982e105a17 Jean*0647 c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
50603588ee Jean*0648         ENDIF
                0649         IF ( diagAdvec ) THEN
                0650           CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
                0651      &                           k,Nr, 1, bi,bj, myThid )
982e105a17 Jean*0652 c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
50603588ee Jean*0653         ENDIF
41d43b2d0d Jean*0654 #endif /* ALLOW_DIAGNOSTICS */
                0655 
597720f5c5 Jean*0656 C--   Dissipation term inside the Adams-Bashforth:
                0657         IF ( momViscosity .AND. momDissip_In_AB) THEN
                0658           DO j=jMin,jMax
                0659            DO i=iMin,iMax
                0660              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
                0661            ENDDO
                0662           ENDDO
                0663         ENDIF
                0664 
52c561f327 Jean*0665 C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
                0666 C     and save gW_[n] into gwNm1 for the next time step.
cba4501825 Jean*0667 #ifdef ALLOW_ADAMSBASHFORTH_3
                0668         CALL ADAMS_BASHFORTH3(
94e3f14b29 Jean*0669      I                         bi, bj, k, Nr,
0b6eaf7dba Jean*0670      U                         gW(1-OLx,1-OLy,1,bi,bj), gwNm,
                0671      O                         gw_AB,
cba4501825 Jean*0672      I                         nHydStartAB, myIter, myThid )
                0673 #else /* ALLOW_ADAMSBASHFORTH_3 */
d455cbf76e Jean*0674         CALL ADAMS_BASHFORTH2(
94e3f14b29 Jean*0675      I                         bi, bj, k, Nr,
0b6eaf7dba Jean*0676      U                         gW(1-OLx,1-OLy,1,bi,bj),
                0677      U                         gwNm1(1-OLx,1-OLy,1,bi,bj),
                0678      O                         gw_AB,
117ee419f5 Jean*0679      I                         nHydStartAB, myIter, myThid )
cba4501825 Jean*0680 #endif /* ALLOW_ADAMSBASHFORTH_3 */
b88f9ec5cb Jean*0681 #ifdef ALLOW_DIAGNOSTICS
                0682         IF ( diag_AB ) THEN
67662de2f2 Jean*0683           CALL DIAGNOSTICS_FILL(gw_AB,'AB_gW   ',k,1,2,bi,bj,myThid)
b88f9ec5cb Jean*0684         ENDIF
                0685 #endif /* ALLOW_DIAGNOSTICS */
88830be691 Alis*0686 
597720f5c5 Jean*0687 C--   Dissipation term outside the Adams-Bashforth:
                0688         IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
                0689           DO j=jMin,jMax
                0690            DO i=iMin,iMax
                0691              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
                0692            ENDDO
                0693           ENDDO
                0694         ENDIF
                0695 
52c561f327 Jean*0696 C-    end of the k loop
                0697       ENDDO
ffb313c34a Alis*0698 
a25a2def2e Jean*0699 #ifdef ALLOW_DIAGNOSTICS
                0700       IF (useDiagnostics) THEN
                0701         CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid)
                0702         CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid)
                0703       ENDIF
                0704 #endif /* ALLOW_DIAGNOSTICS */
                0705 
0bd9d32119 Jean*0706 #endif /* ALLOW_MOM_COMMON */
88830be691 Alis*0707 #endif /* ALLOW_NONHYDROSTATIC */
                0708 
                0709       RETURN
                0710       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!