<p><b>dwj07@fsu.edu</b> 2012-02-03 14:33:51 -0700 (Fri, 03 Feb 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
More clean up for monotonic advection.<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-03 19:28:49 UTC (rev 1464)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-02-03 21:33:51 UTC (rev 1465)
@@ -23,49 +23,40 @@
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: uh, w, h, tend_h
- real (kind=RKIND), intent(in) :: dt
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: current tracer values
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh !< Input: Thichness weighted velocitiy
+ real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical velocitiy
+ real (kind=RKIND), dimension(:,:), intent(in) :: h !< Input: Thickness
+ real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Tendency for thickness field
+ real (kind=RKIND), intent(in) :: dt !< Input: Timestep
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance !< Input: Distance between vertical interfaces
+ real (kind=RKIND), dimension(:), intent(in) :: zWeightK !< Input: Weight for tracers to map to edge, weight from K cell
+ real (kind=RKIND), dimension(:), intent(in) :: zWeightKm1 !< Input: Weight for tracers to map to edge, weight from K-1 cell
type (mesh_type), intent(in) :: grid
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
+ 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), dimension(:), pointer :: dvEdge, areaCell
- integer, dimension(:,:), pointer :: cellsOnEdge
-
- integer, dimension(:,:), pointer :: advCellsForEdge
- integer, dimension(:), pointer :: nAdvCellsForEdge
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( grid % nVertLevels, grid % nCells ) :: tracer_cur, tendency_temp, h_new
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: s_max, s_min
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: scale_in, scale_out
-
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nEdges ) :: horiz_flux
- real (kind=RKIND), dimension( grid % nVertLevels + 1, grid % nCells ) :: vert_flux
-
- integer :: nVertLevels, num_tracers, icellmax, kmax
-
- real (kind=RKIND) :: coef_3rd_order
-
- real (kind=RKIND) :: tracer_turb_flux, z1,z2,z3,z4,zm,z0,zp
-
- real (kind=RKIND) :: flux_upwind
- real (kind=RKIND) :: s_min_update, s_max_update, s_upwind, scale_factor
-
- integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop
-
real (kind=RKIND), parameter :: eps = 1.e-20
coef_3rd_order = config_coef_3rd_order
dvEdge => grid % dvEdge % array
cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnCell => grid % cellsOnCell % array
areaCell => grid % areaCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
@@ -73,22 +64,45 @@
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
-
num_tracers = size(tracers,dim=1)
+ ! 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 nEdges arrays
+ allocate(horiz_flux(nVertLevels, nEdges))
+
+ ! allocate nVertLevels+1 and nCells arrays
+ allocate(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)
+ 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
- do iCell = 1, grid%nCells
+ do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
tendency_temp(k, iCell) = 0.0
- h_new(k, iCell) = h(k, iCell) + dt * tend_h(k, iCell)
s_min(k, iCell) = tracer_cur(k, iCell)
s_max(k, iCell) = tracer_cur(k, iCell)
- end do
- end do
+ end do ! k loop
+ end do ! iCell loop
vert_flux = 0.0
horiz_flux = 0.0
@@ -96,7 +110,7 @@
!
! Compute high order vertical flux. Also determine bounds on tracer_cur.
!
- do iCell=1,grid % nCells
+ 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))
@@ -123,29 +137,28 @@
vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
! pull s_min and s_max from the (horizontal) surrounding cells
- do i=1, grid % nEdgesOnCell % array(iCell)
- do k=1, maxLevelCell( grid % cellsOnCell % array(i, iCell))
- s_max(k,iCell) = max(s_max(k,iCell),tracer_cur(k, grid % CellsOnCell % array(i,iCell)))
- s_min(k,iCell) = min(s_min(k,iCell),tracer_cur(k, grid % CellsOnCell % array(i,iCell)))
- end do
- end do
- end do
+ 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)))
+ end do ! k loop
+ end do ! i loop over nEdgesOnCell
+ end do ! iCell Loop
!
! high order horizontal flux
!
- horiz_flux(:,:) = 0.
- do iEdge=1,grid%nEdges
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do i=1,nAdvCellsForEdge(iEdge)
+ do i = 1, nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
- do k=1,maxLevelCell(iCell)
+ 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)
- end do
- end do
- end do
+ end do ! k loop
+ end do ! i loop over nAdvCellsForEdge
+ end do ! iEdge loop
!
! low order upwind vertical flux (monotonic and diffused)
@@ -153,7 +166,7 @@
! Store left over high order flux in vert_flux array.
! Upwind fluxes are accumulated in tendency_temp
!
- do iCell = 1, grid % nCells
+ 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)
@@ -161,21 +174,21 @@
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
- end do
+ end do ! k loop
! scale_in contains the total remaining high order flux into iCell
! it is positive.
! scale_out contains the total remaining high order flux out of iCell
! it is negative
- do k=1,maxLevelCell(iCell)
+ 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)))
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))
- end do
- end do
+ end do ! k Loop
+ end do ! iCell Loop
!
! low order upwind horizontal flux (monotinc and diffused)
@@ -183,10 +196,10 @@
! Store left over high order flux in horiz_flux array
! Upwind fluxes are accumulated in tendency_temp
!
- do iEdge=1,grid%nEdges
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+ 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)
@@ -197,11 +210,11 @@
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)
- end do
- end do
+ end do ! k loop
+ end do ! iEdge loop
! Build the factors for the FCT
- do iCell = 1, grid % nCells
+ 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)
@@ -212,24 +225,23 @@
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) )
+ end do ! k loop
+ end do ! iCell loop
- end do
- end do
-
! rescale the high order horizontal fluxes
- do iEdge = 1, grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge = 1, nEdges
+ 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
- end do
- end do
+ end do ! k loop
+ end do ! iEdge loop
! rescale the high order vertical flux
- do iCell=1,grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k = 2, maxLevelCell(iCell)
flux = vert_flux(k,iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical.
@@ -238,33 +250,42 @@
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
- end do
- end do
+ end do ! k loop
+ end do ! iCell loop
! Accumulate the scaled high order horizontal tendencies
- do iEdge=1,grid%nEdges
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
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)
- end do
- end do
+ end do ! k loop
+ end do ! iEdge loop
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
- do iCell=1,grid % nCellsSolve
- do k=1,maxLevelCell(iCell)
+ 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)
- end do
- end do
+ end do ! k loop
+ end do ! iCell loop
! do iCell = 1, grid%nCells
! do k=1, grid%maxLevelCell(iCell)
! tracer_next_in(iTracer,k,iCell) = max(0.,tracer_next(k,iCell))
-! end do
-! end do
+! end do ! k loop
+! end do ! iCell loop
+ end do ! iTracer loop
- end do ! loop over tracers
+ 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)
end subroutine mpas_tracer_advection_mono_tend!}}}
end module mpas_tracer_advection_mono
</font>
</pre>