pkg/mom_common/mom_w_sidedrag.F

7d6d4ca1be Jean*0001 #include "MOM_COMMON_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: MOM_W_SIDEDRAG
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE MOM_W_SIDEDRAG(
                0008      I        bi,bj,k,
                0009      I        wFld, del2w,
                0010      I        rThickC, recip_rThickC,
                0011      I        viscAh_W, viscA4_W,
                0012      O        gwSideDrag,
                0013      I        myThid)
                0014
                0015 C !DESCRIPTION:
                0016 C Calculates the drag terms due to the no-slip condition on viscous stresses:
                0017 C \begin{equation*}
                0018 C G^w_{drag} = - \frac{2}{\Delta x_w} (A_h w - A_4 \nabla^2 w)
                0019 C \end{equation*}
                0020
                0021 C !USES: ===============================================================
                0022       IMPLICIT NONE
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "SURFACE.h"
                0028
                0029 C !INPUT PARAMETERS: ===================================================
                0030 C  bi,bj                :: tile indices
                0031 C  k                    :: vertical level
                0032 C  wFld                 :: vertical velocity
                0033 C  del2w                :: Laplacian of vertical velocity
                0034 C  rThickC              :: thickness of W-cell
                0035 C  recip_rThickC        :: reciprol of W-cell thickness
                0036 C  viscAh_W             :: harmonic   horizontal viscosity (at W-Cell center)
                0037 C  viscA4_W             :: biharmonic horizontal viscosity (at W-Cell center)
                0038 C  myThid               :: my Thread Id number
                0039       INTEGER bi,bj,k
                0040       _RL  wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0041       _RL  del2w   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL  rThickC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0043       _RL recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL  viscAh_W(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0045       _RL  viscA4_W(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0046       INTEGER myThid
                0047
                0048 C !OUTPUT PARAMETERS: ==================================================
                0049 C  gwSideDrag           :: side drag tendency
                0050       _RL gwSideDrag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0051
                0052 #ifdef ALLOW_NONHYDROSTATIC
                0053 C !LOCAL VARIABLES: ====================================================
                0054 C  i,j                  :: loop indices
                0055 C  heightStepW          :: height of topographic step to West
                0056 C  heightStepE          :: height of topographic step to East
                0057 C  heightStepS          :: height of topographic step to South
                0058 C  heightStepN          :: height of topographic step to North
                0059       INTEGER i,j
                0060       _RL heightStepW, heightStepE
                0061       _RL heightStepS, heightStepN
                0062 CEOP
                0063
                0064 C--   Laplacian and bi-harmonic terms: using variable-Viscosity coeff.
                0065 C     from MOM_CALC_VISC, consistent with dissipation in the interior
                0066
                0067       DO j=2-Oly,sNy+Oly-1
                0068        DO i=2-Olx,sNx+Olx-1
                0069 C-      this will not give any side-drag along thin wall.
                0070 C       (but this might just be what we want ...)
                0071         heightStepW = MAX( 0. _d 0, rThickC(i,j) - rThickC(i-1,j) )
                0072         heightStepE = MAX( 0. _d 0, rThickC(i,j) - rThickC(i+1,j) )
                0073         heightStepS = MAX( 0. _d 0, rThickC(i,j) - rThickC(i,j-1) )
                0074         heightStepN = MAX( 0. _d 0, rThickC(i,j) - rThickC(i,j+1) )
                0075 C- jmc: take directly viscAh_W & viscA4_W (no bar_i, bar_j) since
                0076 C       they are already at the center (and not above u,v point).
                0077         gwSideDrag(i,j) =
                0078      &   -sideDragFactor
                0079      &    *( heightStepW*_dyG( i ,j,bi,bj)*_recip_dxC( i ,j,bi,bj)
                0080      &         *( viscAh_W(i,j,k,bi,bj)*wFld(i,j,k,bi,bj)
                0081      &                                 *cosFacV(j,bi,bj)
                0082      &           -viscA4_W(i,j,k,bi,bj)*del2w(i,j)
                0083 #ifdef COSINEMETH_III
                0084      &                                 *sqCosFacU(j,bi,bj)
                0085 #else
                0086      &                                 *cosFacU(j,bi,bj)
                0087 #endif
                0088      &          )
                0089      &     +heightStepE*_dyG(i+1,j,bi,bj)*_recip_dxC(i+1,j,bi,bj)
                0090      &         *( viscAh_W(i,j,k,bi,bj)*wFld(i,j,k,bi,bj)
                0091      &                                 *cosFacV(j,bi,bj)
                0092      &           -viscA4_W(i,j,k,bi,bj)*del2w(i,j)
                0093 #ifdef COSINEMETH_III
                0094      &                                 *sqCosFacU(j,bi,bj)
                0095 #else
                0096      &                                 *cosFacU(j,bi,bj)
                0097 #endif
                0098      &          )
                0099      &     +heightStepS*_dxG(i,j,bi,bj)*_recip_dyC(i,j,bi,bj)
                0100      &         *( viscAh_W(i,j,k,bi,bj)*wFld(i,j,k,bi,bj)
                0101 #ifdef ISOTROPIC_COS_SCALING
                0102      &                                 *cosFacV(j,bi,bj)
                0103 #endif
                0104      &           -viscA4_W(i,j,k,bi,bj)*del2w(i,j)
                0105 #ifdef ISOTROPIC_COS_SCALING
                0106 # ifdef COSINEMETH_III
                0107      &                                 *sqCosFacV(j,bi,bj)
                0108 # else
                0109      &                                 *cosFacV(j,bi,bj)
                0110 # endif
                0111 #endif
                0112      &          )
                0113      &     +heightStepN*_dxG(i,j+1,bi,bj)*_recip_dyC(i,j+1,bi,bj)
                0114      &         *( viscAh_W(i,j,k,bi,bj)*wFld(i,j,k,bi,bj)
                0115 #ifdef ISOTROPIC_COS_SCALING
                0116      &                                 *cosFacV(j+1,bi,bj)
                0117 #endif
                0118      &           -viscA4_W(i,j,k,bi,bj)*del2w(i,j)
                0119 #ifdef ISOTROPIC_COS_SCALING
                0120 # ifdef COSINEMETH_III
                0121      &                                 *sqCosFacV(j+1,bi,bj)
                0122 # else
                0123      &                                 *cosFacV(j+1,bi,bj)
                0124 # endif
                0125 #endif
                0126      &          )
                0127      &    )*recip_rThickC(i,j)*recip_rA(i,j,bi,bj)
                0128        ENDDO
                0129       ENDDO
                0130
                0131 #ifdef ALLOW_DIAGNOSTICS
                0132       IF (useDiagnostics) THEN
                0133         CALL DIAGNOSTICS_FILL(gwSideDrag,'WSidDrag',k,1,2,bi,bj,myThid)
                0134       ENDIF
                0135 #endif /* ALLOW_DIAGNOSTICS */
                0136
                0137 #endif /* ALLOW_NONHYDROSTATIC */
                0138
                0139       RETURN
                0140       END