pkg/mom_vecinv/mom_vi_u_coriolis_c4.F

cec2469d72 Alis*0001 #include "MOM_VECINV_OPTIONS.h"
aea29c8517 Alis*0002
2be0e061d5 Jean*0003 CBOP
                0004 C     !ROUTINE: MOM_VI_U_CORIOLIS_C4
                0005 C     !INTERFACE:
                0006       SUBROUTINE MOM_VI_U_CORIOLIS_C4(
                0007      I        bi,bj,k,
3370e71df9 Mart*0008      I        selectVortScheme,highOrderVorticity,upwindVorticity,
aea29c8517 Alis*0009      I        vFld,omega3,r_hFacZ,
                0010      O        uCoriolisTerm,
                0011      I        myThid)
2be0e061d5 Jean*0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | S/R MOM_VI_U_CORIOLIS_C4
                0015 C     |==========================================================*
7fe6343684 Jean*0016 C     | o Calculate flux (in Y-dir.) of vorticity at U point
                0017 C     |   using 4th order (or 1rst order) interpolation
2be0e061d5 Jean*0018 C     *==========================================================*
                0019 C     \ev
                0020
                0021 C     !USES:
aea29c8517 Alis*0022       IMPLICIT NONE
                0023
                0024 C     == Global variables ==
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
b487ac3d03 Jean*0027 #include "GRID.h"
                0028 #ifdef ALLOW_EXCH2
f9f661930b Jean*0029 #include "W2_EXCH2_SIZE.h"
b487ac3d03 Jean*0030 #include "W2_EXCH2_TOPOLOGY.h"
                0031 #endif /* ALLOW_EXCH2 */
aea29c8517 Alis*0032
2be0e061d5 Jean*0033 C     !INPUT/OUTPUT PARAMETERS:
aea29c8517 Alis*0034 C     == Routine arguments ==
2be0e061d5 Jean*0035       INTEGER bi,bj,k
aea29c8517 Alis*0036       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0037       _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0038       _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0039       _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
3370e71df9 Mart*0040       INTEGER selectVortScheme
                0041       LOGICAL highOrderVorticity,upwindVorticity
aea29c8517 Alis*0042       INTEGER myThid
2be0e061d5 Jean*0043 CEOP
aea29c8517 Alis*0044
                0045 C     == Local variables ==
7fe6343684 Jean*0046 C     msgBuf :: Informational/error meesage buffer
                0047       CHARACTER*(MAX_LEN_MBUF) msgBuf
2be0e061d5 Jean*0048       INTEGER i,j
b487ac3d03 Jean*0049       _RL vort3r(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df5ebd9340 Jean*0050       _RL vBarXY,vort3u,Rjp,Rjm,Rj
b487ac3d03 Jean*0051       _RL vBarXm,vBarXp
                0052
                0053       LOGICAL northWestCorner, northEastCorner,
                0054      &        southWestCorner, southEastCorner
                0055       INTEGER myFace
                0056 #ifdef ALLOW_EXCH2
                0057       INTEGER myTile
                0058 #endif /* ALLOW_EXCH2 */
                0059       _RL oneSixth, oneTwelve
aea29c8517 Alis*0060       LOGICAL fourthVort3
df5ebd9340 Jean*0061 C jmc: not sure about these 1/6 & 1/12 factors (should use the same)
b487ac3d03 Jean*0062       PARAMETER(oneSixth=1.D0/6.D0 , oneTwelve=1.D0/12.D0)
2be0e061d5 Jean*0063       PARAMETER(fourthVort3=.TRUE. )
                0064
b487ac3d03 Jean*0065 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0066
                0067       DO j=1-Oly,sNy+Oly
                0068        DO i=1-Olx,sNx+Olx
                0069          vort3r(i,j) = r_hFacZ(i,j)*omega3(i,j)
                0070        ENDDO
                0071       ENDDO
                0072
                0073 C--   Special stuff for Cubed Sphere
7fe6343684 Jean*0074       IF ( useCubedSphereExchange.AND.highOrderVorticity ) THEN
b487ac3d03 Jean*0075
                0076 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0077        myTile = W2_myTileList(bi,bj)
b487ac3d03 Jean*0078        myFace = exch2_myFace(myTile)
                0079        southWestCorner = exch2_isWedge(myTile).EQ.1
                0080      &             .AND. exch2_isSedge(myTile).EQ.1
                0081        southEastCorner = exch2_isEedge(myTile).EQ.1
                0082      &             .AND. exch2_isSedge(myTile).EQ.1
                0083        northEastCorner = exch2_isEedge(myTile).EQ.1
                0084      &             .AND. exch2_isNedge(myTile).EQ.1
                0085        northWestCorner = exch2_isWedge(myTile).EQ.1
                0086      &             .AND. exch2_isNedge(myTile).EQ.1
                0087 #else
                0088        myFace = bi
                0089        southWestCorner = .TRUE.
                0090        southEastCorner = .TRUE.
                0091        northWestCorner = .TRUE.
                0092        northEastCorner = .TRUE.
                0093 #endif /* ALLOW_EXCH2 */
                0094        IF ( southWestCorner ) THEN
                0095          i = 1
                0096          j = 1
                0097          vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i+1,j) )*0.5 _d 0
                0098        ENDIF
                0099        IF ( southEastCorner ) THEN
                0100          i = sNx+1
                0101          j = 1
                0102          vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i-1,j) )*0.5 _d 0
                0103        ENDIF
                0104        IF ( northWestCorner ) THEN
                0105          i = 1
                0106          j = sNy+1
                0107          vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i+1,j) )*0.5 _d 0
                0108        ENDIF
                0109        IF ( northEastCorner ) THEN
                0110          i = sNx+1
                0111          j = sNy+1
                0112          vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i-1,j) )*0.5 _d 0
                0113        ENDIF
                0114
                0115 C--   End of special stuff for Cubed Sphere.
                0116       ENDIF
                0117
                0118 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2be0e061d5 Jean*0119
