Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:42:20 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
91672e10e3 Alis*0001 #include "MONITOR_OPTIONS.h"
03fb26d88c Alis*0002 
7633b97660 Ed H*0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: MON_KE
                0006 
                0007 C     !INTERFACE:
03fb26d88c Alis*0008       SUBROUTINE MON_KE(
7633b97660 Ed H*0009      I     myIter, myThid )
                0010 
                0011 C     !DESCRIPTION:
eeca444c33 Jean*0012 C     Calculates stats for Kinetic Energy, (barotropic) Potential Energy
                0013 C                      and total Angular Momentum
03fb26d88c Alis*0014 
7633b97660 Ed H*0015 C     !USES:
                0016       IMPLICIT NONE
03fb26d88c Alis*0017 #include "SIZE.h"
                0018 #include "EEPARAMS.h"
33b48b30ea Jean*0019 #include "PARAMS.h"
03fb26d88c Alis*0020 #include "DYNVARS.h"
1389d71047 Chri*0021 #include "MONITOR.h"
3eb69c3147 Jean*0022 #include "GRID.h"
a2558899b3 Jean*0023 #include "SURFACE.h"
03fb26d88c Alis*0024 
7633b97660 Ed H*0025 C     !INPUT PARAMETERS:
a2558899b3 Jean*0026       INTEGER myIter, myThid
7633b97660 Ed H*0027 CEOP
03fb26d88c Alis*0028 
7633b97660 Ed H*0029 C     !LOCAL VARIABLES:
eeca444c33 Jean*0030       INTEGER bi, bj
                0031       INTEGER i,j,k
                0032       INTEGER ks, kp1
9426193690 Jean*0033       _RL numPnts,theVol,tmpVal, mskp1, msk_1
831ca981ca Jean*0034       _RL abFac1, abFac2, R_drK, cosLat
a2558899b3 Jean*0035       _RL theMax,theMean,theVolMean,potEnMean
831ca981ca Jean*0036       _RL totAMu, totAMs
33b48b30ea Jean*0037       _RL tileMean(nSx,nSy)
                0038       _RL tileVlAv(nSx,nSy)
                0039       _RL tilePEav(nSx,nSy)
                0040       _RL tileVol (nSx,nSy)
eeca444c33 Jean*0041       _RL tileAMu (nSx,nSy)
                0042       _RL tileAMs (nSx,nSy)
831ca981ca Jean*0043       _RL tmpFld(1:sNx,1:sNy)
                0044       _RS cos2LatG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
604f9e5dc9 Jean*0045 #ifdef ALLOW_NONHYDROSTATIC
                0046       _RL tmpWke
                0047 #endif
831ca981ca Jean*0048 #ifdef ALLOW_ADAMSBASHFORTH_3
                0049       INTEGER m1, m2
                0050 #endif
03fb26d88c Alis*0051 
e32308643d Jean*0052 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0053 
a2558899b3 Jean*0054       numPnts=0.
                0055       theVol=0.
03fb26d88c Alis*0056       theMax=0.
3eb69c3147 Jean*0057       theMean=0.
                0058       theVolMean=0.
a2558899b3 Jean*0059       potEnMean =0.
03fb26d88c Alis*0060 
                0061       DO bj=myByLo(myThid),myByHi(myThid)
                0062        DO bi=myBxLo(myThid),myBxHi(myThid)
33b48b30ea Jean*0063         tileVol(bi,bj)  = 0. _d 0
                0064         tileMean(bi,bj) = 0. _d 0
                0065         tileVlAv(bi,bj) = 0. _d 0
                0066         tilePEav(bi,bj) = 0. _d 0
                0067         DO k=1,Nr
604f9e5dc9 Jean*0068          kp1 = MIN(k+1,Nr)
                0069          mskp1 = 1.
9426193690 Jean*0070          IF ( k.GE.Nr ) mskp1 = 0.
                0071 C- Note: Present NH implementation does not account for D.w/dt at k=1.
                0072 C        Consequently, wVel(k=1) does not contribute to NH KE (msk_1=0).
                0073          msk_1 = 1.
