pkg/mom_common/mom_visc_qgl_limit.F

f59d76b0dd Ed D*0001 #include "MOM_COMMON_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: MOM_VISC_QGL_LIMIT
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE MOM_VISC_QGL_LIMIT(
                0008      I                           bi, bj, k,
                0009      O                           stretching,
                0010      I                           Nsquare,
                0011      I                           uFld, vFld, vort3,
                0012      I                           myTime, myIter, myThid )
                0013
                0014 C     !DESCRIPTION:
                0015 C     *==========================================================*
                0016 C     | SUBROUTINE MOM_VISC_QGL_LIMIT
                0017 C     | Limit the contribution of the vortex stretching term.
                0018 C     | When the flow is unstratified, the algorithm tries to
                0019 C     | divide by zero. This subroutine ensures that the result
                0020 C     | remains finite.
                0021 C     *==========================================================*
                0022
                0023 C !USES: ===============================================================
                0024       IMPLICIT NONE
                0025
                0026 C     == Global variables ==
                0027 #include "SIZE.h"
                0028 #include "EEPARAMS.h"
                0029 #include "PARAMS.h"
                0030 #include "GRID.h"
                0031 c#include "DYNVARS.h"
                0032 c#include "MOM_VISC.h"
                0033
                0034 C !INPUT PARAMETERS: ===================================================
                0035 C  bi,bj                :: tile indices
                0036 C  k                    :: vertical level
                0037 C  Nsquare              :: buoyancy frequency
                0038 C  uFld                 :: U velocity (east flow on spherical grid)
                0039 C  vFld                 :: V velocity (north flow on spherical grid)
                0040 C  vort3                :: Relative vorticity
                0041 C  myTime               :: current time of simulation ( s )
                0042 C  myIter               :: current iteration number of simulation
                0043 C  myThid               :: my Thread Id number
                0044       INTEGER bi,bj,k
                0045       _RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0048       _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049
                0050       _RL     myTime
                0051       INTEGER myIter
                0052       INTEGER myThid
                0053
                0054 C !OUTPUT PARAMETERS: ==================================================
                0055 C  stretching           :: vortex stretching term (to be limited)
                0056       _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057
                0058 #ifdef ALLOW_LEITH_QG
                0059
                0060 C !LOCAL VARIABLES: ====================================================
                0061       INTEGER i,j
                0062       _RL U_scale_sq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       _RL Ro_g_sq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064 C      _RL Bu_g(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0065       _RL Fr_g_sq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0066       _RL stretching_hold(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL QGL_lim_epsil
                0068       _RL vort3C
                0069
                0070 CEOP
                0071
                0072 C Initialise local variables
                0073       DO j=1-OLy,sNy+OLy
                0074        DO i=1-OLx,sNx+OLx
                0075 C       square of the velocity scale
                0076         U_scale_sq(i,j)         = 0. _d 0
                0077 C       square of the gridscale Rossby number
                0078         Ro_g_sq(i,j)           = 0. _d 0
                0079 C       Gridscale Burger number
                0080 C        Bu_g(i,j)            = 0. _d 0
                0081 C       Square of the gridscale Froude number
                0082         Fr_g_sq(i,j)            = 0. _d 0
                0083         stretching_hold(i,j) = 0. _d 0
                0084        ENDDO
                0085       ENDDO
                0086
                0087       QGL_lim_epsil = 1. _d -24
                0088
                0089 C Put a cap on the stretching term.
                0090
                0091        DO j=2-OLy,sNy+OLy-1
                0092         DO I=2-OLx,sNx+OLx-1
                0093           U_scale_sq(i,j) = 0.5* (
                0094      &             ( uFld(i,j)*uFld(i,j)
                0095      &              + uFld(i+1,j)*uFld(i+1,j) )
                0096      &           + ( vFld(i,j)*vFld(i,j)
                0097      &              + vFld(i,j+1)*vFld(i,j+1) )
                0098      &                             )
                0099
                0100 C        Grid scale Rossby number, squared: U^2 / (fL)^2
                0101           Ro_g_sq(i,j) = U_scale_sq(i,j) * recip_rA(i,j,bi,bj) /
                0102      &                  MAX(QGL_lim_epsil, fCori(i,j,bi,bj)**2)
                0103 C        Grid scale Burger number: N^2 H^2 / (f^2 L^2)
                0104 C         Bu_g(i,j) = Nsquare(i,j) * drf(k)**2 * recip_rA(i,j,bi,bj) /
                0105 C     &                 (fCori(i,j,bi,bj)**2 * PI)
                0106 C        Grid scale Froude number, squared: U^2 pi^2 /(N^2 H)^2
                0107 C          Include a small number to prevent division by zero
                0108           Fr_g_sq(i,j) = U_scale_sq(i,j) * PI * PI /
                0109      &                 (MAX((Nsquare(i,j) * drF(k))**2,
                0110      &                      QGL_lim_epsil))
                0111
                0112 C     Make the scheme gracefully transition to something else
                0113 C     as stratification goes to zero
                0114 C
                0115 C        Implement eqn. (55) from Bachman et al. (2017) JGR-Oceans
                0116 C         stretching_hold(i,j) = MIN( ABS(stretching(i,j)),
                0117 C     &    ABS(vort3(i,j) / MAX(Bu_g(i,j),Ro_g_sq(i,j)) ) )
                0118
                0119 C        Implement eqn. (56) from Bachman et al. (2017) JGR-Oceans
                0120 C          This limiter goes to 2D Leith as stratification -> 0
                0121           vort3C = halfRL*halfRL*(vort3(i,j) + vort3(i+1,j)
                0122      &                            + vort3(i,j+1) + vort3(i+1,j+1))
                0123
                0124           stretching_hold(i,j) = min( ABS(stretching(i,j)),
                0125      &                  ABS(vort3C * Fr_g_sq(i,j) /
                0126      &                       (Ro_g_sq(i,j) + Fr_g_sq(i,j)**2
                0127      &                         + QGL_lim_epsil) )
                0128      &                               )
                0129
                0130           stretching(i,j) = SIGN(stretching_hold(i,j), stretching(i,j))
                0131
                0132         ENDDO
                0133        ENDDO
                0134
                0135 #endif /* ALLOW_LEITH_QG */
                0136
                0137       RETURN
                0138       END