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

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

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


Revision 1.24 - (hide annotations) (download)
Sun Feb 17 04:07:30 2013 UTC (11 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, HEAD
Changes since 1.23: +121 -41 lines
Change diagnostics of total angular momentum:
- do not add Eta contribution if using rigid-lid;
- correction for Adams-Bashforth done in zonal-wind part (instead of
  in Eta contribution);
- horizontal discretisation: compute Zonal wind contribution separately for
  each component (instead of from cell centered averaged) and, regarding
  Eta contribution, use the 4 grid-cell corner averaged value of Omega*r^2;
In the linear dynamics case (momAdvection=F and linear FreeSurf),
 get exact conservation of AM if using vectorInvariant momentum.

1 jmc 1.24 C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_ke.F,v 1.23 2012/12/26 23:57:46 jmc Exp $
2 adcroft 1.1 C $Name: $
3    
4 adcroft 1.11 #include "MONITOR_OPTIONS.h"
5 adcroft 1.1
6 edhill 1.12 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7     CBOP
8     C !ROUTINE: MON_KE
9    
10     C !INTERFACE:
11 adcroft 1.1 SUBROUTINE MON_KE(
12 edhill 1.12 I myIter, myThid )
13    
14     C !DESCRIPTION:
15 jmc 1.22 C Calculates stats for Kinetic Energy, (barotropic) Potential Energy
16     C and total Angular Momentum
17 edhill 1.12
18     C !USES:
19 adcroft 1.1 IMPLICIT NONE
20     #include "SIZE.h"
21     #include "EEPARAMS.h"
22 jmc 1.16 #include "PARAMS.h"
23 adcroft 1.1 #include "DYNVARS.h"
24 cnh 1.5 #include "MONITOR.h"
25 jmc 1.8 #include "GRID.h"
26 jmc 1.10 #include "SURFACE.h"
27 adcroft 1.1
28 edhill 1.12 C !INPUT PARAMETERS:
29 jmc 1.10 INTEGER myIter, myThid
30 edhill 1.12 CEOP
31 adcroft 1.1
32 edhill 1.12 C !LOCAL VARIABLES:
33 jmc 1.22 INTEGER bi, bj
34     INTEGER i,j,k
35     INTEGER ks, kp1
36 jmc 1.19 _RL numPnts,theVol,tmpVal, mskp1, msk_1
37 jmc 1.24 _RL abFac1, abFac2, R_drK, cosLat
38 jmc 1.10 _RL theMax,theMean,theVolMean,potEnMean
39 jmc 1.24 _RL totAMu, totAMs
40 jmc 1.16 _RL tileMean(nSx,nSy)
41     _RL tileVlAv(nSx,nSy)
42     _RL tilePEav(nSx,nSy)
43     _RL tileVol (nSx,nSy)
44 jmc 1.22 _RL tileAMu (nSx,nSy)
45     _RL tileAMs (nSx,nSy)
46 jmc 1.24 _RL tmpFld(1:sNx,1:sNy)
47     _RS cos2LatG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48 jmc 1.18 #ifdef ALLOW_NONHYDROSTATIC
49     _RL tmpWke
50     #endif
51 jmc 1.24 #ifdef ALLOW_ADAMSBASHFORTH_3
52     INTEGER m1, m2
53     #endif
54 adcroft 1.1
55 jmc 1.21 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
56    
57 jmc 1.10 numPnts=0.
58     theVol=0.
59 adcroft 1.1 theMax=0.
60     theMean=0.
61 jmc 1.8 theVolMean=0.
62 jmc 1.10 potEnMean =0.
63 adcroft 1.1
64     DO bj=myByLo(myThid),myByHi(myThid)
65     DO bi=myBxLo(myThid),myBxHi(myThid)
66 jmc 1.16 tileVol(bi,bj) = 0. _d 0
67     tileMean(bi,bj) = 0. _d 0
68     tileVlAv(bi,bj) = 0. _d 0
69     tilePEav(bi,bj) = 0. _d 0
70     DO k=1,Nr
71 jmc 1.18 kp1 = MIN(k+1,Nr)
72     mskp1 = 1.
73 jmc 1.19 IF ( k.GE.Nr ) mskp1 = 0.
74     C- Note: Present NH implementation does not account for D.w/dt at k=1.
75     C Consequently, wVel(k=1) does not contribute to NH KE (msk_1=0).
76     msk_1 = 1.
77 jmc 1.22 IF ( k.EQ.1 .AND. selectNHfreeSurf.LE.0 ) msk_1 = 0.
78 jmc 1.16 DO j=1,sNy
79     DO i=1,sNx
80     tileVol(bi,bj) = tileVol(bi,bj)
81     & + rA(i,j,bi,bj)*deepFac2C(k)
82     & *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
83 jmc 1.21 & *maskInC(i,j,bi,bj)
84 jmc 1.9
85     C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F)
86 jmc 1.16 c tmpVal=0.25*( uVel( i , j ,k,bi,bj)*uVel( i , j ,k,bi,bj)
87     c & +uVel(i+1, j ,k,bi,bj)*uVel(i+1, j ,k,bi,bj)
88     c & +vVel( i , j ,k,bi,bj)*vVel( i , j ,k,bi,bj)
89     c & +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) )
90     c tileVlAv(bi,bj) = tileVlAv(bi,bj)
91     c & +tmpVal*rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
92 jmc 1.9
93     C- Energy conservative form (like in pkg/mom_fluxform/mom_calc_ke.F)
94     C this is the safe way to check the energy conservation
95     C with no assumption on how grid spacing & area are defined.
96     tmpVal=0.25*(
97     & uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj)
98 heimbach 1.15 & *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)
99 jmc 1.9 & +uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
100 heimbach 1.15 & *dyG(i+1,j,bi,bj)*dxC(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)
101 jmc 1.9 & +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj)
102 heimbach 1.15 & *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)
103 jmc 1.9 & +vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
104 heimbach 1.15 & *dxG(i,j+1,bi,bj)*dyC(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)
105 jmc 1.21 & )*maskInC(i,j,bi,bj)
106 jmc 1.16 tileVlAv(bi,bj) = tileVlAv(bi,bj)
107     & + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
108 heimbach 1.15 tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
109 jmc 1.9
110 jmc 1.18 #ifdef ALLOW_NONHYDROSTATIC
111     IF ( nonHydrostatic ) THEN
112     tmpWke = 0.25*
113 jmc 1.19 & ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)*msk_1
114 jmc 1.18 & *deepFac2F( k )*rhoFacF( k )
115     & +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1
116     & *deepFac2F(kp1)*rhoFacF(kp1)
117 jmc 1.21 & )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
118 jmc 1.18 tileVlAv(bi,bj) = tileVlAv(bi,bj)
119     & + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj)
120     tmpVal = tmpVal
121     & + tmpWke*recip_deepFac2C(k)*recip_rhoFacC(k)
122     ENDIF
123     #endif
124    
125 jmc 1.16 theMax=MAX(theMax,tmpVal)
126 adcroft 1.1 IF (tmpVal.NE.0.) THEN
127 jmc 1.16 tileMean(bi,bj)=tileMean(bi,bj)+tmpVal
128 jmc 1.10 numPnts=numPnts+1.
129 adcroft 1.1 ENDIF
130 jmc 1.9
131 adcroft 1.1 ENDDO
132     ENDDO
133     ENDDO
134 jmc 1.10 C- Potential Energy (external mode):
135 jmc 1.16 DO j=1,sNy
136     DO i=1,sNx
137 jmc 1.10 tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj)
138     & *etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
139     C- jmc: if geoid not flat (phi0surf), needs to add this term.
140     C not sure for atmos/ocean in P ; or atmos. loading in ocean-Z
141     tmpVal = tmpVal
142     & + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)
143 jmc 1.16 tilePEav(bi,bj) = tilePEav(bi,bj)
144 jmc 1.20 & + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)
145     & *maskInC(i,j,bi,bj)
146 jmc 1.10 c tmpVal = etaN(i,j,bi,bj)
147     c & + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
148 jmc 1.16 c tilePEav(bi,bj) = tilePEav(bi,bj)
149 jmc 1.10 c & + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal
150 jmc 1.20 c & *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
151 jmc 1.10 ENDDO
152     ENDDO
153     C- end bi,bj loops
154 adcroft 1.1 ENDDO
155     ENDDO
156 jmc 1.17 _GLOBAL_SUM_RL(numPnts,myThid)
157     _GLOBAL_MAX_RL(theMax,myThid)
158 jmc 1.16 CALL GLOBAL_SUM_TILE_RL( tileMean, theMean , myThid )
159     CALL GLOBAL_SUM_TILE_RL( tileVol , theVol , myThid )
160     CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )
161     CALL GLOBAL_SUM_TILE_RL( tilePEav, potEnMean , myThid )
162 jmc 1.10 IF (numPnts.NE.0.) theMean=theMean/numPnts
163     IF (theVol.NE.0.) THEN
164     theVolMean=theVolMean/theVol
165     potEnMean = potEnMean/theVol
166     ENDIF
167    
168 jmc 1.22 C-- Compute total angular momentum
169     IF ( mon_output_AM ) THEN
170     DO bj=myByLo(myThid),myByHi(myThid)
171     DO bi=myBxLo(myThid),myBxHi(myThid)
172 jmc 1.24 C- Calculate contribution from zonal velocity
173     abFac1 = 0. _d 0
174     abFac2 = 0. _d 0
175     #ifdef ALLOW_ADAMSBASHFORTH_3
176     m1 = 1 + mod(myIter+1,2)
177     m2 = 1 + mod( myIter ,2)
178     IF ( myIter.GE.2 ) abFac2 = beta_AB
179     IF ( myIter.GE.1 ) abFac1 = -( alph_AB + abFac2 )
180     #else
181     IF ( myIter.GE.1 ) abFac1 = -( 0.5 _d 0 + abEps )
182     #endif
183     C- contribution from uVel component: 1rst integrate vertically
184 jmc 1.22 DO j=1,sNy
185     DO i=1,sNx
186 jmc 1.24 tmpFld(i,j) = 0. _d 0
187 jmc 1.22 ENDDO
188     ENDDO
189     DO k=1,Nr
190 jmc 1.24 R_drK = rSphere*deepFacC(k)*deepFac2C(k)
191     & *rhoFacC(k)*drF(k)
192 jmc 1.22 DO j=1,sNy
193     DO i=1,sNx
194 jmc 1.24 #ifdef ALLOW_ADAMSBASHFORTH_3
195     tmpVal = abFac1*guNm(i,j,k,bi,bj,m1)
196     & + abFac2*guNm(i,j,k,bi,bj,m2)
197     #else
198     tmpVal = abFac1*guNm1(i,j,k,bi,bj)
199     #endif
200     tmpVal = tmpVal*deltaTMom + uVel(i,j,k,bi,bj)
201     tmpFld(i,j) = tmpFld(i,j)
202     & + R_drK*tmpVal*_hFacW(i,j,k,bi,bj)
203 jmc 1.22 ENDDO
204     ENDDO
205     ENDDO
206 jmc 1.24 C- and then integrate horizontally over this tile
207     DO j=1,sNy
208     DO i=1,sNx
209     cosLat = COS( deg2rad*
210     & ( yG(i,j,bi,bj) + yG(i,j+1,bi,bj) )*halfRL )
211     tmpFld(i,j) = tmpFld(i,j)*u2zonDir(i,j,bi,bj)
212     & *cosLat*rAw(i,j,bi,bj)
213     & *maskInW(i,j,bi,bj)
214     ENDDO
215     ENDDO
216     tileAMu(bi,bj) = 0. _d 0
217     DO j=1,sNy
218     DO i=1,sNx
219     tileAMu(bi,bj) = tileAMu(bi,bj) + tmpFld(i,j)
220     ENDDO
221     ENDDO
222     C- contribution from vVel component: 1rst integrate vertically
223     DO j=1,sNy
224     DO i=1,sNx
225     tmpFld(i,j) = 0. _d 0
226     ENDDO
227     ENDDO
228     DO k=1,Nr
229     R_drK = rSphere*deepFacC(k)*deepFac2C(k)
230     & *rhoFacC(k)*drF(k)
231     DO j=1,sNy
232     DO i=1,sNx
233 jmc 1.23 #ifdef ALLOW_ADAMSBASHFORTH_3
234 jmc 1.24 tmpVal = abFac1*gvNm(i,j,k,bi,bj,m1)
235     & + abFac2*gvNm(i,j,k,bi,bj,m2)
236 jmc 1.23 #else
237 jmc 1.24 tmpVal = abFac1*gvNm1(i,j,k,bi,bj)
238 jmc 1.23 #endif
239 jmc 1.24 tmpVal = tmpVal*deltaTMom + vVel(i,j,k,bi,bj)
240     tmpFld(i,j) = tmpFld(i,j)
241     & + R_drK*tmpVal*_hFacS(i,j,k,bi,bj)
242     ENDDO
243     ENDDO
244     ENDDO
245     C- and then integrate horizontally over this tile
246     DO j=1,sNy
247     DO i=1,sNx
248     cosLat = COS( deg2rad*
249     & ( yG(i,j,bi,bj) + yG(i+1,j,bi,bj) )*halfRL )
250     tmpFld(i,j) = tmpFld(i,j)*v2zonDir(i,j,bi,bj)
251     & *cosLat*rAs(i,j,bi,bj)
252     & *maskInS(i,j,bi,bj)
253     ENDDO
254     ENDDO
255     DO j=1,sNy
256     DO i=1,sNx
257     tileAMu(bi,bj) = tileAMu(bi,bj) + tmpFld(i,j)
258     ENDDO
259     ENDDO
260     C- Calculate contribution from mass distribution anomaly (i.e., free-surface)
261     IF ( exactConserv ) THEN
262 jmc 1.23 DO j=1,sNy
263     DO i=1,sNx
264     #ifdef EXACT_CONSERV
265 jmc 1.24 tmpFld(i,j) = etaHnm1(i,j,bi,bj)
266     #else
267     tmpFld(i,j) = 0.
268 jmc 1.23 #endif
269     ENDDO
270     ENDDO
271     ELSE
272 jmc 1.22 DO j=1,sNy
273     DO i=1,sNx
274 jmc 1.24 tmpFld(i,j) = etaN(i,j,bi,bj)
275 jmc 1.22 ENDDO
276     ENDDO
277 jmc 1.23 ENDIF
278 jmc 1.24 C- calculate angular momentum from mass-distribution anomaly
279     C using square of radial distance (averaged @ center point)
280     DO j=1-OLy,sNy+OLy
281     DO i=1-OLx,sNx+OLx
282     cosLat = COS( deg2rad*yG(i,j,bi,bj) )
283     cos2LatG(i,j) = cosLat*cosLat
284     ENDDO
285     ENDDO
286     DO j=1,sNy
287     DO i=1,sNx
288     tmpFld(i,j) = tmpFld(i,j)
289     & *omega*rSphere*rSphere
290     & *( ( cos2LatG(i,j) + cos2LatG(i+1,j+1) )
291     & + ( cos2LatG(i+1,j) + cos2LatG(i,j+1) )
292     & )*0.25 _d 0
293     ENDDO
294     ENDDO
295     DO j=1,sNy
296     DO i=1,sNx
297     ks = kSurfC(i,j,bi,bj)
298     tmpFld(i,j) = tmpFld(i,j)
299     & *maskInC(i,j,bi,bj)*deepFac2F(ks)
300     & *rA(i,j,bi,bj)*deepFac2F(ks)*rhoFacF(ks)
301     ENDDO
302     ENDDO
303     tileAMs(bi,bj) = 0. _d 0
304     DO j=1,sNy
305     DO i=1,sNx
306     tileAMs(bi,bj) = tileAMs(bi,bj) + tmpFld(i,j)
307     ENDDO
308     ENDDO
309 jmc 1.22 C- end bi,bj loops
310     ENDDO
311     ENDDO
312     CALL GLOBAL_SUM_TILE_RL( tileAMu , totAMu, myThid )
313     CALL GLOBAL_SUM_TILE_RL( tileAMs , totAMs, myThid )
314    
315 jmc 1.23 C-- Print stats for total Angular Momentum (per unit area, in kg/s):
316 jmc 1.22 CALL MON_SET_PREF('am',myThid)
317     totAMu = totAMu*rUnit2mass
318     totAMs = totAMs*rUnit2mass
319     IF ( globalArea.GT.0. ) totAMu = totAMu/globalArea
320     IF ( globalArea.GT.0. ) totAMs = totAMs/globalArea
321     CALL MON_OUT_RL( mon_string_none, totAMs,
322     & '_eta_mean', myThid )
323     CALL MON_OUT_RL( mon_string_none, totAMu,
324     & '_uZo_mean', myThid )
325 jmc 1.24 totAMu = totAMu + freeSurfFac*totAMs
326 jmc 1.22 CALL MON_OUT_RL( mon_string_none, totAMu,
327     & '_tot_mean', myThid )
328    
329     ENDIF
330    
331 jmc 1.10 C-- Print stats for (barotropic) Potential Energy:
332     CALL MON_SET_PREF('pe_b',myThid)
333     CALL MON_OUT_RL(mon_string_none,potEnMean,
334     & mon_foot_mean,myThid)
335 adcroft 1.1
336 jmc 1.10 C-- Print stats for KE
337     CALL MON_SET_PREF('ke',myThid)
338 jmc 1.8 CALL MON_OUT_RL(mon_string_none,theMax,mon_foot_max,myThid)
339 jmc 1.10 c CALL MON_OUT_RL(mon_string_none,theMean,mon_foot_mean,myThid)
340 jmc 1.8 CALL MON_OUT_RL(mon_string_none,theVolMean,
341 jmc 1.10 & mon_foot_mean,myThid)
342 jmc 1.9 CALL MON_OUT_RL(mon_string_none,theVol,
343     & mon_foot_vol,myThid)
344 adcroft 1.1
345     RETURN
346     END

  ViewVC Help
Powered by ViewVC 1.1.22