/[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.11 by jmc, Sun Mar 30 21:46:48 2008 UTC revision 1.12 by jmc, Wed Apr 2 21:18:45 2008 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "MOM_VECINV_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    
6    CBOP
7    C     !ROUTINE: MOM_VI_U_CORIOLIS
8    C     !INTERFACE:
9        SUBROUTINE MOM_VI_U_CORIOLIS(        SUBROUTINE MOM_VI_U_CORIOLIS(
10       I        bi,bj,k,       I                     bi, bj, k,
11       I        vFld,omega3,hFacZ,r_hFacZ,       I                     vFld, omega3, hFacZ, r_hFacZ,
12       O        uCoriolisTerm,       O                     uCoriolisTerm,
13       I        myThid)       I                     myThid )
14        IMPLICIT NONE  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          _RS     epsil
51          PARAMETER( upwindVort3 = .FALSE. )
52    
53        epsil = 1. _d -9        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          IF ( use_original_hFac ) THEN  
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-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))
63       &     +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,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 ( selectVortScheme.EQ.2 ) 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            ELSE           vBarXY= 0.5*(
87             vort3u=Zp*r_hFacZ(i,j+1)*omega3(i,j+1)       &      (vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i, j )
88            ENDIF       &      +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i, j ))
          ELSE  
           Zm=Zm*r_hFacZ(i, j )*omega3(i, j )  
           Zp=Zp*r_hFacZ(i,j+1)*omega3(i,j+1)  
           vort3u=0.5*( Zm + Zp )  
          ENDIF  
          vBarXY=1.  
         ELSE  
 c--      test a different formulation (relatively to hFac)  
          vBarXY=0.5*(  
      &      (vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i,j)  
      &      +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)       &     +(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))       &      +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) )       &               )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) )
92           IF (upwindVorticity) THEN           IF (upwindVort3) THEN
93            IF (vBarXY.GT.0.) THEN            IF (vBarXY.GT.0.) THEN
94             vort3u=omega3(i,j)             vort3u=omega3(i,j)
95            ELSE            ELSE
# Line 92  c--      test a different formulation (r Line 98  c--      test a different formulation (r
98           ELSE           ELSE
99             vort3u=0.5*(omega3(i,j)+omega3(i,j+1))             vort3u=0.5*(omega3(i,j)+omega3(i,j+1))
100           ENDIF           ENDIF
101          ENDIF           uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
102         &                              * _maskW(i,j,k,bi,bj)
103            ENDDO
104           ENDDO
105    
106          ELSEIF ( selectVortScheme.EQ.2 ) THEN
107    C--   using energy conserving scheme (used by Sadourny in JAS 75 paper)
108    
109           DO j=1-Oly,sNy+Oly-1
110            DO i=2-Olx,sNx+Olx
111             vBarXm=0.5*(
112         &       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
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          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(i,j+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.11  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22