/[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.8 - (hide annotations) (download)
Sat Apr 3 21:17:10 2004 UTC (20 years, 2 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint52n_post, checkpoint53d_post, checkpoint54a_pre, checkpoint55c_post, checkpoint54e_post, checkpoint54a_post, checkpoint53c_post, checkpoint55d_pre, checkpoint57b_post, checkpoint57c_pre, checkpoint55j_post, checkpoint56b_post, checkpoint55h_post, checkpoint53b_post, checkpoint54b_post, checkpoint53b_pre, checkpoint55b_post, checkpoint54d_post, checkpoint56c_post, checkpoint52m_post, checkpoint55, checkpoint53a_post, checkpoint57a_post, checkpoint54, checkpoint54f_post, checkpoint55g_post, checkpoint55f_post, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, checkpoint53, checkpoint53g_post, checkpoint55e_post, checkpoint53f_post, checkpoint55a_post, checkpoint53d_pre, checkpoint54c_post, checkpoint56a_post, checkpoint55d_post
Changes since 1.7: +15 -10 lines
 o removing duplicate code in mon_out.F in preparation for MNC output
 o convert all monitor files to protex-style comments

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

  ViewVC Help
Powered by ViewVC 1.1.22