<p><b>dwj07@fsu.edu</b> 2010-09-02 11:14:50 -0600 (Thu, 02 Sep 2010)</p><p><br>
        namelist.input.sw: Added config_stats_interval<br>
<br>
        module_global_diagnostics.F: modified variables for more readable names.<br>
                                                                 modified and checked for accurate computations.<br>
<br>
        mpas_interface.F: config_stats_interval = 0 causes no global diagnostics.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/namelist.input.sw
===================================================================
--- trunk/mpas/namelist.input.sw        2010-09-01 23:21:01 UTC (rev 492)
+++ trunk/mpas/namelist.input.sw        2010-09-02 17:14:50 UTC (rev 493)
@@ -4,6 +4,7 @@
    config_dt = 172.8
    config_ntimesteps = 7500
    config_output_interval = 500
+   config_stats_interval = 0
    config_visc  = 0.0
 /
 

Modified: trunk/mpas/src/core_sw/module_global_diagnostics.F
===================================================================
--- trunk/mpas/src/core_sw/module_global_diagnostics.F        2010-09-01 23:21:01 UTC (rev 492)
+++ trunk/mpas/src/core_sw/module_global_diagnostics.F        2010-09-02 17:14:50 UTC (rev 493)
@@ -28,7 +28,7 @@
       ! 1. Define the array to integrate, and the variable for the value above.
       ! 2. Allocate the array with the correct dimensions.
       ! 3. Fill the array with the data to be integrated.
-      !     eg. G_h = Sum(h dA)/Sum(dA), See below for array filling
+      !     eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
       ! 4. Call Function to compute Global Stat that you want.
       ! 5. Finish computing the global stat/integral
       ! 6. Write out your global stat to the file
@@ -46,26 +46,36 @@
       integer :: nCells
 
       ! Step 1
-      ! 1. Define the array to integrate, and the variable for the value above.
-      real (kind=RKIND) ::  areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
-      real (kind=RKIND), dimension(:), pointer ::  areaCell, dcEdge, dvEdge, areaTriangle, areaEdge, h_s
-      real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex
+      ! 1. Define the array to integrate, and the variable for the value to be stored in after the integration
+      real (kind=RKIND), dimension(:), pointer ::  areaCell, dcEdge, dvEdge, areaTriangle, h_s, fCell, fEdge
+      real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex, weightsOnEdge
 
-      real (kind=RKIND), dimension(:), allocatable :: eres, h_ref
-      real (kind=RKIND), dimension(:,:), allocatable :: pe_vertex, u2, h2, hb, pv_temp, pve_bot, h_top, h_bot, cor_energy, sshr, ke_sink, pe_sink
       real (kind=RKIND), dimension(:,:,:), pointer :: tracers
 
-      real (kind=RKIND) ::  global_temp, gh, gpv, gpe, ge, ger, e_cor, sshr_val, h_ref_val, ke_sink_val, pe_sink_val, ke
+      real (kind=RKIND), dimension(:), allocatable :: volumeWeightedPotentialEnergyReservoir, averageThickness
+      real (kind=RKIND), dimension(:), allocatable :: potentialEnstrophyReservior, areaEdge, h_s_edge
+
+      real (kind=RKIND), dimension(:,:), allocatable :: cellVolume, cellArea, volumeWeightedPotentialVorticity
+      real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnstrophy, vertexVolume, volumeWeightedKineticEnergy 
+      real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnergy, volumeWeightedPotentialEnergyTopography 
+      real (kind=RKIND), dimension(:,:), allocatable :: keTend_CoriolisForce, keTend_PressureGradient 
+      real (kind=RKIND), dimension(:,:), allocatable ::peTend_DivThickness, refAreaWeightedSurfaceHeight, refAreaWeightedSurfaceHeight_edge
+
+      real (kind=RKIND) :: sumCellVolume, sumCellArea, sumVertexVolume, sumrefAreaWeightedSurfaceHeight
+
+      real (kind=RKIND) :: globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, globalEnergy 
+      real (kind=RKIND) :: globalCoriolisEnergyTendency, globalKEPETendency, globalPotentialEnstrophyReservoir 
+      real (kind=RKIND) :: globalKineticEnergy, globalPotentialEnergy, globalPotentialEnergyReservoir
+      real (kind=RKIND) :: globalKineticEnergyTendency, globalPotentialEnergyTendency
+      real (kind=RKIND) ::  global_temp, workpv, q
       real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
