<p><b>mhecht@lanl.gov</b> 2010-06-17 16:28:08 -0600 (Thu, 17 Jun 2010)</p><p>Made new scalar_tend code the only code (no more use of dual routines<br>
for tracer numbers le/gt 2). O(3) and O(4) code both drop down to O(2)<br>
when near_land is hardwired true, producing 4output.vtk file which<br>
matches that produced with config_tracer_adv_order of 2. Should still<br>
be cautious, as neither routine has actually been run with land. The<br>
O(4) code still needs visual inspection of results.<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 20:14:56 UTC (rev 360)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F 2010-06-17 22:28:08 UTC (rev 361)
@@ -124,9 +124,7 @@
block => domain % blocklist
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 compute_scalar_tend_le2(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call compute_scalar_tend_gt2(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_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
@@ -385,7 +383,6 @@
end subroutine compute_tend
-
subroutine compute_scalar_tend(tend, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
@@ -401,94 +398,6 @@
type (grid_state), intent(in) :: s
type (grid_meta), intent(in) :: grid
- integer :: iCell, iEdge, k, iTracer, cell1, cell2
- real (kind=RKIND) :: flux, tracer_edge
-
- tend % tracers % array(:,:,:) = 0.0
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
- flux = s % u % array(k,iEdge) * grid % dvEdge % array(iEdge) * s % h_edge % array(k,iEdge) * tracer_edge
- tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) - flux
- tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) + flux
- end do
- end do
- end if
- end do
-
- do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevelsSolve
- do iTracer=1,grid % nTracers
- tend % tracers % array(iTracer,k,iCell) = tend % tracers % array(iTracer,k,iCell) / grid % areaCell % array(iCell)
- end do
- end do
- end do
-
- end subroutine compute_scalar_tend
-
- subroutine compute_scalar_tend_le2(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (grid_state), intent(inout) :: tend
- type (grid_state), intent(in) :: s
- type (grid_meta), intent(in) :: grid
-
- integer :: iCell, iEdge, k, iTracer, cell1, cell2
- real (kind=RKIND) :: flux, tracer_edge
-
- tend % tracers % array(:,:,:) = 0.0
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- do iTracer=1,2
- tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
- flux = s % u % array(k,iEdge) * grid % dvEdge % array(iEdge) * s % h_edge % array(k,iEdge) * tracer_edge
- tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) - flux
- tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) + flux
- end do
- end do
- end if
- end do
-
- do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevelsSolve
- do iTracer=1,2
- tend % tracers % array(iTracer,k,iCell) = tend % tracers % array(iTracer,k,iCell) / grid % areaCell % array(iCell)
- end do
- end do
- end do
-
- end subroutine compute_scalar_tend_le2
-
- subroutine compute_scalar_tend_gt2(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (grid_state), intent(inout) :: tend
- type (grid_state), intent(in) :: s
- type (grid_meta), intent(in) :: grid
-
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
real (kind=RKIND) :: flux, tracer_edge
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
@@ -518,34 +427,8 @@
tracer_tend => tend % tracers % array
!! ----------------------------
-!!$! commenting this out, relying on tend's being zeroed in _le2 routine:
-!!$ tracer_tend(:,:,:) = 0.0
+ tracer_tend(:,:,:) = 0.0
-!!$! this was the original code:
-!!$ do iEdge=1,grid % nEdges
-!!$ cell1 = cellsOnEdge(1,iEdge)
-!!$ cell2 = cellsOnEdge(2,iEdge)
-!!$ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-!!$ do k=1,grid % nVertLevels
-!!$ do iTracer=3,grid % nTracers
-!!$ tracer_edge = 0.5 * (tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2))
-!!$ flux = s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * dvEdge(iEdge) * tracer_edge
-!!$ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux
-!!$ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux
-!!$ end do
-!!$ end do
-!!$ end if
-!!$ end do
-!!$
-!!$ do iCell=1,grid % nCellsSolve
-!!$ do k=1,grid % nVertLevelsSolve
-!!$ do iTracer=3,grid % nTracers
-!!$ tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) / areaCell(iCell)
-!!$ end do
-!!$ end do
-!!$ end do
-
-!!$! here is new (equivalent) code:
if (config_tracer_adv_order == 2) then
do iEdge=1,grid%nEdges
@@ -553,7 +436,7 @@
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
- do iTracer=3,grid % nTracers
+ do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2))
flux = s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * dvEdge(iEdge) * tracer_edge
tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
@@ -591,7 +474,7 @@
d2fdx2_cell1 = 0.0
d2fdx2_cell2 = 0.0
- do iTracer=3,grid % nTracers
+ do iTracer=1,grid % nTracers
! 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
@@ -613,15 +496,15 @@
d2fdx2_cell2 = d2fdx2_cell2 + &
deriv_two(i+1,2,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
end do
- end if
+ end if ! end conditional on near_land
- !-- look one way
+ !-- if u > 0:
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) * ( &
0.5*(tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- look the other way
+ !-- else u <= 0:
else
flux = dvEdge(iEdge) * s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * ( &
0.5*(tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2)) &
@@ -647,19 +530,59 @@
do k=1,grid % nVertLevels
- do iTracer=3,grid % nTracers
+ !-- 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
+
+!!$ !-- I think the O(4) code should work correctly with land,
+!!$ !-- but I haven't verified this. --Matthew Hecht
+!!$ if (near_land) then
+!!$ print *, 'for now, have a stop here,'
+!!$ print *, 'as fourth order code needs a bit more inspection'
+!!$ print *, 'before it is used with land'
+!!$ exit
+!!$ end if
+!!$
+!!$ !-- temporary kluge, which was used to verify that
+!!$ !-- when near_land is fixed to be true
+!!$ !-- this code correctly steps down to O(2) fluxes:
+!!$ near_land = .true.
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+ do iTracer=1,grid % nTracers
+ ! if land was found in higher order stencil then leave d2fdx2 terms as zero:
+ 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)
- if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ !-- 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 ! end conditional on near_land
flux = dvEdge(iEdge) * s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * ( &
0.5*(tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2)) &
@@ -675,7 +598,7 @@
endif
- end subroutine compute_scalar_tend_gt2
+ end subroutine compute_scalar_tend
subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
</font>
</pre>