Back to home page

MITgcm
 
 

    

File indexing completed on 2021-05-05 05:11:54 UTC

view on githubraw file Latest commit 70ceddf2 on 2021-05-04 21:30:49 UTC
f59d76b0dd Ed D*0001 #include "GMREDI_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C !ROUTINE: GMREDI_CALC_QGLEITH
                0005 
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE GMREDI_CALC_QGLEITH(
                0008      O             leithQG_K,
                0009      I             bi, bj, myTime, myIter, myThid )
                0010 
                0011 C     !DESCRIPTION:
                0012 C     *==========================================================*
                0013 C     | SUBROUTINE GMREDI_CALC_QGLEITH
                0014 C     | Calculate QG Leith contribution to GMRedi tensor.
                0015 C     | leithQG_K is located at the cell centre.
                0016 C     *==========================================================*
                0017 
                0018 C !USES: ===============================================================
                0019       IMPLICIT NONE
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
                0024 #include "DYNVARS.h"
                0025 #ifdef ALLOW_MOM_COMMON
                0026 # include "MOM_VISC.h"
                0027 #endif
                0028 
                0029 C !INPUT PARAMETERS: ===================================================
                0030 C     bi, bj    :: tile indices
                0031 C     myTime    :: Current time in simulation
                0032 C     myIter    :: Current iteration number in simulation
                0033 C     myThid    :: My Thread Id. number
                0034       INTEGER bi, bj
                0035       _RL     myTime
                0036       INTEGER myIter
                0037       INTEGER myThid
                0038 
                0039 C !OUTPUT PARAMETERS: ==================================================
                0040 C  leithQG_K    :: Horizontal LeithQG viscosity, to add to GM coefficient
                0041       _RL leithQG_K(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0042 
                0043 #ifdef ALLOW_GM_LEITH_QG
                0044 #ifdef ALLOW_MOM_COMMON
                0045 
                0046 C !LOCAL VARIABLES: ====================================================
                0047       _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0048       _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049       _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0050       _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0051       _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0052       _RL Nsquare (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0053       _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0054       _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0055       _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0056       _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0058       _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0059       _RL grdVrt, grdDiv
                0060       _RL leithQG2fac
70ceddf25f Timo*0061       _RL deepFac3, L3, sqarg
f59d76b0dd Ed D*0062 C     i,j,k    :: Loop counters
                0063       INTEGER i,j,k
                0064 CEOP
                0065 
                0066 C--   Initialise terms
                0067 
70ceddf25f Timo*0068       IF (useFullLeith) THEN
                0069 C     Uses correct calculation for gradients, but might not work on cube sphere
                0070         leithQG2fac = (viscC2LeithQG/pi)**6
                0071       ELSE
                0072 C     Uses approximate gradients, but works on cube sphere. No reason to use this
                0073 C      unless `useFullLeith` fails for your setup.
                0074         leithQG2fac = (viscC2LeithQG/pi)**3
                0075       ENDIF
f59d76b0dd Ed D*0076 
                0077       DO j=1-OLy,sNy+OLy
                0078         DO i=1-OLx,sNx+OLx
                0079           hFacZ(i,j)         = 0. _d 0
                0080           r_hFacZ(i,j)       = 0. _d 0
                0081           uFld(i,j)          = 0. _d 0
                0082           vFld(i,j)          = 0. _d 0
                0083           stretching(i,j)    = 0. _d 0
                0084           Nsquare(i,j)       = 0. _d 0
                0085           vort3(i,j)         = 0. _d 0
                0086           hDiv(i,j)          = 0. _d 0
                0087           divDx(i,j)         = 0. _d 0
                0088           divDy(i,j)         = 0. _d 0
                0089           vrtDx(i,j)         = 0. _d 0
                0090           vrtDy(i,j)         = 0. _d 0
                0091         ENDDO
                0092       ENDDO
                0093 
                0094 C     start k loop since momentum code is inside one
                0095       DO k=1,Nr
                0096 
                0097         deepFac3 = deepFac2C(k)*deepFacC(k)
                0098 
                0099 C--     Calculate open water fraction at vorticity points
                0100         CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
                0101 
                0102 C       Make local copies of horizontal flow field
                0103         DO j=1-OLy,sNy+OLy
                0104          DO i=1-OLx,sNx+OLx
                0105           uFld(i,j) = uVel(i,j,k,bi,bj)
                0106           vFld(i,j) = vVel(i,j,k,bi,bj)
                0107          ENDDO
                0108         ENDDO
                0109 
                0110         CALL MOM_CALC_RELVORT3( bi,bj,k,uFld,vFld,hFacZ,vort3,myThid )
                0111         CALL MOM_CALC_HDIV( bi,bj,k,2,uFld,vFld,hDiv,myThid )
                0112 
                0113 C vorticity needs to be masked if using vector invariant momentum
                0114 #ifdef ALLOW_MOM_VECINV
                0115 C-      mask vort3
                0116         DO j=1-OLy,sNy+OLy
                0117           DO i=1-OLx,sNx+OLx
                0118             IF ( hFacZ(i,j).EQ.zeroRS ) THEN
                0119               vort3(i,j)   = 0.
                0120             ENDIF
                0121           ENDDO
                0122         ENDDO
                0123 #endif /* ALLOW_MOM_VECINV */
                0124 
                0125 C     Having calculated the quantitites, use them to calculate
                0126 C       LeithQG coefficient
                0127 
                0128         CALL MOM_VISC_QGL_STRETCH(bi,bj,k,
                0129      O                                   stretching, Nsquare,
                0130      I                                   myTime, myIter, myThid)
                0131         CALL MOM_VISC_QGL_LIMIT(bi,bj,k,
                0132      O                                 stretching,
                0133      I                                 Nsquare, uFld,vFld,vort3,
                0134      I                                 myTime, myIter, myThid)
                0135 
                0136 C--     horizontal gradient of horizontal divergence:
                0137 C-       gradient in x direction:
                0138         IF (useCubedSphereExchange) THEN
                0139 C        to compute d/dx(hDiv), fill corners with appropriate values:
                0140           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
                0141      &                                hDiv, bi,bj, myThid )
                0142         ENDIF
                0143         DO j=2-OLy,sNy+OLy-1
                0144           DO i=2-OLx,sNx+OLx-1
                0145             divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))
                0146      &                  *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
                0147           ENDDO
                0148         ENDDO
                0149 
                0150 C-      gradient in y direction:
                0151         IF (useCubedSphereExchange) THEN
                0152 C        to compute d/dy(hDiv), fill corners with appropriate values:
                0153           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
                0154      &                                hDiv, bi,bj, myThid )
                0155         ENDIF
                0156         DO j=2-OLy,sNy+OLy-1
                0157           DO i=2-OLx,sNx+OLx-1
                0158             divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))
                0159      &                  *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
                0160           ENDDO
                0161         ENDDO
                0162 
                0163 C       horizontal gradient of vorticity and vortex stretching:
                0164 C        In the case of using QG Leith, we want to add a term
                0165 C        before calculating vector magnitude.
                0166 C        gradient in x direction:
                0167         DO j=2-OLy,sNy+OLy
                0168           DO i=2-OLx,sNx+OLx-1
                0169             vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
                0170      &                  *recip_dxG(i,j,bi,bj)*recip_deepFacC(k)
                0171      &                  *maskS(i,j,k,bi,bj)
                0172 #ifdef ALLOW_OBCS
                0173      &                  *maskInS(i,j,bi,bj)
                0174 #endif
                0175 C        Average d/dx of stretching onto V-points to match vrtDX
