/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_ad_moc/mon_vort3.F
ViewVC logotype

Annotation of /MITgcm_contrib/heimbach/OpenAD/code_ad_moc/mon_vort3.F

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


Revision 1.1 - (hide annotations) (download)
Mon Sep 15 22:16:01 2008 UTC (16 years, 10 months ago) by utke
Branch: MAIN
CVS Tags: HEAD
keep this version around

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

  ViewVC Help
Powered by ViewVC 1.1.22