eeca444c33 Jean*0074          IF ( k.EQ.1 .AND. selectNHfreeSurf.LE.0 ) msk_1 = 0.
33b48b30ea Jean*0075          DO j=1,sNy
                0076           DO i=1,sNx
                0077            tileVol(bi,bj) = tileVol(bi,bj)
                0078      &                    + rA(i,j,bi,bj)*deepFac2C(k)
                0079      &                     *rhoFacC(k)*drF(k)*_hFacC(i,j,k,bi,bj)
e32308643d Jean*0080      &                     *maskInC(i,j,bi,bj)
ec5f8c35c4 Jean*0081 
                0082 C- Vector Invariant form (like in pkg/mom_vecinv/mom_vi_calc_ke.F)
33b48b30ea Jean*0083 c          tmpVal=0.25*( uVel( i , j ,k,bi,bj)*uVel( i , j ,k,bi,bj)
                0084 c    &                  +uVel(i+1, j ,k,bi,bj)*uVel(i+1, j ,k,bi,bj)
                0085 c    &                  +vVel( i , j ,k,bi,bj)*vVel( i , j ,k,bi,bj)
                0086 c    &                  +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) )
                0087 c          tileVlAv(bi,bj) = tileVlAv(bi,bj)
                0088 c    &              +tmpVal*rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
ec5f8c35c4 Jean*0089 
                0090 C- Energy conservative form (like in pkg/mom_fluxform/mom_calc_ke.F)
                0091 C    this is the safe way to check the energy conservation
                0092 C    with no assumption on how grid spacing & area are defined.
                0093            tmpVal=0.25*(
                0094      &       uVel( i ,j,k,bi,bj)*uVel( i ,j,k,bi,bj)
616600b8d2 Patr*0095      &         *dyG( i ,j,bi,bj)*dxC( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)
ec5f8c35c4 Jean*0096      &      +uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
616600b8d2 Patr*0097      &         *dyG(i+1,j,bi,bj)*dxC(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)
ec5f8c35c4 Jean*0098      &      +vVel(i, j ,k,bi,bj)*vVel(i, j ,k,bi,bj)
616600b8d2 Patr*0099      &         *dxG(i, j ,bi,bj)*dyC(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)
ec5f8c35c4 Jean*0100      &      +vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
616600b8d2 Patr*0101      &         *dxG(i,j+1,bi,bj)*dyC(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)
e32308643d Jean*0102      &        )*maskInC(i,j,bi,bj)
33b48b30ea Jean*0103            tileVlAv(bi,bj) = tileVlAv(bi,bj)
                0104      &                     + tmpVal*deepFac2C(k)*rhoFacC(k)*drF(k)
616600b8d2 Patr*0105            tmpVal= tmpVal*_recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
ec5f8c35c4 Jean*0106 
604f9e5dc9 Jean*0107 #ifdef ALLOW_NONHYDROSTATIC
                0108            IF ( nonHydrostatic ) THEN
                0109             tmpWke = 0.25*
