<p><b>gaw06e@fsu.edu</b> 2011-05-28 15:48:34 -0600 (Sat, 28 May 2011)</p><p>- changed compute_tend() to use cellArea approximation at boundaryCells<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F        2011-05-26 21:09:12 UTC (rev 856)
+++ branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F        2011-05-28 21:48:34 UTC (rev 857)
@@ -263,6 +263,7 @@
circulation, vorticity, ke, pv_edge, divergence, h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: boundaryCell
real (kind=RKIND) :: r, u_diffusion
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
@@ -311,6 +312,8 @@
nVertLevels = grid % nVertLevels
u_src => grid % u_src % array
+
+ boundaryCell => grid % boundaryCell % array
!
! Compute height tendency for each cell
@@ -326,8 +329,13 @@
end do
end do
do iCell=1,grid % nCellsSolve
+         if (boundaryCell(1,iCell).eq.1) then
+ r = 1.0 / ( areaCell(iCell) * 2.0 )
+ else
+ r = 1.0 / areaCell(iCell)
+ end if
do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ tend_h(k,iCell) = tend_h(k,iCell) * r
end do
end do
@@ -466,7 +474,11 @@
end do
end do
do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
+ if (boundaryCell(1,iCell).eq.1) then
+ r = 1.0 / ( areaCell(iCell) * 2.0 )
+ else
+ r = 1.0 / areaCell(iCell)
+ end if
do k = 1,nVertLevels
delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
end do
@@ -498,7 +510,7 @@
deallocate(delsq_circulation)
deallocate(delsq_vorticity)
- end if
+ end if ! if (config_h_mom_eddy_visc4 > 0.0)
! Compute u (velocity) tendency from wind stress (u_src)
if(config_wind_stress) then
@@ -508,6 +520,7 @@
end do
endif
+ ! Modify tend_u for bottom drag
if (config_bottom_drag) then
do iEdge=1,grid % nEdges
! bottom drag is the same as POP:
@@ -544,6 +557,7 @@
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
integer, dimension(:,:), pointer :: boundaryEdge
+ integer, dimension(:,:), pointer :: boundaryVertex
real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
@@ -564,6 +578,7 @@
cellsOnEdge => grid % cellsOnEdge % array
boundaryCell=> grid % boundaryCell % array
boundaryEdge=> grid % boundaryEdge % array
+ boundaryVertex => grid % boundaryVertex % array
areaCell => grid % areaCell % array
tracer_tend => tend % tracers % array
@@ -829,6 +844,7 @@
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
h_vertex, vorticity_cell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
+ integer, dimension(:,:), pointer :: boundaryVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
@@ -880,22 +896,25 @@
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
-
+ boundaryVertex => grid % boundaryVertex % array
+
!
! Find those cells that have an edge on the boundary
!
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- enddo
- enddo
+
+! boundaryCell(:,:) = 0
+! do iEdge=1,nEdges
+! do k=1,nVertLevels
+! if(boundaryEdge(k,iEdge).eq.1) then
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! boundaryCell(k,cell1) = 1
+! boundaryCell(k,cell2) = 1
+! endif
+! enddo
+! enddo
+
!
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
@@ -1261,7 +1280,7 @@
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.2) then
tend_u(k,iEdge) = 0.0
endif
</font>
</pre>