C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vi_u_coriolis_c4.F,v 1.15 2017/04/28 17:17:14 mlosch Exp $ C $Name: $ #include "MOM_VECINV_OPTIONS.h" CBOP C !ROUTINE: MOM_VI_U_CORIOLIS_C4 C !INTERFACE: SUBROUTINE MOM_VI_U_CORIOLIS_C4( I bi,bj,k, I selectVortScheme,highOrderVorticity,upwindVorticity, I vFld,omega3,r_hFacZ, O uCoriolisTerm, I myThid) C !DESCRIPTION: \bv C *==========================================================* C | S/R MOM_VI_U_CORIOLIS_C4 C |==========================================================* C | o Calculate flux (in Y-dir.) of vorticity at U point C | using 4th order (or 1rst order) interpolation C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #ifdef ALLOW_EXCH2 #include "W2_EXCH2_SIZE.h" #include "W2_EXCH2_TOPOLOGY.h" #endif /* ALLOW_EXCH2 */ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == INTEGER bi,bj,k _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER selectVortScheme LOGICAL highOrderVorticity,upwindVorticity INTEGER myThid CEOP C == Local variables == C msgBuf :: Informational/error meesage buffer CHARACTER*(MAX_LEN_MBUF) msgBuf INTEGER i,j _RL vort3r(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vBarXY,vort3u,Rjp,Rjm,Rj _RL vBarXm,vBarXp LOGICAL northWestCorner, northEastCorner, & southWestCorner, southEastCorner INTEGER myFace #ifdef ALLOW_EXCH2 INTEGER myTile #endif /* ALLOW_EXCH2 */ _RL oneSixth, oneTwelve LOGICAL fourthVort3 C jmc: not sure about these 1/6 & 1/12 factors (should use the same) PARAMETER(oneSixth=1.D0/6.D0 , oneTwelve=1.D0/12.D0) PARAMETER(fourthVort3=.TRUE. ) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx vort3r(i,j) = r_hFacZ(i,j)*omega3(i,j) ENDDO ENDDO C-- Special stuff for Cubed Sphere IF ( useCubedSphereExchange.AND.highOrderVorticity ) THEN #ifdef ALLOW_EXCH2 myTile = W2_myTileList(bi,bj) myFace = exch2_myFace(myTile) southWestCorner = exch2_isWedge(myTile).EQ.1 & .AND. exch2_isSedge(myTile).EQ.1 southEastCorner = exch2_isEedge(myTile).EQ.1 & .AND. exch2_isSedge(myTile).EQ.1 northEastCorner = exch2_isEedge(myTile).EQ.1 & .AND. exch2_isNedge(myTile).EQ.1 northWestCorner = exch2_isWedge(myTile).EQ.1 & .AND. exch2_isNedge(myTile).EQ.1 #else myFace = bi southWestCorner = .TRUE. southEastCorner = .TRUE. northWestCorner = .TRUE. northEastCorner = .TRUE. #endif /* ALLOW_EXCH2 */ IF ( southWestCorner ) THEN i = 1 j = 1 vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i+1,j) )*0.5 _d 0 ENDIF IF ( southEastCorner ) THEN i = sNx+1 j = 1 vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i-1,j) )*0.5 _d 0 ENDIF IF ( northWestCorner ) THEN i = 1 j = sNy+1 vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i+1,j) )*0.5 _d 0 ENDIF IF ( northEastCorner ) THEN i = sNx+1 j = sNy+1 vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i-1,j) )*0.5 _d 0 ENDIF C-- End of special stuff for Cubed Sphere. ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( selectVortScheme.EQ.0 ) THEN C-- using Sadourny Enstrophy conserving discretization: c DO j=2-Oly,sNy+Oly-2 c DO i=2-Olx,sNx+Olx DO j=1,sNy DO i=1,sNx+1 vBarXY=0.25*( & (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj)) & +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj)) & ) IF (upwindVorticity) THEN IF (vBarXY.GT.0.) THEN vort3u=vort3r(i,j) ELSE vort3u=vort3r(i,j+1) ENDIF ELSEIF (fourthVort3) THEN #ifdef ALLOW_OBCS Rjp = ( vort3r(i,j+2) - vort3r(i,j+1) )*maskInW(i,j+1,bi,bj) Rjm = ( vort3r(i, j ) - vort3r(i,j-1) )*maskInW(i,j-1,bi,bj) #else Rjp = vort3r(i,j+2) - vort3r(i,j+1) Rjm = vort3r(i, j ) - vort3r(i,j-1) #endif vort3u=0.5*( (vort3r(i,j) + vort3r(i,j+1)) & -oneTwelve*(Rjp-Rjm) & ) ELSE vort3u=0.5*( vort3r(i,j) + vort3r(i,j+1) ) ENDIF uCoriolisTerm(i,j) = vort3u*vBarXY*recip_dxC(i,j,bi,bj) & * _maskW(i,j,k,bi,bj) ENDDO ENDDO ELSEIF ( selectVortScheme.EQ.2 ) THEN C-- using Energy conserving discretization: c DO j=2-Oly,sNy+Oly-2 c DO i=2-Olx,sNx+Olx DO j=1,sNy DO i=1,sNx+1 vBarXm=0.5*( & vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) ) vBarXp=0.5*( & vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) ) IF (upwindVorticity) THEN IF ( (vBarXm+vBarXp) .GT.0.) THEN vort3u=vBarXm*vort3r(i, j ) ELSE vort3u=vBarXp*vort3r(i,j+1) ENDIF ELSEIF (fourthVort3) THEN #ifdef ALLOW_OBCS Rjp = ( vort3r(i,j+2) - vort3r(i,j+1) )*maskInW(i,j+1,bi,bj) Rjm = ( vort3r(i, j ) - vort3r(i,j-1) )*maskInW(i,j-1,bi,bj) #else Rjp = vort3r(i,j+2) - vort3r(i,j+1) Rjm = vort3r(i, j ) - vort3r(i,j-1) #endif Rj = vort3r(i,j+1) - vort3r(i, j ) Rjp = vort3r(i,j+1) -oneSixth*( Rjp-Rj ) Rjm = vort3r(i, j ) -oneSixth*( Rj-Rjm ) c Rjp = vort3r(i,j+1) -oneSixth*( vort3r(i,j+2)-vort3r(i, j ) ) c Rjm = vort3r(i, j ) +oneSixth*( vort3r(i,j+1)-vort3r(i,j-1) ) vort3u=0.5*( vBarXm*Rjm + vBarXp*Rjp ) ELSE vort3u=0.5*( vBarXm*vort3r(i, j ) + vBarXp*vort3r(i,j+1) ) ENDIF uCoriolisTerm(i,j) = vort3u*recip_dxC(i,j,bi,bj) & * _maskW(i,j,k,bi,bj) ENDDO ENDDO ELSE WRITE(msgBuf,'(A,I5,A)') & 'MOM_VI_U_CORIOLIS_C4: selectVortScheme=', selectVortScheme, & ' not implemented' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS_C4' ENDIF RETURN END