<p><b>duda</b> 2009-08-18 10:24:30 -0600 (Tue, 18 Aug 2009)</p><p>Move computation of diagnostic fields used in tendency computation for<br>
h and u out of compute_tend and into a new subroutine<br>
compute_solve_diagnostics. Now, all diagnostic fields used in the<br>
solver should be computed in compute_solve_diagnostics, and all fields<br>
needed only in the model output can be computed in<br>
compute_output_diagnostics, which is called from the main solver loop<br>
in module_sw_solver.<br>
<br>
This change permits the diagnostic fields to appear in the output.nc<br>
file, including at the initial time, while only defining the<br>
computation of these fields in once.<br>
<br>
M module_time_integration.F<br>
M module_sw_solver.F<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/module_sw_solver.F
===================================================================
--- trunk/swmodel/module_sw_solver.F 2009-08-18 00:10:26 UTC (rev 31)
+++ trunk/swmodel/module_sw_solver.F 2009-08-18 16:24:30 UTC (rev 32)
@@ -37,9 +37,11 @@
ntimesteps = config_ntimesteps
- ! Initialize simulation time to 0 for all blocks
+ ! Compute diagnostic fields needed in solve loop, and initialize
+ ! simulation time to 0 for all blocks
block_ptr => domain % blocklist
do while (associated(block_ptr))
+ call compute_solve_diagnostics(block_ptr % time_levs(1) % state, block_ptr % mesh)
block_ptr % time_levs(1) % state % time = 0.0
block_ptr => block_ptr % next
end do
@@ -90,7 +92,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
- call compute_diagnostics(block_ptr % time_levs(1) % state, block_ptr % mesh)
+ call compute_output_diagnostics(block_ptr % time_levs(1) % state, block_ptr % mesh)
block_ptr => block_ptr % next
end do
@@ -100,7 +102,7 @@
end subroutine write_output_frame
- subroutine compute_diagnostics(state, grid)
+ subroutine compute_output_diagnostics(state, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields for a domain
!
@@ -118,23 +120,7 @@
integer :: i, eoe
integer :: iEdge, k
- !
- ! Compute v (tangential) velocities
- !
- state % v % array(:,:) = 0.0
- do iEdge = 1, grid % nEdgesSolve
- do i=1,grid % nEdgesOnEdge % array(iEdge)
- eoe = grid % edgesOnEdge % array(i,iEdge)
- do k = 1,grid % nVertLevelsSolve
+ end subroutine compute_output_diagnostics
- state % v % array(k,iEdge) = state % v % array(k,iEdge) + &
- grid % weightsOnEdge % array(i,iEdge) * state % u % array(k, eoe)
- end do
- end do
- end do
-
- end subroutine compute_diagnostics
-
-
end module sw_solver
Modified: trunk/swmodel/module_time_integration.F
===================================================================
--- trunk/swmodel/module_time_integration.F 2009-08-18 00:10:26 UTC (rev 31)
+++ trunk/swmodel/module_time_integration.F 2009-08-18 16:24:30 UTC (rev 32)
@@ -72,14 +72,12 @@
! Take first RK step
block => domain % blocklist
do while (associated(block))
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:)
- block % intermediate_step(TEMP) % tracers % array(:,:,:) = block % time_levs(1) % state % tracers % array(:,:,:)
+ ! Copy fEdge, fVertex, and h_s into provisional array for use in 2nd, 3rd, and 4th RK intermediate steps
block % intermediate_step(TEMP) % fEdge % array(:) = block % time_levs(1) % state % fEdge % array(:)
block % intermediate_step(TEMP) % fVertex % array(:) = block % time_levs(1) % state % fVertex % array(:)
block % intermediate_step(TEMP) % h_s % array(:) = block % time_levs(1) % state % h_s % array(:)
- call compute_tend(block % intermediate_step(Q1), block % intermediate_step(TEMP), block % mesh)
- call compute_scalar_tend(block % intermediate_step(Q1), block % intermediate_step(TEMP), block % mesh)
+ call compute_tend(block % intermediate_step(Q1), block % time_levs(1) % state, block % mesh)
+ call compute_scalar_tend(block % intermediate_step(Q1), block % time_levs(1) % state, block % mesh)
block % intermediate_step(Q1) % u % array(:,:) = dt * block % intermediate_step(Q1) % u % array(:,:)
block % intermediate_step(Q1) % h % array(:,:) = dt * block % intermediate_step(Q1) % h % array(:,:)
block % intermediate_step(Q1) % tracers % array(:,:,:) = dt * block % intermediate_step(Q1) % tracers % array(:,:,:)
@@ -121,6 +119,7 @@
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
+ call compute_solve_diagnostics(block % intermediate_step(TEMP), block % mesh)
call compute_tend(block % intermediate_step(Q2), block % intermediate_step(TEMP), block % mesh)
call compute_scalar_tend(block % intermediate_step(Q2), block % intermediate_step(TEMP), block % mesh)
block % intermediate_step(Q2) % u % array(:,:) = dt * block % intermediate_step(Q2) % u % array(:,:)
@@ -164,6 +163,7 @@
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
+ call compute_solve_diagnostics(block % intermediate_step(TEMP), block % mesh)
call compute_tend(block % intermediate_step(Q3), block % intermediate_step(TEMP), block % mesh)
call compute_scalar_tend(block % intermediate_step(Q3), block % intermediate_step(TEMP), block % mesh)
block % intermediate_step(Q3) % u % array(:,:) = dt * block % intermediate_step(Q3) % u % array(:,:)
@@ -207,6 +207,7 @@
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
+ call compute_solve_diagnostics(block % intermediate_step(TEMP), block % mesh)
call compute_tend(block % intermediate_step(Q4), block % intermediate_step(TEMP), block % mesh)
call compute_scalar_tend(block % intermediate_step(Q4), block % intermediate_step(TEMP), block % mesh)
block % intermediate_step(Q4) % u % array(:,:) = dt * block % intermediate_step(Q4) % u % array(:,:)
@@ -264,6 +265,8 @@
block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
+ call compute_solve_diagnostics(block % time_levs(2) % state, block % mesh)
+
block => block % next
end do
@@ -293,17 +296,22 @@
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, tend_h, tend_u, &
+ real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- h_s => s % h_s % array
- fVertex => s % fVertex % array
- fEdge => s % fEdge % array
- h => s % h % array
- u => s % u % array
+ h_s => s % h_s % array
+ fVertex => s % fVertex % array
+ fEdge => s % fEdge % array
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ ke => s % ke % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -321,12 +329,8 @@
areaTriangle => grid % areaTriangle % array
vh => tend % vh % array
- h_edge => tend % h_edge % array
tend_h => tend % h % array
tend_u => tend % u % array
- circulation => tend % circulation % array
- vorticity => tend % vorticity % array
- ke => tend % ke % array
pv_edge => tend % pv_edge % array
nCells = grid % nCells
@@ -336,19 +340,6 @@
!
- ! Compute height on cell edges at velocity locations
- !
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
- end do
-
- !
! Compute height tendency for each cell
!
tend_h(:,:) = 0.0
@@ -374,44 +365,6 @@
end do
end do
- !
- ! Compute circulation and relative vorticity at each vertex
- !
- circulation(:,:) = 0.0
- do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
- if (verticesOnEdge(2,iEdge) > 0) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
- end do
- do iVertex=1,nVertices
- do k=1,nVertLevels
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
- end do
-
- !
- ! Compute kinetic energy in each cell
- !
- ke(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- do k=1,nVertLevels
- ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- end do
- end do
- do k=1,nVertLevels
- ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- end do
- end do
-
#ifdef LANL_FORMULATION
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
@@ -550,4 +503,132 @@
end subroutine compute_scalar_tend
+
+ subroutine compute_solve_diagnostics(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(inout) :: s
+ 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, pv_vertex, workpv, q
+
+ 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
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+
+
+ h_s => s % h_s % array
+ fVertex => s % fVertex % array
+ fEdge => s % fEdge % array
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ vh => s % vh % array
+ h_edge => s % h_edge % array
+ tend_h => s % h % array
+ tend_u => s % u % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ ke => s % ke % array
+ pv_edge => s % pv_edge % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ !
+ ! Compute height on cell edges at velocity locations
+ !
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
+ end if
+ end do
+
+ !
+ ! Compute circulation and relative vorticity at each vertex
+ !
+ circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ if (verticesOnEdge(1,iEdge) > 0) then
+ do k=1,nVertLevels
+ circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
+ end do
+ end if
+ if (verticesOnEdge(2,iEdge) > 0) then
+ do k=1,nVertLevels
+ circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
+ end do
+ end if
+ end do
+ do iVertex=1,nVertices
+ do k=1,nVertLevels
+ vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+ end do
+ end do
+
+ !
+ ! Compute kinetic energy in each cell
+ !
+ ke(:,:) = 0.0
+ do iCell=1,nCells
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ do k=1,nVertLevels
+ ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ end do
+ end do
+ do k=1,nVertLevels
+ ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+ !
+ ! Compute v (tangential) velocities
+ !
+ v(:,:) = 0.0
+ do iEdge = 1,nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ do k = 1,nVertLevels
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
+ end do
+ end do
+
+ end subroutine compute_solve_diagnostics
+
end module time_integration
</font>
</pre>