+
       integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
-      integer :: timeLevel,k,i
+      integer :: timeLevel, eoe, iLevel, iCell, iEdge, iVertex
+      integer :: fileID, iCell1, iCell2, j
 
-      integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
-
-      integer, parameter :: kMaxVariables = 1024 ! this must be a little more than double the number of variables to be reduced
-      real (kind=RKIND), dimension(kMaxVariables) :: sums, mins, maxes, averages, verticalSumMins, verticalSumMaxes, reductions
-
-      integer :: fileID, iCell1, iCell2, iEdge
+      integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgesOnEdge
+      integer, dimension(:), pointer :: nEdgesOnEdge
       
       cellsOnEdge =&gt; grid % cellsOnEdge % array
       edgesOnCell =&gt; grid % edgesOnCell % array
@@ -81,8 +91,14 @@
       dcEdge =&gt; grid % dcEdge % array
       dvEdge =&gt; grid % dvEdge % array
       areaTriangle =&gt; grid % areaTriangle % array
+      fCell =&gt; grid % fCell % array
+      fEdge =&gt; grid % fEdge % array
+      edgesOnEdge       =&gt; grid % edgesOnEdge % array
+      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
+
       allocate(areaEdge(1:nEdgesSolve))
       areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
+      weightsOnEdge     =&gt; grid % weightsOnEdge % array
 
       h =&gt; state % h % array
       u =&gt; state % u % array
@@ -96,114 +112,171 @@
 
       ! Step 2
       ! 2. Allocate the array with the correct dimensions.
-      allocate(h_top(nVertLevels,nCellsSolve))
-      allocate(h_bot(nVertLevels,nCellsSolve))
-      allocate(h2(nVertLevels,nCellsSolve))
-      allocate(hb(nVertLevels,nCellsSolve))
-      allocate(sshr(nVertLevels,nCellsSolve))
+      allocate(cellVolume(nVertLevels,nCellsSolve))
+      allocate(cellArea(nVertLevels,nCellsSolve))
+      allocate(refAreaWeightedSurfaceHeight(nVertLevels,nCellsSolve))
+      allocate(refAreaWeightedSurfaceHeight_edge(nVertLevels,nEdgesSolve))
+      allocate(volumeWeightedPotentialVorticity(nVertLevels,nVerticesSolve))
+      allocate(volumeWeightedPotentialEnstrophy(nVertLevels,nVerticesSolve))
+      allocate(potentialEnstrophyReservior(nCellsSolve))
+      allocate(vertexVolume(nVertLevels,nVerticesSolve))
+      allocate(volumeWeightedKineticEnergy(nVertLevels,nEdgesSolve))
+      allocate(volumeWeightedPotentialEnergy(nVertLevels,nCellsSolve))
+      allocate(volumeWeightedPotentialEnergyTopography(nVertLevels,nCellsSolve))
+      allocate(volumeWeightedPotentialEnergyReservoir(nCellsSolve))
+      allocate(keTend_CoriolisForce(nVertLevels,nEdgesSolve))
+      allocate(keTend_PressureGradient(nVertLevels,nEdgesSolve))
+      allocate(peTend_DivThickness(nVertLevels,nCells))
 
-      allocate(h_ref(nCellsSolve))
+      allocate(averageThickness(nCellsSolve))
 
-      allocate(pv_temp(nVertLevels,nVerticesSolve))
-      allocate(pe_vertex(nVertLevels,nVerticesSolve))
-      allocate(pve_bot(nVertLevels,nVerticesSolve))
+      allocate(h_s_edge(nEdgesSOlve))
 
-      allocate(u2(nVertLevels,nEdgesSolve))
-      allocate(cor_energy(nVertLevels,nEdgesSolve))
 
