<p><b>mhecht@lanl.gov</b> 2010-06-17 14:14:56 -0600 (Thu, 17 Jun 2010)</p><p>Third order advection should now drop down to second order if land is<br>
found in the higher order stencil. Note that this has only been tested<br>
under test problem 6, which has no land, so have only verified that<br>
code has not been broken.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F        2010-06-17 19:47:05 UTC (rev 359)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F        2010-06-17 20:14:56 UTC (rev 360)
@@ -493,6 +493,10 @@
real (kind=RKIND) :: flux, tracer_edge
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ ! this could be a precomputed 3-d array, trading memory use
+ ! for elimination of conditionals:
+ logical :: near_land
+
!! --- adding pointers here ---
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
! can delete "_new"?
@@ -570,26 +574,47 @@
do k=1,grid % nVertLevels
+ !-- check if land points are in the higher order stencil:
+ near_land = .false.
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ !-- check neighbor, across edge i of cell 1
+ if ( grid % CellsOnCell % array (i,cell1) <= 0) &
+ near_land = .true.
+ end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ !-- check neighbor, across edge i of cell 2
+ if ( grid % CellsOnCell % array (i,cell2) <= 0) &
+ near_land = .true.
+ end do
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
do iTracer=3,grid % nTracers
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracer_new(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracer_new(iTracer,k,cell2)
+ ! if land was found in higher order stencil then leave d2fdx2 terms as zero
+ ! (this is sufficient to cause routine to revert to O(2) flux):
+ if (.not. near_land) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracer_new(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracer_new(iTracer,k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ !-- check neighbor, across edge i of cell 1
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- !-- check neighbor, across edge i of cell 1
- if ( grid % CellsOnCell % array (i,cell1) > 0) &
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ !-- check neighbor, across edge i of cell 2
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+ end if
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- !-- check neighbor, across edge i of cell 2?
- if ( grid % CellsOnCell % array (i,cell2) > 0) &
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
!-- look one way
if (s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) > 0) then
flux = dvEdge(iEdge) * s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * ( &
</font>
</pre>