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

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

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

revision 1.3 by heimbach, Wed Sep 5 17:46:03 2001 UTC revision 1.13 by jmc, Mon May 5 22:45:00 2008 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    
6        SUBROUTINE MOM_VI_U_CORIOLIS(  CBOP
7       I        bi,bj,K,  C     !ROUTINE: MOM_VI_U_CORIOLIS
8       I        vFld,omega3,r_hFacZ,  C     !INTERFACE:
9       O        uCoriolisTerm,        SUBROUTINE MOM_VI_U_CORIOLIS(
10       I        myThid)       I                     bi, bj, k,
11         I                     vFld, omega3, hFacZ, r_hFacZ,
12         O                     uCoriolisTerm,
13         I                     myThid )
14    C     !DESCRIPTION: \bv
15    C     *==========================================================*
16    C     | S/R MOM_VI_U_CORIOLIS
17    C     |==========================================================*
18    C     | o Calculate flux (in Y-dir.) of vorticity at U point
19    C     |   using 2nd 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"
# Line 20  C     == Global variables == Line 29  C     == Global variables ==
29  #include "GRID.h"  #include "GRID.h"
30  #include "PARAMS.h"  #include "PARAMS.h"
31    
32    C     !INPUT/OUTPUT PARAMETERS:
33  C     == Routine arguments ==  C     == Routine arguments ==
34        INTEGER bi,bj,K        INTEGER bi, bj, k
35        _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     vFld     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36        _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     omega3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
37        _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS     hFacZ    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38          _RS     r_hFacZ  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39        _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40        INTEGER myThid        INTEGER myThid
41    CEOP
42    
43  C     == Local variables ==  C     == Local variables ==
44        INTEGER I,J  C     msgBuf :: Informational/error meesage buffer
45        _RL vBarXY,vort3u        CHARACTER*(MAX_LEN_MBUF) msgBuf
46        LOGICAL upwindVort3        LOGICAL upwindVort3
47          INTEGER i, j
48          _RL     vBarXY, vBarXm, vBarXp
49          _RL     vort3u
50          _RL     vort3mj, vort3ij, vort3mp, vort3ip
51          _RL     oneThird, tmpFac
52          _RS     epsil
53          PARAMETER( upwindVort3 = .FALSE. )
54    
55          epsil = 1. _d -9
56          tmpFac = 1. _d 0
57    c     oneThird = 1. _d 0 / ( 1. _d 0 + 2.*tmpFac )
58          oneThird = 1. _d 0 / 3. _d 0
59    
60        upwindVort3=.FALSE.        IF ( selectVortScheme.EQ.0 ) THEN
61    C--   using enstrophy conserving scheme (Shallow-Water Eq.) by Sadourny, JAS 75
62    
63         DO J=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
64          DO I=2-Olx,sNx+Olx          DO i=2-Olx,sNx+Olx
65           vBarXY=0.25*(           vBarXY=0.25*(
66       &       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)
67       &      +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))
68       &      +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)
69       &      +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))
70  c        vBarXY=0.25*(vFld( i ,j)+vFld( i ,j+1)       &               )
 c    &               +vFld(i-1,j)+vFld(i-1,j+1))  
