Back to home page

MITgcm
 
 

    

File indexing completed on 2024-02-01 06:10:04 UTC

view on githubraw file Latest commit 427e24e1 on 2024-01-31 16:50: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
b88f9ec5cb Jean*0125       LOGICAL diagDiss, diagAdvec, 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 )
67662de2f2 Jean*0136         diag_AB   = DIAGNOSTICS_IS_ON( 'AB_gW   ', myThid )
41d43b2d0d Jean*0137       ELSE
                0138         diagDiss  = .FALSE.
                0139         diagAdvec = .FALSE.
b88f9ec5cb Jean*0140         diag_AB   = .FALSE.
41d43b2d0d Jean*0141       ENDIF
                0142 #endif /* ALLOW_DIAGNOSTICS */
                0143 
0bd9d32119 Jean*0144       biharmonicVisc = viscA4W.NE.zeroRL
                0145      &           .OR. ( useVariableVisc .AND. useBiharmonicVisc )
                0146 
d455cbf76e Jean*0147 C--   Initialise gW to zero
                0148       DO k=1,Nr
                0149         DO j=1-OLy,sNy+OLy
                0150          DO i=1-OLx,sNx+OLx
88830be691 Alis*0151            gW(i,j,k,bi,bj) = 0.
                0152          ENDDO
                0153         ENDDO
d455cbf76e Jean*0154       ENDDO
597720f5c5 Jean*0155 C-    Initialise gwDiss to zero
                0156       DO j=1-OLy,sNy+OLy
                0157        DO i=1-OLx,sNx+OLx
                0158          gwDiss(i,j) = 0.
                0159        ENDDO
                0160       ENDDO
50603588ee Jean*0161       IF (momViscosity) THEN
                0162 C-    Initialize del2w to zero:
935bd32a21 Jean*0163         DO j=1-OLy,sNy+OLy
                0164          DO i=1-OLx,sNx+OLx
50603588ee Jean*0165            del2w(i,j) = 0. _d 0
                0166          ENDDO
                0167         ENDDO
                0168       ENDIF
88830be691 Alis*0169 
9091e25e07 Jean*0170 C--   Boundaries condition at top (vertical advection of vertical momentum):
982e105a17 Jean*0171       DO j=1-OLy,sNy+OLy
                0172        DO i=1-OLx,sNx+OLx
                0173          flxAdvUp(i,j) = 0.
                0174 c        flxDisUp(i,j) = 0.
                0175        ENDDO
                0176       ENDDO
                0177 
d455cbf76e Jean*0178 C---  Sweep down column
982e105a17 Jean*0179       DO k=1,Nr
                0180         km1 = MAX( k-1, 1 )
                0181         kp1 = MIN( k+1,Nr )
                0182         mskM1 = 1.
                0183         mskP1 = 1.
                0184         IF ( k.EQ. 1 ) mskM1 = 0.
                0185         IF ( k.EQ.Nr ) mskP1 = 0.
                0186         IF ( k.GT.1 ) THEN
597720f5c5 Jean*0187 C--   Compute grid factor arround a W-point:
53169f2ab0 Jean*0188 #ifdef CALC_GW_NEW_THICK
935bd32a21 Jean*0189          DO j=1-OLy,sNy+OLy
                0190           DO i=1-OLx,sNx+OLx
53169f2ab0 Jean*0191            IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
                0192      &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
                0193              recip_rThickC(i,j) = 0.
                0194            ELSE
                0195 C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
                0196              recip_rThickC(i,j) = 1. _d 0 /
                0197      &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
                0198      &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
                0199      &        )
                0200            ENDIF
982e105a17 Jean*0201           ENDDO
53169f2ab0 Jean*0202          ENDDO
982e105a17 Jean*0203          IF (momViscosity) THEN
935bd32a21 Jean*0204          DO j=1-OLy,sNy+OLy
                0205           DO i=1-OLx,sNx+OLx
53169f2ab0 Jean*0206            rThickC_C(i,j) = MAX( zeroRS,
                0207      &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
                0208      &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
                0209      &                         )
                0210           ENDDO
                0211          ENDDO
935bd32a21 Jean*0212          DO j=1-OLy,sNy+OLy
                0213           DO i=1-OLx+1,sNx+OLx
