/[MITgcm]/MITgcm/pkg/shap_filt/shap_filt_relvort3.F
ViewVC logotype

Diff of /MITgcm/pkg/shap_filt/shap_filt_relvort3.F

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

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