/[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.3 by adcroft, Mon Jun 4 14:25:53 2001 UTC revision 1.24 by jmc, Sun Feb 17 04:07:30 2013 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_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                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"
25    #include "GRID.h"
26    #include "SURFACE.h"
27    
28    C     !INPUT PARAMETERS:
29          INTEGER myIter, myThid
30    CEOP
31    
32  C     === Routine arguments ===  C     !LOCAL VARIABLES:
33        INTEGER myThid        INTEGER bi, bj
34          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
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     === Local variables ====  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       INTEGER bi,bj,I,J,K  
       _RL tmpVal,theMax,theMean  
       INTEGER numPnts  
56    
57          numPnts=0.
58          theVol=0.
59        theMax=0.        theMax=0.
60        theMean=0.        theMean=0.
61        numPnts=0        theVolMean=0.
62          potEnMean =0.
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             tmpVal=0.25*( uVel( I , J ,K,bi,bj)*uVel( I , J ,K,bi,bj)          tilePEav(bi,bj) = 0. _d 0
70       &                  +uVel(I+1, J ,K,bi,bj)*uVel(I+1, J ,K,bi,bj)          DO k=1,Nr
71       &                  +vVel( I , J ,K,bi,bj)*vVel( I , J ,K,bi,bj)           kp1 = MIN(k+1,Nr)
72       &                  +vVel( I ,J+1,K,bi,bj)*vVel( I ,J+1,K,bi,bj) )           mskp1 = 1.
73             theMax=max(theMax,tmpVal)           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)
86    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    
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         &         *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)
100         &         *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)
102         &         *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)
104         &         *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               tileVlAv(bi,bj) = tileVlAv(bi,bj)
107         &                     + 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)
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    
131            ENDDO            ENDDO
132           ENDDO           ENDDO
133          ENDDO          ENDDO
134    C- Potential Energy (external mode):
135             DO j=1,sNy
136              DO i=1,sNx
137               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               tilePEav(bi,bj) = tilePEav(bi,bj)
144         &            + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)
145         &                    *maskInC(i,j,bi,bj)
146    c          tmpVal = etaN(i,j,bi,bj)
147    c    &            + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
148    c          tilePEav(bi,bj) = tilePEav(bi,bj)
149    c    &        + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal
150    c    &                  *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
151              ENDDO
152             ENDDO
153    C- end bi,bj loops
154         ENDDO         ENDDO
155        ENDDO        ENDDO
156        _GLOBAL_MAX_R8(theMax,myThid)        _GLOBAL_SUM_RL(numPnts,myThid)
157        _GLOBAL_SUM_R8(theMean,myThid)        _GLOBAL_MAX_RL(theMax,myThid)
158        tmpVal=float(numPnts)        CALL GLOBAL_SUM_TILE_RL( tileMean, theMean   , myThid )
159        _GLOBAL_SUM_R8(tmpVal,myThid)        CALL GLOBAL_SUM_TILE_RL( tileVol , theVol    , myThid )
160        IF (tmpVal.NE.0.) theMean=theMean*tmpVal        CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )
161          CALL GLOBAL_SUM_TILE_RL( tilePEav, potEnMean , myThid )
162        _BEGIN_MASTER( myThid )        IF (numPnts.NE.0.) theMean=theMean/numPnts
163        WRITE(*,'(A,24x,A,1PE22.14)')        IF (theVol.NE.0.) THEN
164       &      'MON_KE: ','  max=',theMax          theVolMean=theVolMean/theVol
165        WRITE(*,'(A,24x,A,1PE22.14)')          potEnMean = potEnMean/theVol
166       &      'MON_KE: ',' mean=',theMean        ENDIF
167        _END_MASTER( )  
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:
332          CALL MON_SET_PREF('pe_b',myThid)
333          CALL MON_OUT_RL(mon_string_none,potEnMean,
334         &         mon_foot_mean,myThid)
335    
336    C--   Print stats for KE
337          CALL MON_SET_PREF('ke',myThid)
338          CALL MON_OUT_RL(mon_string_none,theMax,mon_foot_max,myThid)
339    c     CALL MON_OUT_RL(mon_string_none,theMean,mon_foot_mean,myThid)
340          CALL MON_OUT_RL(mon_string_none,theVolMean,
341         &         mon_foot_mean,myThid)
342          CALL MON_OUT_RL(mon_string_none,theVol,
343         &         mon_foot_vol,myThid)
344    
345        RETURN        RETURN
346        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22