/[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.5 by adcroft, Wed Jun 2 13:23:55 2004 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,hFacZ,r_hFacZ,  C     !INTERFACE:
9       O        uCoriolisTerm,        SUBROUTINE MOM_VI_U_CORIOLIS(
10       I        myThid)       I                     bi, bj, k,
11        IMPLICIT NONE       I                     vFld, omega3, hFacZ, r_hFacZ,
12         O                     uCoriolisTerm,
13         I                     myThid )
14    C     !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
16  C     | S/R MOM_VI_U_CORIOLIS  C     | S/R MOM_VI_U_CORIOLIS
17  C     | o Calculate meridional flux of vorticity at U point  C     |==========================================================*
18    C     | o Calculate flux (in Y-dir.) of vorticity at U point
19    C     |   using 2nd order interpolation
20  C     *==========================================================*  C     *==========================================================*
21    C     \ev
22    
23    C     !USES:
24          IMPLICIT NONE
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 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)        _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        LOGICAL use_original_hFac  C     msgBuf :: Informational/error meesage buffer
45        INTEGER I,J        CHARACTER*(MAX_LEN_MBUF) msgBuf
46        _RL vBarXY,vort3u,Zp,Zm        LOGICAL upwindVort3
47        _RS epsil        INTEGER i, j
48        PARAMETER ( use_original_hFac=.FALSE. )        _RL     vBarXY, vBarXm, vBarXp
49          _RL     vort3u
50        epsil = 1. _d -9        _RL     vort3mj, vort3ij, vort3mp, vort3ip
51          _RL     oneThird, tmpFac
52        DO J=1-Oly,sNy+Oly-1        _RS     epsil
53         DO I=2-Olx,sNx+Olx        PARAMETER( upwindVort3 = .FALSE. )
54          IF ( use_original_hFac ) THEN  
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          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
64            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           IF (upwindVorticity) THEN       &               )
71             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            ELSE
75             vort3u=omega3(I,J+1)*r_hFacZ(i,j+1)             vort3u=omega3(i,j+1)*r_hFacZ(i,j+1)
76            ENDIF            ENDIF
77           ELSE           ELSE
78             vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)             vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)
79       &                +omega3(i,j+1)*r_hFacZ(i,j+1))       &                +omega3(i,j+1)*r_hFacZ(i,j+1))
80           ENDIF           ENDIF
81          ELSEIF ( SadournyCoriolis ) THEN           uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
82           Zm=0.5*(       &                              * _maskW(i,j,k,bi,bj)
83       &       vFld( i , j )*dxG( i , j ,bi,bj)*hFacS( i , j ,k,bi,bj)          ENDDO
84       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacS(i-1, j ,k,bi,bj) )         ENDDO
85           Zp=0.5*(  
86       &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacS( i ,j+1,k,bi,bj)        ELSEIF ( selectVortScheme.EQ.1 ) THEN
87       &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacS(i-1,j+1,k,bi,bj) )  C--   same as above, with different formulation (relatively to hFac)
88           IF (upwindVorticity) THEN  
89            IF ( (Zm+Zp) .GT.0.) THEN         DO j=1-Oly,sNy+Oly-1
90             vort3u=Zm*r_hFacZ(i, j )*omega3(i, j )          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            ELSE
101             vort3u=Zp*r_hFacZ(i,j+1)*omega3(i,j+1)             vort3u=omega3(i,j+1)
102            ENDIF            ENDIF
103           ELSE           ELSE
104            Zm=Zm*r_hFacZ(i, j )*omega3(i, j )             vort3u=0.5*(omega3(i,j)+omega3(i,j+1))
           Zp=Zp*r_hFacZ(i,j+1)*omega3(i,j+1)  
           vort3u=0.5*( Zm + Zp )  