53169f2ab0 Jean*0214            rThickC_W(i,j) = MAX( zeroRS,
50603588ee Jean*0215      &                           MIN( rSurfW(i,j,bi,bj), rC(k-1) )
                0216      &                          -MAX(  rLowW(i,j,bi,bj), rC(k)   )
53169f2ab0 Jean*0217      &                         )
                0218 C     W-Cell Western face area:
                0219            xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
                0220 c    &                              *deepFacF(k)
                0221           ENDDO
                0222          ENDDO
935bd32a21 Jean*0223          DO j=1-OLy+1,sNy+OLy
                0224           DO i=1-OLx,sNx+OLx
53169f2ab0 Jean*0225            rThickC_S(i,j) = MAX( zeroRS,
50603588ee Jean*0226      &                           MIN( rSurfS(i,j,bi,bj), rC(k-1) )
                0227      &                          -MAX(  rLowS(i,j,bi,bj), rC(k)   )
53169f2ab0 Jean*0228      &                         )
                0229 C     W-Cell Southern face area:
                0230            yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
                0231 c    &                              *deepFacF(k)
                0232 C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
                0233 C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
                0234           ENDDO
                0235          ENDDO
982e105a17 Jean*0236          ENDIF
53169f2ab0 Jean*0237 #else /* CALC_GW_NEW_THICK */
935bd32a21 Jean*0238          DO j=1-OLy,sNy+OLy
                0239           DO i=1-OLx,sNx+OLx
597720f5c5 Jean*0240 C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
                0241            IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
                0242              recip_rThickC(i,j) = 0.
                0243            ELSE
                0244              recip_rThickC(i,j) = 1. _d 0 /
85162e9502 Jean*0245      &        ( drF(k-1)*halfRS
597720f5c5 Jean*0246      &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
                0247      &        )
                0248            ENDIF
1487584115 Jean*0249 c          IF (momViscosity) THEN
                0250 #ifdef NONLIN_FRSURF
4606c28752 Jean*0251            rThickC_C(i,j) =
1487584115 Jean*0252      &          drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0253      &        + drF( k )*MIN( h0FacC(i,j,k  ,bi,bj), halfRS )
                0254 #else
4606c28752 Jean*0255            rThickC_C(i,j) =
1487584115 Jean*0256      &          drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0257      &        + drF( k )*MIN( _hFacC(i,j,k  ,bi,bj), halfRS )
                0258 #endif
                0259            rThickC_W(i,j) =
                0260      &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0261      &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
                0262            rThickC_S(i,j) =
                0263      &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
                0264      &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
597720f5c5 Jean*0265 C     W-Cell Western face area:
                0266            xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
4606c28752 Jean*0267 c    &                              *deepFacF(k)
597720f5c5 Jean*0268 C     W-Cell Southern face area:
                0269            yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
4606c28752 Jean*0270 c    &                              *deepFacF(k)
                0271 C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
                0272 C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
1487584115 Jean*0273 c          ENDIF
982e105a17 Jean*0274           ENDDO
597720f5c5 Jean*0275          ENDDO
53169f2ab0 Jean*0276 #endif /* CALC_GW_NEW_THICK */
982e105a17 Jean*0277         ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
935bd32a21 Jean*0278          DO j=1-OLy,sNy+OLy
                0279           DO i=1-OLx,sNx+OLx
982e105a17 Jean*0280            recip_rThickC(i,j) = recip_drC(k)
                0281 c          rThickC_C(i,j) = drC(k)
                0282 c          rThickC_W(i,j) = drC(k)
                0283 c          rThickC_S(i,j) = drC(k)
                0284 c          xA(i,j) = _dyG(i,j,bi,bj)*drC(k)
                0285 c          yA(i,j) = _dxG(i,j,bi,bj)*drC(k)
                0286           ENDDO
                0287          ENDDO
                0288         ENDIF
                0289 
                0290 C--   local copy of wVel:
935bd32a21 Jean*0291         DO j=1-OLy,sNy+OLy
                0292           DO i=1-OLx,sNx+OLx
