<p><b>mhecht@lanl.gov</b> 2010-04-21 09:47:27 -0600 (Wed, 21 Apr 2010)</p><p>Successful port of O(2) advection code from hyd_atmos core to sw, in<br>
preparation for less trivial porting of higher order routines. This<br>
version corresponds to version 1.10 in Matthew's local cvs repo.<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-04-21 15:42:44 UTC (rev 194)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F 2010-04-21 15:47:27 UTC (rev 195)
@@ -124,7 +124,10 @@
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(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 enforce_uBC(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -426,7 +429,148 @@
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 > 0 .and. cell2 > 0) 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 :: iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_edge
+
+ !! adding pointers here
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell ! for areaCell, dvEdge
+ ! can delete "_new"?
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracer_new, tracer_tend ! for tracer_tend
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ dvEdge => grid % dvEdge % array
+ tracer_new => s % tracers % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ areaCell => grid % areaCell % array
+ tracer_tend => tend % tracers % array
+ !! --------------------
+
+
+!! okay to comment this out, relying on tend's being zeroed in _le2 routine?
+!!$ tracer_tend(:,:,:) = 0.0
+
+!!$! this was the original code:
+!!$ do iEdge=1,grid % nEdges
+!!$ cell1 = cellsOnEdge(1,iEdge)
+!!$ cell2 = cellsOnEdge(2,iEdge)
+!!$ if (cell1 > 0 .and. cell2 > 0) 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 code:
+!!$ if (config_tracer_adv_order == 2) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) 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/areaCell(cell1)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ end do
+ end do
+ end if
+ end do
+
+!!$ do iEdge=1,grid%nEdges
+!!$ cell1 = cellsOnEdge(1,iEdge)
+!!$ cell2 = cellsOnEdge(2,iEdge)
+!!$ if (cell1 > 0 .and. cell2 > 0) 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/areaCell(cell1)
+!!$ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+!!$ end do
+!!$ end do
+!!$ end if
+!!$ end do
+
+!!$ else if (config_tracer_adv_order == 3) then
+!!$ endif
+
+ end subroutine compute_scalar_tend_gt2
+
+
subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
</font>
</pre>