71           IF (upwindVort3) THEN           IF (upwindVort3) THEN
72            IF (vBarXY.GT.0.) THEN            IF (vBarXY.GT.0.) THEN
73             vort3u=omega3(I,J)*r_hFacZ(i,j)             vort3u=omega3(i,j)*r_hFacZ(i,j)
74              ELSE
75               vort3u=omega3(i,j+1)*r_hFacZ(i,j+1)
76              ENDIF
77             ELSE
78               vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)
79         &                +omega3(i,j+1)*r_hFacZ(i,j+1))
80             ENDIF
81             uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
82         &                              * _maskW(i,j,k,bi,bj)
83            ENDDO
84           ENDDO
85    
86          ELSEIF ( selectVortScheme.EQ.1 ) THEN
87    C--   same as above, with different formulation (relatively to hFac)
88    
89           DO j=1-Oly,sNy+Oly-1
90            DO i=2-Olx,sNx+Olx
91             vBarXY= 0.5*(
92         &      (vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i, j )
93         &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i, j ))
94         &     +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacZ(i,j+1)
95         &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1))
96         &               )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) )
97             IF (upwindVort3) THEN
98              IF (vBarXY.GT.0.) THEN
99               vort3u=omega3(i,j)
100              ELSE
101               vort3u=omega3(i,j+1)
102              ENDIF
103             ELSE
104               vort3u=0.5*(omega3(i,j)+omega3(i,j+1))
105             ENDIF
106             uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
107         &                              * _maskW(i,j,k,bi,bj)
108            ENDDO
109           ENDDO
110    
111          ELSEIF ( selectVortScheme.EQ.2 ) THEN
112    C--   using energy conserving scheme (used by Sadourny in JAS 75 paper)
113    
114           DO j=1-Oly,sNy+Oly-1
115            DO i=2-Olx,sNx+Olx
116             vBarXm=0.5*(
117         &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
118         &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
119             vBarXp=0.5*(
120         &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
121         &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
122             IF (upwindVort3) THEN
123              IF ( (vBarXm+vBarXp) .GT.0.) THEN
124               vort3u=vBarXm*r_hFacZ(i, j )*omega3(i, j )
125            ELSE            ELSE
126             vort3u=omega3(I,J+1)*r_hFacZ(i,j+1)             vort3u=vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
127            ENDIF            ENDIF
128           ELSE           ELSE
129            vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)             vort3u = ( vBarXm*r_hFacZ(i, j )*omega3(i, j )
130       &               +omega3(i,j+1)*r_hFacZ(i,j+1))       &               +vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
131         &              )*0.5 _d 0
132           ENDIF           ENDIF
133           uCoriolisTerm(i,j)=           uCoriolisTerm(i,j)= +vort3u*recip_dxC(i,j,bi,bj)
134       &   +vort3u*vBarXY*recip_dxc(i,j,bi,bj)       &                              * _maskW(i,j,k,bi,bj)
135       &   *_maskW(I,J,K,bi,bj)          ENDDO
136  cph *note* put these comments after end of continued line         ENDDO
137  cph        to ensure TAMC compatibility  
138  C high order vorticity advection term        ELSEIF ( selectVortScheme.EQ.3 ) THEN
139  c    &   ...  C--   using energy & enstrophy conserving scheme
140  C linear Coriolis term  C     (from Sadourny, described by Burridge & Haseler, ECMWF Rep.4, 1977)
141  c    &   +0.5*(fCoriG(I,J,bi,bj)+fCoriG(I,J+1,bi,bj))*vBarXY  
142  C full nonlinear Coriolis term  C     domain where uCoriolisTerm is valid :
143  c    &   +0.5*(omega3(I,J)+omega3(I,J+1))*vBarXY  C     [ 3-Olx : sNx+Olx-1 ] x [ 2-Oly : sNy+Oly-1 ]
144  C correct energy conserving form of Coriolis term  C     (=> might need overlap of 3 if using CD-scheme)
145  c    &   +0.5*( fCori(I  ,J,bi,bj)*vBarY(I  ,J,K,bi,bj) +         DO j=1-Oly,sNy+Oly-1
146  c    &          fCori(I-1,J,bi,bj)*vBarY(I-1,J,K,bi,bj)  )          DO i=2-Olx,sNx+Olx-1
147  C original form of Coriolis term (copied from calc_mom_rhs)           vort3mj= ( r_hFacZ(i, j )*omega3(i, j )
148  c    &   +0.5*(fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj))*vBarXY       &            +(r_hFacZ(i,j+1)*omega3(i,j+1)
149         &             +r_hFacZ(i-1,j)*omega3(i-1,j)
150         &            ))*oneThird
151    c    &            )*tmpFac)*oneThird
152         &      *vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj)
153             vort3ij= ( r_hFacZ(i, j )*omega3(i, j )
154         &            +(r_hFacZ(i,j+1)*omega3(i,j+1)
155         &             +r_hFacZ(i+1,j)*omega3(i+1,j)
156         &            ))*oneThird
157    c    &            )*tmpFac)*oneThird
158         &      *vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
159             vort3mp= ( r_hFacZ(i,j+1)*omega3(i,j+1)
160         &            +(r_hFacZ(i, j )*omega3(i, j )
161         &             +r_hFacZ(i-1,j+1)*omega3(i-1,j+1)
162         &            ))*oneThird
163    c    &            )*tmpFac)*oneThird
164         &      *vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj)
165             vort3ip= ( r_hFacZ(i,j+1)*omega3(i,j+1)
166         &            +(r_hFacZ(i, j )*omega3(i, j )
167         &             +r_hFacZ(i+1,j+1)*omega3(i+1,j+1)
168         &            ))*oneThird
169    c    &            )*tmpFac)*oneThird
170         &      *vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
171    C---
172             uCoriolisTerm(i,j)= +( (vort3mj+vort3ij)+(vort3mp+vort3ip) )
173         &                     *0.25 _d 0 *recip_dxC(i,j,bi,bj)
174         &                                * _maskW(i,j,k,bi,bj)
175            ENDDO
176           ENDDO
177    
178          ELSE
179            WRITE(msgBuf,'(A,I5,A)')
180         &   'MOM_VI_U_CORIOLIS: selectVortScheme=', selectVortScheme,
181         &   ' not implemented'
182            CALL PRINT_ERROR( msgBuf, myThid )
183            STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS'
184    
185          ENDIF
186    
187          IF ( useJamartMomAdv ) THEN
188           DO j=1-Oly,sNy+Oly-1
189            DO i=2-Olx,sNx+Olx-1
190             uCoriolisTerm(i,j) = uCoriolisTerm(i,j)
191         &           * 4. _d 0 * _hFacW(i,j,k,bi,bj)
192         &           / MAX( epsil,
193         &                  (_hFacS(i, j ,k,bi,bj)+_hFacS(i-1, j ,k,bi,bj))
194         &                 +(_hFacS(i,j+1,k,bi,bj)+_hFacS(i-1,j+1,k,bi,bj))
195         &                )
196          ENDDO          ENDDO
197         ENDDO         ENDDO
198          ENDIF
199    
200        RETURN        RETURN
201        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22