982e105a17 Jean*0293             wFld(i,j) = wVel(i,j,k,bi,bj)
                0294           ENDDO
                0295         ENDDO
597720f5c5 Jean*0296 
d455cbf76e Jean*0297 C--   horizontal bi-harmonic dissipation
0bd9d32119 Jean*0298         IF ( momViscosity .AND. k.GT.1 .AND. biharmonicVisc ) THEN
50603588ee Jean*0299 
d455cbf76e Jean*0300 C-    calculate the horizontal Laplacian of vertical flow
0b6cbae535 Mart*0301 C     Zonal flux d/dx W
50603588ee Jean*0302           IF ( useCubedSphereExchange ) THEN
                0303 C     to compute d/dx(W), fill corners with appropriate values:
                0304             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
                0305      &                                 wFld, bi,bj, myThid )
                0306           ENDIF
935bd32a21 Jean*0307           DO j=1-OLy,sNy+OLy
                0308            flx_EW(1-OLx,j)=0.
                0309            DO i=1-OLx+1,sNx+OLx
597720f5c5 Jean*0310             flx_EW(i,j) =
50603588ee Jean*0311      &               ( wFld(i,j) - wFld(i-1,j) )
597720f5c5 Jean*0312      &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
0b6cbae535 Mart*0313 #ifdef COSINEMETH_III
4606c28752 Jean*0314      &              *sqCosFacU(j,bi,bj)
0b6cbae535 Mart*0315 #endif
b09ee35a04 Jean*0316 #ifdef ALLOW_OBCS
                0317      &              *maskInW(i,j,bi,bj)
                0318 #endif
0b6cbae535 Mart*0319            ENDDO
d455cbf76e Jean*0320           ENDDO
50603588ee Jean*0321 
0b6cbae535 Mart*0322 C     Meridional flux d/dy W
50603588ee Jean*0323           IF ( useCubedSphereExchange ) THEN
                0324 C     to compute d/dy(W), fill corners with appropriate values:
                0325             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
                0326      &                                 wFld, bi,bj, myThid )
                0327           ENDIF
935bd32a21 Jean*0328           DO i=1-OLx,sNx+OLx
                0329            flx_NS(i,1-OLy)=0.
0b6cbae535 Mart*0330           ENDDO
935bd32a21 Jean*0331           DO j=1-OLy+1,sNy+OLy
                0332            DO i=1-OLx,sNx+OLx
597720f5c5 Jean*0333             flx_NS(i,j) =
50603588ee Jean*0334      &               ( wFld(i,j) - wFld(i,j-1) )
597720f5c5 Jean*0335      &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
0b6cbae535 Mart*0336 #ifdef ISOTROPIC_COS_SCALING
                0337 #ifdef COSINEMETH_III
597720f5c5 Jean*0338      &              *sqCosFacV(j,bi,bj)
0b6cbae535 Mart*0339 #endif
                0340 #endif
b09ee35a04 Jean*0341 #ifdef ALLOW_OBCS
                0342      &              *maskInS(i,j,bi,bj)
                0343 #endif
0b6cbae535 Mart*0344            ENDDO
                0345           ENDDO
52c561f327 Jean*0346 
0b6cbae535 Mart*0347 C     del^2 W
50603588ee Jean*0348 C     Divergence of horizontal fluxes
935bd32a21 Jean*0349           DO j=1-OLy,sNy+OLy-1
                0350            DO i=1-OLx,sNx+OLx-1
50603588ee Jean*0351             del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
                0352      &                    +( flx_NS(i,j+1)-flx_NS(i,j) )
597720f5c5 Jean*0353      &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
4606c28752 Jean*0354      &                    *recip_deepFac2F(k)
0b6cbae535 Mart*0355            ENDDO
                0356           ENDDO
50603588ee Jean*0357 C end if biharmonic viscosity
d455cbf76e Jean*0358         ENDIF
0b6cbae535 Mart*0359 
982e105a17 Jean*0360         IF ( momViscosity .AND. k.GT.1 ) THEN
597720f5c5 Jean*0361 C Viscous Flux on Western face
                0362           DO j=jMin,jMax
                0363            DO i=iMin,iMax+1
                0364              flx_EW(i,j)=
                0365      &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
                0366      &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
