/[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.21 - (show annotations) (download)
Mon May 2 15:00:51 2011 UTC (13 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62z, checkpoint62y, checkpoint62x, checkpoint63g, checkpoint64, checkpoint63, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64a
Changes since 1.20: +7 -6 lines 
use maskInC in KE statistics.
1 C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_ke.F,v 1.20 2009/12/21 00:06:03 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, msk_1
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 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
44
45 numPnts=0.
46 theVol=0.
47 theMax=0.
48 theMean=0.
49 theVolMean=0.
50 potEnMean =0.
51
52 DO bj=myByLo(myThid),myByHi(myThid)
53 DO bi=myBxLo(myThid),myBxHi(myThid)
54 tileVol(bi,bj) = 0. _d 0
55 tileMean(bi,bj) = 0. _d 0
56 tileVlAv(bi,bj) = 0. _d 0
57 tilePEav(bi,bj) = 0. _d 0
58 DO k=1,Nr
59 kp1 = MIN(k+1,Nr)
60 mskp1 = 1.
61 IF ( k.GE.Nr ) mskp1 = 0.
62 C- Note: Present NH implementation does not account for D.w/dt at k=1.
63 C Consequently, wVel(k=1) does not contribute to NH KE (msk_1=0).
64 msk_1 = 1.
65 IF ( k.EQ. 1 ) msk_1 = 0.
66 DO j=1,sNy
67 DO i=1,sNx
68 tileVol(bi,bj) = tileVol(bi,bj)
69 & + rA(i,j,bi,bj)*deepFac2C(k)
70 & *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
71 & *maskInC(i,j,bi,bj)
72
73 C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F)
74 c tmpVal=0.25*( uVel( i , j ,k,bi,bj)*uVel( i , j ,k,bi,bj)
75 c & +uVel(i+1, j ,k,bi,bj)*uVel(i+1, j ,k,bi,bj)
76 c & +vVel( i , j ,k,bi,bj)*vVel( i , j ,k,bi,bj)
77 c & +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) )
78 c tileVlAv(bi,bj) = tileVlAv(bi,bj)
79 c & +tmpVal*rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
80
81 C- Energy conservative form (like in pkg/mom_fluxform/mom_calc_ke.F)
82 C this is the safe way to check the energy conservation
83 C with no assumption on how grid spacing & area are defined.
84 tmpVal=0.25*(
85 & uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj)
86 & *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)
87 & +uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
88 & *dyG(i+1,j,bi,bj)*dxC(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)
89 & +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj)
90 & *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)
91 & +vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
92 & *dxG(i,j+1,bi,bj)*dyC(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)
93 & )*maskInC(i,j,bi,bj)
94 tileVlAv(bi,bj) = tileVlAv(bi,bj)
95 & + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
96 tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
97
98 #ifdef ALLOW_NONHYDROSTATIC
99 IF ( nonHydrostatic ) THEN
100 tmpWke = 0.25*
101 & ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)*msk_1
102 & *deepFac2F( k )*rhoFacF( k )
103 & +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1
104 & *deepFac2F(kp1)*rhoFacF(kp1)
105 & )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
106 tileVlAv(bi,bj) = tileVlAv(bi,bj)
107 & + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj)
108 tmpVal = tmpVal
109 & + tmpWke*recip_deepFac2C(k)*recip_rhoFacC(k)
110 ENDIF
111 #endif
112
113 theMax=MAX(theMax,tmpVal)
114 IF (tmpVal.NE.0.) THEN
115 tileMean(bi,bj)=tileMean(bi,bj)+tmpVal
116 numPnts=numPnts+1.
117 ENDIF
118
119 ENDDO
120 ENDDO
121 ENDDO
122 C- Potential Energy (external mode):
123 DO j=1,sNy
124 DO i=1,sNx
125 tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj)
126 & *etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
127 C- jmc: if geoid not flat (phi0surf), needs to add this term.
128 C not sure for atmos/ocean in P ; or atmos. loading in ocean-Z
129 tmpVal = tmpVal
130 & + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)
131 tilePEav(bi,bj) = tilePEav(bi,bj)
132 & + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)
133 & *maskInC(i,j,bi,bj)
134 c tmpVal = etaN(i,j,bi,bj)
135 c & + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
136 c tilePEav(bi,bj) = tilePEav(bi,bj)
137 c & + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal
138 c & *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
139 ENDDO
140 ENDDO
141 C- end bi,bj loops
142 ENDDO
143 ENDDO
144 _GLOBAL_SUM_RL(numPnts,myThid)
145 _GLOBAL_MAX_RL(theMax,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
  ViewVC Help
Powered by ViewVC 1.1.22