/[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.20 by jmc, Mon Dec 21 00:06:03 2009 UTC revision 1.23 by jmc, Wed Dec 26 23:57:46 2012 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          INTEGER i,j,k
35          INTEGER ks, kp1
36        _RL numPnts,theVol,tmpVal, mskp1, msk_1        _RL numPnts,theVol,tmpVal, mskp1, msk_1
37          _RL weight0, weight1
38        _RL theMax,theMean,theVolMean,potEnMean        _RL theMax,theMean,theVolMean,potEnMean
39          _RL uBarC, vBarC, 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 radDist(1:sNx,1:sNy)
47  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
48        _RL tmpWke        _RL tmpWke
49  #endif  #endif
50    
51    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
52    
53        numPnts=0.        numPnts=0.
54        theVol=0.        theVol=0.
55        theMax=0.        theMax=0.
# Line 60  C     !LOCAL VARIABLES: Line 70  C     !LOCAL VARIABLES:
70  C- Note: Present NH implementation does not account for D.w/dt at k=1.  C- Note: Present NH implementation does not account for D.w/dt at k=1.
71  C        Consequently, wVel(k=1) does not contribute to NH KE (msk_1=0).  C        Consequently, wVel(k=1) does not contribute to NH KE (msk_1=0).
72           msk_1 = 1.           msk_1 = 1.
73           IF ( k.EQ. 1 ) msk_1 = 0.           IF ( k.EQ.1 .AND. selectNHfreeSurf.LE.0 ) msk_1 = 0.
74           DO j=1,sNy           DO j=1,sNy
75            DO i=1,sNx            DO i=1,sNx
76             tileVol(bi,bj) = tileVol(bi,bj)             tileVol(bi,bj) = tileVol(bi,bj)
77       &                    + rA(i,j,bi,bj)*deepFac2C(k)       &                    + rA(i,j,bi,bj)*deepFac2C(k)
78       &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)       &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
79         &                     *maskInC(i,j,bi,bj)
80    
81  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)
82  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 87  C    with no assumption on how grid spac Line 98  C    with no assumption on how grid spac
98       &         *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)
99       &      +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)
100       &         *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)
101       &        )       &        )*maskInC(i,j,bi,bj)
102             tileVlAv(bi,bj) = tileVlAv(bi,bj)             tileVlAv(bi,bj) = tileVlAv(bi,bj)
103       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
104             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 99  C    with no assumption on how grid spac Line 110  C    with no assumption on how grid spac
110       &                             *deepFac2F( k )*rhoFacF( k )       &                             *deepFac2F( k )*rhoFacF( k )
111       &         +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
112       &                             *deepFac2F(kp1)*rhoFacF(kp1)       &                             *deepFac2F(kp1)*rhoFacF(kp1)
113       &        )*maskC(i,j,k,bi,bj)       &        )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
114              tileVlAv(bi,bj) = tileVlAv(bi,bj)              tileVlAv(bi,bj) = tileVlAv(bi,bj)
115       &             + 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)
116              tmpVal = tmpVal              tmpVal = tmpVal
# Line 150  C- end bi,bj loops Line 161  C- end bi,bj loops
161          potEnMean = potEnMean/theVol          potEnMean = potEnMean/theVol
162        ENDIF        ENDIF
163    
164    C--   Compute total angular momentum
165          IF ( mon_output_AM ) THEN
166           DO bj=myByLo(myThid),myByHi(myThid)
167            DO bi=myBxLo(myThid),myBxHi(myThid)
168    C-    calculate radial distance
169             DO j=1,sNy
170              DO i=1,sNx
171                radDist(i,j) = rSphere*COS( deg2rad*yC(i,j,bi,bj) )
172         &                            *maskInC(i,j,bi,bj)
173              ENDDO
174             ENDDO
175    C-    calculate contribution from zonal velocity
176             tileAMu(bi,bj) = 0. _d 0
177             tileAMs(bi,bj) = 0. _d 0
178             DO k=1,Nr
179              DO j=1,sNy
180               DO i=1,sNx
181                uBarC = (uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))*0.5 _d 0
182                vBarC = (vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj))*0.5 _d 0
183                tmpVal = ( angleCosC(i,j,bi,bj)*uBarC
184         &                -angleSinC(i,j,bi,bj)*vBarC
185         &               )*radDist(i,j)*deepFacC(k)
186                tileAMu(bi,bj) = tileAMu(bi,bj)
187         &             + tmpVal*rA(i,j,bi,bj)*deepFac2C(k)
188         &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
189               ENDDO
190              ENDDO
191             ENDDO
192    C-    and contribution from mass distribution anomaly (i.e., free-surface)
193    c        IF ( .FALSE. ) THEN
194             IF ( exactConserv ) THEN
195    C-    To improve balance between zonal-wind (AMu) and mass (AMs) AM, need
196    C     a consistent time-stepping of meridional mass transport in Coriolis
197    C     (zonal momentum, go through AB) and mass transport divergence;
198    C     try to account for AB-2 in AMs the same way it applies to f*v:
199    #ifdef ALLOW_ADAMSBASHFORTH_3
200              weight1 = alph_AB
201    #else
202              weight1 = 0.5 _d 0 + abEps
203    #endif
204              weight0 = 1.0 _d 0 - weight1
205              DO j=1,sNy
206               DO i=1,sNx
207                ks = kSurfC(i,j,bi,bj)
208                tmpVal = weight1*etaH(i,j,bi,bj)
209    #ifdef EXACT_CONSERV
210         &             + weight0*etaHnm1(i,j,bi,bj)
211    #endif
212                tmpVal = omega*tmpVal
213         &             * radDist(i,j)*radDist(i,j)*deepFac2F(ks)
214                tileAMs(bi,bj) = tileAMs(bi,bj)
215         &             + tmpVal*rA(i,j,bi,bj)*deepFac2F(ks)*rhoFacF(ks)
216               ENDDO
217              ENDDO
218             ELSE
219              DO j=1,sNy
220               DO i=1,sNx
221                ks = kSurfC(i,j,bi,bj)
222                tmpVal = omega*etaN(i,j,bi,bj)
223         &             * radDist(i,j)*radDist(i,j)*deepFac2F(ks)
224                tileAMs(bi,bj) = tileAMs(bi,bj)
225         &             + tmpVal*rA(i,j,bi,bj)*deepFac2F(ks)*rhoFacF(ks)
226               ENDDO
227              ENDDO
228             ENDIF
229    C- end bi,bj loops
230            ENDDO
231           ENDDO
232           CALL GLOBAL_SUM_TILE_RL( tileAMu , totAMu, myThid )
233           CALL GLOBAL_SUM_TILE_RL( tileAMs , totAMs, myThid )
234    
235    C--   Print stats for total Angular Momentum (per unit area, in kg/s):
236           CALL MON_SET_PREF('am',myThid)
237           totAMu = totAMu*rUnit2mass
238           totAMs = totAMs*rUnit2mass
239           IF ( globalArea.GT.0. ) totAMu = totAMu/globalArea
240           IF ( globalArea.GT.0. ) totAMs = totAMs/globalArea
241           CALL MON_OUT_RL( mon_string_none, totAMs,
242         &         '_eta_mean', myThid )
243           CALL MON_OUT_RL( mon_string_none, totAMu,
244         &         '_uZo_mean', myThid )
245           totAMu = totAMu + totAMs
246           CALL MON_OUT_RL( mon_string_none, totAMu,
247         &         '_tot_mean', myThid )
248    
249          ENDIF
250    
251  C--   Print stats for (barotropic) Potential Energy:  C--   Print stats for (barotropic) Potential Energy:
252        CALL MON_SET_PREF('pe_b',myThid)        CALL MON_SET_PREF('pe_b',myThid)
253        CALL MON_OUT_RL(mon_string_none,potEnMean,        CALL MON_OUT_RL(mon_string_none,potEnMean,
# Line 166  c     CALL MON_OUT_RL(mon_string_none,th Line 264  c     CALL MON_OUT_RL(mon_string_none,th
264    
265        RETURN        RETURN
266        END        END
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22