/[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.18 by jmc, Sat Nov 28 20:52:54 2009 UTC revision 1.24 by jmc, Sun Feb 17 04:07:30 2013 UTC
# Line 12  C     !INTERFACE: Line 12  C     !INTERFACE:
12       I     myIter, myThid )       I     myIter, myThid )
13    
14  C     !DESCRIPTION:  C     !DESCRIPTION:
15  C     Calculates stats for Kinetic energy  C     Calculates stats for Kinetic Energy, (barotropic) Potential Energy
16    C                      and total Angular Momentum
17    
18  C     !USES:  C     !USES:
19        IMPLICIT NONE        IMPLICIT NONE
# Line 29  C     !INPUT PARAMETERS: Line 30  C     !INPUT PARAMETERS:
30  CEOP  CEOP
31    
32  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
33        INTEGER bi,bj,i,j,k,kp1        INTEGER bi, bj
34        _RL numPnts,theVol,tmpVal,mskp1        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)        _RL tileMean(nSx,nSy)
41        _RL tileVlAv(nSx,nSy)        _RL tileVlAv(nSx,nSy)
42        _RL tilePEav(nSx,nSy)        _RL tilePEav(nSx,nSy)
43        _RL tileVol (nSx,nSy)        _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  #ifdef ALLOW_NONHYDROSTATIC
49        _RL tmpWke        _RL tmpWke
50  #endif  #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 56  C     !LOCAL VARIABLES: Line 70  C     !LOCAL VARIABLES:
70          DO k=1,Nr          DO k=1,Nr
71           kp1 = MIN(k+1,Nr)           kp1 = MIN(k+1,Nr)
72           mskp1 = 1.           mskp1 = 1.
73           IF ( k.GE.Nr) mskp1 = 0.           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           DO j=1,sNy
79            DO i=1,sNx            DO i=1,sNx
80             tileVol(bi,bj) = tileVol(bi,bj)             tileVol(bi,bj) = tileVol(bi,bj)
81       &                    + rA(i,j,bi,bj)*deepFac2C(k)       &                    + rA(i,j,bi,bj)*deepFac2C(k)
82       &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)       &                     *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)
# Line 83  C    with no assumption on how grid spac Line 102  C    with no assumption on how grid spac
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             tileVlAv(bi,bj) = tileVlAv(bi,bj)             tileVlAv(bi,bj) = tileVlAv(bi,bj)
107       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
108             tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)             tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
# Line 91  C    with no assumption on how grid spac Line 110  C    with no assumption on how grid spac
110  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
111             IF ( nonHydrostatic ) THEN             IF ( nonHydrostatic ) THEN
112              tmpWke = 0.25*              tmpWke = 0.25*
113       &        ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)       &        ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)*msk_1
114       &                             *deepFac2F( k )*rhoFacF( k )       &                             *deepFac2F( k )*rhoFacF( k )
115       &         +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1       &         +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1
116       &                             *deepFac2F(kp1)*rhoFacF(kp1)       &                             *deepFac2F(kp1)*rhoFacF(kp1)
117       &        )*maskC(i,j,k,bi,bj)       &        )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
118              tileVlAv(bi,bj) = tileVlAv(bi,bj)              tileVlAv(bi,bj) = tileVlAv(bi,bj)
119       &             + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj)       &             + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj)
120              tmpVal = tmpVal              tmpVal = tmpVal
# Line 122  C       not sure for atmos/ocean in P ; Line 141  C       not sure for atmos/ocean in P ;
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             tilePEav(bi,bj) = tilePEav(bi,bj)             tilePEav(bi,bj) = tilePEav(bi,bj)
144       &            + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)*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          tilePEav(bi,bj) = tilePEav(bi,bj)  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
 c        theMean    = theMean    + tileMean(bi,bj)  
 c        theVol     = theVol     + tileVol(bi,bj)  
 c        theVolMean = theVolMean + tileVlAv(bi,bj)  
 c        potEnMean  = potEnMean  + tilePEav(bi,bj)  
153  C- end bi,bj loops  C- end bi,bj loops
154         ENDDO         ENDDO
155        ENDDO        ENDDO
156        _GLOBAL_SUM_RL(numPnts,myThid)        _GLOBAL_SUM_RL(numPnts,myThid)
157        _GLOBAL_MAX_RL(theMax,myThid)        _GLOBAL_MAX_RL(theMax,myThid)
 c     _GLOBAL_SUM_RL(theMean,myThid)  
 c     _GLOBAL_SUM_RL(theVol,myThid)  
 c     _GLOBAL_SUM_RL(theVolMean,myThid)  
 c     _GLOBAL_SUM_RL(potEnMean, myThid)  
158        CALL GLOBAL_SUM_TILE_RL( tileMean, theMean   , myThid )        CALL GLOBAL_SUM_TILE_RL( tileMean, theMean   , myThid )
159        CALL GLOBAL_SUM_TILE_RL( tileVol , theVol    , myThid )        CALL GLOBAL_SUM_TILE_RL( tileVol , theVol    , myThid )
160        CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )        CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )
# Line 153  c     _GLOBAL_SUM_RL(potEnMean, myThid) Line 165  c     _GLOBAL_SUM_RL(potEnMean, myThid)
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,
# Line 169  c     CALL MON_OUT_RL(mon_string_none,th Line 344  c     CALL MON_OUT_RL(mon_string_none,th
344    
345        RETURN        RETURN
346        END        END
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22