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

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

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


Revision 1.1 - (show annotations) (download)
Fri May 2 22:27:01 2003 UTC (21 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint50d_pre, checkpoint50c_post, checkpoint50d_post
monitor global statistic of vorticity. note: contribution of singular points
 (cube-sphere corner & N+S pole of Lat-Lon grid) not yet tested.

1 C $Header: $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #undef MONITOR_TEST_HFACZ
6
7 SUBROUTINE MON_VORT3(
8 I myIter, myThid )
9 C *==========================================================*
10 C | SUBROUTINE MON_VORT3
11 C | o Calculates stats for Vorticity (z-component)
12 C *==========================================================*
13 IMPLICIT NONE
14
15 C === Global data ===
16 #include "SIZE.h"
17 #include "EEPARAMS.h"
18 #include "PARAMS.h"
19 #include "DYNVARS.h"
20 #include "MONITOR.h"
21 #include "GRID.h"
22
23 C === Routine arguments ===
24 INTEGER myIter, myThid
25
26 C === Local variables ====
27 INTEGER bi,bj,i,j,k
28 INTEGER iMax,jMax
29 _RL numPnts, theVol, theArea, tmpVal, tmpAre, tmpVol
30 _RL theMin, theMax, theMean, theVar, volMean, volVar, theSD
31 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33 _RL AZcorner
34 #ifdef MONITOR_TEST_HFACZ
35 _RL tmpFac
36 _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
37 #endif
38
39 theMin = 1. _d 20
40 theMax =-1. _d 20
41 theArea= 0. _d 0
42 theMean= 0. _d 0
43 theVar = 0. _d 0
44 theVol = 0. _d 0
45 volMean= 0. _d 0
46 volVar = 0. _d 0
47 theSD = 0. _d 0
48 AZcorner = 1. _d 0
49
50 DO bj=myByLo(myThid),myByHi(myThid)
51 DO bi=myBxLo(myThid),myBxHi(myThid)
52 #ifdef MONITOR_TEST_HFACZ
53 tmpFac = 0.
54 IF( implicDiv2Dflow.GT.0 .AND. abEps.GT.-0.5 )
55 & tmpFac = (0.5+abEps) / implicDiv2Dflow
56 DO j=1-Oly,sNy+Oly
57 DO i=1-Olx,sNx+Olx
58 etaFld(i,j) = etaH(i,j,bi,bj)
59 & + tmpFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
60 ENDDO
61 ENDDO
62 #endif
63 DO k=1,Nr
64
65 iMax = sNx
66 jMax = sNy
67 DO j=1,sNy
68 DO i=1,sNx
69 #ifdef MONITOR_TEST_HFACZ
70 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
71 C- Test various definitions of hFacZ (for 1 layer, flat bottom ocean):
72 c hFacZ(i,j) = 1. +
73 c & 0.25 _d 0*( etaFld(i-1,j-1)
74 c & + etaFld( i ,j-1)
75 c & + etaFld(i-1, j )
76 c & + etaFld( i , j )
77 c & )*recip_drF(k)
78 c hFacZ(i,j) = 1. +
79 c & 0.25 _d 0*( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
80 c & + etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
81 c & + etaFld(i-1, j )*rA(i-1, j ,bi,bj)
82 c & + etaFld( i , j )*rA( i , j ,bi,bj)
83 c & )*recip_drF(k)*recip_rAz(i,j,bi,bj)
84 hFacZ(i,j) = 1. + 0.125 _d 0*
85 & ( ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
86 & +etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
87 & )*recip_rAw(i,j-1,bi,bj)
88 & + ( etaFld(i-1, j )*rA(i-1, j ,bi,bj)
89 & +etaFld( i , j )*rA( i , j ,bi,bj)
90 & )*recip_rAw(i, j ,bi,bj)
91 & + ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
92 & +etaFld(i-1, j )*rA(i-1, j ,bi,bj)
93 & )*recip_rAs(i-1,j,bi,bj)
94 & + ( etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
95 & + etaFld( i , j )*rA( i , j ,bi,bj)
96 & )*recip_rAs( i ,j,bi,bj)
97 & )*recip_drF(k)
98 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
99 #else
100 C- Standard definition of hFac at vorticity point:
101 hFacZ(i,j) =
102 & 0.25 _d 0*( _hFacW(i,j-1,k,bi,bj)
103 & + _hFacW(i, j ,k,bi,bj)
104 & + _hFacS(i-1,j,k,bi,bj)
105 & + _hFacS( i ,j,k,bi,bj)
106 & )
107 #endif /* MONITOR_TEST_HFACZ */
108 vort3(i,j) = recip_rAz(i,j,bi,bj)*(
109 & vVel( i ,j,k,bi,bj)*dyC( i ,j,bi,bj)
110 & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
111 & -uVel(i, j ,k,bi,bj)*dxC(i, j ,bi,bj)
112 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
113 & )
114 ENDDO
115 ENDDO
116
117 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
118 C Special stuff for Cubed Sphere:
119 IF (useCubedSphereExchange) THEN
120 c AZcorner = 0.75 _d 0
121 iMax = sNx+1
122 jMax = sNy+1
123 DO i=1,iMax
124 hFacZ(i,jMax)=0.
125 vort3(i,jMax)=0.
126 ENDDO
127 DO j=1,jMax
128 hFacZ(iMax,j)=0.
129 vort3(iMax,j)=0.
130 ENDDO
131 C-- S.W. corner:
132 i=1
133 j=1
134 vort3(i,j)=
135 & +recip_rAz(i,j,bi,bj)/AZcorner*(
136 & vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
137 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
138 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
139 & )
140 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
141 & + _hFacW(i, j ,k,bi,bj)
142 & + _hFacS( i ,j,k,bi,bj)
143 & )/3. _d 0
144 C-- S.E. corner:
145 i=iMax
146 j=1
147 vort3(I,J)=
148 & +recip_rAz(I,J,bi,bj)/AZcorner*(
149 & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
150 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
151 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
152 & )
153 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
154 & + _hFacW(i, j ,k,bi,bj)
155 & + _hFacS(i-1,j,k,bi,bj)
156 & )/3. _d 0
157 C-- N.W. corner:
158 i=1
159 j=jMax
160 vort3(i,j)=
161 & +recip_rAz(i,j,bi,bj)/AZcorner*(
162 & vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
163 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
164 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
165 & )
166 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
167 & + _hFacW(i, j ,k,bi,bj)
168 & + _hFacS( i ,j,k,bi,bj)
169 & )/3. _d 0
170 C-- N.E. corner:
171 i=iMax
172 j=jMax
173 vort3(i,j)=
174 & +recip_rAz(i,j,bi,bj)/AZcorner*(
175 & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
176 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
177 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
178 & )
179 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
180 & + _hFacW(i, j ,k,bi,bj)
181 & + _hFacS(i-1,j,k,bi,bj)
182 & )/3. _d 0
183 ENDIF
184
185 C- Special stuff for North & South Poles, LatLon grid
186 IF ( usingSphericalPolarGrid ) THEN
187 IF (yG(1,sNy+1,bi,bj).EQ.90.) THEN
188 jMax = sNy+1
189 j = jMax
190 DO i=1,sNx
191 vort3(i,j) = 0.
192 vort3(1,j) = vort3(1,j)
193 & + uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
194 hFacZ(i,j) = 0.
195 #ifndef MONITOR_TEST_HFACZ
196 hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j-1,k,bi,bj)
197 ENDDO
198 #else
199 hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j-1)
200 ENDDO
201 hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
202 #endif
203 hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
204 vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
205 ENDIF
206 IF (yG(1,1,bi,bj).EQ.-90.) THEN
207 j = 1
208 DO i=1,sNx
209 vort3(i,j) = 0.
210 vort3(1,j) = vort3(1,j)
211 & - uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
212 hFacZ(i,j) = 0.
213 #ifndef MONITOR_TEST_HFACZ
214 hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j,k,bi,bj)
215 ENDDO
216 #else
217 hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j)
218 ENDDO
219 hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
220 #endif
221 hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
222 vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
223 ENDIF
224 ENDIF
225
226 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
227
228 DO J=1,jMax
229 DO I=1,iMax
230 IF (hFacZ(i,j).GT.0. _d 0) THEN
231 tmpVal = vort3(i,j)
232 tmpAre = rAz(i,j,bi,bj)*drF(k)
233 tmpVol = rAz(i,j,bi,bj)*drF(k)*hFacZ(i,j)
234 theArea= theArea + tmpAre
235 C- min,max of relative vorticity ("r")
236 theMin = MIN(theMin,tmpVal)
237 theMax = MAX(theMax,tmpVal)
238 C- average & std.dev of absolute vorticity ("a")
239 tmpVal = tmpVal + fCoriG(i,j,bi,bj)
240 theMean= theMean+ tmpAre*tmpVal
241 theVar = theVar + tmpAre*tmpVal*tmpVal
242 C- average & std.dev of potential vorticity ("p")
243 tmpVal = tmpVal / hFacZ(i,j)
244 theVol = theVol + tmpVol
245 volMean= volMean+ tmpVol*tmpVal
246 volVar = volVar + tmpVol*tmpVal*tmpVal
247 ENDIF
248 ENDDO
249 ENDDO
250
251 ENDDO
252 ENDDO
253 ENDDO
254
255 theMin = -theMin
256 _GLOBAL_MAX_R8(theMin, myThid)
257 _GLOBAL_MAX_R8(theMax, myThid)
258 _GLOBAL_SUM_R8(theArea,myThid)
259 _GLOBAL_SUM_R8(theVol, myThid)
260 _GLOBAL_SUM_R8(theMean,myThid)
261 _GLOBAL_SUM_R8(theVar, myThid)
262 _GLOBAL_SUM_R8(volMean,myThid)
263 _GLOBAL_SUM_R8(volVar ,myThid)
264 theMin = -theMin
265 IF (theArea.GT.0.) THEN
266 theMean= theMean/theArea
267 theVar = theVar /theArea
268 theVar = theVar - theMean*theMean
269 c IF (theVar.GT.0.) theSD = SQRT(theVar)
270 IF (theVar.GT.0.) theVar = SQRT(theVar)
271 ENDIF
272 IF (theVol.GT.0.) THEN
273 volMean= volMean/theVol
274 volVar = volVar /theVol
275 volVar = volVar - volMean*volMean
276 IF (volVar.GT.0.) theSD = SQRT(volVar)
277 ENDIF
278
279 C- Print stats for (relative/absolute) Vorticity AND Pot.Vort.
280 CALL MON_SET_PREF('vort',myThid)
281 CALL MON_OUT_RL(mon_string_none,theMin, '_r_min', myThid)
282 CALL MON_OUT_RL(mon_string_none,theMax, '_r_max', myThid)
283 CALL MON_OUT_RL(mon_string_none,theMean,'_a_mean', myThid)
284 CALL MON_OUT_RL(mon_string_none,theVar, '_a_sd', myThid)
285 CALL MON_OUT_RL(mon_string_none,volMean,'_p_mean', myThid)
286 CALL MON_OUT_RL(mon_string_none,theSD, '_p_sd', myThid)
287 c CALL MON_OUT_RL(mon_string_none,theVol,mon_foot_vol,myThid)
288
289 RETURN
290 END

  ViewVC Help
Powered by ViewVC 1.1.22