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

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

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

revision 1.7.4.3 by heimbach, Tue Jun 24 23:10:27 2003 UTC revision 1.24 by jmc, Sun Feb 17 04:07:30 2013 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "MONITOR_OPTIONS.h"  #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(        SUBROUTINE MON_KE(
12       I                  myIter, myThid )       I     myIter, myThid )
13  C     /==========================================================\  
14  C     | SUBROUTINE MON_KE                                        |  C     !DESCRIPTION:
15  C     | o Calculates stats for Kinetic energy                    |  C     Calculates stats for Kinetic Energy, (barotropic) Potential Energy
16  C     |==========================================================|  C                      and total Angular Momentum
 C     \==========================================================/  
       IMPLICIT NONE  
17    
18  C     === Global data ===  C     !USES:
19          IMPLICIT NONE
20  #include "SIZE.h"  #include "SIZE.h"
21  #include "EEPARAMS.h"  #include "EEPARAMS.h"
22    #include "PARAMS.h"
23  #include "DYNVARS.h"  #include "DYNVARS.h"
24  #include "MONITOR.h"  #include "MONITOR.h"
25  #include "GRID.h"  #include "GRID.h"
26  #include "SURFACE.h"  #include "SURFACE.h"
27    
28  C     === Routine arguments ===  C     !INPUT PARAMETERS:
29        INTEGER myIter, myThid        INTEGER myIter, myThid
30    CEOP
31    
32  C     === Local variables ====  C     !LOCAL VARIABLES:
33        INTEGER bi,bj,I,J,K        INTEGER bi, bj
34        _RL numPnts,theVol,tmpVal,tmpVol        INTEGER i,j,k
35          INTEGER ks, kp1
36          _RL numPnts,theVol,tmpVal, mskp1, msk_1
37          _RL abFac1, abFac2, R_drK, cosLat
38        _RL theMax,theMean,theVolMean,potEnMean        _RL theMax,theMean,theVolMean,potEnMean
39          _RL totAMu, totAMs
40          _RL tileMean(nSx,nSy)
41          _RL tileVlAv(nSx,nSy)
42          _RL tilePEav(nSx,nSy)
43          _RL tileVol (nSx,nSy)
44          _RL tileAMu (nSx,nSy)
45          _RL tileAMs (nSx,nSy)
46          _RL tmpFld(1:sNx,1:sNy)
47          _RS cos2LatG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48    #ifdef ALLOW_NONHYDROSTATIC
49          _RL tmpWke
50    #endif
51    #ifdef ALLOW_ADAMSBASHFORTH_3
52          INTEGER m1, m2
53    #endif
54    
55    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
56    
57        numPnts=0.        numPnts=0.
58        theVol=0.        theVol=0.
# Line 37  C     === Local variables ==== Line 63  C     === Local variables ====
63    
64        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
65         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
66          DO K=1,Nr          tileVol(bi,bj)  = 0. _d 0
67           DO J=1,sNy          tileMean(bi,bj) = 0. _d 0
68            DO I=1,sNx          tileVlAv(bi,bj) = 0. _d 0
69             theVol=theVol+rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)          tilePEav(bi,bj) = 0. _d 0
70            DO k=1,Nr
71             kp1 = MIN(k+1,Nr)
72             mskp1 = 1.
73             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             IF ( k.EQ.1 .AND. selectNHfreeSurf.LE.0 ) msk_1 = 0.
78             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         &                     *maskInC(i,j,bi,bj)
84    
85  C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F)  C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F)
86  c          tmpVal=0.25*( uVel( I , J ,K,bi,bj)*uVel( I , J ,K,bi,bj)  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)  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)  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) )  c    &                  +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) )
90  c          theVolMean=theVolMean+tmpVal  c          tileVlAv(bi,bj) = tileVlAv(bi,bj)
91  c    &           *ra(i,j,bi,bj)*drf(k)*hFacC(i,j,k,bi,bj)  c    &              +tmpVal*rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
92    
93  C- Energy conservative form (like in pkg/mom_fluxform/mom_calc_ke.F)  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  C    this is the safe way to check the energy conservation
95  C    with no assumption on how grid spacing & area are defined.  C    with no assumption on how grid spacing & area are defined.
96             tmpVal=0.25*(             tmpVal=0.25*(
97       &       uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj)       &       uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj)
98       &         *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)       &         *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)
99       &      +uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)       &      +uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
100       &         *dyG(i+1,j,bi,bj)*dxC(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)       &         *dyG(i+1,j,bi,bj)*dxC(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)
101       &      +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj)       &      +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj)
102       &         *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)       &         *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)
103       &      +vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)       &      +vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
104       &         *dxG(i,j+1,bi,bj)*dyC(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)       &         *dxG(i,j+1,bi,bj)*dyC(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)
105       &        )       &        )*maskInC(i,j,bi,bj)
106             theVolMean= theVolMean + tmpVal*drF(k)             tileVlAv(bi,bj) = tileVlAv(bi,bj)
107             tmpVal= tmpVal*recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
108               tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
109    
110    #ifdef ALLOW_NONHYDROSTATIC
111               IF ( nonHydrostatic ) THEN
112                tmpWke = 0.25*
113         &        ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)*msk_1
114         &                             *deepFac2F( k )*rhoFacF( k )
115         &         +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1
116         &                             *deepFac2F(kp1)*rhoFacF(kp1)
117         &        )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
118                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             theMax=max(theMax,tmpVal)             theMax=MAX(theMax,tmpVal)
126             IF (tmpVal.NE.0.) THEN             IF (tmpVal.NE.0.) THEN
127              theMean=theMean+tmpVal              tileMean(bi,bj)=tileMean(bi,bj)+tmpVal
128              numPnts=numPnts+1.              numPnts=numPnts+1.
129             ENDIF             ENDIF
130    
# Line 76  C    with no assumption on how grid spac Line 132  C    with no assumption on how grid spac
132           ENDDO           ENDDO
133          ENDDO          ENDDO
134  C- Potential Energy (external mode):  C- Potential Energy (external mode):
135           DO J=1,sNy           DO j=1,sNy
136            DO I=1,sNx            DO i=1,sNx
137             tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj)             tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj)
138       &                      *etaN(i,j,bi,bj)*etaN(i,j,bi,bj)       &                      *etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
139  C- jmc: if geoid not flat (phi0surf), needs to add this term.  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  C       not sure for atmos/ocean in P ; or atmos. loading in ocean-Z
141             tmpVal = tmpVal             tmpVal = tmpVal
142       &            + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)       &            + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)
143             potEnMean = potEnMean             tilePEav(bi,bj) = tilePEav(bi,bj)
144       &               + tmpVal*rA(i,j,bi,bj)*maskH(i,j,bi,bj)       &            + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)
145         &                    *maskInC(i,j,bi,bj)
146  c          tmpVal = etaN(i,j,bi,bj)  c          tmpVal = etaN(i,j,bi,bj)
147  c    &            + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)  c    &            + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
148  c          potEnMean = potEnMean  c          tilePEav(bi,bj) = tilePEav(bi,bj)
149  c    &        + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal  c    &        + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal
150  c    &                  *rA(i,j,bi,bj)*maskH(i,j,bi,bj)  c    &                  *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
151            ENDDO            ENDDO
152           ENDDO           ENDDO
153  C- end bi,bj loops  C- end bi,bj loops
154         ENDDO         ENDDO
155        ENDDO        ENDDO
156        _GLOBAL_SUM_R8(numPnts,myThid)        _GLOBAL_SUM_RL(numPnts,myThid)
157        _GLOBAL_MAX_R8(theMax,myThid)        _GLOBAL_MAX_RL(theMax,myThid)
158        _GLOBAL_SUM_R8(theMean,myThid)        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        IF (numPnts.NE.0.) theMean=theMean/numPnts        IF (numPnts.NE.0.) theMean=theMean/numPnts
       _GLOBAL_SUM_R8(theVol,myThid)  
       _GLOBAL_SUM_R8(theVolMean,myThid)  
       _GLOBAL_SUM_R8(potEnMean, myThid)  
