<p><b>ringler@lanl.gov</b> 2009-10-27 14:45:53 -0600 (Tue, 27 Oct 2009)</p><p><br>
added option to specify no-slip boundary condition at edge locations<br>
        uBC is added to registry<br>
                uBC == 0 denotes regular velocity points<br>
                uBC == 1 deontes locations where u and v are set equal to zero<br>
                we can extend uBC to have other values in order to enforce other boundary conditions<br>
<br>
        uBC is "iro" and must be included in grid.nc files that contain initial condition data (config_test_case == 0)<br>
<br>
        boundary conditions are enforced at each stage in the RK4 procedure and after the last stage<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/Registry/Registry
===================================================================
--- trunk/swmodel/Registry/Registry        2009-10-27 20:39:38 UTC (rev 66)
+++ trunk/swmodel/Registry/Registry        2009-10-27 20:45:53 UTC (rev 67)
@@ -73,6 +73,9 @@
var real fVertex ( nVertices ) iro fVertex
var real h_s ( nCells ) iro h_s
+# Boundary conditions: read from input, saved in restart and written to output
+var integer uBC ( nVertLevels nEdges ) iro uBC
+
# Prognostic variables: read from input, saved in restart, and written to output
var real u ( nVertLevels nEdges Time ) iro u
var real h ( nVertLevels nCells Time ) iro h
Modified: trunk/swmodel/src/module_time_integration.F
===================================================================
--- trunk/swmodel/src/module_time_integration.F        2009-10-27 20:39:38 UTC (rev 66)
+++ trunk/swmodel/src/module_time_integration.F        2009-10-27 20:45:53 UTC (rev 67)
@@ -124,6 +124,7 @@
do while (associated(block))
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)
+ call enforce_uBC(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
block => block % next
end do
@@ -211,6 +212,7 @@
end if
call compute_solve_diagnostics(dt, block % time_levs(2) % state, block % mesh)
+ call enforce_uBC(block % time_levs(2) % state, block % time_levs(1) % state, block % mesh)
block => block % next
end do
@@ -461,9 +463,9 @@
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
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r
+ real (kind=RKIND) :: r, h1, h2
h => s % h % array
@@ -506,6 +508,8 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ uBC => grid % uBC % array
+
!
! Compute height on cell edges at velocity locations
!
@@ -705,6 +709,76 @@
enddo
enddo
+ !
+ ! set pv_edge = fEdge / h_edge at boundary points
+ !
+ if (maxval(uBC).gt.0) then
+ do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k = 1,nVertLevels
+ if(uBC(k,iEdge).eq.1) then
+ if(cell1.gt.0) h1 = h(k,cell1)
+ if(cell2.gt.0) h2 = h(k,cell2)
+ pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
+ endif
+ enddo
+ enddo
+ endif
+
+
end subroutine compute_solve_diagnostics
+ subroutine enforce_uBC(tend, s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Enforce any boundary conditions on the normal velocity at each edge
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: tend_u set to zero at uBC == 1 locations
+ ! u set to zero at uBC == 1 locations
+ ! v set to zero at uBC == 1 locations
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ implicit none
+
+ type (grid_state), intent(inout) :: tend
+ type (grid_state), intent(inout) :: s
+ type (grid_meta), intent(in) :: grid
+
+ integer, dimension(:,:), pointer :: uBC
+ real (kind=RKIND), dimension(:,:), pointer :: u, tend_u, v, tend_v
+ integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: iEdge, k
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ uBC => grid % uBC % array
+ u => s % u % array
+ v => s % v % array
+ tend_u => tend % u % array
+ tend_v => tend % v % array
+
+ if(maxval(uBC).le.0) return
+
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+
+ if(uBC(k,iEdge).eq.1) then
+ u(k,iEdge) = 0.0
+ v(k,iEdge) = 0.0
+ tend_u(k,iEdge) = 0.0
+ tend_v(k,iEdge) = 0.0
+ endif
+
+ enddo
+ enddo
+
+ end subroutine enforce_uBC
+
+
end module time_integration
</font>
</pre>