105           ENDIF           ENDIF
106           vBarXY=1.           uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
107          ELSE       &                              * _maskW(i,j,k,bi,bj)
108  c--      test a different formulation (relatively to hFac)          ENDDO
109           vBarXY=0.5*(         ENDDO
110       &       vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i,j)  
111       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i,j)        ELSEIF ( selectVortScheme.EQ.2 ) THEN
112       &      +vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacZ(i,j+1)  C--   using energy conserving scheme (used by Sadourny in JAS 75 paper)
113       &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1)  
114       &              )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) )         DO j=1-Oly,sNy+Oly-1
115           IF (upwindVorticity) THEN          DO i=2-Olx,sNx+Olx
116            IF (vBarXY.GT.0.) THEN           vBarXm=0.5*(
117             vort3u=omega3(i,j)       &       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)             vort3u=vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
127            ENDIF            ENDIF
128           ELSE           ELSE
129             vort3u=0.5*(omega3(i,j)+omega3(i,j+1))             vort3u = ( vBarXm*r_hFacZ(i, j )*omega3(i, j )
130         &               +vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
131         &              )*0.5 _d 0
132           ENDIF           ENDIF
133          ENDIF           uCoriolisTerm(i,j)= +vort3u*recip_dxC(i,j,bi,bj)
134         &                              * _maskW(i,j,k,bi,bj)
135            ENDDO
136           ENDDO
137    
138          ELSEIF ( selectVortScheme.EQ.3 ) THEN
139    C--   using energy & enstrophy conserving scheme
140    C     (from Sadourny, described by Burridge & Haseler, ECMWF Rep.4, 1977)
141    
142    C     domain where uCoriolisTerm is valid :
143    C     [ 3-Olx : sNx+Olx-1 ] x [ 2-Oly : sNy+Oly-1 ]
144    C     (=> might need overlap of 3 if using CD-scheme)
145           DO j=1-Oly,sNy+Oly-1
146            DO i=2-Olx,sNx+Olx-1
147             vort3mj= ( r_hFacZ(i, j )*omega3(i, j )
148         &            +(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          IF (useJamartWetpoints)        ELSE
179       &   vBarXY = vBarXY * 4. _d 0 * hFacW(i,j,k,bi,bj)          WRITE(msgBuf,'(A,I5,A)')
180       &   * MAX( epsil, hFacS( i , j ,k,bi,bj)+hFacS(i-1, j ,k,bi,bj)       &   'MOM_VI_U_CORIOLIS: selectVortScheme=', selectVortScheme,
181       &                +hFacS( j ,i+1,k,bi,bj)+hFacS(i-1,j+1,k,bi,bj) )       &   ' not implemented'
182            CALL PRINT_ERROR( msgBuf, myThid )
183          uCoriolisTerm(i,j)=          STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS'
184       &   +vort3u*vBarXY*recip_dxC(i,j,bi,bj)*_maskW(i,j,k,bi,bj)  
185  cph *note* put these comments after end of continued line        ENDIF
186  cph        to ensure TAMC compatibility  
187  C high order vorticity advection term        IF ( useJamartMomAdv ) THEN
188  c    &   ...         DO j=1-Oly,sNy+Oly-1
189  C linear Coriolis term          DO i=2-Olx,sNx+Olx-1
190  c    &   +0.5*(fCoriG(I,J,bi,bj)+fCoriG(I,J+1,bi,bj))*vBarXY           uCoriolisTerm(i,j) = uCoriolisTerm(i,j)
191  C full nonlinear Coriolis term       &           * 4. _d 0 * _hFacW(i,j,k,bi,bj)
192  c    &   +0.5*(omega3(I,J)+omega3(I,J+1))*vBarXY       &           / MAX( epsil,
193  C correct energy conserving form of Coriolis term       &                  (_hFacS(i, j ,k,bi,bj)+_hFacS(i-1, j ,k,bi,bj))
194  c    &   +0.5*( fCori(I  ,J,bi,bj)*vBarY(I  ,J,K,bi,bj) +       &                 +(_hFacS(i,j+1,k,bi,bj)+_hFacS(i-1,j+1,k,bi,bj))
195  c    &          fCori(I-1,J,bi,bj)*vBarY(I-1,J,K,bi,bj)  )       &                )
196  C original form of Coriolis term (copied from calc_mom_rhs)          ENDDO
 c    &   +0.5*(fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj))*vBarXY  
197         ENDDO         ENDDO
198        ENDDO        ENDIF
199    
200        RETURN        RETURN
201        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22