41a8c7589b Jean*0367      &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
1487584115 Jean*0368      &              *cosFacU(j,bi,bj)
597720f5c5 Jean*0369      &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
                0370      &              *(del2w(i,j)-del2w(i-1,j))
41a8c7589b Jean*0371      &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
0b6cbae535 Mart*0372 #ifdef COSINEMETH_III
597720f5c5 Jean*0373      &              *sqCosFacU(j,bi,bj)
0b6cbae535 Mart*0374 #else
1487584115 Jean*0375      &              *cosFacU(j,bi,bj)
0b6cbae535 Mart*0376 #endif
597720f5c5 Jean*0377            ENDDO
                0378           ENDDO
                0379 C Viscous Flux on Southern face
                0380           DO j=jMin,jMax+1
                0381            DO i=iMin,iMax
                0382              flx_NS(i,j)=
                0383      &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
                0384      &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
41a8c7589b Jean*0385      &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
1487584115 Jean*0386 #ifdef ISOTROPIC_COS_SCALING
                0387      &              *cosFacV(j,bi,bj)
                0388 #endif
597720f5c5 Jean*0389      &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
                0390      &              *(del2w(i,j)-del2w(i,j-1))
41a8c7589b Jean*0391      &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
597720f5c5 Jean*0392 #ifdef ISOTROPIC_COS_SCALING
0b6cbae535 Mart*0393 #ifdef COSINEMETH_III
597720f5c5 Jean*0394      &              *sqCosFacV(j,bi,bj)
52c561f327 Jean*0395 #else
1487584115 Jean*0396      &              *cosFacV(j,bi,bj)
0b6cbae535 Mart*0397 #endif
597720f5c5 Jean*0398 #endif
                0399            ENDDO
                0400           ENDDO
                0401 C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
                0402           DO j=jMin,jMax
                0403            DO i=iMin,iMax
                0404 C     Interpolate vert viscosity to center of tracer-cell (level k):
efd2c352d2 Jean*0405              viscLoc = ( kappaRU(i,j,k)  +kappaRU(i+1,j,k)
                0406      &                  +kappaRU(i,j,k+1)+kappaRU(i+1,j,k+1)
                0407      &                  +kappaRV(i,j,k)  +kappaRV(i,j+1,k)
                0408      &                  +kappaRV(i,j,k+1)+kappaRV(i,j+1,k+1)
597720f5c5 Jean*0409      &                 )*0.125 _d 0
                0410              flx_Dn(i,j) =
982e105a17 Jean*0411      &          - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
597720f5c5 Jean*0412      &                     -wVel(i,j, k ,bi,bj) )*rkSign
41a8c7589b Jean*0413      &                   *recip_drF(k)*rA(i,j,bi,bj)
4606c28752 Jean*0414      &                   *deepFac2C(k)*rhoFacC(k)
597720f5c5 Jean*0415            ENDDO
                0416           ENDDO
9091e25e07 Jean*0417           IF ( k.EQ.2 ) THEN
                0418 C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)
                0419            DO j=jMin,jMax
                0420             DO i=iMin,iMax
                0421 C     Interpolate horizontally (but not vertically) vert viscosity to center:
                0422 C     Although background visc. might be defined at k=1, this is not
                0423 C     generally true when using variable visc. (from vertical mixing scheme).
                0424 C     Therefore, no vert. interp. and only horizontal interpolation.
efd2c352d2 Jean*0425              viscLoc = ( kappaRU(i,j,k) + kappaRU(i+1,j,k)
                0426      &                  +kappaRV(i,j,k) + kappaRV(i,j+1,k)
9091e25e07 Jean*0427      &                 )*0.25 _d 0
                0428              flxDisUp(i,j) =
                0429      &          - viscLoc*( wVel(i,j, k ,bi,bj)
                0430      &                     -wVel(i,j,k-1,bi,bj) )*rkSign
                0431      &                   *recip_drF(k-1)*rA(i,j,bi,bj)
                0432      &                   *deepFac2C(k-1)*rhoFacC(k-1)
                0433 C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero)
                0434 c            flxDisUp(i,j) = 0.
                0435             ENDDO
                0436            ENDDO
                0437           ENDIF
