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

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

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

Revision 1.18 - (show annotations) (download)
Sat Nov 28 20:52:54 2009 UTC (14 years, 6 months ago) by jmc
Branch: MAIN
Changes since 1.17: +24 -3 lines 
add wVel contribution if nonHydrostatic
1 C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_ke.F,v 1.17 2009/04/28 18:16:53 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_KE
9
10 C !INTERFACE:
11 SUBROUTINE MON_KE(
12 I myIter, myThid )
13
14 C !DESCRIPTION:
15 C Calculates stats for Kinetic energy
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 "MONITOR.h"
24 #include "GRID.h"
25 #include "SURFACE.h"
26
27 C !INPUT PARAMETERS:
28 INTEGER myIter, myThid
29 CEOP
30
31 C !LOCAL VARIABLES:
32 INTEGER bi,bj,i,j,k,kp1
33 _RL numPnts,theVol,tmpVal,mskp1
34 _RL theMax,theMean,theVolMean,potEnMean
35 _RL tileMean(nSx,nSy)
36 _RL tileVlAv(nSx,nSy)
37 _RL tilePEav(nSx,nSy)
38 _RL tileVol (nSx,nSy)
39 #ifdef ALLOW_NONHYDROSTATIC
40 _RL tmpWke
41 #endif
42
43 numPnts=0.
44 theVol=0.
45 theMax=0.
46 theMean=0.
47 theVolMean=0.
48 potEnMean =0.
49
50 DO bj=myByLo(myThid),myByHi(myThid)
51 DO bi=myBxLo(myThid),myBxHi(myThid)
52 tileVol(bi,bj) = 0. _d 0
53 tileMean(bi,bj) = 0. _d 0
54 tileVlAv(bi,bj) = 0. _d 0
55 tilePEav(bi,bj) = 0. _d 0
56 DO k=1,Nr
57 kp1 = MIN(k+1,Nr)
58 mskp1 = 1.
59 IF ( k.GE.Nr) mskp1 = 0.
60 DO j=1,sNy
61 DO i=1,sNx
62 tileVol(bi,bj) = tileVol(bi,bj)
63 & + rA(i,j,bi,bj)*deepFac2C(k)
64 & *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
65
66 C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F)
67 c tmpVal=0.25*( uVel( i , j ,k,bi,bj)*uVel( i , j ,k,bi,bj)
68 c & +uVel(i+1, j ,k,bi,bj)*uVel(i+1, j ,k,bi,bj)
69 c & +vVel( i , j ,k,bi,bj)*vVel( i , j ,k,bi,bj)
70 c & +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) )
71 c tileVlAv(bi,bj) = tileVlAv(bi,bj)
72 c & +tmpVal*rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
73
74 C- Energy conservative form (like in pkg/mom_fluxform/mom_calc_ke.F)
75 C this is the safe way to check the energy conservation
76 C with no assumption on how grid spacing & area are defined.
77 tmpVal=0.25*(
78 & uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj)
79 & *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)
80 & +uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
81 & *dyG(i+1,j,bi,bj)*dxC(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)
82 & +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj)
83 & *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)
84 & +vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
85 & *dxG(i,j+1,bi,bj)*dyC(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)
86 & )
87 tileVlAv(bi,bj) = tileVlAv(bi,bj)
88 & + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
89 tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
90
91 #ifdef ALLOW_NONHYDROSTATIC
92 IF ( nonHydrostatic ) THEN
93 tmpWke = 0.25*
94 & ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)
95 & *deepFac2F( k )*rhoFacF( k )
96 & +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1
97 & *deepFac2F(kp1)*rhoFacF(kp1)
98 & )*maskC(i,j,k,bi,bj)
99 tileVlAv(bi,bj) = tileVlAv(bi,bj)
100 & + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj)
101 tmpVal = tmpVal
102 & + tmpWke*recip_deepFac2C(k)*recip_rhoFacC(k)
103 ENDIF
104 #endif
105
106 theMax=MAX(theMax,tmpVal)
107 IF (tmpVal.NE.0.) THEN
108 tileMean(bi,bj)=tileMean(bi,bj)+tmpVal
109 numPnts=numPnts+1.
110 ENDIF
111
112 ENDDO
113 ENDDO
114 ENDDO
115 C- Potential Energy (external mode):
116 DO j=1,sNy
117 DO i=1,sNx
118 tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj)
119 & *etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
120 C- jmc: if geoid not flat (phi0surf), needs to add this term.
121 C not sure for atmos/ocean in P ; or atmos. loading in ocean-Z
122 tmpVal = tmpVal
123 & + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)
124 tilePEav(bi,bj) = tilePEav(bi,bj)
125 & + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)*maskH(i,j,bi,bj)
126 c tmpVal = etaN(i,j,bi,bj)
127 c & + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
128 c tilePEav(bi,bj) = tilePEav(bi,bj)
129 c & + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal
130 c & *rA(i,j,bi,bj)*maskH(i,j,bi,bj)
131 ENDDO
132 ENDDO
133 c theMean = theMean + tileMean(bi,bj)
134 c theVol = theVol + tileVol(bi,bj)
135 c theVolMean = theVolMean + tileVlAv(bi,bj)
136 c potEnMean = potEnMean + tilePEav(bi,bj)
137 C- end bi,bj loops
138 ENDDO
139 ENDDO
140 _GLOBAL_SUM_RL(numPnts,myThid)
141 _GLOBAL_MAX_RL(theMax,myThid)
142 c _GLOBAL_SUM_RL(theMean,myThid)
143 c _GLOBAL_SUM_RL(theVol,myThid)
144 c _GLOBAL_SUM_RL(theVolMean,myThid)
145 c _GLOBAL_SUM_RL(potEnMean, myThid)
146 CALL GLOBAL_SUM_TILE_RL( tileMean, theMean , myThid )
147 CALL GLOBAL_SUM_TILE_RL( tileVol , theVol , myThid )
148 CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )
149 CALL GLOBAL_SUM_TILE_RL( tilePEav, potEnMean , myThid )
150 IF (numPnts.NE.0.) theMean=theMean/numPnts
151 IF (theVol.NE.0.) THEN
152 theVolMean=theVolMean/theVol
153 potEnMean = potEnMean/theVol
154 ENDIF
155
156 C-- Print stats for (barotropic) Potential Energy:
157 CALL MON_SET_PREF('pe_b',myThid)
158 CALL MON_OUT_RL(mon_string_none,potEnMean,
159 & mon_foot_mean,myThid)
160
161 C-- Print stats for KE
162 CALL MON_SET_PREF('ke',myThid)
163 CALL MON_OUT_RL(mon_string_none,theMax,mon_foot_max,myThid)
164 c CALL MON_OUT_RL(mon_string_none,theMean,mon_foot_mean,myThid)
165 CALL MON_OUT_RL(mon_string_none,theVolMean,
166 & mon_foot_mean,myThid)
167 CALL MON_OUT_RL(mon_string_none,theVol,
168 & mon_foot_vol,myThid)
169
170 RETURN
171 END
172
173 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
  ViewVC Help
Powered by ViewVC 1.1.22