<p><b>dwj07@fsu.edu</b> 2010-07-23 14:23:10 -0600 (Fri, 23 Jul 2010)</p><p><br>
Adding global diagnostics module for shallow water core.<br>
<br>
Modified:<br>
mpas_interface - added call to global_diagnostics<br>
Registry - made h_vertex an array that could be used and saved, and vorticity_cell<br>
module_time_integration - computed h_vertex to be stored in an array, as well as vorticity_cell<br>
Makefile - added module_global_diagnostics to make list<br>
<br>
Added:<br>
module_global_diagnostics - compute actual diagnostics and write them to a file.<br>
<br>
A mpas/src/core_sw/module_global_diagnostics.F<br>
M mpas/src/core_sw/mpas_interface.F<br>
M mpas/src/core_sw/Registry<br>
M mpas/src/core_sw/module_time_integration.F<br>
M mpas/src/core_sw/Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/src/core_sw/Makefile
===================================================================
--- trunk/mpas/src/core_sw/Makefile 2010-07-23 18:28:16 UTC (rev 394)
+++ trunk/mpas/src/core_sw/Makefile 2010-07-23 20:23:10 UTC (rev 395)
@@ -2,7 +2,8 @@
OBJS = module_test_cases.o \
module_time_integration.o \
- mpas_interface.o
+ module_global_diagnostics.o \
+ mpas_interface.o
all: core_sw
@@ -13,7 +14,7 @@
module_time_integration.o:
-mpas_interface.o: module_test_cases.o module_time_integration.o
+mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: trunk/mpas/src/core_sw/Registry
===================================================================
--- trunk/mpas/src/core_sw/Registry 2010-07-23 18:28:16 UTC (rev 394)
+++ trunk/mpas/src/core_sw/Registry 2010-07-23 20:23:10 UTC (rev 395)
@@ -6,6 +6,7 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
+namelist integer sw_model config_stats_interval 100
namelist real sw_model config_visc 0.0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -79,6 +80,7 @@
var real kiteAreasOnVertex ( vertexDegree nVertices ) iro kiteAreasOnVertex - -
var real fEdge ( nEdges ) iro fEdge - -
var real fVertex ( nVertices ) iro fVertex - -
+var real fCell ( nCells ) iro fCell - -
var real h_s ( nCells ) iro h_s - -
# Arrays required for reconstruction of velocity field
@@ -97,6 +99,7 @@
var real v ( nVertLevels nEdges Time ) o v - -
var real divergence ( nVertLevels nCells Time ) o divergence - -
var real vorticity ( nVertLevels nVertices Time ) o vorticity - -
+var real vorticity_cell ( nVertLevels nCells Time ) o vorticity_cell - -
var real pv_edge ( nVertLevels nEdges Time ) o pv_edge - -
var real h_edge ( nVertLevels nEdges Time ) o h_edge - -
var real ke ( nVertLevels nCells Time ) o ke - -
@@ -113,3 +116,4 @@
var real circulation ( nVertLevels nVertices Time ) - circulation - -
var real gradPVt ( nVertLevels nEdges Time ) - gradPVt - -
var real gradPVn ( nVertLevels nEdges Time ) - gradPVn - -
+var real h_vertex ( nVertLevels nVertices Time ) - h_vertex - -
Added: trunk/mpas/src/core_sw/module_global_diagnostics.F
===================================================================
--- trunk/mpas/src/core_sw/module_global_diagnostics.F (rev 0)
+++ trunk/mpas/src/core_sw/module_global_diagnostics.F 2010-07-23 20:23:10 UTC (rev 395)
@@ -0,0 +1,311 @@
+module global_diagnostics
+
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+
+ implicit none
+ save
+ public
+
+ contains
+
+ subroutine computeGlobalDiagnostics(dminfo, state, grid, timeIndex, dt)
+
+ ! Note: this routine assumes that there is only one block per processor. No looping
+ ! is preformed over blocks.
+ ! dminfo is the domain info needed for global communication
+ ! state contains the state variables needed to compute global diagnostics
+ ! grid conains the meta data about the grid
+ ! timeIndex is the current time step counter
+ ! dt is the duration of each time step
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! INSTRUCTIONS !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! To add a new Diagnostic as a Global Stat, follow these steps.
+ ! 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
+ ! 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
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ type (grid_state), intent(inout) :: state
+ type (grid_meta), intent(in) :: grid
+ integer, intent(in) :: timeIndex
+ real (kind=RKIND), intent(in) :: dt
+
+ integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, nCellsGlobal, nEdgesGlobal, nVerticesGlobal, iTracer
+ 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
+
+ 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) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
+ integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
+ integer :: timeLevel,k,i
+
+ 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
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+
+ nVertLevels = grid % nVertLevels
+ nCellsSolve = grid % nCellsSolve
+ nEdgesSolve = grid % nEdgesSolve
+ nVerticesSolve = grid % nVerticesSolve
+ nCells = grid % nCells
+
+ h_s => grid % h_s % array
+ areaCell => grid % areaCell % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaTriangle => grid % areaTriangle % array
+ allocate(areaEdge(1:nEdgesSolve))
+ areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
+
+ h => state % h % array
+ u => state % u % array
+ v => state % v % array
+ tracers => state % tracers % array
+ h_edge => state % h_edge % array
+ h_vertex => state % h_vertex % array
+ pv_edge => state % pv_edge % array
+ pv_vertex => state % pv_vertex % array
+ pv_cell => state % pv_cell % array
+
+ ! 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(h_ref(nCellsSolve))
+
+ allocate(pv_temp(nVertLevels,nVerticesSolve))
+ allocate(pe_vertex(nVertLevels,nVerticesSolve))
+ allocate(pve_bot(nVertLevels,nVerticesSolve))
+
+ allocate(u2(nVertLevels,nEdgesSolve))
+ allocate(cor_energy(nVertLevels,nEdgesSolve))
+
+ allocate(ke_sink(nVertLevels,nEdgesSolve))
+ allocate(pe_sink(nVertLevels,nCells))
+
+ 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))
+
+ 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)
+ end do
+ pe_sink(i,:) = pe_sink(i,1:nCellsSolve)*gravity*(h(i,1:nCellsSolve)+h_s(1:nCellsSolve))
+ 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)
+
+ call computeGlobalSum(dminfo, 1, nCellsSolve, h_ref, h_ref_val)
+
+ ! Global Area Average of h
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, h_top, gh)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, h_bot, global_temp)
+
+ ! Step 5
+ ! 5. Finish computing the global stat/integral
+ gh = gh/global_temp
+
+ ! Mean SSH for PE Res
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, sshr, sshr_val)
+
+ h_ref(:) = (sshr_val/global_temp)-h_s(1:nCellsSolve)
+
+
+ ! 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)
+
+ gpv = gpv/global_temp
+ gpe = gpe/global_temp
+
+ ! 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
+
+ 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
+
+ ! Global Integral of spurious Coriolis energy for Time to KE doubling
+ call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, cor_energy, e_cor)
+
+ ! 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')
+ 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))
+ close(fileID)
+ end if
+
+ deallocate(areaEdge)
+ end subroutine computeGlobalDiagnostics
+
+ integer function getFreeUnit()
+ implicit none
+
+ integer :: index
+ logical :: isOpened
+
+ getFreeUnit = 0
+ do index = 1,99
+ if((index /= 5) .and. (index /= 6)) then
+ inquire(unit = index, opened = isOpened)
+ if( .not. isOpened) then
+ getFreeUnit = index
+ return
+ end if
+ end if
+ end do
+ end function getFreeUnit
+
+ subroutine computeGlobalSum(dminfo, nVertLevels, nElements, field, globalSum)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalSum
+
+ real (kind=RKIND) :: localSum
+
+ localSum = sum(field)
+ call dmpar_sum_real(dminfo, localSum, globalSum)
+
+ end subroutine computeGlobalSum
+
+ subroutine computeGlobalMin(dminfo, nVertLevels, nElements, field, globalMin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMin
+
+ real (kind=RKIND) :: localMin
+
+ localMin = minval(field)
+ call dmpar_min_real(dminfo, localMin, globalMin)
+
+ end subroutine computeGlobalMin
+
+ subroutine computeGlobalMax(dminfo, nVertLevels, nElements, field, globalMax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMax
+
+ real (kind=RKIND) :: localMax
+
+ localMax = maxval(field)
+ call dmpar_max_real(dminfo, localMax, globalMax)
+
+ end subroutine computeGlobalMax
+
+ subroutine computeGlobalVertSumHorizMin(dminfo, nVertLevels, nElements, field, globalMin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMin
+
+ real (kind=RKIND) :: localMin
+
+ localMin = minval(sum(field,1))
+ call dmpar_min_real(dminfo, localMin, globalMin)
+
+ end subroutine computeGlobalVertSumHorizMin
+
+ subroutine computeGlobalVertSumHorizMax(dminfo, nVertLevels, nElements, field, globalMax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMax
+
+ real (kind=RKIND) :: localMax
+
+ localMax = maxval(sum(field,1))
+ call dmpar_max_real(dminfo, localMax, globalMax)
+
+ end subroutine computeGlobalVertSumHorizMax
+
+end module global_diagnostics
Modified: trunk/mpas/src/core_sw/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_sw/module_time_integration.F 2010-07-23 18:28:16 UTC (rev 394)
+++ trunk/mpas/src/core_sw/module_time_integration.F 2010-07-23 20:23:10 UTC (rev 395)
@@ -239,12 +239,12 @@
type (grid_meta), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
+ real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, divergence
+ circulation, vorticity, ke, pv_edge, divergence, h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: u_diffusion, visc
@@ -445,12 +445,13 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv
+ real (kind=RKIND) :: flux, vorticity_abs, workpv
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
+ circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
+ h_vertex, vorticity_cell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -461,6 +462,7 @@
v => s % v % array
vh => s % vh % array
h_edge => s % h_edge % array
+ h_vertex => s % h_vertex % array
tend_h => s % h % array
tend_u => s % u % array
circulation => s % circulation % array
@@ -470,6 +472,7 @@
pv_edge => s % pv_edge % array
pv_vertex => s % pv_vertex % array
pv_cell => s % pv_cell % array
+ vorticity_cell => s % vorticity_cell % array
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
@@ -621,13 +624,13 @@
if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
end do
do k=1,nVertLevels
- h_vertex = 0.0
+ h_vertex(k,iVertex) = 0.0
do i=1,grid % vertexDegree
- h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+ h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
- h_vertex = h_vertex / areaTriangle(iVertex)
+ h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
- pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
end do
end do VTX_LOOP
@@ -674,12 +677,14 @@
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
!
pv_cell(:,:) = 0.0
+ vorticity_cell(:,:) = 0.0
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
+ vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
Modified: trunk/mpas/src/core_sw/mpas_interface.F
===================================================================
--- trunk/mpas/src/core_sw/mpas_interface.F 2010-07-23 18:28:16 UTC (rev 394)
+++ trunk/mpas/src/core_sw/mpas_interface.F 2010-07-23 20:23:10 UTC (rev 395)
@@ -52,15 +52,32 @@
use grid_types
use time_integration
+ use timer
+ use global_diagnostics
implicit none
type (domain_type), intent(inout) :: domain
integer, intent(in) :: itimestep
real (kind=RKIND), intent(in) :: dt
+ type (block_type), pointer :: block_ptr
call timestep(domain, dt)
+ if (mod(itimestep, config_stats_interval) == 0) then
+ block_ptr => domain % blocklist
+ if(associated(block_ptr % next)) then
+ write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
+ 'that there is only one block per processor.'
+ end if
+
+ call timer_start("global_diagnostics")
+ call computeGlobalDiagnostics(domain % dminfo, &
+ block_ptr % time_levs(2) % state, block_ptr % mesh, &
+ itimestep, dt)
+ call timer_stop("global_diagnostics")
+ end if
+
end subroutine mpas_timestep
</font>
</pre>