/[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.19 by jmc, Mon Nov 30 03:58:22 2009 UTC revision 1.22 by jmc, Fri Dec 21 01:00:09 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 theMax,theMean,theVolMean,potEnMean        _RL theMax,theMean,theVolMean,potEnMean
38          _RL uBarC, vBarC, totAMu, totAMs
39        _RL tileMean(nSx,nSy)        _RL tileMean(nSx,nSy)
40        _RL tileVlAv(nSx,nSy)        _RL tileVlAv(nSx,nSy)
41        _RL tilePEav(nSx,nSy)        _RL tilePEav(nSx,nSy)
42        _RL tileVol (nSx,nSy)        _RL tileVol (nSx,nSy)
43          _RL tileAMu (nSx,nSy)
44          _RL tileAMs (nSx,nSy)
45          _RL radDist(1:sNx,1:sNy)
46  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
47        _RL tmpWke        _RL tmpWke
48  #endif  #endif
49    
50    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
51    
52        numPnts=0.        numPnts=0.
53        theVol=0.        theVol=0.
54        theMax=0.        theMax=0.
# Line 60  C     !LOCAL VARIABLES: Line 69  C     !LOCAL VARIABLES:
69  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.
70  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).
71           msk_1 = 1.           msk_1 = 1.
72           IF ( k.EQ. 1 ) msk_1 = 0.           IF ( k.EQ.1 .AND. selectNHfreeSurf.LE.0 ) msk_1 = 0.
73           DO j=1,sNy           DO j=1,sNy
74            DO i=1,sNx            DO i=1,sNx
75             tileVol(bi,bj) = tileVol(bi,bj)             tileVol(bi,bj) = tileVol(bi,bj)
76       &                    + rA(i,j,bi,bj)*deepFac2C(k)       &                    + rA(i,j,bi,bj)*deepFac2C(k)
77       &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)       &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
78         &                     *maskInC(i,j,bi,bj)
79    
80  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)
81  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 97  C    with no assumption on how grid spac
97       &         *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)
98       &      +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)
99       &         *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)
100       &        )       &        )*maskInC(i,j,bi,bj)
101             tileVlAv(bi,bj) = tileVlAv(bi,bj)             tileVlAv(bi,bj) = tileVlAv(bi,bj)
102       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)       &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
103             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 109  C    with no assumption on how grid spac
109       &                             *deepFac2F( k )*rhoFacF( k )       &                             *deepFac2F( k )*rhoFacF( k )
110       &         +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
111       &                             *deepFac2F(kp1)*rhoFacF(kp1)       &                             *deepFac2F(kp1)*rhoFacF(kp1)
112       &        )*maskC(i,j,k,bi,bj)       &        )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
113              tileVlAv(bi,bj) = tileVlAv(bi,bj)              tileVlAv(bi,bj) = tileVlAv(bi,bj)
114       &             + 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)
115              tmpVal = tmpVal              tmpVal = tmpVal
# Line 126  C       not sure for atmos/ocean in P ; Line 136  C       not sure for atmos/ocean in P ;
136             tmpVal = tmpVal             tmpVal = tmpVal
137       &            + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)       &            + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)
138             tilePEav(bi,bj) = tilePEav(bi,bj)             tilePEav(bi,bj) = tilePEav(bi,bj)
139       &            + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)*maskH(i,j,bi,bj)       &            + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)
140         &                    *maskInC(i,j,bi,bj)
141  c          tmpVal = etaN(i,j,bi,bj)  c          tmpVal = etaN(i,j,bi,bj)
142  c    &            + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)  c    &            + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
143  c          tilePEav(bi,bj) = tilePEav(bi,bj)  c          tilePEav(bi,bj) = tilePEav(bi,bj)
144  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
145  c    &                  *rA(i,j,bi,bj)*maskH(i,j,bi,bj)  c    &                  *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
146            ENDDO            ENDDO
147           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)  
148  C- end bi,bj loops  C- end bi,bj loops
149         ENDDO         ENDDO
150        ENDDO        ENDDO
151        _GLOBAL_SUM_RL(numPnts,myThid)        _GLOBAL_SUM_RL(numPnts,myThid)
152        _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)  
153        CALL GLOBAL_SUM_TILE_RL( tileMean, theMean   , myThid )        CALL GLOBAL_SUM_TILE_RL( tileMean, theMean   , myThid )
154        CALL GLOBAL_SUM_TILE_RL( tileVol , theVol    , myThid )        CALL GLOBAL_SUM_TILE_RL( tileVol , theVol    , myThid )
155        CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )        CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )
# Line 157  c     _GLOBAL_SUM_RL(potEnMean, myThid) Line 160  c     _GLOBAL_SUM_RL(potEnMean, myThid)
160          potEnMean = potEnMean/theVol          potEnMean = potEnMean/theVol
161        ENDIF        ENDIF
162    
163    C--   Compute total angular momentum
164          IF ( mon_output_AM ) THEN
165           DO bj=myByLo(myThid),myByHi(myThid)
166            DO bi=myBxLo(myThid),myBxHi(myThid)
167    C-    calculate radial distance
168             DO j=1,sNy
169              DO i=1,sNx
170                radDist(i,j) = rSphere*COS( deg2rad*yC(i,j,bi,bj) )
171         &                            *maskInC(i,j,bi,bj)
172              ENDDO
173             ENDDO
174    C-    calculate contribution from zonal velocity
175             tileAMu(bi,bj) = 0. _d 0
176             tileAMs(bi,bj) = 0. _d 0
177             DO k=1,Nr
178              DO j=1,sNy
179               DO i=1,sNx
180                uBarC = (uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))*0.5 _d 0
181                vBarC = (vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj))*0.5 _d 0
182                tmpVal = ( angleCosC(i,j,bi,bj)*uBarC
183         &                -angleSinC(i,j,bi,bj)*vBarC
184         &               )*radDist(i,j)*deepFacC(k)
185                tileAMu(bi,bj) = tileAMu(bi,bj)
186         &             + tmpVal*rA(i,j,bi,bj)*deepFac2C(k)
187         &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
188               ENDDO
189              ENDDO
190             ENDDO
191    C-    add contribution from mass distribution anomaly (i.e., free-surface)
192    c        IF ( nonlinFreeSurf.GT.0 ) THEN
193              DO j=1,sNy
194               DO i=1,sNx
195                ks = kSurfC(i,j,bi,bj)
196                tmpVal = omega*etaN(i,j,bi,bj)
197         &             * radDist(i,j)*radDist(i,j)*deepFac2F(ks)
198                tileAMs(bi,bj) = tileAMs(bi,bj)
199         &             + tmpVal*rA(i,j,bi,bj)*deepFac2F(ks)*rhoFacF(ks)
200               ENDDO
201              ENDDO
202    c        ENDIF
203    C- end bi,bj loops
204            ENDDO
205           ENDDO
206           CALL GLOBAL_SUM_TILE_RL( tileAMu , totAMu, myThid )
207           CALL GLOBAL_SUM_TILE_RL( tileAMs , totAMs, myThid )
208    
209    C--   Print stats for total Angular Momentum:
210           CALL MON_SET_PREF('am',myThid)
211           totAMu = totAMu*rUnit2mass
212           totAMs = totAMs*rUnit2mass
213           IF ( globalArea.GT.0. ) totAMu = totAMu/globalArea
214           IF ( globalArea.GT.0. ) totAMs = totAMs/globalArea
215           CALL MON_OUT_RL( mon_string_none, totAMs,
216         &         '_eta_mean', myThid )
217           CALL MON_OUT_RL( mon_string_none, totAMu,
218         &         '_uZo_mean', myThid )
219           totAMu = totAMu + totAMs
220           CALL MON_OUT_RL( mon_string_none, totAMu,
221         &         '_tot_mean', myThid )
222    
223          ENDIF
224    
225  C--   Print stats for (barotropic) Potential Energy:  C--   Print stats for (barotropic) Potential Energy:
226        CALL MON_SET_PREF('pe_b',myThid)        CALL MON_SET_PREF('pe_b',myThid)
227        CALL MON_OUT_RL(mon_string_none,potEnMean,        CALL MON_OUT_RL(mon_string_none,potEnMean,
# Line 173  c     CALL MON_OUT_RL(mon_string_none,th Line 238  c     CALL MON_OUT_RL(mon_string_none,th
238    
239        RETURN        RETURN
240        END        END
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22