C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_ke.F,v 1.24 2013/02/17 04:07:30 jmc Exp $ C $Name: $ #include "MONITOR_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: MON_KE C !INTERFACE: SUBROUTINE MON_KE( I myIter, myThid ) C !DESCRIPTION: C Calculates stats for Kinetic Energy, (barotropic) Potential Energy C and total Angular Momentum C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "MONITOR.h" #include "GRID.h" #include "SURFACE.h" C !INPUT PARAMETERS: INTEGER myIter, myThid CEOP C !LOCAL VARIABLES: INTEGER bi, bj INTEGER i,j,k INTEGER ks, kp1 _RL numPnts,theVol,tmpVal, mskp1, msk_1 _RL abFac1, abFac2, R_drK, cosLat _RL theMax,theMean,theVolMean,potEnMean _RL totAMu, totAMs _RL tileMean(nSx,nSy) _RL tileVlAv(nSx,nSy) _RL tilePEav(nSx,nSy) _RL tileVol (nSx,nSy) _RL tileAMu (nSx,nSy) _RL tileAMs (nSx,nSy) _RL tmpFld(1:sNx,1:sNy) _RS cos2LatG(1-OLx:sNx+OLx,1-OLy:sNy+OLy) #ifdef ALLOW_NONHYDROSTATIC _RL tmpWke #endif #ifdef ALLOW_ADAMSBASHFORTH_3 INTEGER m1, m2 #endif C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| numPnts=0. theVol=0. theMax=0. theMean=0. theVolMean=0. potEnMean =0. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) tileVol(bi,bj) = 0. _d 0 tileMean(bi,bj) = 0. _d 0 tileVlAv(bi,bj) = 0. _d 0 tilePEav(bi,bj) = 0. _d 0 DO k=1,Nr kp1 = MIN(k+1,Nr) mskp1 = 1. IF ( k.GE.Nr ) mskp1 = 0. C- Note: Present NH implementation does not account for D.w/dt at k=1. C Consequently, wVel(k=1) does not contribute to NH KE (msk_1=0). msk_1 = 1. IF ( k.EQ.1 .AND. selectNHfreeSurf.LE.0 ) msk_1 = 0. DO j=1,sNy DO i=1,sNx tileVol(bi,bj) = tileVol(bi,bj) & + rA(i,j,bi,bj)*deepFac2C(k) & *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj) & *maskInC(i,j,bi,bj) C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F) c tmpVal=0.25*( uVel( i , j ,k,bi,bj)*uVel( i , j ,k,bi,bj) c & +uVel(i+1, j ,k,bi,bj)*uVel(i+1, j ,k,bi,bj) c & +vVel( i , j ,k,bi,bj)*vVel( i , j ,k,bi,bj) c & +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) ) c tileVlAv(bi,bj) = tileVlAv(bi,bj) c & +tmpVal*rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj) C- Energy conservative form (like in pkg/mom_fluxform/mom_calc_ke.F) C this is the safe way to check the energy conservation C with no assumption on how grid spacing & area are defined. tmpVal=0.25*( & uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj) & *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj) & +uVel(i+1,j,k,bi,bj)*uVel(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) & +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj) & *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj) & +vVel(i,j+1,k,bi,bj)*vVel(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) & )*maskInC(i,j,bi,bj) tileVlAv(bi,bj) = tileVlAv(bi,bj) & + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k) tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj) #ifdef ALLOW_NONHYDROSTATIC IF ( nonHydrostatic ) THEN tmpWke = 0.25* & ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)*msk_1 & *deepFac2F( k )*rhoFacF( k ) & +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1 & *deepFac2F(kp1)*rhoFacF(kp1) & )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj) tileVlAv(bi,bj) = tileVlAv(bi,bj) & + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj) tmpVal = tmpVal & + tmpWke*recip_deepFac2C(k)*recip_rhoFacC(k) ENDIF #endif theMax=MAX(theMax,tmpVal) IF (tmpVal.NE.0.) THEN tileMean(bi,bj)=tileMean(bi,bj)+tmpVal numPnts=numPnts+1. ENDIF ENDDO ENDDO ENDDO C- Potential Energy (external mode): DO j=1,sNy DO i=1,sNx tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj) & *etaN(i,j,bi,bj)*etaN(i,j,bi,bj) C- jmc: if geoid not flat (phi0surf), needs to add this term. C not sure for atmos/ocean in P ; or atmos. loading in ocean-Z tmpVal = tmpVal & + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj) tilePEav(bi,bj) = tilePEav(bi,bj) & + tmpVal*rA(i,j,bi,bj)*deepFac2F(1) & *maskInC(i,j,bi,bj) c tmpVal = etaN(i,j,bi,bj) c & + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj) c tilePEav(bi,bj) = tilePEav(bi,bj) c & + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal c & *rA(i,j,bi,bj)*maskInC(i,j,bi,bj) ENDDO ENDDO C- end bi,bj loops ENDDO ENDDO _GLOBAL_SUM_RL(numPnts,myThid) _GLOBAL_MAX_RL(theMax,myThid) CALL GLOBAL_SUM_TILE_RL( tileMean, theMean , myThid ) CALL GLOBAL_SUM_TILE_RL( tileVol , theVol , myThid ) CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid ) CALL GLOBAL_SUM_TILE_RL( tilePEav, potEnMean , myThid ) IF (numPnts.NE.0.) theMean=theMean/numPnts IF (theVol.NE.0.) THEN theVolMean=theVolMean/theVol potEnMean = potEnMean/theVol ENDIF C-- Compute total angular momentum IF ( mon_output_AM ) THEN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C- Calculate contribution from zonal velocity abFac1 = 0. _d 0 abFac2 = 0. _d 0 #ifdef ALLOW_ADAMSBASHFORTH_3 m1 = 1 + mod(myIter+1,2) m2 = 1 + mod( myIter ,2) IF ( myIter.GE.2 ) abFac2 = beta_AB IF ( myIter.GE.1 ) abFac1 = -( alph_AB + abFac2 ) #else IF ( myIter.GE.1 ) abFac1 = -( 0.5 _d 0 + abEps ) #endif C- contribution from uVel component: 1rst integrate vertically DO j=1,sNy DO i=1,sNx tmpFld(i,j) = 0. _d 0 ENDDO ENDDO DO k=1,Nr R_drK = rSphere*deepFacC(k)*deepFac2C(k) & *rhoFacC(k)*drF(k) DO j=1,sNy DO i=1,sNx #ifdef ALLOW_ADAMSBASHFORTH_3 tmpVal = abFac1*guNm(i,j,k,bi,bj,m1) & + abFac2*guNm(i,j,k,bi,bj,m2) #else tmpVal = abFac1*guNm1(i,j,k,bi,bj) #endif tmpVal = tmpVal*deltaTMom + uVel(i,j,k,bi,bj) tmpFld(i,j) = tmpFld(i,j) & + R_drK*tmpVal*_hFacW(i,j,k,bi,bj) ENDDO ENDDO ENDDO C- and then integrate horizontally over this tile DO j=1,sNy DO i=1,sNx cosLat = COS( deg2rad* & ( yG(i,j,bi,bj) + yG(i,j+1,bi,bj) )*halfRL ) tmpFld(i,j) = tmpFld(i,j)*u2zonDir(i,j,bi,bj) & *cosLat*rAw(i,j,bi,bj) & *maskInW(i,j,bi,bj) ENDDO ENDDO tileAMu(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx tileAMu(bi,bj) = tileAMu(bi,bj) + tmpFld(i,j) ENDDO ENDDO C- contribution from vVel component: 1rst integrate vertically DO j=1,sNy DO i=1,sNx tmpFld(i,j) = 0. _d 0 ENDDO ENDDO DO k=1,Nr R_drK = rSphere*deepFacC(k)*deepFac2C(k) & *rhoFacC(k)*drF(k) DO j=1,sNy DO i=1,sNx #ifdef ALLOW_ADAMSBASHFORTH_3 tmpVal = abFac1*gvNm(i,j,k,bi,bj,m1) & + abFac2*gvNm(i,j,k,bi,bj,m2) #else tmpVal = abFac1*gvNm1(i,j,k,bi,bj) #endif tmpVal = tmpVal*deltaTMom + vVel(i,j,k,bi,bj) tmpFld(i,j) = tmpFld(i,j) & + R_drK*tmpVal*_hFacS(i,j,k,bi,bj) ENDDO ENDDO ENDDO C- and then integrate horizontally over this tile DO j=1,sNy DO i=1,sNx cosLat = COS( deg2rad* & ( yG(i,j,bi,bj) + yG(i+1,j,bi,bj) )*halfRL ) tmpFld(i,j) = tmpFld(i,j)*v2zonDir(i,j,bi,bj) & *cosLat*rAs(i,j,bi,bj) & *maskInS(i,j,bi,bj) ENDDO ENDDO DO j=1,sNy DO i=1,sNx tileAMu(bi,bj) = tileAMu(bi,bj) + tmpFld(i,j) ENDDO ENDDO C- Calculate contribution from mass distribution anomaly (i.e., free-surface) IF ( exactConserv ) THEN DO j=1,sNy DO i=1,sNx #ifdef EXACT_CONSERV tmpFld(i,j) = etaHnm1(i,j,bi,bj) #else tmpFld(i,j) = 0. #endif ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx tmpFld(i,j) = etaN(i,j,bi,bj) ENDDO ENDDO ENDIF C- calculate angular momentum from mass-distribution anomaly C using square of radial distance (averaged @ center point) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx cosLat = COS( deg2rad*yG(i,j,bi,bj) ) cos2LatG(i,j) = cosLat*cosLat ENDDO ENDDO DO j=1,sNy DO i=1,sNx tmpFld(i,j) = tmpFld(i,j) & *omega*rSphere*rSphere & *( ( cos2LatG(i,j) + cos2LatG(i+1,j+1) ) & + ( cos2LatG(i+1,j) + cos2LatG(i,j+1) ) & )*0.25 _d 0 ENDDO ENDDO DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) tmpFld(i,j) = tmpFld(i,j) & *maskInC(i,j,bi,bj)*deepFac2F(ks) & *rA(i,j,bi,bj)*deepFac2F(ks)*rhoFacF(ks) ENDDO ENDDO tileAMs(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx tileAMs(bi,bj) = tileAMs(bi,bj) + tmpFld(i,j) ENDDO ENDDO C- end bi,bj loops ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( tileAMu , totAMu, myThid ) CALL GLOBAL_SUM_TILE_RL( tileAMs , totAMs, myThid ) C-- Print stats for total Angular Momentum (per unit area, in kg/s): CALL MON_SET_PREF('am',myThid) totAMu = totAMu*rUnit2mass totAMs = totAMs*rUnit2mass IF ( globalArea.GT.0. ) totAMu = totAMu/globalArea IF ( globalArea.GT.0. ) totAMs = totAMs/globalArea CALL MON_OUT_RL( mon_string_none, totAMs, & '_eta_mean', myThid ) CALL MON_OUT_RL( mon_string_none, totAMu, & '_uZo_mean', myThid ) totAMu = totAMu + freeSurfFac*totAMs CALL MON_OUT_RL( mon_string_none, totAMu, & '_tot_mean', myThid ) ENDIF C-- Print stats for (barotropic) Potential Energy: CALL MON_SET_PREF('pe_b',myThid) CALL MON_OUT_RL(mon_string_none,potEnMean, & mon_foot_mean,myThid) C-- Print stats for KE CALL MON_SET_PREF('ke',myThid) CALL MON_OUT_RL(mon_string_none,theMax,mon_foot_max,myThid) c CALL MON_OUT_RL(mon_string_none,theMean,mon_foot_mean,myThid) CALL MON_OUT_RL(mon_string_none,theVolMean, & mon_foot_mean,myThid) CALL MON_OUT_RL(mon_string_none,theVol, & mon_foot_vol,myThid) RETURN END