-      allocate(ke_sink(nVertLevels,nEdgesSolve))
-      allocate(pe_sink(nVertLevels,nCells))
+      cellVolume = 0
+      refAreaWeightedSurfaceHeight = 0
+      refAreaWeightedSurfaceHeight_edge = 0
+      vertexVolume = 0
+      cellArea = 0
+      averageThickness = 0
+      volumeWeightedPotentialVorticity = 0
+      volumeWeightedPotentialEnstrophy = 0
+      volumeWeightedKineticEnergy = 0
+      volumeWeightedPotentialEnergy = 0
+      volumeWeightedPotentialEnergyTopography = 0
+      volumeWeightedPotentialEnergyReservoir = 0
+      keTend_PressureGradient = 0
+      peTend_DivThickness = 0
+      keTend_CoriolisForce = 0
+      h_s_edge = 0
 
-      allocate(eres(nCellsSolve))
-
-      pe_sink = 0
-
       ! Build Arrays for Global Integrals
       ! Step 3
       ! 3. Fill the array with the data to be integrated.
-      !     eg. G_h = Sum(h dA)/Sum(dA), See below for array filling
-      do i = 1,nVertLevels
-        ! G_h Top = Sum(h dA)
-        h_top(i,:) = h(i,1:nCellsSolve)*areaCell(1:nCellsSolve)
-        ! G_h Bot = Sum(dA)
-        h_bot(i,:) = areaCell(1:nCellsSolve)
-        pv_temp(i,:) = pv_vertex(i,1:nVerticesSolve)*h_vertex(i,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve) 
-        pe_vertex(i,:) = pv_vertex(i,1:nVerticesSolve)*pv_vertex(i,1:nVerticesSolve)*h_vertex(i,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
-        pve_bot(i,:) = h_vertex(i,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
-        u2(i,:) = u(i,1:nEdgesSolve)*u(i,1:nEdgesSolve)*h_edge(i,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
-        h2(i,:) = gravity*h(i,1:nCellsSolve)*h(i,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
-        hb(i,:) = gravity*h(i,1:nCellsSolve)*h_s(1:nCellsSolve)*areaCell(1:nCellsSolve)
-        cor_energy(i,:) = h_edge(i,1:nEdgesSolve)*u(i,1:nEdgesSolve)*pv_edge(i,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)
-        sshr(i,:) = areaCell(1:nCellsSolve)*(h(i,1:nCellsSolve)+h_s(1:nCellsSolve))
+      !     eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
+      do iLevel = 1,nVertLevels
+        ! eg. GlobalFluidThickness top (Sum( h dA)) = Sum(cellVolume)
+        cellVolume(iLevel,:) = h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
+        ! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
+        cellArea(iLevel,:) = areaCell(1:nCellsSolve)
+        volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &amp;
+                *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve) 
+        volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &amp; 
+                *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+        vertexVolume(iLevel,:) = h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+        volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &amp;
+                *h_edge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
+        volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
+        volumeWeightedPotentialEnergyTopography(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h_s(1:nCellsSolve)*areaCell(1:nCellsSolve)
+        refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(h(iLevel,1:nCellsSolve)+h_s(1:nCellsSolve))
 
-        do k = 1,nEdgesSolve
-            iCell1 = cellsOnEdge(k,1)
-            iCell2 = cellsOnEdge(k,2)
-            ke_sink(i,k) = areaEdge(k)*h_edge(i,k)*u(i,k)*gravity*(h(i,iCell2)+h_s(iCell2) - h(i,iCell1)-h_s(iCell1))/dcEdge(k)
-            pe_sink(i,iCell1) = pe_sink(i,iCell1) + h_edge(i,k)*u(i,k)*dvEdge(k)
-            pe_sink(i,iCell2) = pe_sink(i,iCell2) - h_edge(i,k)*u(i,k)*dvEdge(k)
+        do iEdge = 1,nEdgesSolve
+            q = 0.0
+            do j = 1,nEdgesOnEdge(iEdge)
+               eoe = edgesOnEdge(j,iEdge)
+               workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
+               q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe) 
+            end do
+            keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
+
+            iCell1 = cellsOnEdge(iEdge,1)
+            iCell2 = cellsOnEdge(iEdge,2)
+
+            refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(h_edge(iLevel,iEdge) + 0.5*(h_s(iCell1) + h_s(iCell2)))
+
+            keTend_PressureGradient(iLevel,iEdge) = areaEdge(iEdge)*h_edge(iLevel,iEdge)*u(iLevel,iEdge) &amp;
+                        *gravity*(h(iLevel,iCell2)+h_s(iCell2) - h(iLevel,iCell1)-h_s(iCell1))/dcEdge(iEdge)
+            peTend_DivThickness(iLevel,iCell1) = peTend_DivThickness(iLevel,iCell1) &amp;
+                        + h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
+            peTend_DivThickness(iLevel,iCell2) = peTend_DivThickness(iLevel,iCell2) &amp;
+                        - h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
         end do
-        pe_sink(i,:) = pe_sink(i,1:nCellsSolve)*gravity*(h(i,1:nCellsSolve)+h_s(1:nCellsSolve))
+
+        peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &amp;
+                   *(h(iLevel,1:nCells)+h_s(1:nCells))
       end do
 
+      do iEdge = 1,nEdgesSolve
+          iCell1 = cellsOnEdge(iEdge,1)
+          iCell2 = cellsOnEdge(iEdge,2)
+          
+          h_s_edge(iEdge) = 0.5*(h_s(iCell1) + h_s(iCell2))
+      end do
+
       ! Step 4
       ! 4. Call Function to compute Global Stat that you want.
-      ! KE and PE Sink Terms
-      call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, ke_sink, ke_sink_val)
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, pe_sink, pe_sink_val)
+      ! Computing Kinetic and Potential Energy Tendency Terms
+      call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, keTend_PressureGradient, globalKineticEnergyTendency)
+      call computeGlobalSum(dminfo, nVertLevels, nCells, peTend_DivThickness, globalPotentialEnergyTendency)
 