70ceddf25f Timo*0176      &                 + halfRL*halfRL*
f59d76b0dd Ed D*0177      &                    ((stretching(i+1,j)-stretching(i,j))
                0178      &                      *recip_dxC(i+1,j,bi,bj)*recip_deepFacC(k)
                0179      &                    + (stretching(i,j)-stretching(i-1,j))
                0180      &                      *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
                0181      &                    + (stretching(i+1,j-1)-stretching(i,j-1))
                0182      &                      *recip_dxC(i,j-1,bi,bj)*recip_deepFacC(k)
                0183      &                    + (stretching(i,j-1)-stretching(i-1,j-1))
                0184      &                      *recip_dxC(i-1,j-1,bi,bj)*recip_deepFacC(k)
70ceddf25f Timo*0185      &                    )*maskS(i,j,k,bi,bj)
f59d76b0dd Ed D*0186 #ifdef ALLOW_OBCS
70ceddf25f Timo*0187      &                     *maskInS(i,j,bi,bj)
f59d76b0dd Ed D*0188 #endif
                0189           ENDDO
                0190         ENDDO
                0191 C-       gradient in y direction:
                0192         DO j=2-OLy,sNy+OLy-1
                0193           DO i=2-OLx,sNx+OLx
                0194             vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
                0195      &                  *recip_dyG(i,j,bi,bj)*recip_deepFacC(k)
                0196      &                  *maskW(i,j,k,bi,bj)
                0197 #ifdef ALLOW_OBCS
                0198      &                  *maskInW(i,j,bi,bj)
                0199 #endif
                0200 C        Average d/dy of stretching onto U-points to match vrtDy
