<p><b>ringler@lanl.gov</b> 2009-10-02 22:43:12 -0600 (Fri, 02 Oct 2009)</p><p><br>
put in code to implement the anticipated potential vorticity method<br>
<br>
note: the code related to "LANL_FORMULATION" is likely not correct<br>
because it requires pv_edge to be valid along the edges of halo cells.<br>
since vorticity is not updated with a halo exhange, pv_vertex along<br>
halo cells in not correct leading to incorrect values of pv_edge.<br>
<br>
work around: run "LANL_FORMULATION" on one proc so there is<br>
no halo. code is correct in this case.<br>
<br>
this error can be corrected by adding a halo exchange<br>
on pv_edge.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/src/module_time_integration.F
===================================================================
--- trunk/swmodel/src/module_time_integration.F 2009-10-02 23:21:51 UTC (rev 55)
+++ trunk/swmodel/src/module_time_integration.F 2009-10-03 04:43:12 UTC (rev 56)
@@ -112,7 +112,7 @@
block => domain % blocklist
do while (associated(block))
- call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
+ call compute_tend(dt, 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
@@ -208,7 +208,7 @@
end subroutine rk4
- subroutine compute_tend(tend, s, grid)
+ subroutine compute_tend(dt, tend, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -222,17 +222,18 @@
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
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, pv_vertex, workpv, q
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, 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
+ circulation, vorticity, ke, pv_edge, gradPVn, gradPVt, pv_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -244,6 +245,9 @@
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
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -302,33 +306,32 @@
#ifdef LANL_FORMULATION
!
- ! Compute height at vertices, pv at vertices, and average pv to edge locations
+ ! Compute pv at the edges
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
- if (cellsOnVertex(1,iVertex) > 0 .and. &
- cellsOnVertex(2,iVertex) > 0 .and. &
- cellsOnVertex(3,iVertex) > 0 .and. &
- edgesOnVertex(1,iVertex) > 0 .and. &
- edgesOnVertex(2,iVertex) > 0 .and. &
- edgesOnVertex(3,iVertex) > 0 &
- ) then
+ do i=1,3
+ iEdge = edgesOnVertex(i,iVertex)
+ if(iEdge > 0) then
do k=1,nVertLevels
- h_vertex = 0.0
- do i=1,3
- h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex = h_vertex / areaTriangle(iVertex)
-
- pv_vertex = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
- do i=1,3
- pv_edge(k,edgesOnVertex(i,iVertex)) = pv_edge(k,edgesOnVertex(i,iVertex)) + 0.5 * pv_vertex
- end do
- end do
- end if
+ 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
@@ -454,13 +457,13 @@
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 :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, 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
+ circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -476,6 +479,10 @@
vorticity => s % vorticity % array
ke => s % ke % array
pv_edge => s % pv_edge % array
+ pv_vertex => s % pv_vertex % array
+ pv_cell => s % pv_cell % array
+ gradPVn => s % gradPVn % array
+ gradPVt => s % gradPVt % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -566,6 +573,77 @@
end do
end do
+
+ ! tdr
+ !
+ ! Compute height at vertices, pv at vertices, and average pv to edge locations
+ ! ( this computes pv_vertex at all vertices bounding real cells )
+ !
+ do iVertex = 1,nVertices
+ if (cellsOnVertex(1,iVertex) > 0 .and. &
+ cellsOnVertex(2,iVertex) > 0 .and. &
+ cellsOnVertex(3,iVertex) > 0 &
+ ) then
+ do k=1,nVertLevels
+ h_vertex = 0.0
+ do i=1,3
+ h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+ end do
+ h_vertex = h_vertex / areaTriangle(iVertex)
+
+ pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ end do
+ end if
+ end do
+ ! tdr
+
+
+ ! tdr
+ !
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells )
+ !
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
+ dvEdge(iEdge)
+ enddo
+ enddo
+
+ ! tdr
+ !
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells )
+ !
+ pv_cell(:,:) = 0.0
+ do iVertex = 1, nVertices
+ do i=1,3
+ iCell = cellsOnVertex(i,iVertex)
+ if( iCell > 0) then
+ do k = 1,nVertLevels
+ pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
+ enddo
+ endif
+ enddo
+ enddo
+ ! tdr
+
+ ! tdr
+ !
+ ! Compute gradient of PV in normal direction
+ ! (tdr: 2009-10-02: this is not correct because the pv_cell in the halo is not correct)
+ !
+ gradPVn(:,:) = 0.0
+ do iEdge = 1,nEdges
+ if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ do k = 1,nVertLevels
+ gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
+ dcEdge(iEdge)
+ enddo
+ endif
+ enddo
+ ! tdr
+
end subroutine compute_solve_diagnostics
end module time_integration
</font>
</pre>