163        IF (theVol.NE.0.) THEN        IF (theVol.NE.0.) THEN
164          theVolMean=theVolMean/theVol          theVolMean=theVolMean/theVol
165          potEnMean = potEnMean/theVol          potEnMean = potEnMean/theVol
166        ENDIF        ENDIF
167    
168    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    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             DO j=1,sNy
185              DO i=1,sNx
186                tmpFld(i,j) = 0. _d 0
187              ENDDO
188             ENDDO
189             DO k=1,Nr
190              R_drK = rSphere*deepFacC(k)*deepFac2C(k)
191         &                   *rhoFacC(k)*drF(k)
192              DO j=1,sNy
193               DO i=1,sNx
194    #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               ENDDO
204              ENDDO
205             ENDDO
206    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    #ifdef ALLOW_ADAMSBASHFORTH_3
234                tmpVal = abFac1*gvNm(i,j,k,bi,bj,m1)
235         &             + abFac2*gvNm(i,j,k,bi,bj,m2)
236    #else
237                tmpVal = abFac1*gvNm1(i,j,k,bi,bj)
238    #endif
239                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              DO j=1,sNy
263               DO i=1,sNx
264    #ifdef EXACT_CONSERV
265                tmpFld(i,j) = etaHnm1(i,j,bi,bj)
266    #else
267                tmpFld(i,j) = 0.
268    #endif
269               ENDDO
270              ENDDO
271             ELSE
272              DO j=1,sNy
273               DO i=1,sNx
274                tmpFld(i,j) = etaN(i,j,bi,bj)
275               ENDDO
276              ENDDO
277             ENDIF
278    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    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    C--   Print stats for total Angular Momentum (per unit area, in kg/s):
316           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           totAMu = totAMu + freeSurfFac*totAMs
326           CALL MON_OUT_RL( mon_string_none, totAMu,
327         &         '_tot_mean', myThid )
328    
329          ENDIF
330    
331  C--   Print stats for (barotropic) Potential Energy:  C--   Print stats for (barotropic) Potential Energy:
332        CALL MON_SET_PREF('pe_b',myThid)        CALL MON_SET_PREF('pe_b',myThid)
333        CALL MON_OUT_RL(mon_string_none,potEnMean,        CALL MON_OUT_RL(mon_string_none,potEnMean,
Legend:
Removed from v.1.7.4.3  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22