<p><b>mhoffman@lanl.gov</b> 2013-02-05 13:29:17 -0700 (Tue, 05 Feb 2013)</p><p>Branch commit - land ice<br>
Adding output of global diagnostics for the initial time level. By default this is turned on, but can be disabled with the config_write_initial_stats namelist option.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry        2013-02-05 17:35:26 UTC (rev 2435)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry        2013-02-05 20:29:17 UTC (rev 2436)
@@ -10,8 +10,10 @@
namelist character land_ice_model config_start_time 0000-01-01_00:00:00
namelist character land_ice_model config_stop_time none
namelist character land_ice_model config_run_duration none
-namelist integer land_ice_model config_stats_interval 100
+namelist integer land_ice_model config_stats_interval 1
% config_stats_interval is how many time steps before calculating & outputting stats via the global_diagnostics module (which has not yet been setup to do much)
+namelist logical land_ice_model config_write_initial_stats true
+% config_write_initial_stats determines if stats are written at the initial time
% Land ice core-specific options:
namelist character land_ice_model config_dycore SIA
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F        2013-02-05 17:35:26 UTC (rev 2435)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F        2013-02-05 20:29:17 UTC (rev 2436)
@@ -17,7 +17,7 @@
subroutine land_ice_compute_global_diagnostics(dminfo, state, grid, timeIndex, dt)
! Note: this routine assumes that there is only one block per processor. No looping
- ! is performed over blocks.
+ ! is performed over blocks. TODO
! 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
@@ -128,12 +128,12 @@
allocate(cellVolume(nVertLevels,nCellsSolve))
allocate(iceArea(nCells+1))
allocate(cellArea(nVertLevels,nCellsSolve))
- allocate(refAreaWeightedSurfaceHeight(nVertLevels,nCellsSolve))
- allocate(refAreaWeightedSurfaceHeight_edge(nVertLevels,nEdgesSolve))
+ !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(vertexVolume(nVertLevels,nVerticesSolve))
!allocate(volumeWeightedKineticEnergy(nVertLevels,nEdgesSolve))
!allocate(volumeWeightedPotentialEnergy(nVertLevels,nCellsSolve))
!allocate(volumeWeightedPotentialEnergyTopography(nVertLevels,nCellsSolve))
@@ -142,18 +142,18 @@
!allocate(keTend_PressureGradient(nVertLevels,nEdgesSolve))
!allocate(peTend_DivThickness(nVertLevels,nCells))
- allocate(averageThickness(nCellsSolve))
+ !allocate(averageThickness(nCellsSolve))
!!! allocate(h_s_edge(nEdgesSOlve))
cellVolume = 0.0_RKIND
iceArea = 0.0_RKIND
- refAreaWeightedSurfaceHeight = 0.0_RKIND
- refAreaWeightedSurfaceHeight_edge = 0.0_RKIND
- vertexVolume = 0.0_RKIND
+ !refAreaWeightedSurfaceHeight = 0.0_RKIND
+ !refAreaWeightedSurfaceHeight_edge = 0.0_RKIND
+ !vertexVolume = 0.0_RKIND
cellArea = 0.0_RKIND
- averageThickness = 0.0_RKIND
+ !averageThickness = 0.0_RKIND
!volumeWeightedPotentialVorticity = 0
!volumeWeightedPotentialEnstrophy = 0
!volumeWeightedKineticEnergy = 0
@@ -184,7 +184,7 @@
!!! *layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
!!! volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
!!! *pv_vertex(iLevel,1:nVerticesSolve)*layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- vertexVolume(iLevel,:) = layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+ !vertexVolume(iLevel,:) = layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
!!! volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &
!!! *layerThicknessEdge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
!!! volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
@@ -246,14 +246,14 @@
! Compute Average Sea Surface Height for Potential Energy and Enstrophy
! Reservoir computations
- call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, refAreaWeightedSurfaceHeight, sumrefAreaWeightedSurfaceHeight)
+ !call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, refAreaWeightedSurfaceHeight, sumrefAreaWeightedSurfaceHeight)
!!! averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
- call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
+ !call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
!!! globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
!!! globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
@@ -293,8 +293,10 @@
if (dminfo % my_proc_id == IO_NODE) then
fileID = land_ice_get_free_unit()
- if (timeIndex/config_stats_interval == 1) then
+ if (config_write_initial_stats .and. (timeIndex == 0)) then
open(fileID, file='GlobalIntegrals.txt',STATUS='unknown')
+ elseif ( .not.(config_write_initial_stats) .and. (timeIndex/config_stats_interval == 1) ) then
+ open(fileID, file='GlobalIntegrals.txt',STATUS='unknown')
else
open(fileID, file='GlobalIntegrals.txt',POSITION='append')
endif
@@ -308,7 +310,12 @@
close(fileID)
end if
+
+ deallocate(cellVolume)
+ deallocate(iceArea)
+ deallocate(cellArea)
deallocate(areaEdge)
+
end subroutine land_ice_compute_global_diagnostics
integer function land_ice_get_free_unit()
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F        2013-02-05 17:35:26 UTC (rev 2435)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F        2013-02-05 20:29:17 UTC (rev 2436)
@@ -23,6 +23,7 @@
use mpas_constants
use land_ice_vel
use mpas_ocn_tracer_advection
+ use land_ice_global_diagnostics
implicit none
@@ -39,6 +40,10 @@
!
! Initialize core
!
+
+ write(0,*) 'Initial time ', startTimeStamp
+ write(6,*) 'Initial time ', startTimeStamp
+
! If dt in years is supplied, use it. Otherwise use dt in seconds
if (config_dt_years .eq. 0.0_RKIND) then
dt = config_dt_seconds
@@ -75,6 +80,11 @@
call mpas_copy_state(block % state % time_levs(i) % state, block % state % time_levs(1) % state)
end do
+ ! TODO Global diagnostics needs to be modified to work on blocks! Until then, just assuming one block
+ if (config_write_initial_stats) then
+ call land_ice_compute_global_diagnostics(domain % dminfo, block % state % time_levs(1) % state, block % mesh, 0, dt)
+ endif
+
block => block % next
end do
@@ -278,8 +288,6 @@
call mpas_timer_start("land ice core run")
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Initial timestep ', timeStamp
- write(6,*) 'Initial timestep ', timeStamp
call mpas_timer_start("write output frame")
call write_output_frame(output_obj, output_frame, domain)
</font>
</pre>