70ceddf25f Timo*0201      &                 + halfRL*halfRL*
f59d76b0dd Ed D*0202      &                    ((stretching(i,j+1)-stretching(i,j))
                0203      &                      *recip_dyC(i,j+1,bi,bj)*recip_deepFacC(k)
                0204      &                    + (stretching(i,j)-stretching(i,j-1))
                0205      &                      *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
                0206      &                    + (stretching(i-1,j+1)-stretching(i-1,j))
                0207      &                      *recip_dyC(i-1,j+1,bi,bj)*recip_deepFacC(k)
                0208      &                    + (stretching(i-1,j)-stretching(i-1,j-1))
                0209      &                      *recip_dyC(i-1,j,bi,bj)*recip_deepFacC(k)
70ceddf25f Timo*0210      &                    )*maskW(i,j,k,bi,bj)
f59d76b0dd Ed D*0211 #ifdef ALLOW_OBCS
70ceddf25f Timo*0212      &                     *maskInW(i,j,bi,bj)
f59d76b0dd Ed D*0213 #endif
                0214           ENDDO
                0215         ENDDO
                0216 
                0217         DO j=2-OLy,sNy+OLy-1
                0218           DO i=2-OLx,sNx+OLx-1
                0219 
                0220 C These are (powers of) length scales
                0221             L3 = L3_D(i,j,bi,bj)*deepFac3
                0222 
                0223             IF (useFullLeith) THEN
                0224 C This is the vector magnitude of the vorticity gradient squared
                0225               grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
                0226      &                        + vrtDx(i,j)*vrtDx(i,j) )
                0227      &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
                0228      &                        + vrtDy(i,j)*vrtDy(i,j) )  )
                0229 
                0230 C This is the vector magnitude of grad (div.v) squared
                0231 C Using it in Leith serves to damp instabilities in w.
                0232               grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
                0233      &                        + divDx(i,j)*divDx(i,j) )
                0234      &                     + (divDy(i,j+1)*divDy(i,j+1)
                0235      &                        + divDy(i,j)*divDy(i,j) )  )
                0236 
70ceddf25f Timo*0237               sqarg = leithQG2fac*(grdVrt+grdDiv)
                0238 #ifdef ALLOW_AUTODIFF
                0239 C avoid derivative of SQRT(0)
                0240               IF (sqarg .GT. 0. _d 0) sqarg = SQRT(sqarg)
                0241 #else
                0242               sqarg = SQRT(sqarg)
                0243 #endif
                0244               LeithQG_K(i,j,k) = sqarg*L3
f59d76b0dd Ed D*0245 
                0246             ELSE
                0247 C but this approximation will work on cube (and differs by as much as 4X)
                0248               grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )
                0249               grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
                0250               grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
                0251 
                0252 C This approximation is good to the same order as above...
                0253               grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
                0254               grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
                0255               grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
                0256 
                0257               LeithQG_K(i,j,k) = leithQG2fac*(grdVrt + grdDiv)*L3
                0258 
                0259             ENDIF
                0260           ENDDO
                0261         ENDDO
                0262 
                0263       ENDDO /* k loop */
                0264 
                0265 #ifdef ALLOW_DIAGNOSTICS
                0266       IF ( useDiagnostics ) THEN
                0267         CALL DIAGNOSTICS_FILL( LeithQG_K, 'GM_LTHQG',
                0268      &                         0, Nr, 2, bi, bj, myThid )
                0269       ENDIF
                0270 #endif /* ALLOW_DIAGNOSTICS */
                0271 
                0272 #endif /* ALLOW_MOM_COMMON */
                0273 #endif /* ALLOW_GM_LEITH_QG */
                0274 
                0275       RETURN
                0276       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!