/[MITgcm]/MITgcm/pkg/mom_vecinv/mom_vi_u_coriolis_c4.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_vecinv/mom_vi_u_coriolis_c4.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by adcroft, Tue Jul 20 17:46:38 2004 UTC revision 1.10 by jmc, Tue Apr 1 01:27:33 2008 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "MOM_VECINV_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    
6        SUBROUTINE MOM_VI_U_CORIOLIS_C4(  CBOP
7       I        bi,bj,K,  C     !ROUTINE: MOM_VI_U_CORIOLIS_C4
8    C     !INTERFACE:
9          SUBROUTINE MOM_VI_U_CORIOLIS_C4(
10         I        bi,bj,k,
11       I        vFld,omega3,r_hFacZ,       I        vFld,omega3,r_hFacZ,
12       O        uCoriolisTerm,       O        uCoriolisTerm,
13       I        myThid)       I        myThid)
14    C     !DESCRIPTION: \bv
15    C     *==========================================================*
16    C     | S/R MOM_VI_U_CORIOLIS_C4
17    C     |==========================================================*
18    C     | o Calculate flux (in Y-dir.) of vorticity at U point
19    C     |   using 4th order (or 1rst order) interpolation
20    C     *==========================================================*
21    C     \ev
22    
23    C     !USES:
24        IMPLICIT NONE        IMPLICIT NONE
 C     /==========================================================\  
 C     | S/R MOM_VI_U_CORIOLIS                                    |  
 C     |==========================================================|  
 C     \==========================================================/  
25    
26  C     == Global variables ==  C     == Global variables ==
27  #include "SIZE.h"  #include "SIZE.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "GRID.h"  
29  #include "PARAMS.h"  #include "PARAMS.h"
30    #include "GRID.h"
31    #ifdef ALLOW_EXCH2
32    #include "W2_EXCH2_TOPOLOGY.h"
33    #include "W2_EXCH2_PARAMS.h"
34    #endif /* ALLOW_EXCH2 */
35    
36    C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38        INTEGER bi,bj,K        INTEGER bi,bj,k
39        _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40        _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41        _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42        _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43        INTEGER myThid        INTEGER myThid
44    CEOP
45    
46  C     == Local variables ==  C     == Local variables ==
47        INTEGER I,J  C     msgBuf :: Informational/error meesage buffer
48          CHARACTER*(MAX_LEN_MBUF) msgBuf
49          INTEGER i,j
50          _RL vort3r(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51        _RL vBarXY,vort3u,Rjp,Rjm        _RL vBarXY,vort3u,Rjp,Rjm
52        LOGICAL upwindVort3        _RL vBarXm,vBarXp
53    
54          LOGICAL northWestCorner, northEastCorner,
55         &        southWestCorner, southEastCorner
56          INTEGER myFace
57    #ifdef ALLOW_EXCH2
58          INTEGER myTile
59    #endif /* ALLOW_EXCH2 */
60          _RL oneSixth, oneTwelve
61        LOGICAL fourthVort3        LOGICAL fourthVort3
62          PARAMETER(oneSixth=1.D0/6.D0 , oneTwelve=1.D0/12.D0)
63          PARAMETER(fourthVort3=.TRUE. )
64    
65    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
66    
67          DO j=1-Oly,sNy+Oly
68           DO i=1-Olx,sNx+Olx
69             vort3r(i,j) = r_hFacZ(i,j)*omega3(i,j)
70           ENDDO
71          ENDDO
72    
73    C--   Special stuff for Cubed Sphere
74          IF ( useCubedSphereExchange.AND.highOrderVorticity ) THEN
75    
76        upwindVort3=.FALSE.  #ifdef ALLOW_EXCH2
77        fourthVort3=.TRUE.         myTile = W2_myTileList(bi)
78           myFace = exch2_myFace(myTile)
79           southWestCorner = exch2_isWedge(myTile).EQ.1
80         &             .AND. exch2_isSedge(myTile).EQ.1
81           southEastCorner = exch2_isEedge(myTile).EQ.1
82         &             .AND. exch2_isSedge(myTile).EQ.1
83           northEastCorner = exch2_isEedge(myTile).EQ.1
84         &             .AND. exch2_isNedge(myTile).EQ.1
85           northWestCorner = exch2_isWedge(myTile).EQ.1
86         &             .AND. exch2_isNedge(myTile).EQ.1
87    #else
88           myFace = bi
89           southWestCorner = .TRUE.
90           southEastCorner = .TRUE.
91           northWestCorner = .TRUE.
92           northEastCorner = .TRUE.
93    #endif /* ALLOW_EXCH2 */
94           IF ( southWestCorner ) THEN
95             i = 1
96             j = 1
97             vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i+1,j) )*0.5 _d 0
98           ENDIF
99           IF ( southEastCorner ) THEN
100             i = sNx+1
101             j = 1
102             vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i-1,j) )*0.5 _d 0
103           ENDIF
104           IF ( northWestCorner ) THEN
105             i = 1
106             j = sNy+1
107             vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i+1,j) )*0.5 _d 0
108           ENDIF
109           IF ( northEastCorner ) THEN
110             i = sNx+1
111             j = sNy+1
112             vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i-1,j) )*0.5 _d 0
113           ENDIF
114    
115    C--   End of special stuff for Cubed Sphere.
116          ENDIF
117    
118    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
119    
120          IF ( selectVortScheme.EQ.0 ) THEN
121    C--   using Sadourny Enstrophy conserving discretization:
122    
123    c      DO j=2-Oly,sNy+Oly-2
124    c       DO i=2-Olx,sNx+Olx
125           DO j=1,sNy
126            DO i=1,sNx+1
127    
        DO J=1-Oly,sNy+Oly-1  
         DO I=2-Olx,sNx+Olx  
