<p><b>dwj07@fsu.edu</b> 2012-02-21 15:49:18 -0700 (Tue, 21 Feb 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding a temporary bounds check. It's commented out by default. <br>
To enable it, uncomment the block below the comment<br>
<br>
do bounds check<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-02-21 21:49:24 UTC (rev 1525)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-02-21 22:49:18 UTC (rev 1526)
@@ -42,9 +42,10 @@
real (kind=RKIND) :: coef_3rd_order, flux_upwind, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor
real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
+ real (kind=RKIND) :: cur_max, cur_min, new_max, new_min
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
- real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, upwind_tendency, inv_h_new, tracer_max, tracer_min
+ real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, tracer_new, upwind_tendency, inv_h_new, tracer_max, tracer_min
real (kind=RKIND), dimension(:,:), pointer :: flux_incoming, flux_outgoing, high_order_horiz_flux, high_order_vert_flux
real (kind=RKIND), parameter :: eps = 1.e-20
@@ -73,6 +74,7 @@
! allocate nCells arrays
allocate(tracer_cur(nVertLevels, nCells))
+ allocate(tracer_new(nVertLevels, nCells))
allocate(upwind_tendency(nVertLevels, nCells))
allocate(inv_h_new(nVertLevels, nCells))
allocate(tracer_max(nVertLevels, nCells))
@@ -101,6 +103,9 @@
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
upwind_tendency(k, iCell) = 0.0
+ !dwj 02/21/12 tracer_new is supposed to be the "new" tracer state. This allows bounds checks.
+ tracer_new(k,iCell) = 0.0
+
! tracer_min(k, iCell) = tracer_cur(k, iCell)
! tracer_max(k, iCell) = tracer_cur(k, iCell)
end do ! k loop
@@ -262,6 +267,10 @@
do k=1,maxLevelEdgeTop(iEdge)
tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+
+ !dwj 02/21/12 tracer_new holds a tendency for now.
+ tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
+ tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
end do ! k loop
end do ! iEdge loop
@@ -269,9 +278,30 @@
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
+
+ !dwj 02/21/12 tracer_new holds a tendency for now.
+ tracer_new(k, iCell) = tracer_new(k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
+
+ !dwj 02/21/12 now, updated tracer_new to be the new state of the tracer.
+ tracer_new(k, iCell) = (tracer_cur(k, iCell)*h(k, iCell) + dt * tracer_new(k, iCell)) * inv_h_new(k, iCell)
end do ! k loop
end do ! iCell loop
+ !dwj 02/21/12 build min and max of cur and new tracer states to check bounds
+ cur_min = minval(tracer_cur)
+ cur_max = maxval(tracer_cur)
+ new_min = minval(tracer_new)
+ new_max = maxval(tracer_new)
+
+ !dwj 02/21/12 do bounds check
+! if(new_min < cur_min ) then
+! write(*,*) 'Minimum out of bounds on tracer ', iTracer, cur_min, new_min
+! end if
+
+! if(new_max > cur_max) then
+! write(*,*) 'Maximum out of bounds on tracer ', iTracer, cur_max, new_max
+! end if
+
! do iCell = 1, grid%nCells
! do k=1, grid%maxLevelCell(iCell)
! tracer_next_in(iTracer,k,iCell) = max(0.,tracer_next(k,iCell))
@@ -280,6 +310,7 @@
end do ! iTracer loop
deallocate(tracer_cur)
+ deallocate(tracer_new)
deallocate(upwind_tendency)
deallocate(inv_h_new)
deallocate(tracer_max)
</font>
</pre>