597720f5c5 Jean*0438 C     Tendency is minus divergence of viscous fluxes:
4606c28752 Jean*0439 C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
597720f5c5 Jean*0440           DO j=jMin,jMax
                0441            DO i=iMin,iMax
                0442              gwDiss(i,j) =
41a8c7589b Jean*0443      &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
                0444      &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
                0445      &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
4606c28752 Jean*0446      &                                          *recip_rhoFacF(k)
41a8c7589b Jean*0447      &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
4606c28752 Jean*0448      &          *recip_deepFac2F(k)
d455cbf76e Jean*0449 C--        prepare for next level (k+1)
597720f5c5 Jean*0450              flxDisUp(i,j)=flx_Dn(i,j)
                0451            ENDDO
                0452           ENDDO
                0453         ENDIF
                0454 
982e105a17 Jean*0455         IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
597720f5c5 Jean*0456 C-     No-slip BCs impose a drag at walls...
1487584115 Jean*0457           CALL MOM_W_SIDEDRAG(
                0458      I               bi,bj,k,
                0459      I               wVel, del2w,
                0460      I               rThickC_C, recip_rThickC,
                0461      I               viscAh_W, viscA4_W,
                0462      O               gwAdd,
                0463      I               myThid )
                0464           DO j=jMin,jMax
                0465            DO i=iMin,iMax
                0466             gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
                0467            ENDDO
                0468           ENDDO
597720f5c5 Jean*0469         ENDIF
                0470 
efc81681f0 Jean*0471 #ifdef ALLOW_SMAG_3D
                0472         IF ( useSmag3D .AND. k.GT.1 ) THEN
                0473              CALL MOM_W_SMAG_3D(
                0474      I         str13, str23, str33,
                0475      I         viscAh3d_00, viscAh3d_13, viscAh3d_23,
                0476      I         rThickC_W, rThickC_S, rThickC_C, recip_rThickC,
                0477      O         gwAdd,
                0478      I         iMin,iMax,jMin,jMax, k, bi, bj, myThid )
                0479           DO j = jMin,jMax
                0480            DO i = iMin,iMax
                0481             gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
                0482            ENDDO
                0483           ENDDO
                0484         ENDIF
                0485 #endif /* ALLOW_SMAG_3D */
                0486 
597720f5c5 Jean*0487 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0488 
                0489         IF ( momAdvection ) THEN
982e105a17 Jean*0490 
                0491          IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
597720f5c5 Jean*0492 C Advective Flux on Western face
                0493           DO j=jMin,jMax
                0494            DO i=iMin,iMax+1
                0495 C     transport through Western face area:
                0496              uTrans = (
982e105a17 Jean*0497      &          drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
                0498      &                  *rhoFacC(km1)*mskM1
597720f5c5 Jean*0499      &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
4606c28752 Jean*0500      &                  *rhoFacC(k)
                0501      &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
982e105a17 Jean*0502              flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
                0503 c            flx_EW(i,j)=
                0504 c    &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
597720f5c5 Jean*0505            ENDDO
                0506           ENDDO
                0507 C Advective Flux on Southern face
                0508           DO j=jMin,jMax+1
                0509            DO i=iMin,iMax
                0510 C     transport through Southern face area:
                0511              vTrans = (
982e105a17 Jean*0512      &          drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
                0513      &                  *rhoFacC(km1)*mskM1
597720f5c5 Jean*0514      &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
4606c28752 Jean*0515      &                  *rhoFacC(k)
                0516      &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
982e105a17 Jean*0517              flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
                0518 c            flx_NS(i,j)=
                0519 c    &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
597720f5c5 Jean*0520            ENDDO
                0521           ENDDO
982e105a17 Jean*0522          ENDIF
597720f5c5 Jean*0523 C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
982e105a17 Jean*0524 c        IF (.TRUE.) THEN
597720f5c5 Jean*0525           DO j=jMin,jMax
                0526            DO i=iMin,iMax
53169f2ab0 Jean*0527 C     NH in p-coord.: advect wSpeed [m/s] with rTrans
                0528              tmp_WbarZ = halfRL*
982e105a17 Jean*0529      &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k )
                0530      &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 )