9426193690 Jean*0110      &        ( wVel(i,j, k, bi,bj)*wVel(i,j, k, bi,bj)*msk_1
604f9e5dc9 Jean*0111      &                             *deepFac2F( k )*rhoFacF( k )
                0112      &         +wVel(i,j,kp1,bi,bj)*wVel(i,j,kp1,bi,bj)*mskp1
                0113      &                             *deepFac2F(kp1)*rhoFacF(kp1)
e32308643d Jean*0114      &        )*maskC(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
604f9e5dc9 Jean*0115             tileVlAv(bi,bj) = tileVlAv(bi,bj)
                0116      &             + tmpWke*rA(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj)
                0117             tmpVal = tmpVal
                0118      &             + tmpWke*recip_deepFac2C(k)*recip_rhoFacC(k)
                0119            ENDIF
                0120 #endif
                0121 
33b48b30ea Jean*0122            theMax=MAX(theMax,tmpVal)
03fb26d88c Alis*0123            IF (tmpVal.NE.0.) THEN
33b48b30ea Jean*0124             tileMean(bi,bj)=tileMean(bi,bj)+tmpVal
a2558899b3 Jean*0125             numPnts=numPnts+1.
03fb26d88c Alis*0126            ENDIF
ec5f8c35c4 Jean*0127 
03fb26d88c Alis*0128           ENDDO
                0129          ENDDO
                0130         ENDDO
a2558899b3 Jean*0131 C- Potential Energy (external mode):
33b48b30ea Jean*0132          DO j=1,sNy
                0133           DO i=1,sNx
a2558899b3 Jean*0134            tmpVal = 0.5 _d 0*Bo_surf(i,j,bi,bj)
                0135      &                      *etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
                0136 C- jmc: if geoid not flat (phi0surf), needs to add this term.
                0137 C       not sure for atmos/ocean in P ; or atmos. loading in ocean-Z
                0138            tmpVal = tmpVal
                0139      &            + phi0surf(i,j,bi,bj)*etaN(i,j,bi,bj)
33b48b30ea Jean*0140            tilePEav(bi,bj) = tilePEav(bi,bj)
d4c49a9bec Jean*0141      &            + tmpVal*rA(i,j,bi,bj)*deepFac2F(1)
                0142      &                    *maskInC(i,j,bi,bj)
a2558899b3 Jean*0143 c          tmpVal = etaN(i,j,bi,bj)
                0144 c    &            + phi0surf(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
33b48b30ea Jean*0145 c          tilePEav(bi,bj) = tilePEav(bi,bj)
a2558899b3 Jean*0146 c    &        + 0.5 _d 0*Bo_surf(i,j,bi,bj)*tmpVal*tmpVal
d4c49a9bec Jean*0147 c    &                  *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
a2558899b3 Jean*0148           ENDDO
                0149          ENDDO
                0150 C- end bi,bj loops
03fb26d88c Alis*0151        ENDDO
                0152       ENDDO
7163a40534 Jean*0153       _GLOBAL_SUM_RL(numPnts,myThid)
                0154       _GLOBAL_MAX_RL(theMax,myThid)
33b48b30ea Jean*0155       CALL GLOBAL_SUM_TILE_RL( tileMean, theMean   , myThid )
                0156       CALL GLOBAL_SUM_TILE_RL( tileVol , theVol    , myThid )
                0157       CALL GLOBAL_SUM_TILE_RL( tileVlAv, theVolMean, myThid )
                0158       CALL GLOBAL_SUM_TILE_RL( tilePEav, potEnMean , myThid )
a2558899b3 Jean*0159       IF (numPnts.NE.0.) theMean=theMean/numPnts
                0160       IF (theVol.NE.0.) THEN
                0161         theVolMean=theVolMean/theVol
                0162         potEnMean = potEnMean/theVol
                0163       ENDIF
03fb26d88c Alis*0164 
eeca444c33 Jean*0165 C--   Compute total angular momentum
                0166       IF ( mon_output_AM ) THEN
                0167        DO bj=myByLo(myThid),myByHi(myThid)
                0168         DO bi=myBxLo(myThid),myBxHi(myThid)
831ca981ca Jean*0169 C-    Calculate contribution from zonal velocity
                0170          abFac1 = 0. _d 0
                0171          abFac2 = 0. _d 0
                0172 #ifdef ALLOW_ADAMSBASHFORTH_3
                0173           m1 = 1 + mod(myIter+1,2)
                0174           m2 = 1 + mod( myIter ,2)
                0175           IF ( myIter.GE.2 ) abFac2 = beta_AB
                0176           IF ( myIter.GE.1 ) abFac1 = -( alph_AB + abFac2 )
                0177 #else
                0178           IF ( myIter.GE.1 ) abFac1 = -( 0.5 _d 0 + abEps )
                0179 #endif
                0180 C-    contribution from uVel component: 1rst integrate vertically
eeca444c33 Jean*0181          DO j=1,sNy
                0182           DO i=1,sNx
831ca981ca Jean*0183             tmpFld(i,j) = 0. _d 0
eeca444c33 Jean*0184           ENDDO
                0185          ENDDO
                0186          DO k=1,Nr
831ca981ca Jean*0187           R_drK = rSphere*deepFacC(k)*deepFac2C(k)
                0188      &                   *rhoFacC(k)*drF(k)
eeca444c33 Jean*0189           DO j=1,sNy
                0190            DO i=1,sNx
831ca981ca Jean*0191 #ifdef ALLOW_ADAMSBASHFORTH_3
                0192             tmpVal = abFac1*guNm(i,j,k,bi,bj,m1)
                0193      &             + abFac2*guNm(i,j,k,bi,bj,m2)
                0194 #else
                0195             tmpVal = abFac1*guNm1(i,j,k,bi,bj)
                0196 #endif
                0197             tmpVal = tmpVal*deltaTMom + uVel(i,j,k,bi,bj)
                0198             tmpFld(i,j) = tmpFld(i,j)
                0199      &             + R_drK*tmpVal*_hFacW(i,j,k,bi,bj)
eeca444c33 Jean*0200            ENDDO
                0201           ENDDO
                0202          ENDDO
831ca981ca Jean*0203 C-    and then integrate horizontally over this tile
                0204          DO j=1,sNy
                0205           DO i=1,sNx
                0206             cosLat = COS( deg2rad*
                0207      &             ( yG(i,j,bi,bj) + yG(i,j+1,bi,bj) )*halfRL )
                0208             tmpFld(i,j) = tmpFld(i,j)*u2zonDir(i,j,bi,bj)
                0209      &                   *cosLat*rAw(i,j,bi,bj)
                0210      &                   *maskInW(i,j,bi,bj)
                0211           ENDDO
                0212          ENDDO
                0213          tileAMu(bi,bj) = 0. _d 0
                0214          DO j=1,sNy
                0215           DO i=1,sNx
                0216             tileAMu(bi,bj) = tileAMu(bi,bj) + tmpFld(i,j)
                0217           ENDDO
                0218          ENDDO
                0219 C-    contribution from vVel component: 1rst integrate vertically
                0220          DO j=1,sNy
                0221           DO i=1,sNx
                0222             tmpFld(i,j) = 0. _d 0
                0223           ENDDO
                0224          ENDDO
                0225          DO k=1,Nr
                0226           R_drK = rSphere*deepFacC(k)*deepFac2C(k)
                0227      &                   *rhoFacC(k)*drF(k)
                0228           DO j=1,sNy
                0229            DO i=1,sNx
ba692aabde Jean*0230 #ifdef ALLOW_ADAMSBASHFORTH_3
831ca981ca Jean*0231             tmpVal = abFac1*gvNm(i,j,k,bi,bj,m1)
                0232      &             + abFac2*gvNm(i,j,k,bi,bj,m2)
ba692aabde Jean*0233 #else
831ca981ca Jean*0234             tmpVal = abFac1*gvNm1(i,j,k,bi,bj)
ba692aabde Jean*0235 #endif
831ca981ca Jean*0236             tmpVal = tmpVal*deltaTMom + vVel(i,j,k,bi,bj)
                0237             tmpFld(i,j) = tmpFld(i,j)
                0238      &             + R_drK*tmpVal*_hFacS(i,j,k,bi,bj)
                0239            ENDDO
                0240           ENDDO
                0241          ENDDO
                0242 C-    and then integrate horizontally over this tile
                0243          DO j=1,sNy
                0244           DO i=1,sNx
                0245             cosLat = COS( deg2rad*
                0246      &             ( yG(i,j,bi,bj) + yG(i+1,j,bi,bj) )*halfRL )
                0247             tmpFld(i,j) = tmpFld(i,j)*v2zonDir(i,j,bi,bj)
                0248      &                   *cosLat*rAs(i,j,bi,bj)
                0249      &                   *maskInS(i,j,bi,bj)
                0250           ENDDO
                0251          ENDDO
                0252          DO j=1,sNy
                0253           DO i=1,sNx
                0254             tileAMu(bi,bj) = tileAMu(bi,bj) + tmpFld(i,j)
                0255           ENDDO
                0256          ENDDO
                0257 C-    Calculate contribution from mass distribution anomaly (i.e., free-surface)
                0258          IF ( exactConserv ) THEN
ba692aabde Jean*0259           DO j=1,sNy
                0260            DO i=1,sNx
                0261 #ifdef EXACT_CONSERV
831ca981ca Jean*0262             tmpFld(i,j) = etaHnm1(i,j,bi,bj)
                0263 #else
                0264             tmpFld(i,j) = 0.
ba692aabde Jean*0265 #endif
                0266            ENDDO
                0267           ENDDO
                0268          ELSE
eeca444c33 Jean*0269           DO j=1,sNy
                0270            DO i=1,sNx
831ca981ca Jean*0271             tmpFld(i,j) = etaN(i,j,bi,bj)
eeca444c33 Jean*0272            ENDDO
                0273           ENDDO
ba692aabde Jean*0274          ENDIF
831ca981ca Jean*0275 C-    calculate angular momentum from mass-distribution anomaly
                0276 C     using square of radial distance (averaged @ center point)
                0277          DO j=1-OLy,sNy+OLy
                0278           DO i=1-OLx,sNx+OLx
                0279             cosLat = COS( deg2rad*yG(i,j,bi,bj) )
                0280             cos2LatG(i,j) = cosLat*cosLat
                0281           ENDDO
                0282          ENDDO
                0283          DO j=1,sNy
                0284           DO i=1,sNx
                0285             tmpFld(i,j) = tmpFld(i,j)
                0286      &        *omega*rSphere*rSphere
                0287      &        *( ( cos2LatG(i,j) + cos2LatG(i+1,j+1) )
                0288      &         + ( cos2LatG(i+1,j) + cos2LatG(i,j+1) )
                0289      &         )*0.25 _d 0
                0290           ENDDO
                0291          ENDDO
                0292          DO j=1,sNy
                0293           DO i=1,sNx
                0294             ks = kSurfC(i,j,bi,bj)
                0295             tmpFld(i,j) = tmpFld(i,j)
                0296      &             *maskInC(i,j,bi,bj)*deepFac2F(ks)
                0297      &             *rA(i,j,bi,bj)*deepFac2F(ks)*rhoFacF(ks)
                0298           ENDDO
                0299          ENDDO
                0300          tileAMs(bi,bj) = 0. _d 0
                0301          DO j=1,sNy
                0302           DO i=1,sNx
                0303             tileAMs(bi,bj) = tileAMs(bi,bj) + tmpFld(i,j)
                0304           ENDDO
                0305          ENDDO
eeca444c33 Jean*0306 C- end bi,bj loops
                0307         ENDDO
                0308        ENDDO
                0309        CALL GLOBAL_SUM_TILE_RL( tileAMu , totAMu, myThid )
                0310        CALL GLOBAL_SUM_TILE_RL( tileAMs , totAMs, myThid )
                0311 
ba692aabde Jean*0312 C--   Print stats for total Angular Momentum (per unit area, in kg/s):
eeca444c33 Jean*0313        CALL MON_SET_PREF('am',myThid)
                0314        totAMu = totAMu*rUnit2mass
                0315        totAMs = totAMs*rUnit2mass
                0316        IF ( globalArea.GT.0. ) totAMu = totAMu/globalArea
                0317        IF ( globalArea.GT.0. ) totAMs = totAMs/globalArea
                0318        CALL MON_OUT_RL( mon_string_none, totAMs,
                0319      &         '_eta_mean', myThid )
                0320        CALL MON_OUT_RL( mon_string_none, totAMu,
                0321      &         '_uZo_mean', myThid )
831ca981ca Jean*0322        totAMu = totAMu + freeSurfFac*totAMs
eeca444c33 Jean*0323        CALL MON_OUT_RL( mon_string_none, totAMu,
                0324      &         '_tot_mean', myThid )
                0325 
                0326       ENDIF
                0327 
a2558899b3 Jean*0328 C--   Print stats for (barotropic) Potential Energy:
                0329       CALL MON_SET_PREF('pe_b',myThid)
                0330       CALL MON_OUT_RL(mon_string_none,potEnMean,
                0331      &         mon_foot_mean,myThid)
                0332 
                0333 C--   Print stats for KE
                0334       CALL MON_SET_PREF('ke',myThid)
3eb69c3147 Jean*0335       CALL MON_OUT_RL(mon_string_none,theMax,mon_foot_max,myThid)
a2558899b3 Jean*0336 c     CALL MON_OUT_RL(mon_string_none,theMean,mon_foot_mean,myThid)
3eb69c3147 Jean*0337       CALL MON_OUT_RL(mon_string_none,theVolMean,
a2558899b3 Jean*0338      &         mon_foot_mean,myThid)
ec5f8c35c4 Jean*0339       CALL MON_OUT_RL(mon_string_none,theVol,
                0340      &         mon_foot_vol,myThid)
03fb26d88c Alis*0341 
                0342       RETURN
                0343       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!