-      call computeGlobalSum(dminfo, 1, nCellsSolve, h_ref, h_ref_val)
+      ! Computing top and bottom of global mass integral
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellVolume, sumCellVolume)
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellArea, sumCellArea)
 
-      ! Global Area Average of h
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, h_top, gh)
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, h_bot, global_temp)
+      globalKineticEnergyTendency = globalKineticEnergyTendency / sumCellVolume
+      globalPotentialEnergyTendency = globalPotentialEnergyTendency / sumCellVolume
 
       ! Step 5
       ! 5. Finish computing the global stat/integral
-      gh = gh/global_temp
+      globalFluidThickness = sumCellVolume/sumCellArea
 
-      ! Mean SSH for PE Res
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, sshr, sshr_val)
+      ! Compute Average Sea Surface Height for Potential Energy and Enstrophy
+      ! Reservoir computations
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, refAreaWeightedSurfaceHeight, sumrefAreaWeightedSurfaceHeight)
 
-      h_ref(:) = (sshr_val/global_temp)-h_s(1:nCellsSolve)
+      averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
 
+      ! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
+      call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
+      call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
+      call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
 
-      ! Gloabl Area-Thickness Average of pv and pe
-      call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, pv_temp, gpv)
-      call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, pe_vertex, gpe)
-      call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, pve_bot, global_temp)
+      globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
+      globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
 
-      gpv = gpv/global_temp
-      gpe = gpe/global_temp
+      ! Compte Potential Enstrophy Reservior
+      potentialEnstrophyReservior(:) = areaCell(:)*fCell(:)*fCell(:)/averageThickness
+      call computeGlobalSum(dminfo, 1, nCellsSolve, potentialEnstrophyReservior, globalPotentialEnstrophyReservoir)
+      globalPotentialEnstrophyReservoir = globalPotentialEnstrophyReservoir/sumCellVolume
 
-      ! Global Integral of Total Energy with subtracting reference res
-      call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, u2, ke)
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, h2, global_temp)
-      ge = ke + global_temp
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, hb, global_temp)
-      ge = ge + global_temp
+      globalPotentialEnstrophy = globalPotentialEnstrophy - globalPotentialEnstrophyReservoir
 