597720f5c5 Jean*0531 C     transport through Lower face area:
4606c28752 Jean*0532              rTrans = halfRL*
                0533      &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
                0534      &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
982e105a17 Jean*0535      &                                   *mskP1
4606c28752 Jean*0536      &              )*rA(i,j,bi,bj)
41a8c7589b Jean*0537              flx_Dn(i,j) = rTrans*tmp_WbarZ
597720f5c5 Jean*0538            ENDDO
                0539           ENDDO
982e105a17 Jean*0540 c        ENDIF
                0541          IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
                0542 C Advective Flux on Upper face of W-Cell (= at surface)
9091e25e07 Jean*0543            DO j=jMin,jMax
                0544             DO i=iMin,iMax
982e105a17 Jean*0545              tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
                0546              rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
                0547      &               *rA(i,j,bi,bj)
9091e25e07 Jean*0548              flxAdvUp(i,j) = rTrans*tmp_WbarZ
                0549 c            flxAdvUp(i,j) = 0.
                0550             ENDDO
                0551            ENDDO
982e105a17 Jean*0552          ENDIF
935bd32a21 Jean*0553 
982e105a17 Jean*0554          IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
597720f5c5 Jean*0555 C     Tendency is minus divergence of advective fluxes:
4606c28752 Jean*0556 C     anelastic: all transports & advect. fluxes are scaled by rhoFac
597720f5c5 Jean*0557           DO j=jMin,jMax
                0558            DO i=iMin,iMax
89a3c0c31d Jean*0559 C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)
                0560 c            IF (k.EQ.2) flxAdvUp(i,j) = 0.
597720f5c5 Jean*0561              gW(i,j,k,bi,bj) =
41a8c7589b Jean*0562      &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
                0563      &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
53169f2ab0 Jean*0564      &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
41a8c7589b Jean*0565      &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
4606c28752 Jean*0566      &          *recip_deepFac2F(k)*recip_rhoFacF(k)
982e105a17 Jean*0567            ENDDO
                0568           ENDDO
935bd32a21 Jean*0569 #ifdef ALLOW_ADDFLUID
                0570           IF ( selectAddFluid.GE.1 ) THEN
                0571            DO j=jMin,jMax
                0572             DO i=iMin,iMax
                0573              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
                0574      &        + wVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
                0575      &          *( addMass(i,j,k,bi,bj)
                0576      &            +addMass(i,j,km1,bi,bj)*mskM1 )
                0577      &          *recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
                0578      &          *recip_deepFac2F(k)*recip_rhoFacF(k)
                0579             ENDDO
                0580            ENDDO
                0581           ENDIF
                0582 #endif /* ALLOW_ADDFLUID */
982e105a17 Jean*0583          ENDIF
                0584 
                0585          DO j=jMin,jMax
                0586            DO i=iMin,iMax
597720f5c5 Jean*0587 C--          prepare for next level (k+1)
                0588              flxAdvUp(i,j)=flx_Dn(i,j)
                0589            ENDDO
982e105a17 Jean*0590          ENDDO
                0591 
                0592 c       ELSE
                0593 C-    if momAdvection / else
                0594 c         DO j=1-OLy,sNy+OLy
                0595 c          DO i=1-OLx,sNx+OLx
                0596 c            gW(i,j,k,bi,bj) = 0. _d 0
                0597 c          ENDDO
                0598 c         ENDDO
                0599 
                0600 C-    endif momAdvection.
597720f5c5 Jean*0601         ENDIF
                0602 
982e105a17 Jean*0603 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0604 
                0605         IF ( useNHMTerms .AND. k.GT.1 ) THEN
597720f5c5 Jean*0606           CALL MOM_W_METRIC_NH(
                0607      I               bi,bj,k,
                0608      I               uVel, vVel,
                0609      O               gwAdd,
                0610      I               myThid )
                0611           DO j=jMin,jMax
                0612            DO i=iMin,iMax
                0613              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
                0614            ENDDO
                0615           ENDDO
                0616         ENDIF