128           vBarXY=0.25*(           vBarXY=0.25*(
129       &       vFld( i , j )*dxG( i , j ,bi,bj)*hFacS( i , j ,k,bi,bj)       &      (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
130       &      +vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacS( i ,j+1,k,bi,bj)       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))
131       &      +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)
132       &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacS(i-1,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))
133  c        vBarXY=0.25*(vFld( i ,j)+vFld( i ,j+1)       &               )
134  c    &               +vFld(i-1,j)+vFld(i-1,j+1))           IF (upwindVorticity) THEN
          IF (upwindVort3) THEN  
135            IF (vBarXY.GT.0.) THEN            IF (vBarXY.GT.0.) THEN
136             vort3u=omega3(I,J)*r_hFacZ(i,j)             vort3u=vort3r(i,j)
137            ELSE            ELSE
138             vort3u=omega3(I,J+1)*r_hFacZ(i,j+1)             vort3u=vort3r(i,j+1)
139            ENDIF            ENDIF
140           ELSEIF (fourthVort3) THEN           ELSEIF (fourthVort3) THEN
141            Rjp=omega3(i,j+2)*r_hFacZ(i,j+2)            Rjp = vort3r(i,j+2) - vort3r(i,j+1)
142       &       -omega3(i,j+1)*r_hFacZ(i,j+1)            Rjm = vort3r(i, j ) - vort3r(i,j-1)
143            Rjm=omega3(i,j)*r_hFacZ(i,j)            vort3u=0.5*( vort3r(i,j) + vort3r(i,j+1)
144       &       -omega3(i,j-1)*r_hFacZ(i,j-1)       &                -oneTwelve*(Rjp-Rjm)
           vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)  
      &               +omega3(i,j+1)*r_hFacZ(i,j+1)  
      &               -1./12.*(Rjp-Rjm)  
145       &               )       &               )
146           ELSE           ELSE
147            vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)            vort3u=0.5*( vort3r(i,j) + vort3r(i,j+1) )
148       &               +omega3(i,j+1)*r_hFacZ(i,j+1))           ENDIF
149    
150             uCoriolisTerm(i,j) =  vort3u*vBarXY*recip_dxC(i,j,bi,bj)
151         &                               * _maskW(i,j,k,bi,bj)
152    
153            ENDDO
154           ENDDO
155    
156          ELSEIF ( selectVortScheme.EQ.2 ) THEN
157    C--   using Energy conserving discretization:
158    
159    c      DO j=2-Oly,sNy+Oly-2
160    c       DO i=2-Olx,sNx+Olx
161           DO j=1,sNy
162            DO i=1,sNx+1
163    
164             vBarXm=0.5*(
165         &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
166         &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
167             vBarXp=0.5*(
168         &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
169         &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
170             IF (upwindVorticity) THEN
171              IF ( (vBarXm+vBarXp) .GT.0.) THEN
172               vort3u=vBarXm*vort3r(i, j )
173              ELSE
174               vort3u=vBarXp*vort3r(i,j+1)
175              ENDIF
176             ELSEIF (fourthVort3) THEN
177              Rjp = vort3r(i,j+1) -oneSixth*( vort3r(i,j+2)-vort3r(i, j ) )
178              Rjm = vort3r(i, j ) +oneSixth*( vort3r(i,j+1)-vort3r(i,j-1) )
179              vort3u=0.5*( vBarXm*Rjm + vBarXp*Rjp )
180             ELSE
181              vort3u=0.5*( vBarXm*vort3r(i, j ) + vBarXp*vort3r(i,j+1) )
182           ENDIF           ENDIF
183           uCoriolisTerm(i,j)=  
184  C high order vorticity advection term           uCoriolisTerm(i,j) =  vort3u*recip_dxC(i,j,bi,bj)
185       &   +vort3u*vBarXY*recip_dxc(i,j,bi,bj)       &                               * _maskW(i,j,k,bi,bj)
186  C linear Coriolis term  
 c    &   +0.5*(fCoriG(I,J,bi,bj)+fCoriG(I,J+1,bi,bj))*vBarXY  
 C full nonlinear Coriolis term  
 c    &   +0.5*(omega3(I,J)+omega3(I,J+1))*vBarXY  
 C correct energy conserving form of Coriolis term  
 c    &   +0.5*( fCori(I  ,J,bi,bj)*vBarY(I  ,J,K,bi,bj) +  
 c    &          fCori(I-1,J,bi,bj)*vBarY(I-1,J,K,bi,bj)  )  
 C original form of Coriolis term (copied from calc_mom_rhs)  
 c    &   +0.5*(fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj))*vBarXY  
      &   *_maskW(I,J,K,bi,bj)  
187          ENDDO          ENDDO
188         ENDDO         ENDDO
189    
190          ELSE
191            WRITE(msgBuf,'(A,I5,A)')
192         &   'MOM_VI_U_CORIOLIS_C4: selectVortScheme=', selectVortScheme,
193         &   ' not implemented'
194            CALL PRINT_ERROR( msgBuf, myThid )
195            STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS_C4'
196    
197          ENDIF
198    
199        RETURN        RETURN
200        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22