7fe6343684 Jean*0120       IF ( selectVortScheme.EQ.0 ) THEN
                0121 C--   using Sadourny Enstrophy conserving discretization:
2be0e061d5 Jean*0122
7fe6343684 Jean*0123 c      DO j=2-Oly,sNy+Oly-2
                0124 c       DO i=2-Olx,sNx+Olx
                0125        DO j=1,sNy
                0126         DO i=1,sNx+1
aea29c8517 Alis*0127
                0128          vBarXY=0.25*(
616600b8d2 Patr*0129      &      (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
                0130      &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))
                0131      &     +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
                0132      &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj))
5d7e0a8948 Jean*0133      &               )
7fe6343684 Jean*0134          IF (upwindVorticity) THEN
aea29c8517 Alis*0135           IF (vBarXY.GT.0.) THEN
b487ac3d03 Jean*0136            vort3u=vort3r(i,j)
aea29c8517 Alis*0137           ELSE
b487ac3d03 Jean*0138            vort3u=vort3r(i,j+1)
aea29c8517 Alis*0139           ENDIF
                0140          ELSEIF (fourthVort3) THEN
df5ebd9340 Jean*0141 #ifdef ALLOW_OBCS
                0142           Rjp = ( vort3r(i,j+2) - vort3r(i,j+1) )*maskInW(i,j+1,bi,bj)
                0143           Rjm = ( vort3r(i, j ) - vort3r(i,j-1) )*maskInW(i,j-1,bi,bj)
                0144 #else
                0145           Rjp =   vort3r(i,j+2) - vort3r(i,j+1)
                0146           Rjm =   vort3r(i, j ) - vort3r(i,j-1)
                0147 #endif
                0148           vort3u=0.5*( (vort3r(i,j) + vort3r(i,j+1))
                0149      &                 -oneTwelve*(Rjp-Rjm)
aea29c8517 Alis*0150      &               )
                0151          ELSE
b487ac3d03 Jean*0152           vort3u=0.5*( vort3r(i,j) + vort3r(i,j+1) )
aea29c8517 Alis*0153          ENDIF
460dc72355 Patr*0154
7fe6343684 Jean*0155          uCoriolisTerm(i,j) =  vort3u*vBarXY*recip_dxC(i,j,bi,bj)
                0156      &                               * _maskW(i,j,k,bi,bj)
                0157
                0158         ENDDO
                0159        ENDDO
                0160
                0161       ELSEIF ( selectVortScheme.EQ.2 ) THEN
                0162 C--   using Energy conserving discretization:
2be0e061d5 Jean*0163
7fe6343684 Jean*0164 c      DO j=2-Oly,sNy+Oly-2
                0165 c       DO i=2-Olx,sNx+Olx
                0166        DO j=1,sNy
                0167         DO i=1,sNx+1
2be0e061d5 Jean*0168
7fe6343684 Jean*0169          vBarXm=0.5*(
                0170      &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
                0171      &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
                0172          vBarXp=0.5*(
                0173      &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
                0174      &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
                0175          IF (upwindVorticity) THEN
                0176           IF ( (vBarXm+vBarXp) .GT.0.) THEN
                0177            vort3u=vBarXm*vort3r(i, j )
                0178           ELSE
                0179            vort3u=vBarXp*vort3r(i,j+1)
                0180           ENDIF
                0181          ELSEIF (fourthVort3) THEN
df5ebd9340 Jean*0182 #ifdef ALLOW_OBCS
                0183           Rjp = ( vort3r(i,j+2) - vort3r(i,j+1) )*maskInW(i,j+1,bi,bj)
                0184           Rjm = ( vort3r(i, j ) - vort3r(i,j-1) )*maskInW(i,j-1,bi,bj)
                0185 #else
                0186           Rjp =   vort3r(i,j+2) - vort3r(i,j+1)
                0187           Rjm =   vort3r(i, j ) - vort3r(i,j-1)
                0188 #endif
                0189           Rj  =   vort3r(i,j+1) - vort3r(i, j )
a70729d476 Jean*0190           Rjp = vort3r(i,j+1) -oneSixth*( Rjp-Rj )
                0191           Rjm = vort3r(i, j ) -oneSixth*( Rj-Rjm )
df5ebd9340 Jean*0192 c         Rjp = vort3r(i,j+1) -oneSixth*( vort3r(i,j+2)-vort3r(i, j ) )
                0193 c         Rjm = vort3r(i, j ) +oneSixth*( vort3r(i,j+1)-vort3r(i,j-1) )
7fe6343684 Jean*0194           vort3u=0.5*( vBarXm*Rjm + vBarXp*Rjp )
                0195          ELSE
                0196           vort3u=0.5*( vBarXm*vort3r(i, j ) + vBarXp*vort3r(i,j+1) )
                0197          ENDIF
                0198
                0199          uCoriolisTerm(i,j) =  vort3u*recip_dxC(i,j,bi,bj)
                0200      &                               * _maskW(i,j,k,bi,bj)
                0201
                0202         ENDDO
aea29c8517 Alis*0203        ENDDO
7fe6343684 Jean*0204
                0205       ELSE
                0206         WRITE(msgBuf,'(A,I5,A)')
                0207      &   'MOM_VI_U_CORIOLIS_C4: selectVortScheme=', selectVortScheme,
                0208      &   ' not implemented'
                0209         CALL PRINT_ERROR( msgBuf, myThid )
                0210         STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS_C4'
                0211
                0212       ENDIF
aea29c8517 Alis*0213
                0214       RETURN
                0215       END