/[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.6 by adcroft, Mon Jun 14 17:48:17 2004 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,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        _RS     epsil
51          PARAMETER( upwindVort3 = .FALSE. )
52        DO J=1-Oly,sNy+Oly-1  
53         DO I=2-Olx,sNx+Olx        epsil = 1. _d -9
54          IF ( use_original_hFac ) THEN  
55          IF ( selectVortScheme.EQ.0 ) THEN
56    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           IF (upwindVorticity) THEN       &               )
66             IF (upwindVort3) THEN
67            IF (vBarXY.GT.0.) THEN            IF (vBarXY.GT.0.) THEN
68             vort3u=omega3(I,J)*r_hFacZ(i,j)             vort3u=omega3(i,j)*r_hFacZ(i,j)
69            ELSE            ELSE
70             vort3u=omega3(I,J+1)*r_hFacZ(i,j+1)             vort3u=omega3(i,j+1)*r_hFacZ(i,j+1)
71            ENDIF            ENDIF
72           ELSE           ELSE
73             vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)             vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j)
74       &                +omega3(i,j+1)*r_hFacZ(i,j+1))       &                +omega3(i,j+1)*r_hFacZ(i,j+1))
75           ENDIF           ENDIF
76          ELSEIF ( SadournyCoriolis ) THEN           uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
77           Zm=0.5*(       &                              * _maskW(i,j,k,bi,bj)
78       &       vFld( i , j )*dxG( i , j ,bi,bj)*hFacS( i , j ,k,bi,bj)          ENDDO
79       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacS(i-1, j ,k,bi,bj) )         ENDDO
80           Zp=0.5*(  
81       &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacS( i ,j+1,k,bi,bj)        ELSEIF ( selectVortScheme.EQ.1 ) THEN
82       &      +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)
83           IF (upwindVorticity) THEN  
84            IF ( (Zm+Zp) .GT.0.) THEN         DO J=1-Oly,sNy+Oly-1
85             vort3u=Zm*r_hFacZ(i, j )*omega3(i, j )          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
93              IF (vBarXY.GT.0.) THEN
94               vort3u=omega3(i,j)
95            ELSE            ELSE
96             vort3u=Zp*r_hFacZ(i,j+1)*omega3(i,j+1)             vort3u=omega3(i,j+1)
97            ENDIF            ENDIF
98           ELSE           ELSE
99            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 )  
100           ENDIF           ENDIF
101           vBarXY=1.           uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
102          ELSE       &                              * _maskW(i,j,k,bi,bj)
103  c--      test a different formulation (relatively to hFac)          ENDDO
104           vBarXY=0.5*(         ENDDO
105       &       vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i,j)  
106       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i,j)        ELSEIF ( selectVortScheme.EQ.2 ) THEN
107       &      +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)
108       &      +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1)  
109       &              )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) )         DO j=1-Oly,sNy+Oly-1
110           IF (upwindVorticity) THEN          DO i=2-Olx,sNx+Olx
111            IF (vBarXY.GT.0.) THEN           vBarXm=0.5*(
112             vort3u=omega3(i,j)       &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
113         &      +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            ELSE
121             vort3u=omega3(i,j+1)             vort3u=vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
122            ENDIF            ENDIF
123           ELSE           ELSE
124             vort3u=0.5*(omega3(i,j)+omega3(i,j+1))             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           ENDIF
128          ENDIF           uCoriolisTerm(i,j)= +vort3u*recip_dxC(i,j,bi,bj)
129         &                              * _maskW(i,j,k,bi,bj)
130            ENDDO
131           ENDDO
132    
133          IF (useJamartMomAdv)        ELSE
134       &   vBarXY = vBarXY * 4. _d 0 * hFacW(i,j,k,bi,bj)          WRITE(msgBuf,'(A,I5,A)')
135       &   * MAX( epsil, hFacS( i , j ,k,bi,bj)+hFacS(i-1, j ,k,bi,bj)       &   'MOM_VI_U_CORIOLIS: selectVortScheme=', selectVortScheme,
136       &                +hFacS( j ,i+1,k,bi,bj)+hFacS(i-1,j+1,k,bi,bj) )       &   ' not implemented'
137            CALL PRINT_ERROR( msgBuf, myThid )
138          uCoriolisTerm(i,j)=          STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS'
139       &   +vort3u*vBarXY*recip_dxC(i,j,bi,bj)*_maskW(i,j,k,bi,bj)  
140  cph *note* put these comments after end of continued line        ENDIF
141  cph        to ensure TAMC compatibility  
142  C high order vorticity advection term        IF ( useJamartMomAdv ) THEN
143  c    &   ...         DO j=1-Oly,sNy+Oly-1
144  C linear Coriolis term          DO i=2-Olx,sNx+Olx
145  c    &   +0.5*(fCoriG(I,J,bi,bj)+fCoriG(I,J+1,bi,bj))*vBarXY           uCoriolisTerm(i,j) = uCoriolisTerm(i,j)
146  C full nonlinear Coriolis term       &           * 4. _d 0 * _hFacW(i,j,k,bi,bj)
147  c    &   +0.5*(omega3(I,J)+omega3(I,J+1))*vBarXY       &           / MAX( epsil,
148  C correct energy conserving form of Coriolis term       &                  (_hFacS(i, j ,k,bi,bj)+_hFacS(i-1, j ,k,bi,bj))
149  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))
150  c    &          fCori(I-1,J,bi,bj)*vBarY(I-1,J,K,bi,bj)  )       &                )
151  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  
152         ENDDO         ENDDO
153        ENDDO        ENDIF
154    
155        RETURN        RETURN
156        END        END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22