/[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.2 by adcroft, Tue May 29 14:01:39 2001 UTC revision 1.12 by jmc, Wed Apr 2 21:18:45 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          _RS     epsil
51          PARAMETER( upwindVort3 = .FALSE. )
52    
53        upwindVort3=.FALSE.        epsil = 1. _d -9
54    
55         DO J=1-Oly,sNy+Oly-1        IF ( selectVortScheme.EQ.0 ) THEN
56          DO I=2-Olx,sNx+Olx  C--   using enstrophy conserving scheme (Shallow-Water Eq.) by Sadourny, JAS 75
57    
58           DO j=1-Oly,sNy+Oly-1
59            DO i=2-Olx,sNx+Olx
60           vBarXY=0.25*(           vBarXY=0.25*(
61       &       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)
62       &      +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))
63       &      +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)
64       &      +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))
65  c        vBarXY=0.25*(vFld( i ,j)+vFld( i ,j+1)       &               )
66  c    &               +vFld(i-1,j)+vFld(i-1,j+1))           IF (upwindVort3) THEN
67              IF (vBarXY.GT.0.) THEN
68               vort3u=omega3(i,j)*r_hFacZ(i,j)
69              ELSE
70               vort3u=omega3(i,j+1)*r_hFacZ(i,j+1)
71              ENDIF
72             ELSE
73               vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)
74         &                +omega3(i,j+1)*r_hFacZ(i,j+1))
75             ENDIF
76             uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
77         &                              * _maskW(i,j,k,bi,bj)
78            ENDDO
79           ENDDO
80    
81          ELSEIF ( selectVortScheme.EQ.1 ) THEN
82    C--   same as above, with different formulation (relatively to hFac)
83    
84           DO J=1-Oly,sNy+Oly-1
85            DO i=2-Olx,sNx+Olx
86             vBarXY= 0.5*(
87         &      (vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i, j )
88         &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i, j ))
89         &     +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacZ(i,j+1)
90         &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1))
91         &               )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) )
92           IF (upwindVort3) THEN           IF (upwindVort3) THEN
93            IF (vBarXY.GT.0.) THEN            IF (vBarXY.GT.0.) THEN
94             vort3u=omega3(I,J)*r_hFacZ(i,j)             vort3u=omega3(i,j)
95            ELSE            ELSE
96             vort3u=omega3(I,J+1)*r_hFacZ(i,j+1)             vort3u=omega3(i,j+1)
97            ENDIF            ENDIF
98           ELSE           ELSE
99            vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)             vort3u=0.5*(omega3(i,j)+omega3(i,j+1))
      &               +omega3(i,j+1)*r_hFacZ(i,j+1))  
100           ENDIF           ENDIF
101           uCoriolisTerm(i,j)=           uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
102  C high order vorticity advection term       &                              * _maskW(i,j,k,bi,bj)
103       &   +vort3u*vBarXY*recip_dxc(i,j,bi,bj)          ENDDO
104  C linear Coriolis term         ENDDO
105  c    &   +0.5*(fCoriG(I,J,bi,bj)+fCoriG(I,J+1,bi,bj))*vBarXY  
106  C full nonlinear Coriolis term        ELSEIF ( selectVortScheme.EQ.2 ) THEN
107  c    &   +0.5*(omega3(I,J)+omega3(I,J+1))*vBarXY  C--   using energy conserving scheme (used by Sadourny in JAS 75 paper)
108  C correct energy conserving form of Coriolis term  
109  c    &   +0.5*( fCori(I  ,J,bi,bj)*vBarY(I  ,J,K,bi,bj) +         DO j=1-Oly,sNy+Oly-1
110  c    &          fCori(I-1,J,bi,bj)*vBarY(I-1,J,K,bi,bj)  )          DO i=2-Olx,sNx+Olx
111  C original form of Coriolis term (copied from calc_mom_rhs)           vBarXm=0.5*(
112  c    &   +0.5*(fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj))*vBarXY       &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
113       &   *_maskW(I,J,K,bi,bj)       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
114             vBarXp=0.5*(
115         &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
116         &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
117             IF (upwindVort3) THEN
118              IF ( (vBarXm+vBarXp) .GT.0.) THEN
119               vort3u=vBarXm*r_hFacZ(i, j )*omega3(i, j )
120              ELSE
121               vort3u=vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
122              ENDIF
123             ELSE
124               vort3u = ( vBarXm*r_hFacZ(i, j )*omega3(i, j )
125         &               +vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
126         &              )*0.5 _d 0
127             ENDIF
128             uCoriolisTerm(i,j)= +vort3u*recip_dxC(i,j,bi,bj)
129         &                              * _maskW(i,j,k,bi,bj)
130            ENDDO
131           ENDDO
132    
133          ELSE
134            WRITE(msgBuf,'(A,I5,A)')
135         &   'MOM_VI_U_CORIOLIS: selectVortScheme=', selectVortScheme,
136         &   ' not implemented'
137            CALL PRINT_ERROR( msgBuf, myThid )
138            STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS'
139    
140          ENDIF
141    
142          IF ( useJamartMomAdv ) THEN
143           DO j=1-Oly,sNy+Oly-1
144            DO i=2-Olx,sNx+Olx
145             uCoriolisTerm(i,j) = uCoriolisTerm(i,j)
146         &           * 4. _d 0 * _hFacW(i,j,k,bi,bj)
147         &           / MAX( epsil,
148         &                  (_hFacS(i, j ,k,bi,bj)+_hFacS(i-1, j ,k,bi,bj))
149         &                 +(_hFacS(i,j+1,k,bi,bj)+_hFacS(i-1,j+1,k,bi,bj))
150         &                )
151          ENDDO          ENDDO
152         ENDDO         ENDDO
153          ENDIF
154    
155        RETURN        RETURN
156        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22