427e24e121 Jean*0617         IF ( select3dCoriScheme.GE.1 .AND. k.GT.1 ) THEN
597720f5c5 Jean*0618           CALL MOM_W_CORIOLIS_NH(
                0619      I               bi,bj,k,
                0620      I               uVel, vVel,
                0621      O               gwAdd,
                0622      I               myThid )
                0623           DO j=jMin,jMax
                0624            DO i=iMin,iMax
                0625              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
                0626            ENDDO
                0627           ENDDO
                0628         ENDIF
88830be691 Alis*0629 
52c561f327 Jean*0630 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ac07cf9a5a Jean*0631 
41d43b2d0d Jean*0632 #ifdef ALLOW_DIAGNOSTICS
50603588ee Jean*0633         IF ( diagDiss )  THEN
                0634           CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
                0635      &                           k, 1, 2, bi,bj, myThid )
                0636 C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
                0637 C        does it only if k=1 (never the case here)
982e105a17 Jean*0638 c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
50603588ee Jean*0639         ENDIF
                0640         IF ( diagAdvec ) THEN
                0641           CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
                0642      &                           k,Nr, 1, bi,bj, myThid )
982e105a17 Jean*0643 c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
50603588ee Jean*0644         ENDIF
41d43b2d0d Jean*0645 #endif /* ALLOW_DIAGNOSTICS */
                0646 
597720f5c5 Jean*0647 C--   Dissipation term inside the Adams-Bashforth:
                0648         IF ( momViscosity .AND. momDissip_In_AB) THEN
                0649           DO j=jMin,jMax
                0650            DO i=iMin,iMax
                0651              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
                0652            ENDDO
                0653           ENDDO
                0654         ENDIF
                0655 
52c561f327 Jean*0656 C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
                0657 C     and save gW_[n] into gwNm1 for the next time step.
cba4501825 Jean*0658 #ifdef ALLOW_ADAMSBASHFORTH_3
                0659         CALL ADAMS_BASHFORTH3(
94e3f14b29 Jean*0660      I                         bi, bj, k, Nr,
0b6eaf7dba Jean*0661      U                         gW(1-OLx,1-OLy,1,bi,bj), gwNm,
                0662      O                         gw_AB,
cba4501825 Jean*0663      I                         nHydStartAB, myIter, myThid )
                0664 #else /* ALLOW_ADAMSBASHFORTH_3 */
d455cbf76e Jean*0665         CALL ADAMS_BASHFORTH2(
94e3f14b29 Jean*0666      I                         bi, bj, k, Nr,
0b6eaf7dba Jean*0667      U                         gW(1-OLx,1-OLy,1,bi,bj),
                0668      U                         gwNm1(1-OLx,1-OLy,1,bi,bj),
                0669      O                         gw_AB,
117ee419f5 Jean*0670      I                         nHydStartAB, myIter, myThid )
cba4501825 Jean*0671 #endif /* ALLOW_ADAMSBASHFORTH_3 */
b88f9ec5cb Jean*0672 #ifdef ALLOW_DIAGNOSTICS
                0673         IF ( diag_AB ) THEN
67662de2f2 Jean*0674           CALL DIAGNOSTICS_FILL(gw_AB,'AB_gW   ',k,1,2,bi,bj,myThid)
b88f9ec5cb Jean*0675         ENDIF
                0676 #endif /* ALLOW_DIAGNOSTICS */
88830be691 Alis*0677 
597720f5c5 Jean*0678 C--   Dissipation term outside the Adams-Bashforth:
                0679         IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
                0680           DO j=jMin,jMax
                0681            DO i=iMin,iMax
                0682              gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
                0683            ENDDO
                0684           ENDDO
                0685         ENDIF
                0686 
52c561f327 Jean*0687 C-    end of the k loop
                0688       ENDDO
ffb313c34a Alis*0689 
a25a2def2e Jean*0690 #ifdef ALLOW_DIAGNOSTICS
                0691       IF (useDiagnostics) THEN
                0692         CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid)
                0693         CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid)
                0694       ENDIF
                0695 #endif /* ALLOW_DIAGNOSTICS */
                0696 
0bd9d32119 Jean*0697 #endif /* ALLOW_MOM_COMMON */
88830be691 Alis*0698 #endif /* ALLOW_NONHYDROSTATIC */
                0699 
                0700       RETURN
                0701       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!