<p><b>xylar@lanl.gov</b> 2010-05-14 11:06:19 -0600 (Fri, 14 May 2010)</p><p>BRANCH COMMIT<br>
<br>
fixing the vector reconstruction code to ignore (for the time being) the cells that have boundary edges since we don't have enough information to do a reconstruction there.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/z_level_mrp/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/operators/module_vector_reconstruction.F        2010-05-14 14:59:55 UTC (rev 273)
+++ branches/ocean_projects/z_level_mrp/mpas/src/operators/module_vector_reconstruction.F        2010-05-14 17:06:19 UTC (rev 274)
@@ -40,6 +40,7 @@
real (kind=RKIND), allocatable, dimension(:,:) :: mwork
real (kind=RKIND), dimension(:,:), pointer :: matrix, invMatrix
real (kind=RKIND), dimension(:,:), pointer :: normals
+ logical :: hasBoundaryEdges
!========================================================
! arrays filled and saved during init procedure
@@ -85,6 +86,7 @@
rHat = X1
call unit_vector_in_R3(rHat)
+ hasBoundaryEdges = .false.
do j=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(j,iCell)
@@ -101,12 +103,23 @@
X2(3) = zCell(cell2)
normals(:,j) = X1(:) - X2(:)
endif
+ if(cell2 == nCells+1) then
+ ! cell2 is the dummy boundary cell and we don't know where the edge is located
+ ! or what the normal vector should be
+ hasBoundaryEdges = .true.
+ exit
+ end if
call unit_vector_in_R3(normals(:,j))
xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
enddo
+ if(hasBoundaryEdges) then
+ ! we don't have enough information to reconstruct the velocity at this point
+ cycle
+ end if
+
npoints = nEdgesOnCell(iCell) ! only loop over the number of edges for this cell
matrixSize = npoints+2 ! room for 2 vector components for constant flow
! basis functions
</font>
</pre>