-      eres(1:nCellsSolve) = areaCell(1:nCellsSolve)*h_ref*h_ref*gravity*0.5
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, eres, global_temp)
-      ger = global_temp
-      eres(1:nCellsSolve) = areaCell(1:nCellsSolve)*h_ref*h_s(1:nCellsSolve)*gravity
-      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, eres, global_temp)
-      ger = ger + global_temp
-      ge = ge - ger
+      ! Compute Kinetic and Potential Energy terms to be combined into total energy
+      call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, volumeWeightedKineticEnergy, globalKineticEnergy)
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergy, globalPotentialEnergy)
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyTopography, global_temp)
 
-      ! Global Integral of spurious Coriolis energy for Time to KE doubling
-      call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, cor_energy, e_cor)
+      globalKineticEnergy = globalKineticEnergy/sumCellVolume
+      globalPotentialEnergy = (globalPotentialEnergy + global_temp)/sumCellVolume
 
+      ! Compute Potential energy reservoir to be subtracted from potential energy term
+      volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*averageThickness*gravity*0.5
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, globalPotentialEnergyReservoir)
+      volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s(1:nCellsSolve)*gravity
+      call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, global_temp)
+
+      globalPotentialEnergyReservoir = (globalPotentialEnergyReservoir + global_temp)/sumCellVolume
+
+      globalPotentialEnergy = globalPotentialEnergy - globalPotentialEnergyReservoir
+      globalEnergy = globalKineticEnergy + globalPotentialEnergy
+
+      ! Compute Coriolis energy tendency term
+      call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, keTend_CoriolisForce, globalCoriolisEnergyTendency)
+      globalCoriolisEnergyTendency = globalCoriolisEnergyTendency/sumCellVolume
+
       ! Step 6
       ! 6. Write out your global stat to the file
       if (dminfo % my_proc_id == IO_NODE) then
          fileID = getFreeUnit()
+
          if (timeIndex/config_stats_interval == 1) then
              open(fileID, file='GlobalIntegrals.txt',STATUS='unknown')
          else
-             open(fileID, file='GlobalIntegrals.txt',ACCESS='append')
+             open(fileID, file='GlobalIntegrals.txt',POSITION='append')
          endif 
-         write(fileID,'(100es24.16)') gh, gpv, gpe, ge, e_cor, abs(ke/e_cor), ke_sink_val+pe_sink_val,abs(ke/(ke_sink_val+pe_sink_val))
+         write(fileID,'(1i, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &amp;
+                        globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &amp;
+                        globalKineticEnergy, globalPotentialEnergy
          close(fileID)
       end if
 

Modified: trunk/mpas/src/core_sw/mpas_interface.F
===================================================================
--- trunk/mpas/src/core_sw/mpas_interface.F        2010-09-01 23:21:01 UTC (rev 492)
+++ trunk/mpas/src/core_sw/mpas_interface.F        2010-09-02 17:14:50 UTC (rev 493)
@@ -64,18 +64,20 @@
 
    call timestep(domain, dt)
 
-   if (mod(itimestep, config_stats_interval) == 0) then
-       block_ptr =&gt; domain % blocklist
-       if(associated(block_ptr % next)) then
-           write(0,*) 'Error: computeGlobalDiagnostics assumes ',&amp;
-                      'that there is only one block per processor.'
+   if(config_stats_interval .gt. 0) then
+       if(mod(itimestep, config_stats_interval) == 0) then
+           block_ptr =&gt; domain % blocklist
+           if(associated(block_ptr % next)) then
+               write(0,*) 'Error: computeGlobalDiagnostics assumes ',&amp;
+                          'that there is only one block per processor.'
+           end if
+
+           call timer_start(&quot;global_diagnostics&quot;)
+           call computeGlobalDiagnostics(domain % dminfo, &amp;
+                    block_ptr % time_levs(2) % state, block_ptr % mesh, &amp;
+                    itimestep, dt)
+           call timer_stop(&quot;global_diagnostics&quot;)
        end if
-
-       call timer_start(&quot;global_diagnostics&quot;)
-       call computeGlobalDiagnostics(domain % dminfo, &amp;
-                block_ptr % time_levs(2) % state, block_ptr % mesh, &amp;
-                itimestep, dt)
-       call timer_stop(&quot;global_diagnostics&quot;)
    end if
 
 end subroutine mpas_timestep

</font>
</pre>