/[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.2 - (show annotations) (download)
Tue May 13 18:18:05 2003 UTC (21 years, 1 month ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint51e_post, checkpoint51k_post, checkpoint52e_pre, checkpoint51o_pre, checkpoint50e_post, checkpoint52e_post, checkpoint51n_pre, checkpoint51l_post, checkpoint51q_post, checkpoint51j_post, checkpoint50g_post, hrcube_1, branch-netcdf, checkpoint52d_pre, checkpoint52b_pre, checkpoint51a_post, checkpoint51c_post, checkpoint51f_pre, checkpoint51, checkpoint51o_post, checkpoint51p_post, checkpoint52a_pre, checkpoint51i_post, checkpoint52, checkpoint51f_post, checkpoint52d_post, checkpoint51b_post, checkpoint51r_post, checkpoint51b_pre, checkpoint52a_post, checkpoint52b_post, checkpoint52f_post, branchpoint-genmake2, checkpoint50h_post, checkpoint52c_post, checkpoint51h_pre, checkpoint51l_pre, checkpoint51g_post, ecco_c52_e35, checkpoint50f_post, checkpoint50f_pre, checkpoint51d_post, checkpoint52i_post, checkpoint51t_post, checkpoint51n_post, checkpoint51i_pre, checkpoint50i_post, checkpoint52i_pre, checkpoint51u_post, checkpoint52h_pre, checkpoint50e_pre, checkpoint52f_pre, checkpoint51m_post, checkpoint51s_post
Branch point for: branch-nonh, branch-genmake2, tg2-branch, checkpoint51n_branch, ecco-branch, netcdf-sm0
Changes since 1.1: +2 -3 lines 
 o split mon_set.F into mon_set_iounit.F and mon_set_pref.F
 o replaced ref's to CPP_OPTIONS with MONITOR_OPTIONS
 o added new s/r monitor_solution.F that checks that model state
   and if unlikely lets the model die cleanly
   - this is to reduce the number of hanging processes we encounter
     if the model dies due to FPEs
1 C $Header: /u/gcmpack/models/MITgcmUV/pkg/monitor/mon_vort3.F,v 1.1 2003/05/02 22:27:01 jmc Exp $
2 C $Name: $
3
4 #include "MONITOR_OPTIONS.h"
5
6 SUBROUTINE MON_VORT3(
7 I myIter, myThid )
8 C *==========================================================*
9 C | SUBROUTINE MON_VORT3
10 C | o Calculates stats for Vorticity (z-component)
11 C *==========================================================*
12 IMPLICIT NONE
13
14 C === Global data ===
15 #include "SIZE.h"
16 #include "EEPARAMS.h"
17 #include "PARAMS.h"
18 #include "DYNVARS.h"
19 #include "MONITOR.h"
20 #include "GRID.h"
21
22 C === Routine arguments ===
23 INTEGER myIter, myThid
24
25 C === Local variables ====
26 INTEGER bi,bj,i,j,k
27 INTEGER iMax,jMax
28 _RL numPnts, theVol, theArea, tmpVal, tmpAre, tmpVol
29 _RL theMin, theMax, theMean, theVar, volMean, volVar, theSD
30 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32 _RL AZcorner
33 #ifdef MONITOR_TEST_HFACZ
34 _RL tmpFac
35 _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36 #endif
37
38 theMin = 1. _d 20
39 theMax =-1. _d 20
40 theArea= 0. _d 0
41 theMean= 0. _d 0
42 theVar = 0. _d 0
43 theVol = 0. _d 0
44 volMean= 0. _d 0
45 volVar = 0. _d 0
46 theSD = 0. _d 0
47 AZcorner = 1. _d 0
48
49 DO bj=myByLo(myThid),myByHi(myThid)
50 DO bi=myBxLo(myThid),myBxHi(myThid)
51 #ifdef MONITOR_TEST_HFACZ
52 tmpFac = 0.
53 IF( implicDiv2Dflow.GT.0 .AND. abEps.GT.-0.5 )
54 & tmpFac = (0.5+abEps) / implicDiv2Dflow
55 DO j=1-Oly,sNy+Oly
56 DO i=1-Olx,sNx+Olx
57 etaFld(i,j) = etaH(i,j,bi,bj)
58 & + tmpFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
59 ENDDO
60 ENDDO
61 #endif
62 DO k=1,Nr
63
64 iMax = sNx
65 jMax = sNy
66 DO j=1,sNy
67 DO i=1,sNx
68 #ifdef MONITOR_TEST_HFACZ
69 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
70 C- Test various definitions of hFacZ (for 1 layer, flat bottom ocean):
71 c hFacZ(i,j) = 1. +
72 c & 0.25 _d 0*( etaFld(i-1,j-1)
73 c & + etaFld( i ,j-1)
74 c & + etaFld(i-1, j )
75 c & + etaFld( i , j )
76 c & )*recip_drF(k)
77 c hFacZ(i,j) = 1. +
78 c & 0.25 _d 0*( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
79 c & + etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
80 c & + etaFld(i-1, j )*rA(i-1, j ,bi,bj)
81 c & + etaFld( i , j )*rA( i , j ,bi,bj)
82 c & )*recip_drF(k)*recip_rAz(i,j,bi,bj)
83 hFacZ(i,j) = 1. + 0.125 _d 0*
84 & ( ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
85 & +etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
86 & )*recip_rAw(i,j-1,bi,bj)
87 & + ( etaFld(i-1, j )*rA(i-1, j ,bi,bj)
88 & +etaFld( i , j )*rA( i , j ,bi,bj)
89 & )*recip_rAw(i, j ,bi,bj)
90 & + ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
91 & +etaFld(i-1, j )*rA(i-1, j ,bi,bj)
92 & )*recip_rAs(i-1,j,bi,bj)
93 & + ( etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
94 & + etaFld( i , j )*rA( i , j ,bi,bj)
95 & )*recip_rAs( i ,j,bi,bj)
96 & )*recip_drF(k)
97 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
98 #else
99 C- Standard definition of hFac at vorticity point:
100 hFacZ(i,j) =
101 & 0.25 _d 0*( _hFacW(i,j-1,k,bi,bj)
102 & + _hFacW(i, j ,k,bi,bj)
103 & + _hFacS(i-1,j,k,bi,bj)
104 & + _hFacS( i ,j,k,bi,bj)
105 & )
106 #endif /* MONITOR_TEST_HFACZ */
107 vort3(i,j) = recip_rAz(i,j,bi,bj)*(
108 & vVel( i ,j,k,bi,bj)*dyC( i ,j,bi,bj)
109 & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
110 & -uVel(i, j ,k,bi,bj)*dxC(i, j ,bi,bj)
111 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
112 & )
113 ENDDO
114 ENDDO
115
116 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
117 C Special stuff for Cubed Sphere:
118 IF (useCubedSphereExchange) THEN
119 c AZcorner = 0.75 _d 0
120 iMax = sNx+1
121 jMax = sNy+1
122 DO i=1,iMax
123 hFacZ(i,jMax)=0.
124 vort3(i,jMax)=0.
125 ENDDO
126 DO j=1,jMax
127 hFacZ(iMax,j)=0.
128 vort3(iMax,j)=0.
129 ENDDO
130 C-- S.W. corner:
131 i=1
132 j=1
133 vort3(i,j)=
134 & +recip_rAz(i,j,bi,bj)/AZcorner*(
135 & vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
136 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
137 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
138 & )
139 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
140 & + _hFacW(i, j ,k,bi,bj)
141 & + _hFacS( i ,j,k,bi,bj)
142 & )/3. _d 0
143 C-- S.E. corner:
144 i=iMax
145 j=1
146 vort3(I,J)=
147 & +recip_rAz(I,J,bi,bj)/AZcorner*(
148 & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
149 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
150 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
151 & )
152 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
153 & + _hFacW(i, j ,k,bi,bj)
154 & + _hFacS(i-1,j,k,bi,bj)
155 & )/3. _d 0
156 C-- N.W. corner:
157 i=1
158 j=jMax
159 vort3(i,j)=
160 & +recip_rAz(i,j,bi,bj)/AZcorner*(
161 & vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
162 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
163 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
164 & )
165 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
166 & + _hFacW(i, j ,k,bi,bj)
167 & + _hFacS( i ,j,k,bi,bj)
168 & )/3. _d 0
169 C-- N.E. corner:
170 i=iMax
171 j=jMax
172 vort3(i,j)=
173 & +recip_rAz(i,j,bi,bj)/AZcorner*(
174 & -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
175 & -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
176 & +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
177 & )
178 hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
179 & + _hFacW(i, j ,k,bi,bj)
180 & + _hFacS(i-1,j,k,bi,bj)
181 & )/3. _d 0
182 ENDIF
183
184 C- Special stuff for North & South Poles, LatLon grid
185 IF ( usingSphericalPolarGrid ) THEN
186 IF (yG(1,sNy+1,bi,bj).EQ.90.) THEN
187 jMax = sNy+1
188 j = jMax
189 DO i=1,sNx
190 vort3(i,j) = 0.
191 vort3(1,j) = vort3(1,j)
192 & + uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
193 hFacZ(i,j) = 0.
194 #ifndef MONITOR_TEST_HFACZ
195 hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j-1,k,bi,bj)
196 ENDDO
197 #else
198 hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j-1)
199 ENDDO
200 hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
201 #endif
202 hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
203 vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
204 ENDIF
205 IF (yG(1,1,bi,bj).EQ.-90.) THEN
206 j = 1
207 DO i=1,sNx
208 vort3(i,j) = 0.
209 vort3(1,j) = vort3(1,j)
210 & - uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
211 hFacZ(i,j) = 0.
212 #ifndef MONITOR_TEST_HFACZ
213 hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j,k,bi,bj)
214 ENDDO
215 #else
216 hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j)
217 ENDDO
218 hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
219 #endif
220 hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
221 vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
222 ENDIF
223 ENDIF
224
225 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
226
227 DO J=1,jMax
228 DO I=1,iMax
229 IF (hFacZ(i,j).GT.0. _d 0) THEN
230 tmpVal = vort3(i,j)
231 tmpAre = rAz(i,j,bi,bj)*drF(k)
232 tmpVol = rAz(i,j,bi,bj)*drF(k)*hFacZ(i,j)
233 theArea= theArea + tmpAre
234 C- min,max of relative vorticity ("r")
235 theMin = MIN(theMin,tmpVal)
236 theMax = MAX(theMax,tmpVal)
237 C- average & std.dev of absolute vorticity ("a")
238 tmpVal = tmpVal + fCoriG(i,j,bi,bj)
239 theMean= theMean+ tmpAre*tmpVal
240 theVar = theVar + tmpAre*tmpVal*tmpVal
241 C- average & std.dev of potential vorticity ("p")
242 tmpVal = tmpVal / hFacZ(i,j)
243 theVol = theVol + tmpVol
244 volMean= volMean+ tmpVol*tmpVal
245 volVar = volVar + tmpVol*tmpVal*tmpVal
246 ENDIF
247 ENDDO
248 ENDDO
249
250 ENDDO
251 ENDDO
252 ENDDO
253
254 theMin = -theMin
255 _GLOBAL_MAX_R8(theMin, myThid)
256 _GLOBAL_MAX_R8(theMax, myThid)
257 _GLOBAL_SUM_R8(theArea,myThid)
258 _GLOBAL_SUM_R8(theVol, myThid)
259 _GLOBAL_SUM_R8(theMean,myThid)
260 _GLOBAL_SUM_R8(theVar, myThid)
261 _GLOBAL_SUM_R8(volMean,myThid)
262 _GLOBAL_SUM_R8(volVar ,myThid)
263 theMin = -theMin
264 IF (theArea.GT.0.) THEN
265 theMean= theMean/theArea
266 theVar = theVar /theArea
267 theVar = theVar - theMean*theMean
268 c IF (theVar.GT.0.) theSD = SQRT(theVar)
269 IF (theVar.GT.0.) theVar = SQRT(theVar)
270 ENDIF
271 IF (theVol.GT.0.) THEN
272 volMean= volMean/theVol
273 volVar = volVar /theVol
274 volVar = volVar - volMean*volMean
275 IF (volVar.GT.0.) theSD = SQRT(volVar)
276 ENDIF
277
278 C- Print stats for (relative/absolute) Vorticity AND Pot.Vort.
279 CALL MON_SET_PREF('vort',myThid)
280 CALL MON_OUT_RL(mon_string_none,theMin, '_r_min', myThid)
281 CALL MON_OUT_RL(mon_string_none,theMax, '_r_max', myThid)
282 CALL MON_OUT_RL(mon_string_none,theMean,'_a_mean', myThid)
283 CALL MON_OUT_RL(mon_string_none,theVar, '_a_sd', myThid)
284 CALL MON_OUT_RL(mon_string_none,volMean,'_p_mean', myThid)
285 CALL MON_OUT_RL(mon_string_none,theSD, '_p_sd', myThid)
286 c CALL MON_OUT_RL(mon_string_none,theVol,mon_foot_vol,myThid)
287
288 RETURN
289 END
  ViewVC Help
Powered by ViewVC 1.1.22