<p><b>ringler@lanl.gov</b> 2009-10-03 00:06:44 -0600 (Sat, 03 Oct 2009)</p><p><br>
remedied the issue with pv_edge not being valid along halo cells<br>
by moving all computations related to pv_edge into compute_solve_diagnostics<br>
then adding a halo exchange for pv_edge.<br>
<br>
note: change to module_sw_solver.F is simply to pass in "dt" so that <br>
upwinding can be completed in diagnostics.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/src/module_sw_solver.F
===================================================================
--- trunk/swmodel/src/module_sw_solver.F 2009-10-03 04:43:12 UTC (rev 56)
+++ trunk/swmodel/src/module_sw_solver.F 2009-10-03 06:06:44 UTC (rev 57)
@@ -43,7 +43,7 @@
! 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)
+ call compute_solve_diagnostics(dt,block_ptr % time_levs(1) % state, block_ptr % mesh)
if (.not. config_do_restart) block_ptr % time_levs(1) % state % xtime % scalar = 0.0
block_ptr => block_ptr % next
end do
Modified: trunk/swmodel/src/module_time_integration.F
===================================================================
--- trunk/swmodel/src/module_time_integration.F 2009-10-03 04:43:12 UTC (rev 56)
+++ trunk/swmodel/src/module_time_integration.F 2009-10-03 06:06:44 UTC (rev 57)
@@ -112,7 +112,7 @@
block => domain % blocklist
do while (associated(block))
- call compute_tend(dt, block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
+ call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
block => block % next
end do
@@ -124,6 +124,9 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % u % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -154,7 +157,7 @@
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
- call compute_solve_diagnostics(block % intermediate_step(PROVIS), block % mesh)
+ call compute_solve_diagnostics(dt, block % intermediate_step(PROVIS), block % mesh)
block => block % next
end do
end if
@@ -200,7 +203,7 @@
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)
+ call compute_solve_diagnostics(dt, block % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -208,7 +211,7 @@
end subroutine rk4
- subroutine compute_tend(dt, tend, s, grid)
+ subroutine compute_tend(tend, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -222,7 +225,6 @@
implicit none
- real (kind=RKIND), intent(in) :: dt
type (grid_state), intent(inout) :: tend
type (grid_state), intent(in) :: s
type (grid_meta), intent(in) :: grid
@@ -233,7 +235,7 @@
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, gradPVn, gradPVt, pv_vertex
+ circulation, vorticity, ke, pv_edge
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -245,9 +247,7 @@
circulation => s % circulation % array
vorticity => s % vorticity % array
ke => s % ke % array
- pv_vertex => s % pv_vertex % array
- gradPVn => s % gradPVn % array
- gradPVt => s % gradPVt % array
+ pv_edge => s % pv_edge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -270,7 +270,6 @@
vh => tend % vh % array
tend_h => tend % h % array
tend_u => tend % u % array
- pv_edge => tend % pv_edge % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -306,32 +305,6 @@
#ifdef LANL_FORMULATION
!
- ! Compute pv at the edges
- !
- pv_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,3
- iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
- do k=1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- enddo
- endif
- end do
- end do
-
- ! tdr
- !
- ! Modify PV edge with upstream bias. Bias has contribution from u and v
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
- ! pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
- enddo
- enddo
-
- !
! Compute u (normal) velocity tendency for each edge (cell face)
!
tend_u(:,:) = 0.0
@@ -442,7 +415,7 @@
end subroutine compute_scalar_tend
- subroutine compute_solve_diagnostics(s, grid)
+ subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
!
@@ -453,6 +426,7 @@
implicit none
+ real (kind=RKIND), intent(in) :: dt
type (grid_state), intent(inout) :: s
type (grid_meta), intent(in) :: grid
@@ -612,6 +586,35 @@
! tdr
!
+ ! Compute pv at the edges
+ ! ( this computes pv_edge at all edges bounding real cells )
+ !
+ pv_edge(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,3
+ iEdge = edgesOnVertex(i,iVertex)
+ if(iEdge > 0) then
+ do k=1,nVertLevels
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+ enddo
+ endif
+ end do
+ end do
+ ! tdr
+
+ ! tdr
+ !
+ ! Modify PV edge with upstream bias.
+ !
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ enddo
+ enddo
+
+
+ ! tdr
+ !
! Compute pv at cell centers
! ( this computes pv_cell for all real cells )
!
@@ -644,6 +647,15 @@
enddo
! tdr
+ ! Modify PV edge with upstream bias. (this can not be done because gradPVn is not correct.)
+ !
+ ! do iEdge = 1,nEdges
+ ! do k = 1,nVertLevels
+ ! pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ ! enddo
+ ! enddo
+
+
end subroutine compute_solve_diagnostics
end module time_integration
</font>
</pre>