/[MITgcm]/MITgcm/pkg/monitor/mon_vort3.F
ViewVC logotype

Annotation of /MITgcm/pkg/monitor/mon_vort3.F

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


Revision 1.10 - (hide annotations) (download)
Wed Apr 27 14:10:06 2005 UTC (19 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57t_post, checkpoint57o_post, checkpoint57m_post, checkpoint57s_post, checkpoint57k_post, checkpoint57i_post, checkpoint57v_post, checkpoint57r_post, checkpoint57h_done, checkpoint57n_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, checkpoint57j_post, checkpoint57h_pre, checkpoint57l_post, checkpoint57h_post
Changes since 1.9: +1 -2 lines
include ${PKG}_OPTIONS.h (if it exists) instead of PACKAGES_CONFIG.h + CPP_OPTIONS.h

1 jmc 1.10 C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_vort3.F,v 1.9 2005/01/27 16:38:22 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4 adcroft 1.2 #include "MONITOR_OPTIONS.h"
5 jmc 1.1
6 edhill 1.8 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7     CBOP
8     C !ROUTINE: MON_VORT3
9    
10     C !INTERFACE:
11 jmc 1.1 SUBROUTINE MON_VORT3(
12 edhill 1.8 I myIter, myThid )
13    
14     C !DESCRIPTION:
15     C Calculates stats for Vorticity (z-component).
16    
17     C !USES:
18 jmc 1.1 IMPLICIT NONE
19     #include "SIZE.h"
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "DYNVARS.h"
23     #include "MONITOR.h"
24     #include "GRID.h"
25 adcroft 1.6 #ifdef ALLOW_EXCH2
26 afe 1.3 #include "W2_EXCH2_TOPOLOGY.h"
27     #include "W2_EXCH2_PARAMS.h"
28 adcroft 1.6 #endif /* ALLOW_EXCH2 */
29 jmc 1.1
30 edhill 1.8 C !INPUT PARAMETERS:
31 jmc 1.1 INTEGER myIter, myThid
32 edhill 1.8 CEOP
33 jmc 1.1
34 edhill 1.8 C !LOCAL VARIABLES:
35 jmc 1.1 INTEGER bi,bj,i,j,k
36     INTEGER iMax,jMax
37 jmc 1.9 _RL theVol, theArea, tmpVal, tmpAre, tmpVol
38 jmc 1.1 _RL theMin, theMax, theMean, theVar, volMean, volVar, theSD
39 jmc 1.9 _RL areaTile, volTile, sumTile, sqsTile, vSumTile, vSqsTile
40 jmc 1.1 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41     _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42     _RL AZcorner
43     #ifdef MONITOR_TEST_HFACZ
44     _RL tmpFac
45     _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46     #endif
47 jmc 1.5 LOGICAL northWestCorner, northEastCorner
48     LOGICAL southWestCorner, southEastCorner
49     INTEGER myTile, iG
50 jmc 1.1
51     theMin = 1. _d 20
52     theMax =-1. _d 20
53     theArea= 0. _d 0
54     theMean= 0. _d 0
55     theVar = 0. _d 0
56     theVol = 0. _d 0
57     volMean= 0. _d 0
58     volVar = 0. _d 0
59     theSD = 0. _d 0
60     AZcorner = 1. _d 0
61    
62     DO bj=myByLo(myThid),myByHi(myThid)
63     DO bi=myBxLo(myThid),myBxHi(myThid)
64 jmc 1.9 areaTile= 0. _d 0
65     volTile = 0. _d 0
66     sumTile = 0. _d 0
67     sqsTile = 0. _d 0
68     vSumTile= 0. _d 0
69     vSqsTile= 0. _d 0
70 jmc 1.1 #ifdef MONITOR_TEST_HFACZ
71     tmpFac = 0.
72     IF( implicDiv2Dflow.GT.0 .AND. abEps.GT.-0.5 )
73     & tmpFac = (0.5+abEps) / implicDiv2Dflow
74     DO j=1-Oly,sNy+Oly
75     DO i=1-Olx,sNx+Olx
76     etaFld(i,j) = etaH(i,j,bi,bj)
77     & + tmpFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
78     ENDDO
79     ENDDO
80     #endif
81     DO k=1,Nr
82    
83     iMax = sNx
84     jMax = sNy
85     DO j=1,sNy
86     DO i=1,sNx
87     #ifdef MONITOR_TEST_HFACZ
88     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89     C- Test various definitions of hFacZ (for 1 layer, flat bottom ocean):
90     c hFacZ(i,j) = 1. +
91     c & 0.25 _d 0*( etaFld(i-1,j-1)
92     c & + etaFld( i ,j-1)
93     c & + etaFld(i-1, j )
94     c & + etaFld( i , j )
95     c & )*recip_drF(k)
96     c hFacZ(i,j) = 1. +
97     c & 0.25 _d 0*( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
98     c & + etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
99     c & + etaFld(i-1, j )*rA(i-1, j ,bi,bj)
100     c & + etaFld( i , j )*rA( i , j ,bi,bj)
101     c & )*recip_drF(k)*recip_rAz(i,j,bi,bj)
102     hFacZ(i,j) = 1. + 0.125 _d 0*
103     & ( ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
104     & +etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
105     & )*recip_rAw(i,j-1,bi,bj)
106     & + ( etaFld(i-1, j )*rA(i-1, j ,bi,bj)
107     & +etaFld( i , j )*rA( i , j ,bi,bj)
108     & )*recip_rAw(i, j ,bi,bj)
109     & + ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
110     & +etaFld(i-1, j )*rA(i-1, j ,bi,bj)
111     & )*recip_rAs(i-1,j,bi,bj)
112     & + ( etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
113     & + etaFld( i , j )*rA( i , j ,bi,bj)
114     & )*recip_rAs( i ,j,bi,bj)
115     & )*recip_drF(k)
116     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
117     #else
118     C- Standard definition of hFac at vorticity point:
119     hFacZ(i,j) =
120     & 0.25 _d 0*( _hFacW(i,j-1,k,bi,bj)
121     & + _hFacW(i, j ,k,bi,bj)
122     & + _hFacS(i-1,j,k,bi,bj)
123     & + _hFacS( i ,j,k,bi,bj)
124     & )
125     #endif /* MONITOR_TEST_HFACZ */
126     vort3(i,j) = recip_rAz(i,j,bi,bj)*(
127     & vVel( i ,j,k,bi,bj)*dyC( i ,j,bi,bj)
128     & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
129     & -uVel(i, j ,k,bi,bj)*dxC(i, j ,bi,bj)
130     & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
131     & )
132     ENDDO
133     ENDDO
134    
135     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
136     C Special stuff for Cubed Sphere:
137     IF (useCubedSphereExchange) THEN
138     c AZcorner = 0.75 _d 0
139     iMax = sNx+1
140     jMax = sNy+1
141     DO i=1,iMax
142     hFacZ(i,jMax)=0.
143     vort3(i,jMax)=0.
144     ENDDO
145     DO j=1,jMax
146     hFacZ(iMax,j)=0.
147     vort3(iMax,j)=0.
148     ENDDO
149 jmc 1.5
150     southWestCorner = .TRUE.
151     southEastCorner = .TRUE.
152     northWestCorner = .TRUE.
153     northEastCorner = .TRUE.
154     iG = bi+(myXGlobalLo-1)/sNx
155 adcroft 1.6 #ifdef ALLOW_EXCH2
156 jmc 1.5 myTile = W2_myTileList(bi)
157     iG = exch2_myFace(myTile)
158     southWestCorner = exch2_isWedge(myTile).EQ.1
159     & .AND. exch2_isSedge(myTile).EQ.1
160     southEastCorner = exch2_isEedge(myTile).EQ.1
161     & .AND. exch2_isSedge(myTile).EQ.1
162     northEastCorner = exch2_isEedge(myTile).EQ.1
163     & .AND. exch2_isNedge(myTile).EQ.1
164     northWestCorner = exch2_isWedge(myTile).EQ.1
165     & .AND. exch2_isNedge(myTile).EQ.1
166 adcroft 1.6 #endif /* ALLOW_EXCH2 */
167 jmc 1.5
168     C-- avoid to count 3 times the same corner:
169     southEastCorner = southEastCorner .AND. iG.EQ.2
170     northWestCorner = northWestCorner .AND. iG.EQ.1
171     northEastCorner = .FALSE.
172    
173 afe 1.3 C-- S.W. corner:
174     IF ( southWestCorner ) THEN
175 jmc 1.1 i=1
176     j=1
177     vort3(i,j)=
178     & +recip_rAz(i,j,bi,bj)/AZcorner*(
179     & vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
180     & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
181     & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
182     & )
183     hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
184     & + _hFacW(i, j ,k,bi,bj)
185     & + _hFacS( i ,j,k,bi,bj)
186     & )/3. _d 0
187 afe 1.3 ENDIF
188     IF ( southEastCorner ) THEN
189 jmc 1.1 C-- S.E. corner:
190     i=iMax
191     j=1
192     vort3(I,J)=
193     & +recip_rAz(I,J,bi,bj)/AZcorner*(
194     & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
195     & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
196     & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
197     & )
198     hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
199     & + _hFacW(i, j ,k,bi,bj)
200     & + _hFacS(i-1,j,k,bi,bj)
201     & )/3. _d 0
202 afe 1.3 ENDIF
203     IF ( northWestCorner ) THEN
204 jmc 1.1 C-- N.W. corner:
205     i=1
206     j=jMax
207     vort3(i,j)=
208     & +recip_rAz(i,j,bi,bj)/AZcorner*(
209     & vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
210     & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
211     & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
212     & )
213     hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
214     & + _hFacW(i, j ,k,bi,bj)
215     & + _hFacS( i ,j,k,bi,bj)
216     & )/3. _d 0
217 afe 1.3 ENDIF
218     IF ( northEastCorner ) THEN
219 jmc 1.1 C-- N.E. corner:
220     i=iMax
221     j=jMax
222     vort3(i,j)=
223     & +recip_rAz(i,j,bi,bj)/AZcorner*(
224     & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
225     & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
226     & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
227     & )
228     hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
229     & + _hFacW(i, j ,k,bi,bj)
230     & + _hFacS(i-1,j,k,bi,bj)
231     & )/3. _d 0
232 afe 1.3 ENDIF
233 jmc 1.1 ENDIF
234    
235     C- Special stuff for North & South Poles, LatLon grid
236     IF ( usingSphericalPolarGrid ) THEN
237     IF (yG(1,sNy+1,bi,bj).EQ.90.) THEN
238     jMax = sNy+1
239     j = jMax
240     DO i=1,sNx
241     vort3(i,j) = 0.
242     vort3(1,j) = vort3(1,j)
243     & + uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
244     hFacZ(i,j) = 0.
245     #ifndef MONITOR_TEST_HFACZ
246     hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j-1,k,bi,bj)
247     ENDDO
248     #else
249     hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j-1)
250     ENDDO
251     hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
252     #endif
253     hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
254     vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
255     ENDIF
256     IF (yG(1,1,bi,bj).EQ.-90.) THEN
257     j = 1
258     DO i=1,sNx
259     vort3(i,j) = 0.
260     vort3(1,j) = vort3(1,j)
261     & - uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
262     hFacZ(i,j) = 0.
263     #ifndef MONITOR_TEST_HFACZ
264     hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j,k,bi,bj)
265     ENDDO
266     #else
267     hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j)
268     ENDDO
269     hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
270     #endif
271     hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
272     vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
273     ENDIF
274     ENDIF
275    
276     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
277    
278     DO J=1,jMax
279     DO I=1,iMax
280     IF (hFacZ(i,j).GT.0. _d 0) THEN
281     tmpVal = vort3(i,j)
282     tmpAre = rAz(i,j,bi,bj)*drF(k)
283     tmpVol = rAz(i,j,bi,bj)*drF(k)*hFacZ(i,j)
284 jmc 1.9 areaTile = areaTile + tmpAre
285 jmc 1.1 C- min,max of relative vorticity ("r")
286     theMin = MIN(theMin,tmpVal)
287     theMax = MAX(theMax,tmpVal)
288     C- average & std.dev of absolute vorticity ("a")
289     tmpVal = tmpVal + fCoriG(i,j,bi,bj)
290 jmc 1.9 sumTile = sumTile + tmpAre*tmpVal
291     sqsTile = sqsTile + tmpAre*tmpVal*tmpVal
292 jmc 1.1 C- average & std.dev of potential vorticity ("p")
293     tmpVal = tmpVal / hFacZ(i,j)
294 jmc 1.9 volTile = volTile + tmpVol
295     vSumTile = vSumTile + tmpVol*tmpVal
296     vSqsTile = vSqsTile + tmpVol*tmpVal*tmpVal
297 jmc 1.1 ENDIF
298     ENDDO
299     ENDDO
300    
301     ENDDO
302 jmc 1.9 theArea= theArea + areaTile
303     theVol = theVol + volTile
304     theMean= theMean + sumTile
305     theVar = theVar + sqsTile
306     volMean= volMean + vSumTile
307     volVar = volVar + vSqsTile
308 jmc 1.1 ENDDO
309     ENDDO
310    
311     theMin = -theMin
312     _GLOBAL_MAX_R8(theMin, myThid)
313     _GLOBAL_MAX_R8(theMax, myThid)
314     _GLOBAL_SUM_R8(theArea,myThid)
315     _GLOBAL_SUM_R8(theVol, myThid)
316     _GLOBAL_SUM_R8(theMean,myThid)
317     _GLOBAL_SUM_R8(theVar, myThid)
318     _GLOBAL_SUM_R8(volMean,myThid)
319     _GLOBAL_SUM_R8(volVar ,myThid)
320     theMin = -theMin
321     IF (theArea.GT.0.) THEN
322     theMean= theMean/theArea
323     theVar = theVar /theArea
324     theVar = theVar - theMean*theMean
325     c IF (theVar.GT.0.) theSD = SQRT(theVar)
326     IF (theVar.GT.0.) theVar = SQRT(theVar)
327     ENDIF
328     IF (theVol.GT.0.) THEN
329     volMean= volMean/theVol
330     volVar = volVar /theVol
331     volVar = volVar - volMean*volMean
332     IF (volVar.GT.0.) theSD = SQRT(volVar)
333     ENDIF
334    
335     C- Print stats for (relative/absolute) Vorticity AND Pot.Vort.
336     CALL MON_SET_PREF('vort',myThid)
337     CALL MON_OUT_RL(mon_string_none,theMin, '_r_min', myThid)
338     CALL MON_OUT_RL(mon_string_none,theMax, '_r_max', myThid)
339     CALL MON_OUT_RL(mon_string_none,theMean,'_a_mean', myThid)
340     CALL MON_OUT_RL(mon_string_none,theVar, '_a_sd', myThid)
341     CALL MON_OUT_RL(mon_string_none,volMean,'_p_mean', myThid)
342     CALL MON_OUT_RL(mon_string_none,theSD, '_p_sd', myThid)
343     c CALL MON_OUT_RL(mon_string_none,theVol,mon_foot_vol,myThid)
344    
345     RETURN
346     END

  ViewVC Help
Powered by ViewVC 1.1.22