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

Contents of /MITgcm/pkg/mom_vecinv/mom_vi_u_coriolis_c4.F

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

Revision 1.13 - (show annotations) (download)
Fri May 13 00:56:30 2011 UTC (14 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint63g, checkpoint64, checkpoint65, checkpoint63, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint65o, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.12: +24 -8 lines 
for OBCS: mask gradient of vorticity with interior mask maskInW,S
1 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vi_u_coriolis_c4.F,v 1.12 2009/06/28 01:08:26 jmc Exp $
2 C $Name: $
3
4 #include "MOM_VECINV_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: MOM_VI_U_CORIOLIS_C4
8 C !INTERFACE:
9 SUBROUTINE MOM_VI_U_CORIOLIS_C4(
10 I bi,bj,k,
11 I vFld,omega3,r_hFacZ,
12 O uCoriolisTerm,
13 I myThid)
14 C !DESCRIPTION: \bv
15 C *==========================================================*
16 C | S/R MOM_VI_U_CORIOLIS_C4
17 C |==========================================================*
18 C | o Calculate flux (in Y-dir.) of vorticity at U point
19 C | using 4th order (or 1rst order) interpolation
20 C *==========================================================*
21 C \ev
22
23 C !USES:
24 IMPLICIT NONE
25
26 C == Global variables ==
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29 #include "PARAMS.h"
30 #include "GRID.h"
31 #ifdef ALLOW_EXCH2
32 #include "W2_EXCH2_SIZE.h"
33 #include "W2_EXCH2_TOPOLOGY.h"
34 #endif /* ALLOW_EXCH2 */
35
36 C !INPUT/OUTPUT PARAMETERS:
37 C == Routine arguments ==
38 INTEGER bi,bj,k
39 _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40 _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41 _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42 _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43 INTEGER myThid
44 CEOP
45
46 C == Local variables ==
47 C msgBuf :: Informational/error meesage buffer
48 CHARACTER*(MAX_LEN_MBUF) msgBuf
49 INTEGER i,j
50 _RL vort3r(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51 _RL vBarXY,vort3u,Rjp,Rjm,Rj
52 _RL vBarXm,vBarXp
53
54 LOGICAL northWestCorner, northEastCorner,
55 & southWestCorner, southEastCorner
56 INTEGER myFace
57 #ifdef ALLOW_EXCH2
58 INTEGER myTile
59 #endif /* ALLOW_EXCH2 */
60 _RL oneSixth, oneTwelve
61 LOGICAL fourthVort3
62 C jmc: not sure about these 1/6 & 1/12 factors (should use the same)
63 PARAMETER(oneSixth=1.D0/6.D0 , oneTwelve=1.D0/12.D0)
64 PARAMETER(fourthVort3=.TRUE. )
65
66 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
67
68 DO j=1-Oly,sNy+Oly
69 DO i=1-Olx,sNx+Olx
70 vort3r(i,j) = r_hFacZ(i,j)*omega3(i,j)
71 ENDDO
72 ENDDO
73
74 C-- Special stuff for Cubed Sphere
75 IF ( useCubedSphereExchange.AND.highOrderVorticity ) THEN
76
77 #ifdef ALLOW_EXCH2
78 myTile = W2_myTileList(bi,bj)
79 myFace = exch2_myFace(myTile)
80 southWestCorner = exch2_isWedge(myTile).EQ.1
81 & .AND. exch2_isSedge(myTile).EQ.1
82 southEastCorner = exch2_isEedge(myTile).EQ.1
83 & .AND. exch2_isSedge(myTile).EQ.1
84 northEastCorner = exch2_isEedge(myTile).EQ.1
85 & .AND. exch2_isNedge(myTile).EQ.1
86 northWestCorner = exch2_isWedge(myTile).EQ.1
87 & .AND. exch2_isNedge(myTile).EQ.1
88 #else
89 myFace = bi
90 southWestCorner = .TRUE.
91 southEastCorner = .TRUE.
92 northWestCorner = .TRUE.
93 northEastCorner = .TRUE.
94 #endif /* ALLOW_EXCH2 */
95 IF ( southWestCorner ) THEN
96 i = 1
97 j = 1
98 vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i+1,j) )*0.5 _d 0
99 ENDIF
100 IF ( southEastCorner ) THEN
101 i = sNx+1
102 j = 1
103 vort3r(i,j-1) = ( vort3r(i,j-1) + vort3r(i-1,j) )*0.5 _d 0
104 ENDIF
105 IF ( northWestCorner ) THEN
106 i = 1
107 j = sNy+1
108 vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i+1,j) )*0.5 _d 0
109 ENDIF
110 IF ( northEastCorner ) THEN
111 i = sNx+1
112 j = sNy+1
113 vort3r(i,j+1) = ( vort3r(i,j+1) + vort3r(i-1,j) )*0.5 _d 0
114 ENDIF
115
116 C-- End of special stuff for Cubed Sphere.
117 ENDIF
118
119 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
120
121 IF ( selectVortScheme.EQ.0 ) THEN
122 C-- using Sadourny Enstrophy conserving discretization:
123
124 c DO j=2-Oly,sNy+Oly-2
125 c DO i=2-Olx,sNx+Olx
126 DO j=1,sNy
127 DO i=1,sNx+1
128
129 vBarXY=0.25*(
130 & (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
131 & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj))
132 & +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
133 & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj))
134 & )
135 IF (upwindVorticity) THEN
136 IF (vBarXY.GT.0.) THEN
137 vort3u=vort3r(i,j)
138 ELSE
139 vort3u=vort3r(i,j+1)
140 ENDIF
141 ELSEIF (fourthVort3) THEN
142 #ifdef ALLOW_OBCS
143 Rjp = ( vort3r(i,j+2) - vort3r(i,j+1) )*maskInW(i,j+1,bi,bj)
144 Rjm = ( vort3r(i, j ) - vort3r(i,j-1) )*maskInW(i,j-1,bi,bj)
145 #else
146 Rjp = vort3r(i,j+2) - vort3r(i,j+1)
147 Rjm = vort3r(i, j ) - vort3r(i,j-1)
148 #endif
149 vort3u=0.5*( (vort3r(i,j) + vort3r(i,j+1))
150 & -oneTwelve*(Rjp-Rjm)
151 & )
152 ELSE
153 vort3u=0.5*( vort3r(i,j) + vort3r(i,j+1) )
154 ENDIF
155
156 uCoriolisTerm(i,j) = vort3u*vBarXY*recip_dxC(i,j,bi,bj)
157 & * _maskW(i,j,k,bi,bj)
158
159 ENDDO
160 ENDDO
161
162 ELSEIF ( selectVortScheme.EQ.2 ) THEN
163 C-- using Energy conserving discretization:
164
165 c DO j=2-Oly,sNy+Oly-2
166 c DO i=2-Olx,sNx+Olx
167 DO j=1,sNy
168 DO i=1,sNx+1
169
170 vBarXm=0.5*(
171 & vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
172 & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
173 vBarXp=0.5*(
174 & vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
175 & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
176 IF (upwindVorticity) THEN
177 IF ( (vBarXm+vBarXp) .GT.0.) THEN
178 vort3u=vBarXm*vort3r(i, j )
179 ELSE
180 vort3u=vBarXp*vort3r(i,j+1)
181 ENDIF
182 ELSEIF (fourthVort3) THEN
183 #ifdef ALLOW_OBCS
184 Rjp = ( vort3r(i,j+2) - vort3r(i,j+1) )*maskInW(i,j+1,bi,bj)
185 Rjm = ( vort3r(i, j ) - vort3r(i,j-1) )*maskInW(i,j-1,bi,bj)
186 #else
187 Rjp = vort3r(i,j+2) - vort3r(i,j+1)
188 Rjm = vort3r(i, j ) - vort3r(i,j-1)
189 #endif
190 Rj = vort3r(i,j+1) - vort3r(i, j )
191 Rjp = vort3r(i,j+1) -oneSixth*( Rjp+Rj )
192 Rjm = vort3r(i, j ) +oneSixth*( Rj+Rjm )
193 c Rjp = vort3r(i,j+1) -oneSixth*( vort3r(i,j+2)-vort3r(i, j ) )
194 c Rjm = vort3r(i, j ) +oneSixth*( vort3r(i,j+1)-vort3r(i,j-1) )
195 vort3u=0.5*( vBarXm*Rjm + vBarXp*Rjp )
196 ELSE
197 vort3u=0.5*( vBarXm*vort3r(i, j ) + vBarXp*vort3r(i,j+1) )
198 ENDIF
199
200 uCoriolisTerm(i,j) = vort3u*recip_dxC(i,j,bi,bj)
201 & * _maskW(i,j,k,bi,bj)
202
203 ENDDO
204 ENDDO
205
206 ELSE
207 WRITE(msgBuf,'(A,I5,A)')
208 & 'MOM_VI_U_CORIOLIS_C4: selectVortScheme=', selectVortScheme,
209 & ' not implemented'
210 CALL PRINT_ERROR( msgBuf, myThid )
211 STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS_C4'
212
213 ENDIF
214
215 RETURN
216 END
  ViewVC Help
Powered by ViewVC 1.1.22