<p><b>dwj07@fsu.edu</b> 2012-02-06 15:03:56 -0700 (Mon, 06 Feb 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Code cleanup.<br>
Renaming several variables to be more understandable.<br>
<br>
Expanding comments.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-02-06 21:41:29 UTC (rev 1469)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-02-06 22:03:56 UTC (rev 1470)
@@ -40,17 +40,18 @@
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge
- real (kind=RKIND) :: coef_3rd_order, flux_upwind, s_min_update, s_max_update, s_upwind, scale_factor
- real (kind=RKIND) :: flux, tracer_weight
+ 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), dimension(:), pointer :: dvEdge, areaCell
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
- real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, tendency_temp, h_new, s_max, s_min
- real (kind=RKIND), dimension(:,:), pointer :: scale_in, scale_out, horiz_flux, vert_flux
+ real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, 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
coef_3rd_order = config_coef_3rd_order
+ ! Initialize pointers
dvEdge => grid % dvEdge % array
cellsOnEdge => grid % cellsOnEdge % array
cellsOnCell => grid % cellsOnCell % array
@@ -72,159 +73,157 @@
! allocate nCells arrays
allocate(tracer_cur(nVertLevels, nCells))
- allocate(tendency_temp(nVertLevels, nCells))
- allocate(h_new(nVertLevels, nCells))
- allocate(s_max(nVertLevels, nCells))
- allocate(s_min(nVertLevels, nCells))
- allocate(scale_in(nVertLevels, nCells))
- allocate(scale_out(nVertLevels, nCells))
+ allocate(upwind_tendency(nVertLevels, nCells))
+ allocate(inv_h_new(nVertLevels, nCells))
+ allocate(tracer_max(nVertLevels, nCells))
+ allocate(tracer_min(nVertLevels, nCells))
+ allocate(flux_incoming(nVertLevels, nCells))
+ allocate(flux_outgoing(nVertLevels, nCells))
! allocate nEdges arrays
- allocate(horiz_flux(nVertLevels, nEdges))
+ allocate(high_order_horiz_flux(nVertLevels, nEdges))
! allocate nVertLevels+1 and nCells arrays
- allocate(vert_flux(nVertLevels+1, nCells))
+ allocate(high_order_vert_flux(nVertLevels+1, nCells))
do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
- h_new(k, iCell) = h(k, iCell) + dt * tend_h(k, iCell)
+ inv_h_new(k, iCell) = 1.0 / (h(k, iCell) + dt * tend_h(k, iCell))
end do
end do
! Loop over tracers. One tracer is advected at a time. It is copied into a temporary array in order to improve locality
do iTracer = 1, num_tracers
+ ! Initialize variables for use in this iTracer iteration
do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
- tendency_temp(k, iCell) = 0.0
+ upwind_tendency(k, iCell) = 0.0
- s_min(k, iCell) = tracer_cur(k, iCell)
- s_max(k, iCell) = tracer_cur(k, iCell)
+! tracer_min(k, iCell) = tracer_cur(k, iCell)
+! tracer_max(k, iCell) = tracer_cur(k, iCell)
end do ! k loop
end do ! iCell loop
- vert_flux = 0.0
- horiz_flux = 0.0
+ high_order_vert_flux = 0.0
+ high_order_horiz_flux = 0.0
- !
- ! Compute high order vertical flux. Also determine bounds on tracer_cur.
- !
+ ! Compute the high order vertical flux. Also determine bounds on tracer_cur.
do iCell = 1, nCells
k = 1
- vert_flux(k,iCell) = 0.
- s_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- s_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+! high_order_vert_flux(k,iCell) = 0.
+ tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+ tracer_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
k = 2
- vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
- s_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- s_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+ high_order_vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
+ tracer_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+ tracer_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
do k=3,maxLevelCell(iCell)-1
- vert_flux(k,iCell) = flux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
+ high_order_vert_flux(k,iCell) = flux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
tracer_cur(k ,iCell),tracer_cur(k+1,iCell), &
w(k,iCell), coef_3rd_order )
- s_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- s_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+ tracer_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+ tracer_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
end do
k = maxLevelCell(iCell)
- vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
- s_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
- s_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
+ high_order_vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
+ tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
+ tracer_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
- vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
+! high_order_vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
- ! pull s_min and s_max from the (horizontal) surrounding cells
+ ! pull tracer_min and tracer_max from the (horizontal) surrounding cells
do i = 1, nEdgesOnCell(iCell)
do k=1, maxLevelCell(cellsOnCell(i, iCell))
- s_max(k,iCell) = max(s_max(k,iCell),tracer_cur(k, cellsOnCell(i,iCell)))
- s_min(k,iCell) = min(s_min(k,iCell),tracer_cur(k, cellsOnCell(i,iCell)))
+ tracer_max(k,iCell) = max(tracer_max(k,iCell),tracer_cur(k, cellsOnCell(i,iCell)))
+ tracer_min(k,iCell) = min(tracer_min(k,iCell),tracer_cur(k, cellsOnCell(i,iCell)))
end do ! k loop
end do ! i loop over nEdgesOnCell
end do ! iCell Loop
- !
- ! high order horizontal flux
- !
+ ! Compute the high order horizontal flux
do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
do i = 1, nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k = 1, maxLevelCell(iCell)
tracer_weight = uh(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
- horiz_flux(k,iEdge) = horiz_flux(k,iEdge) + tracer_weight* tracer_cur(k,iCell)
+ high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) + tracer_weight* tracer_cur(k,iCell)
end do ! k loop
end do ! i loop over nAdvCellsForEdge
end do ! iEdge loop
- !
! low order upwind vertical flux (monotonic and diffused)
! Remove low order flux from the high order flux.
- ! Store left over high order flux in vert_flux array.
- ! Upwind fluxes are accumulated in tendency_temp
- !
+ ! Store left over high order flux in high_order_vert_flux array.
+ ! Upwind fluxes are accumulated in upwind_tendency
do iCell = 1, nCells
do k = 2, maxLevelCell(iCell)
!dwj: wrong for ocean?
! flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
flux_upwind = min(0.,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.,w(k,iCell))*tracer_cur(k,iCell)
- tendency_temp(k-1,iCell) = tendency_temp(k-1,iCell) + flux_upwind
- tendency_temp(k ,iCell) = tendency_temp(k ,iCell) - flux_upwind
- vert_flux(k,iCell) = vert_flux(k,iCell) - flux_upwind
+ upwind_tendency(k-1,iCell) = upwind_tendency(k-1,iCell) + flux_upwind
+ upwind_tendency(k ,iCell) = upwind_tendency(k ,iCell) - flux_upwind
+ high_order_vert_flux(k,iCell) = high_order_vert_flux(k,iCell) - flux_upwind
end do ! k loop
- ! scale_in contains the total remaining high order flux into iCell
+ ! flux_incoming contains the total remaining high order flux into iCell
! it is positive.
- ! scale_out contains the total remaining high order flux out of iCell
+ ! flux_outgoing contains the total remaining high order flux out of iCell
! it is negative
do k = 1, maxLevelCell(iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical
-! scale_in (k,iCell) = -(min(0.,vert_flux(k+1,iCell))-max(0.,vert_flux(k,iCell)))
-! scale_out(k,iCell) = -(max(0.,vert_flux(k+1,iCell))-min(0.,vert_flux(k,iCell)))
+! flux_incoming (k,iCell) = -(min(0.,high_order_vert_flux(k+1,iCell))-max(0.,high_order_vert_flux(k,iCell)))
+! flux_outgoing(k,iCell) = -(max(0.,high_order_vert_flux(k+1,iCell))-min(0.,high_order_vert_flux(k,iCell)))
- scale_in (k, iCell) = max(0.0, vert_flux(k+1, iCell)) - min(0.0, vert_flux(k, iCell))
- scale_out(k, iCell) = min(0.0, vert_flux(k+1, iCell)) - max(0.0, vert_flux(k, iCell))
+ flux_incoming (k, iCell) = max(0.0, high_order_vert_flux(k+1, iCell)) - min(0.0, high_order_vert_flux(k, iCell))
+ flux_outgoing(k, iCell) = min(0.0, high_order_vert_flux(k+1, iCell)) - max(0.0, high_order_vert_flux(k, iCell))
end do ! k Loop
end do ! iCell Loop
- !
! low order upwind horizontal flux (monotinc and diffused)
! Remove low order flux from the high order flux
- ! Store left over high order flux in horiz_flux array
- ! Upwind fluxes are accumulated in tendency_temp
- !
+ ! Store left over high order flux in high_order_horiz_flux array
+ ! Upwind fluxes are accumulated in upwind_tendency
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
do k = 1, maxLevelEdgeTop(iEdge)
flux_upwind = dvEdge(iEdge) * (max(0.,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.,uh(k,iEdge))*tracer_cur(k,cell2))
- tendency_temp(k,cell1) = tendency_temp(k,cell1) - flux_upwind / areaCell(cell1)
- tendency_temp(k,cell2) = tendency_temp(k,cell2) + flux_upwind / areaCell(cell2)
- horiz_flux(k,iEdge) = horiz_flux(k,iEdge) - flux_upwind
+ high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
+ upwind_tendency(k,cell1) = upwind_tendency(k,cell1) - flux_upwind * invAreaCell1
+ upwind_tendency(k,cell2) = upwind_tendency(k,cell2) + flux_upwind * invAreaCell2
+
! Accumulate remaining high order fluxes
- scale_out(k,cell1) = scale_out(k,cell1) - max(0.,horiz_flux(k,iEdge)) / areaCell(cell1)
- scale_in (k,cell1) = scale_in (k,cell1) - min(0.,horiz_flux(k,iEdge)) / areaCell(cell1)
- scale_out(k,cell2) = scale_out(k,cell2) + min(0.,horiz_flux(k,iEdge)) / areaCell(cell2)
- scale_in (k,cell2) = scale_in (k,cell2) + max(0.,horiz_flux(k,iEdge)) / areaCell(cell2)
+ flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
+ flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
+ flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
+ flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
end do ! k loop
end do ! iEdge loop
- ! Build the factors for the FCT
+ ! Build the factors for the FCT
+ ! Computed using the bounds that were computed previously, and the bounds on the newly updated value
+ ! Factors are placed in the flux_incoming and flux_outgoing arrays
do iCell = 1, nCells
do k = 1, maxLevelCell(iCell)
- s_min_update = (tracer_cur(k,iCell)*h(k,iCell) + dt*(tendency_temp(k,iCell)+scale_out(k,iCell)))/h_new(k,iCell)
- s_max_update = (tracer_cur(k,iCell)*h(k,iCell) + dt*(tendency_temp(k,iCell)+scale_in(k,iCell)))/h_new(k,iCell)
- s_upwind = (tracer_cur(k,iCell)*h(k,iCell) + dt*tendency_temp(k,iCell))/h_new(k,iCell)
+ tracer_min_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*(upwind_tendency(k,iCell)+flux_outgoing(k,iCell))) * inv_h_new(k,iCell)
+ tracer_max_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*(upwind_tendency(k,iCell)+flux_incoming(k,iCell))) * inv_h_new(k,iCell)
+ tracer_upwind_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*upwind_tendency(k,iCell)) * inv_h_new(k,iCell)
- scale_factor = (s_max(k,iCell)-s_upwind)/(s_max_update-s_upwind+eps)
- scale_in(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ scale_factor = (tracer_max(k,iCell)-tracer_upwind_new)/(tracer_max_new-tracer_upwind_new+eps)
+ flux_incoming(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
- scale_factor = (s_upwind-s_min(k,iCell))/(s_upwind-s_min_update+eps)
- scale_out(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ scale_factor = (tracer_upwind_new-tracer_min(k,iCell))/(tracer_upwind_new-tracer_min_new+eps)
+ flux_outgoing(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
end do ! k loop
end do ! iCell loop
@@ -233,23 +232,23 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k = 1, maxLevelEdgeTop(iEdge)
- flux = horiz_flux(k,iEdge)
- flux = max(0.,flux) * min(scale_out(k,cell1), scale_in(k,cell2)) &
- + min(0.,flux) * min(scale_in(k,cell1), scale_out(k,cell2))
- horiz_flux(k,iEdge) = flux
+ flux = high_order_horiz_flux(k,iEdge)
+ flux = max(0.,flux) * min(flux_outgoing(k,cell1), flux_incoming(k,cell2)) &
+ + min(0.,flux) * min(flux_incoming(k,cell1), flux_outgoing(k,cell2))
+ high_order_horiz_flux(k,iEdge) = flux
end do ! k loop
end do ! iEdge loop
! rescale the high order vertical flux
do iCell = 1, nCellsSolve
do k = 2, maxLevelCell(iCell)
- flux = vert_flux(k,iCell)
+ flux = high_order_vert_flux(k,iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical.
-! flux = max(0.,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell)) &
-! + min(0.,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell))
- flux = max(0.,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell)) &
- + min(0.,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell))
- vert_flux(k,iCell) = flux
+! flux = max(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
+! + min(0.,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
+ flux = max(0.,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell)) &
+ + min(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell))
+ high_order_vert_flux(k,iCell) = flux
end do ! k loop
end do ! iCell loop
@@ -257,16 +256,19 @@
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
do k=1,maxLevelEdgeTop(iEdge)
- tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - horiz_flux(k, iEdge)/areaCell(cell1)
- tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + horiz_flux(k, iEdge)/areaCell(cell2)
+ 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
end do ! k loop
end do ! iEdge loop
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(k+1, iCell) - vert_flux(k, iCell)) + tendency_temp(k,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)
end do ! k loop
end do ! iCell loop
@@ -278,14 +280,14 @@
end do ! iTracer loop
deallocate(tracer_cur)
- deallocate(tendency_temp)
- deallocate(h_new)
- deallocate(s_max)
- deallocate(s_min)
- deallocate(scale_in)
- deallocate(scale_out)
- deallocate(horiz_flux)
- deallocate(vert_flux)
+ deallocate(upwind_tendency)
+ deallocate(inv_h_new)
+ deallocate(tracer_max)
+ deallocate(tracer_min)
+ deallocate(flux_incoming)
+ deallocate(flux_outgoing)
+ deallocate(high_order_horiz_flux)
+ deallocate(high_order_vert_flux)
end subroutine mpas_tracer_advection_mono_tend!}}}
end module mpas_tracer_advection_mono
</font>
</pre>