1 |
C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_relvort3.F,v 1.8 2009/05/12 19:56:36 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "MOM_COMMON_OPTIONS.h" |
5 |
#undef CALC_CS_CORNER_EXTENDED |
6 |
|
7 |
SUBROUTINE MOM_CALC_RELVORT3( |
8 |
I bi,bj,k, |
9 |
I uFld, vFld, hFacZ, |
10 |
O vort3, |
11 |
I myThid) |
12 |
IMPLICIT NONE |
13 |
C *==========================================================* |
14 |
C | S/R MOM_CALC_RELVORT3 |
15 |
C *==========================================================* |
16 |
C *==========================================================* |
17 |
|
18 |
C == Global variables == |
19 |
#include "SIZE.h" |
20 |
#include "EEPARAMS.h" |
21 |
#include "PARAMS.h" |
22 |
#include "GRID.h" |
23 |
#ifdef ALLOW_EXCH2 |
24 |
#include "W2_EXCH2_SIZE.h" |
25 |
#include "W2_EXCH2_TOPOLOGY.h" |
26 |
#endif /* ALLOW_EXCH2 */ |
27 |
C == Routine arguments == |
28 |
C myThid - Instance number for this innvocation of CALC_MOM_RHS |
29 |
INTEGER bi,bj,k |
30 |
_RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
31 |
_RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
32 |
_RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
33 |
_RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
34 |
INTEGER myThid |
35 |
|
36 |
C == Local variables == |
37 |
INTEGER i,j |
38 |
LOGICAL northWestCorner, northEastCorner, |
39 |
& southWestCorner, southEastCorner |
40 |
INTEGER myFace |
41 |
#ifdef ALLOW_EXCH2 |
42 |
INTEGER myTile |
43 |
#endif /* ALLOW_EXCH2 */ |
44 |
|
45 |
#ifdef ALLOW_AUTODIFF_TAMC |
46 |
DO J=1-Oly,sNy+Oly |
47 |
DO I=1-Olx,sNx+Olx |
48 |
vort3(I,J) = 0. _d 0 |
49 |
ENDDO |
50 |
ENDDO |
51 |
#endif |
52 |
|
53 |
DO J=2-Oly,sNy+Oly |
54 |
DO I=2-Olx,sNx+Olx |
55 |
|
56 |
C Horizontal curl of flow field - ignoring lopping factors |
57 |
vort3(I,J)= |
58 |
& recip_rAz(I,J,bi,bj)*( |
59 |
& ( vFld(I,J)*dyC(I,J,bi,bj) |
60 |
& -vFld(I-1,J)*dyC(I-1,J,bi,bj) ) |
61 |
& -( uFld(I,J)*dxC(I,J,bi,bj) |
62 |
& -uFld(I,J-1)*dxC(I,J-1,bi,bj) ) |
63 |
& ) |
64 |
|
65 |
C Horizontal curl of flow field - including lopping factors |
66 |
c IF (hFacZ(i,j).NE.0.) THEN |
67 |
c vort3(I,J)= |
68 |
c & recip_rAz(I,J,bi,bj)*( |
69 |
c & vFld(I,J)*dyc(I,J,bi,bj)*_hFacW(i,j,k,bi,bj) |
70 |
c & -vFld(I-1,J)*dyc(I-1,J,bi,bj)*_hFacW(i-1,j,k,bi,bj) |
71 |
c & -uFld(I,J)*dxc(I,J,bi,bj)*_hFacS(i,j,k,bi,bj) |
72 |
c & +uFld(I,J-1)*dxc(I,J-1,bi,bj)*_hFacS(i,j-1,k,bi,bj) |
73 |
c & ) |
74 |
c & /hFacZ(i,j) |
75 |
c ELSE |
76 |
c vort3(I,J)=0. |
77 |
c ENDIF |
78 |
|
79 |
ENDDO |
80 |
ENDDO |
81 |
|
82 |
C Special stuff for Cubed Sphere |
83 |
IF (useCubedSphereExchange) THEN |
84 |
#ifdef ALLOW_EXCH2 |
85 |
myTile = W2_myTileList(bi,bj) |
86 |
myFace = exch2_myFace(myTile) |
87 |
southWestCorner = exch2_isWedge(myTile).EQ.1 |
88 |
& .AND. exch2_isSedge(myTile).EQ.1 |
89 |
southEastCorner = exch2_isEedge(myTile).EQ.1 |
90 |
& .AND. exch2_isSedge(myTile).EQ.1 |
91 |
northEastCorner = exch2_isEedge(myTile).EQ.1 |
92 |
& .AND. exch2_isNedge(myTile).EQ.1 |
93 |
northWestCorner = exch2_isWedge(myTile).EQ.1 |
94 |
& .AND. exch2_isNedge(myTile).EQ.1 |
95 |
#else |
96 |
myFace = bi |
97 |
southWestCorner = .TRUE. |
98 |
southEastCorner = .TRUE. |
99 |
northWestCorner = .TRUE. |
100 |
northEastCorner = .TRUE. |
101 |
#endif /* ALLOW_EXCH2 */ |
102 |
|
103 |
IF ( southWestCorner ) THEN |
104 |
C U(0,1) D(0,1) U(1,1) TILE |
105 |
C | | |
106 |
C V(-1,1) --- Z(0,1) --- V(0,1) --- Z(1,1) --- V(1,1) --- |
107 |
C | | |
108 |
C U(0,0) D(0,0) U(1,0) D(1,0) |
109 |
C | | |
110 |
C --- V(0,0) --- Z(1,0) --- V(1,0) --- |
111 |
C | |
112 |
C U(1,-1) |
113 |
I=1 |
114 |
J=1 |
115 |
C- to get the same truncation, independent from the face Nb, |
116 |
C do (1+2)+3, and always in the same order (exch3 convention order): |
117 |
vort3(I,J)= |
118 |
& +recip_rAz(I,J,bi,bj)*( |
119 |
& ( vFld(I,J)*dyC(I,J,bi,bj) |
120 |
& -uFld(I,J)*dxC(I,J,bi,bj) ) |
121 |
& + uFld(I,J-1)*dxC(I,J-1,bi,bj) |
122 |
& ) |
123 |
C- the quick way, but do not get the same truncation on the 3 faces: |
124 |
c vort3(I,J)= |
125 |
c & +recip_rAz(I,J,bi,bj)*( |
126 |
c & vFld(I,J)*dyC(I,J,bi,bj) |
127 |
c & -uFld(I,J)*dxC(I,J,bi,bj) |
128 |
c & +uFld(I,J-1)*dxC(I,J-1,bi,bj) |
129 |
c & ) |
130 |
#ifdef CALC_CS_CORNER_EXTENDED |
131 |
vort3(I-1,J)= |
132 |
& recip_rAz(I-1,J,bi,bj)*( |
133 |
& vFld(I-1,J)*dyC(I-1,J,bi,bj) |
134 |
& -vFld(I-2,J)*dyC(I-2,J,bi,bj) |
135 |
& -uFld(I-1,J)*dxC(I-1,J,bi,bj) |
136 |
& +vFld(I+0,J-1)*dyC(I+0,J-1,bi,bj) |
137 |
& ) |
138 |
& *maskS(i-1,j,k,bi,bj)*maskS(i-2,j,k,bi,bj) |
139 |
& *maskW(i-1,j,k,bi,bj)*maskS(i,j-1,k,bi,bj) |
140 |
vort3(I,J-1)=vort3(I-1,J) |
141 |
#endif |
142 |
ENDIF |
143 |
|
144 |
IF ( southEastCorner ) THEN |
145 |
C TILE U(N+1,1) D(N+1,1) U(N+2,1) |
146 |
C | | |
147 |
C V(N,1) --- Z(N+1,1) --- V(N+1,1) --- Z(N+2,1) --- V(N+3,1) --- |
148 |
C | | |
149 |
C D(N,0) U(N+1,0) D(N+1,0) U(N+2,0) |
150 |
C | | |
151 |
C V(N,0) --- Z(N+1,0) --- V(N+1,0) --- |
152 |
C | | |
153 |
C U(N+1,-1) |
154 |
I=sNx+1 |
155 |
J=1 |
156 |
C- to get the same truncation, independent from the face Nb, |
157 |
C (exch3 convention order): |
158 |
IF ( myFace.EQ.2 ) THEN |
159 |
vort3(I,J)= |
160 |
& +recip_rAz(I,J,bi,bj)*( |
161 |
& (-uFld(I,J)*dxC(I,J,bi,bj) |
162 |
& -vFld(I-1,J)*dyC(I-1,J,bi,bj) ) |
163 |
& + uFld(I,J-1)*dxC(I,J-1,bi,bj) |
164 |
& ) |
165 |
ELSEIF ( myFace.EQ.4 ) THEN |
166 |
vort3(I,J)= |
167 |
& +recip_rAz(I,J,bi,bj)*( |
168 |
& (-vFld(I-1,J)*dyC(I-1,J,bi,bj) |
169 |
& +uFld(I,J-1)*dxC(I,J-1,bi,bj) ) |
170 |
& - uFld(I,J)*dxC(I,J,bi,bj) |
171 |
& ) |
172 |
ELSE |
173 |
vort3(I,J)= |
174 |
& +recip_rAz(I,J,bi,bj)*( |
175 |
& (+uFld(I,J-1)*dxC(I,J-1,bi,bj) |
176 |
& -uFld(I,J)*dxC(I,J,bi,bj) ) |
177 |
& - vFld(I-1,J)*dyC(I-1,J,bi,bj) |
178 |
& ) |
179 |
ENDIF |
180 |
C- the quick way, but do not get the same truncation on the 3 faces: |
181 |
c vort3(I,J)= |
182 |
c & +recip_rAz(I,J,bi,bj)*( |
183 |
c & -vFld(I-1,J)*dyC(I-1,J,bi,bj) |
184 |
c & -uFld(I,J)*dxC(I,J,bi,bj) |
185 |
c & +uFld(I,J-1)*dxC(I,J-1,bi,bj) |
186 |
c & ) |
187 |
#ifdef CALC_CS_CORNER_EXTENDED |
188 |
vort3(I+1,J)= |
189 |
& recip_rAz(I+1,J,bi,bj)*( |
190 |
& vFld(I+1,J)*dyC(I+1,J,bi,bj) |
191 |
& -vFld(I-0,J)*dyC(I-0,J,bi,bj) |
192 |
& -uFld(I+1,J)*dxC(I+1,J,bi,bj) |
193 |
& -vFld(I-1,J-1)*dyC(I-1,J-1,bi,bj) |
194 |
& ) |
195 |
& *maskS(i+1,j,k,bi,bj)*maskS(i-0,j,k,bi,bj) |
196 |
& *maskW(i+1,j,k,bi,bj)*maskS(i-1,j-1,k,bi,bj) |
197 |
vort3(I,J-1)=vort3(I+1,J) |
198 |
#endif |
199 |
ENDIF |
200 |
|
201 |
IF ( northWestCorner ) THEN |
202 |
C U(1,N+2) |
203 |
C | |
204 |
C --- V(0,N+1) --- Z(1,N+2) --- V(1,N+2) --- |
205 |
C | | |
206 |
C U(0,N+1) D(0,N+1) U(1,N+1) D(1,N+1) |
207 |
C | | |
208 |
C V(-1,N+1) --- Z(0,N+1) --- V(0,N+1) --- Z(1,N+1) --- V(1,N+1) --- |
209 |
C | | |
210 |
C U(0,N) D(0,N) U(1,N) TILE |
211 |
I=1 |
212 |
J=sNy+1 |
213 |
C- to get the same truncation, independent from the face Nb, |
214 |
C (exch3 convention order): |
215 |
IF ( myFace.EQ.1 ) THEN |
216 |
vort3(I,J)= |
217 |
& +recip_rAz(I,J,bi,bj)*( |
218 |
& (+uFld(I,J-1)*dxC(I,J-1,bi,bj) |
219 |
& +vFld(I,J)*dyC(I,J,bi,bj) ) |
220 |
& -uFld(I,J)*dxC(I,J,bi,bj) |
221 |
& ) |
222 |
ELSEIF ( myFace.EQ.3 ) THEN |
223 |
vort3(I,J)= |
224 |
& +recip_rAz(I,J,bi,bj)*( |
225 |
& (-uFld(I,J)*dxC(I,J,bi,bj) |
226 |
& +uFld(I,J-1)*dxC(I,J-1,bi,bj) ) |
227 |
& + vFld(I,J)*dyC(I,J,bi,bj) |
228 |
& ) |
229 |
ELSE |
230 |
vort3(I,J)= |
231 |
& +recip_rAz(I,J,bi,bj)*( |
232 |
& (+vFld(I,J)*dyC(I,J,bi,bj) |
233 |
& -uFld(I,J)*dxC(I,J,bi,bj) ) |
234 |
& + uFld(I,J-1)*dxC(I,J-1,bi,bj) |
235 |
& ) |
236 |
ENDIF |
237 |
C- the quick way, but do not get the same truncation on the 3 faces: |
238 |
c vort3(I,J)= |
239 |
c & +recip_rAz(I,J,bi,bj)*( |
240 |
c & vFld(I,J)*dyC(I,J,bi,bj) |
241 |
c & -uFld(I,J)*dxC(I,J,bi,bj) |
242 |
c & +uFld(I,J-1)*dxC(I,J-1,bi,bj) |
243 |
c & ) |
244 |
#ifdef CALC_CS_CORNER_EXTENDED |
245 |
vort3(I-1,J)= |
246 |
& recip_rAz(I-1,J,bi,bj)*( |
247 |
& vFld(I-1,J)*dyC(I-1,J,bi,bj) |
248 |
& -vFld(I-2,J)*dyC(I-2,J,bi,bj) |
249 |
& +vFld(I-0,J+1)*dyC(I-0,J+1,bi,bj) |
250 |
& +uFld(I-1,J-1)*dxC(I-1,J-1,bi,bj) |
251 |
& ) |
252 |
& *maskS(i-1,j,k,bi,bj)*maskS(i-2,j,k,bi,bj) |
253 |
& *maskS(i,j+1,k,bi,bj)*maskW(i-1,j-1,k,bi,bj) |
254 |
vort3(I,J+1)=vort3(I-1,J) |
255 |
#endif |
256 |
ENDIF |
257 |
|
258 |
IF ( northEastCorner ) THEN |
259 |
C U(N+1,N+2) |
260 |
C | | |
261 |
C V(N,N+2) --- Z(N+1,N+2) --- V(N+1,N+2) --- |
262 |
C | | |
263 |
C D(N,N+1) U(N+1,N+1) D(N+1,N+1) U(N+2,N+1) |
264 |
C | | |
265 |
C V(N,N+1) --- Z(N+1,N+1) --- V(N+1,N+1) --- Z(N+2,N+1) --- V(N+3,N+1) --- |
266 |
C | | |
267 |
C TILE U(N+1,N) D(N+1,N) U(N+2,N) |
268 |
I=sNx+1 |
269 |
J=sNy+1 |
270 |
C- to get the same truncation, independent from the face Nb: |
271 |
C (exch3 convention order): |
272 |
IF ( MOD(myFace,2).EQ.1 ) THEN |
273 |
vort3(I,J)= |
274 |
& +recip_rAz(I,J,bi,bj)*( |
275 |
& (-uFld(I,J)*dxC(I,J,bi,bj) |
276 |
& -vFld(I-1,J)*dyC(I-1,J,bi,bj) ) |
277 |
& + uFld(I,J-1)*dxC(I,J-1,bi,bj) |
278 |
& ) |
279 |
ELSE |
280 |
vort3(I,J)= |
281 |
& +recip_rAz(I,J,bi,bj)*( |
282 |
& (+uFld(I,J-1)*dxC(I,J-1,bi,bj) |
283 |
& -uFld(I,J)*dxC(I,J,bi,bj) ) |
284 |
& - vFld(I-1,J)*dyC(I-1,J,bi,bj) |
285 |
& ) |
286 |
ENDIF |
287 |
C- the quick way, but do not get the same truncation on the 3 faces: |
288 |
c vort3(I,J)= |
289 |
c & +recip_rAz(I,J,bi,bj)*( |
290 |
c & -vFld(I-1,J)*dyC(I-1,J,bi,bj) |
291 |
c & -uFld(I,J)*dxC(I,J,bi,bj) |
292 |
c & +uFld(I,J-1)*dxC(I,J-1,bi,bj) |
293 |
c & ) |
294 |
#ifdef CALC_CS_CORNER_EXTENDED |
295 |
vort3(I+1,J)= |
296 |
& recip_rAz(I+1,J,bi,bj)*( |
297 |
& vFld(I+1,J)*dyC(I+1,J,bi,bj) |
298 |
& -vFld(I-0,J)*dyC(I-0,J,bi,bj) |
299 |
& -vFld(I-1,J+1)*dyC(I-1,J+1,bi,bj) |
300 |
& +uFld(I+1,J-1)*dxC(I+1,J-1,bi,bj) |
301 |
& ) |
302 |
& *maskS(i+1,j,k,bi,bj)*maskS(i-0,j,k,bi,bj) |
303 |
& *maskS(i-1,j+1,k,bi,bj)*maskW(i+1,j-1,k,bi,bj) |
304 |
vort3(I,J+1)=vort3(I+1,J) |
305 |
#endif |
306 |
ENDIF |
307 |
ENDIF |
308 |
|
309 |
RETURN |
310 |
END |