<p><b>dwj07@fsu.edu</b> 2012-08-27 11:16:34 -0600 (Mon, 27 Aug 2012)</p><p><br>
        -- BRANCH COMMIT --<br>
<br>
        Committing forgotten shared advection modules.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_helpers.F
===================================================================
--- branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_helpers.F                                (rev 0)
+++ branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_helpers.F        2012-08-27 17:16:34 UTC (rev 2126)
@@ -0,0 +1,279 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  mpas_tracer_advection_helpers
+!
+!&gt; \brief MPAS tracer advection helper functions
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This module contains helper functions tracer advection.
+!
+!-----------------------------------------------------------------------
+module mpas_tracer_advection_helpers
+
+   use mpas_kind_types
+   use mpas_configure
+
+   implicit none
+   save
+
+   contains
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  routine mpas_tracer_advection_coefficients
+!
+!&gt; \brief MPAS tracer advection coefficients
+!&gt; \author Doug Jacobsen, Bill Skamarock
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This routine precomputes the advection coefficients for horizontal
+!&gt;  advection of tracers.
+!
+!-----------------------------------------------------------------------
+   subroutine mpas_tracer_advection_coefficients( grid, err, maxLevelCell_in, highOrderAdvectionMask_in, boundaryCell_in )!{{{
+
+      implicit none
+      type (mesh_type) :: grid !&lt; Input: Grid information
+      integer, intent(out) :: err !&lt; Input/Output: Error flag
+      integer, dimension(:), pointer, optional :: maxLevelCell_in !&lt; Input - optional: Index to last real cell
+      integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !&lt; Input - optional: Mask for high order advection
+      integer, dimension(:,:), pointer, optional :: boundaryCell_in !&lt; Input- optional: Mask for boundary cells
+
+      real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+      real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
+      integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask, boundaryCell
+      integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge, maxLevelCell
+
+      integer, dimension(:), pointer :: cell_indices 
+      integer, dimension(:,:), pointer :: sorted_cell_indices
+      integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell, k, nVertLevels
+      logical :: addcell, highOrderAdvection
+
+      type (hashtable) :: cell_hash
+
+      deriv_two =&gt; grid % deriv_two % array
+      adv_coefs =&gt; grid % adv_coefs % array
+      adv_coefs_2nd =&gt; grid % adv_coefs_2nd % array
+      adv_coefs_3rd =&gt; grid % adv_coefs_3rd % array
+      cellsOnCell =&gt; grid % cellsOnCell % array
+      cellsOnEdge =&gt; grid % cellsOnEdge % array
+      advCellsForEdge =&gt; grid % advCellsForEdge % array
+      nEdgesOnCell =&gt; grid % nEdgesOnCell % array
+      nAdvCellsForEdge =&gt; grid % nAdvCellsForEdge % array
+
+      nVertLevels = grid % nVertLevels
+
+      allocate(cell_indices(grid % maxEdges2 + 2))
+      allocate(sorted_cell_indices(2, grid % maxEdges2 + 2))
+
+      err = 0
+
+      if(present(maxLevelCell_in)) then
+        maxLevelCell =&gt; maxLevelCell_in
+      else
+        allocate(maxLevelCell(grid % nCells+1))
+        maxLevelCell(:) = nVertLevels
+      end if
+
+      if(present(highOrderAdvectionMask_in)) then
+        highOrderAdvectionMask =&gt; highOrderAdvectionMask_in
+        highOrderAdvectionMask = 0
+      end if
+
+      if(present(boundaryCell_in)) then
+        boundaryCell =&gt; boundaryCell_in
+      else
+        allocate(boundaryCell(nVertLevels, grid % nCells+1))
+        boundaryCell(:,:) = 0
+      end if
+
+      do iEdge = 1, grid % nEdges
+        nAdvCellsForEdge(iEdge) = 0
+        cell1 = cellsOnEdge(1,iEdge)
+        cell2 = cellsOnEdge(2,iEdge)
+        
+
+        if(present(highOrderAdvectionMask_in)) then
+          do k = 1, nVertLevels
+            if (boundaryCell(k, cell1) == 1 .or. boundaryCell(k, cell2) == 1) then
+              highOrderAdvectionMask(k, iEdge) = 0
+            else
+              highOrderAdvectionMask(k, iEdge) = 1
+            end if
+          end do
+        end if
+
+        !
+        ! do only if this edge flux is needed to update owned cells
+        !
+        if (cell1 &lt;= grid % nCells .and. cell2 &lt;= grid % nCells) then
+           ! Insert cellsOnEdge to list of advection cells
+           call mpas_hash_init(cell_hash)
+           call mpas_hash_insert(cell_hash, cell1)
+           call mpas_hash_insert(cell_hash, cell2)
+           cell_indices(1) = cell1
+           cell_indices(2) = cell2
+           sorted_cell_indices(1, 1) = grid % indexToCellID % array(cell1)
+           sorted_cell_indices(2, 1) = cell1
+           sorted_cell_indices(1, 2) = grid % indexToCellID % array(cell2)
+           sorted_cell_indices(2, 2) = cell2
+           n = 2
+
+           ! Build unique list of cells used for advection on edge
+           do i = 1, nEdgesOnCell(cell1)
+             if(.not. mpas_hash_search(cell_hash, cellsOnCell(i, cell1))) then
+               n = n + 1
+               cell_indices(n) = cellsOnCell(i, cell1)
+               sorted_cell_indices(1, n) = grid % indexToCellID % array(cellsOnCell(i, cell1))
+               sorted_cell_indices(2, n) = cellsOnCell(i, cell1)
+               call mpas_hash_insert(cell_hash, cellsOnCell(i, cell1))
+             end if
+           end do ! loop over i
+
+           do i = 1, nEdgesOnCell(cell2)
+             if(.not. mpas_hash_search(cell_hash, cellsOnCell(i, cell2))) then
+               n = n + 1
+               cell_indices(n) = cellsOnCell(i, cell2)
+               sorted_cell_indices(1, n) = grid % indexToCellID % array(cellsOnCell(i, cell2))
+               sorted_cell_indices(2, n) = cellsOnCell(i, cell2)
+               call mpas_hash_insert(cell_hash, cellsOnCell(i, cell2))
+             end if
+           end do ! loop over i
+
+           call mpas_hash_destroy(cell_hash)
+
+           call mpas_quicksort(n, sorted_cell_indices)
+
+           nAdvCellsForEdge(iEdge) = n
+           do iCell = 1, nAdvCellsForEdge(iEdge)
+             advCellsForEdge(iCell, iEdge) = sorted_cell_indices(2, iCell)
+           end do ! loop over iCell
+
+           adv_coefs(:,iEdge) = 0.
+!          adv_coefs_2nd(:,iEdge) = 0.
+           adv_coefs_3rd(:,iEdge) = 0.
+
+           k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cell1))
+           if(k &lt;= nAdvCellsForEdge(iEdge)) then
+             adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + deriv_two(1,1,iEdge)
+             adv_coefs_3rd(k, iEdge) = adv_coefs_3rd(k, iEdge) + deriv_two(1,1,iEdge)
+           end if
+
+           do iCell = 1, nEdgesOnCell(cell1)
+             k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cellsOnCell(iCell,cell1)))
+             if(k &lt;= nAdvCellsForEdge(iEdge)) then
+               adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + deriv_two(iCell+1, 1, iEdge)
+               adv_coefs_3rd(k, iEdge) = adv_coefs_3rd(k, iEdge) + deriv_two(iCell+1, 1, iEdge)
+             end if
+           end do ! loop over iCell
+
+           k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cell2))
+           if(k &lt;= nAdvCellsForEdge(iEdge)) then
+             adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + deriv_two(1,2,iEdge)
+             adv_coefs_3rd(k, iEdge) = adv_coefs_3rd(k, iEdge) + deriv_two(1,2,iEdge)
+           end if
+
+           do iCell = 1, nEdgesOnCell(cell2)
+             k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cellsOnCell(iCell,cell2)))
+             if(k &lt;= nAdvCellsForEdge(iEdge)) then
+               adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + deriv_two(iCell+1, 2, iEdge)
+               adv_coefs_3rd(k, iEdge) = adv_coefs_3rd(k, iEdge) + deriv_two(iCell+1, 2, iEdge)
+             end if
+           end do ! loop over iCell
+
+           do iCell = 1,nAdvCellsForEdge(iEdge)
+             adv_coefs    (iCell,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs    (iCell,iEdge) / 12.
+             adv_coefs_3rd(iCell,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs_3rd(iCell,iEdge) / 12.
+           end do ! loop over iCell
+
+           k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cell1))
+           if(k &lt;= nAdvCellsForEdge(iEdge)) then
+             adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + 0.5
+             adv_coefs_2nd(k, iEdge) = adv_coefs_2nd(k, iEdge) + 0.5
+           end if
+
+           k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cell2))
+           if(k &lt;= nAdvCellsForEdge(iEdge)) then
+             adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + 0.5
+             adv_coefs_2nd(k, iEdge) = adv_coefs_2nd(k, iEdge) + 0.5
+           end if
+
+           do iCell=1,nAdvCellsForEdge(iEdge)
+             adv_coefs    (iCell,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs    (iCell,iEdge)
+             adv_coefs_2nd(iCell,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_2nd(iCell,iEdge)
+             adv_coefs_3rd(iCell,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_3rd(iCell,iEdge)
+           end do ! loop over iCell
+        end if
+      end do ! end loop over edges
+
+      deallocate(cell_indices)
+      deallocate(sorted_cell_indices)
+
+      ! If 2nd order advection, set masks appropriately.
+      if(config_horiz_tracer_adv_order == 2 .and. present(highOrderAdvectionMask_in)) then
+        highOrderAdvectionMask = 0
+      end if
+
+      if(.not.present(maxLevelCell_in)) then
+        deallocate(maxLevelCell)
+      end if
+
+      if(.not.present(boundaryCell_in)) then
+        deallocate(boundaryCell)
+      end if
+
+   end subroutine mpas_tracer_advection_coefficients!}}}
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  function mpas_tracer_advection_vflux4
+!
+!&gt; \brief MPAS 4th order vertical tracer advection stencil
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This function provides the stencil for 4th order vertical advection of tracers.
+!
+!-----------------------------------------------------------------------
+   real function mpas_tracer_advection_vflux4(q_im2, q_im1, q_i, q_ip1, w)!{{{
+        real (kind=RKIND), intent(in) :: q_im2 !&lt; Input: Tracer value at index i-2
+        real (kind=RKIND), intent(in) :: q_im1 !&lt; Input: Tracer value at index i-1
+        real (kind=RKIND), intent(in) :: q_i !&lt; Input: Tracer value at index i
+        real (kind=RKIND), intent(in) :: q_ip1 !&lt; Input: Tracer value at index i+1
+        real (kind=RKIND), intent(in) :: w !&lt; Input: vertical veloicity
+        mpas_tracer_advection_vflux4 = w*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+   end function!}}}
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  function mpas_tracer_advection_vflux3
+!
+!&gt; \brief MPAS 3rd order vertical tracer advection stencil
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This function provides the stencil for 3rd order vertical advection of tracers.
+!
+!-----------------------------------------------------------------------
+   real function mpas_tracer_advection_vflux3( q_im2, q_im1, q_i, q_ip1, w, coef)!{{{
+        real (kind=RKIND), intent(in) :: q_im2 !&lt; Input: Tracer value at index i-2
+        real (kind=RKIND), intent(in) :: q_im1 !&lt; Input: Tracer value at index i-1
+        real (kind=RKIND), intent(in) :: q_i !&lt; Input: Tracer value at index i
+        real (kind=RKIND), intent(in) :: q_ip1 !&lt; Input: Tracer value at index i+1
+        real (kind=RKIND), intent(in) :: w !&lt; Input: vertical veloicity
+        real (kind=RKIND), intent(in) :: coef !&lt; Input: Advection coefficient
+
+        !dwj 02/21/12 flux3 is different in and atmosphere
+        if(config_dzdk_positive) then
+          mpas_tracer_advection_vflux3 = (w * (7.0 * (q_i + q_im1) - (q_ip1 + q_im2)) + coef * abs(w) * ((q_ip1 - q_im2) - 3.0*(q_i-q_im1)))/12.0
+        else
+          mpas_tracer_advection_vflux3 = (w * (7.0 * (q_i + q_im1) - (q_ip1 + q_im2)) - coef * abs(w) * ((q_ip1 - q_im2) - 3.0*(q_i-q_im1)))/12.0
+        end if
+   end function!}}}
+
+end module mpas_tracer_advection_helpers

Added: branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_mono.F                                (rev 0)
+++ branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_mono.F        2012-08-27 17:16:34 UTC (rev 2126)
@@ -0,0 +1,475 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  mpas_tracer_advection_mono
+!
+!&gt; \brief MPAS monotonic tracer advection with FCT
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This module contains routines for monotonic advection of tracers using a FCT
+!
+!-----------------------------------------------------------------------
+module mpas_tracer_advection_mono
+
+   use mpas_kind_types
+   use mpas_grid_types
+   use mpas_configure
+   use mpas_dmpar
+
+   use mpas_tracer_advection_helpers
+
+   implicit none
+   private
+   save 
+
+   real (kind=RKIND) :: coef_3rd_order
+   logical :: vert2ndOrder, vert3rdOrder, vert4thOrder
+
+   public :: mpas_tracer_advection_mono_tend, &amp;
+             mpas_tracer_advection_mono_init
+
+   contains
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  routine mpas_tracer_advection_mono_tend
+!
+!&gt; \brief MPAS monotonic tracer advection tendency with FCT
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This routine computes the monotonic tracer advection tendencity using a FCT.
+!&gt;  Both horizontal and vertical.
+!
+!-----------------------------------------------------------------------
+   subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+   !
+   ! Input: s - current model state
+   !        grid - grid metadata
+   !
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+
+      real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !&lt; Input: current tracer values
+      real (kind=RKIND), dimension(:,:), intent(in) :: uh !&lt; Input: Thichness weighted velocitiy
+      real (kind=RKIND), dimension(:,:), intent(in) :: w !&lt; Input: Vertical velocitiy
+      real (kind=RKIND), dimension(:,:), intent(in) :: h !&lt; Input: Thickness
+      real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !&lt; Input: Distance between vertical interfaces of a cell
+      real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !&lt; Input: Tendency for thickness field
+      real (kind=RKIND), intent(in) :: dt !&lt; Input: Timestep
+      type (mesh_type), intent(in) :: grid !&lt; Input: Grid information
+      real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !&lt; Input/Output: Tracer tendency
+      integer, dimension(:), pointer, optional :: maxLevelCell_in !&lt; Input - optional: Index to max level at cell center
+      integer, dimension(:), pointer, optional :: maxLevelEdgeTop_in !&lt; Input - optional: Index to max level at edge with non-land cells on both sides
+      integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !&lt; Input- optional: Mask for high order advection
+
+      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, highOrderAdvectionMask, lowOrderAdvectionMask
+
+      real (kind=RKIND) :: 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) :: cur_max, cur_min, new_max, new_min
+      real (kind=RKIND) :: verticalWeightK, verticalWeightKm1
+      real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+      real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
+      real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, tracer_new, 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-10
+
+      type (field2dReal), pointer :: high_order_horiz_flux_field
+
+      ! Get dimensions
+      nCells = grid % nCells
+      nCellsSolve = grid % nCellsSolve
+      nEdges = grid % nEdges
+      nVertLevels = grid % nVertLevels
+      num_tracers = size(tracers,dim=1)
+
+      ! Initialize pointers
+      dvEdge      =&gt; grid % dvEdge % array
+      cellsOnEdge =&gt; grid % cellsOnEdge % array
+      cellsOnCell =&gt; grid % cellsOnCell % array
+      areaCell    =&gt; grid % areaCell % array
+
+      nEdgesOnCell =&gt; grid % nEdgesOnCell % array
+      nAdvCellsForEdge =&gt; grid % nAdvCellsForEdge % array
+      advCellsForEdge =&gt; grid % advCellsForEdge % array
+      adv_coefs =&gt; grid % adv_coefs % array
+      adv_coefs_3rd =&gt; grid % adv_coefs_3rd % array
+      
+      ! Setup arrays from optional arguments
+      if(present(maxLevelCell_in)) then
+        maxLevelCell =&gt; maxLevelCell_in
+      else
+        allocate(maxLevelCell(nCells+1))
+        maxLevelCell(:) = nVertLevels
+      end if
+
+      if(present(maxLevelEdgeTop_in)) then
+        maxLevelEdgeTop =&gt; maxLevelEdgeTop_in
+      else
+        allocate(maxLevelEdgeTop(nEdges+1))
+        maxLevelEdgeTop(:) = nVertLevels
+      end if
+
+      if(present(highOrderAdvectionMask_in)) then
+        highOrderAdvectionMask =&gt; highOrderAdvectionMask_in
+      else
+        allocate(highOrderAdvectionMask(nVertLevels, nEdges+1))
+        if(config_horiz_tracer_adv_order == 2) then
+          highOrderAdvectionMask = 0
+        else
+          highOrderAdvectionMask = 1
+        end if
+
+      end if
+
+      ! Setup high order horizontal flux field
+      allocate(high_order_horiz_flux_field)
+      nullify(high_order_horiz_flux_field % next)
+      high_order_horiz_flux_field % block =&gt; grid % block
+      high_order_horiz_flux_field % sendList =&gt; grid % xEdge % sendList
+      high_order_horiz_flux_field % recvList =&gt; grid % xEdge % recvList
+      high_order_horiz_flux_field % copyList =&gt; grid % xEdge % copyList
+      high_order_horiz_flux_field % dimSizes(1) = nVertLevels
+      high_order_horiz_flux_field % dimSizes(2) = nEdges+1
+      allocate(high_order_horiz_flux_field % array(high_order_horiz_flux_field % dimSizes(1), high_order_horiz_flux_field % dimSizes(2)))
+      high_order_horiz_flux =&gt; high_order_horiz_flux_field % array
+
+      ! allocate nCells arrays
+      allocate(tracer_new(nVertLevels, nCells+1))
+      allocate(tracer_cur(nVertLevels, nCells+1))
+      allocate(upwind_tendency(nVertLevels, nCells+1))
+      allocate(inv_h_new(nVertLevels, nCells+1))
+      allocate(tracer_max(nVertLevels, nCells+1))
+      allocate(tracer_min(nVertLevels, nCells+1))
+      allocate(flux_incoming(nVertLevels, nCells+1))
+      allocate(flux_outgoing(nVertLevels, nCells+1))
+
+      ! allocate nVertLevels+1 and nCells arrays
+      allocate(high_order_vert_flux(nVertLevels+1, nCells))
+
+      do iCell = 1, nCells
+        do k=1, maxLevelCell(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)
+            upwind_tendency(k, iCell) = 0.0
+
+            !tracer_new is supposed to be the &quot;new&quot; tracer state. This allows bounds checks.
+            if (config_check_monotonicity) then
+              tracer_new(k,iCell) = 0.0
+            end if
+          end do ! k loop
+        end do ! iCell loop
+
+        high_order_vert_flux = 0.0
+        high_order_horiz_flux = 0.0
+
+        !  Compute the high order vertical flux. Also determine bounds on tracer_cur. 
+        do iCell = 1, nCells
+
+          k = 1
+          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 = max(1, min(maxLevelCell(iCell), 2))
+          verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*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
+             if(vert4thOrder) then
+               high_order_vert_flux(k, iCell) = mpas_tracer_advection_vflux4( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell),  &amp;
+                                      tracer_cur(k  ,iCell),tracer_cur(k+1,iCell), w(k,iCell))
+             else if(vert3rdOrder) then
+               high_order_vert_flux(k, iCell) = mpas_tracer_advection_vflux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell),  &amp;
+                                      tracer_cur(k  ,iCell),tracer_cur(k+1,iCell), w(k,iCell), coef_3rd_order )
+             else if (vert2ndOrder) then
+               verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+               verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+               high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
+             end if
+             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 = max(1, maxLevelCell(iCell))
+          verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*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))
+
+          ! pull tracer_min and tracer_max from the (horizontal) surrounding cells
+          do i = 1, nEdgesOnCell(iCell)
+            do k=1, maxLevelCell(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
+
+        !  Compute the high order horizontal flux
+        do iEdge = 1, nEdges
+          cell1 = cellsOnEdge(1, iEdge)
+          cell2 = cellsOnEdge(2, iEdge)
+
+          ! Compute 2nd order fluxes where needed.
+          do k = 1, maxLevelEdgeTop(iEdge)
+            tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5) * uh(k, iEdge)
+
+            high_order_horiz_flux(k, iEdge) = high_order_horiz_flux(k, iedge) + tracer_weight * (tracer_cur(k, cell1) + tracer_cur(k, cell2))
+          end do ! k loop
+
+          ! Compute 3rd or 4th fluxes where requested.
+          do i = 1, nAdvCellsForEdge(iEdge)
+            iCell = advCellsForEdge(i,iEdge)
+            do k = 1, maxLevelCell(iCell)
+              tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+
+              tracer_weight = uh(k,iEdge)*tracer_weight
+              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
+
+        call mpas_dmpar_exch_halo_field(high_order_horiz_flux_field)
+
+        !  low order upwind vertical flux (monotonic and diffused)
+        !  Remove low order flux from the high order flux.
+        !  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 02/03/12 and Atmosphere are different in vertical
+            if(config_dzdk_positive) then
+              flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
+            else
+              flux_upwind = min(0.,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.,w(k,iCell))*tracer_cur(k,iCell)
+            end if
+            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
+
+          ! flux_incoming contains the total remaining high order flux into iCell
+          !          it is positive.
+          ! flux_outgoing contains the total remaining high order flux out of iCell
+          !           it is negative
+          do k = 1, maxLevelCell(iCell)
+            ! dwj 02/03/12 and Atmosphere are different in vertical
+            if(config_dzdk_positive) then
+              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)))
+            else
+              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 if
+          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 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))
+            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
+            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
+        ! 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)
+            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 = (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 = (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
+
+        !  rescale the high order horizontal fluxes
+        do iEdge = 1, nEdges
+          cell1 = cellsOnEdge(1,iEdge)
+          cell2 = cellsOnEdge(2,iEdge)
+          do k = 1, maxLevelEdgeTop(iEdge)
+            flux = high_order_horiz_flux(k,iEdge)
+            flux = max(0.,flux) * min(flux_outgoing(k,cell1), flux_incoming(k,cell2)) &amp;
+                 + 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 =  high_order_vert_flux(k,iCell)
+            ! dwj 02/03/12 and Atmosphere are different in vertical.
+            if(config_dzdk_positive) then
+              flux = max(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k  ,iCell)) &amp;
+                   + min(0.,flux) * min(flux_outgoing(k  ,iCell), flux_incoming(k-1,iCell))
+            else
+              flux = max(0.,flux) * min(flux_outgoing(k  ,iCell), flux_incoming(k-1,iCell)) &amp;
+                   + min(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k  ,iCell))
+            end if
+            high_order_vert_flux(k,iCell) = flux
+          end do ! k loop
+        end do ! iCell loop
+
+        ! Accumulate the scaled high order horizontal tendencies
+        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) - high_order_horiz_flux(k, iEdge) * invAreaCell1
+            tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+
+            if (config_check_monotonicity) then
+              !tracer_new holds a tendency for now.
+              tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
+              tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+            end if
+          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) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
+
+            if (config_check_monotonicity) then
+              !tracer_new holds a tendency for now. Only for a check on monotonicity
+              tracer_new(k, iCell) = tracer_new(k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
+
+              !tracer_new is now the new state of the tracer. Only for a check on monotonicity
+              tracer_new(k, iCell) = (tracer_cur(k, iCell)*h(k, iCell) + dt * tracer_new(k, iCell)) * inv_h_new(k, iCell)
+            end if
+          end do ! k loop
+        end do ! iCell loop
+
+        if (config_check_monotonicity) then
+          !build min and max bounds on old and new tracer for check on monotonicity.
+          do iCell = 1, nCellsSolve
+            do k = 1, maxLevelCell(iCell)
+              if(tracer_new(k,iCell) &lt; tracer_min(k, iCell)-eps) then
+                write(*,*) 'Minimum out of bounds on tracer ', iTracer, tracer_min(k, iCell), tracer_new(k,iCell)
+              end if
+
+              if(tracer_new(k,iCell) &gt; tracer_max(k,iCell)+eps) then
+                write(*,*) 'Maximum out of bounds on tracer ', iTracer, tracer_max(k, iCell), tracer_new(k,iCell)
+              end if
+            end do
+          end do
+        end if
+      end do ! iTracer loop
+
+      deallocate(tracer_new)
+      deallocate(tracer_cur)
+      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)
+      deallocate(high_order_horiz_flux_field)
+
+      if(.not.present(maxLevelCell_in)) then
+        deallocate(maxLevelCell)
+      end if
+
+      if(.not.present(maxLevelEdgeTop_in)) then
+        deallocate(maxLevelEdgeTop)
+      end if
+
+      if(.not.present(highOrderAdvectionMask_in)) then
+        deallocate(highOrderAdvectionMask)
+      end if
+   end subroutine mpas_tracer_advection_mono_tend!}}}
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  routine mpas_tracer_advection_mono_init
+!
+!&gt; \brief MPAS initialize monotonic tracer advection tendency with FCT
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This routine initializes the monotonic tracer advection tendencity using a FCT.
+!
+!-----------------------------------------------------------------------
+   subroutine mpas_tracer_advection_mono_init(err)!{{{
+      integer, intent(inout) :: err !&lt; Input: Error Flags
+
+      integer :: err_tmp
+
+      err = 0
+
+      vert2ndOrder = .false.
+      vert3rdOrder = .false.
+      vert4thOrder = .false.
+
+      if ( config_horiz_tracer_adv_order == 3) then
+          coef_3rd_order = config_coef_3rd_order
+      else if(config_horiz_tracer_adv_order == 2 .or. config_horiz_tracer_adv_order == 4) then
+          coef_3rd_order = 0.0
+      end if
+
+      if (config_vert_tracer_adv_order == 3) then
+          vert3rdOrder = .true.
+      else if (config_vert_tracer_adv_order == 4) then
+          vert4thOrder = .true.
+      else
+          vert2ndOrder = .true.
+          if(config_vert_tracer_adv_order /= 2) then
+            write(6,*) 'Invalid value for config_vert_tracer_adv_order, defaulting to 2nd order'
+          end if
+      end if
+
+
+   end subroutine mpas_tracer_advection_mono_init!}}}
+
+end module mpas_tracer_advection_mono

Added: branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_std.F                                (rev 0)
+++ branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_std.F        2012-08-27 17:16:34 UTC (rev 2126)
@@ -0,0 +1,295 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  mpas_tracer_advection_std
+!
+!&gt; \brief MPAS non-monotonic tracer advection
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This module contains routines for non-monotonic advection of tracers
+!
+!-----------------------------------------------------------------------
+module mpas_tracer_advection_std
+
+   use mpas_kind_types
+   use mpas_grid_types
+   use mpas_configure
+   use mpas_dmpar
+
+   use mpas_tracer_advection_helpers
+
+   implicit none
+   private
+   save 
+
+   real (kind=RKIND) :: coef_3rd_order
+   logical :: vert2ndOrder, vert3rdOrder, vert4thOrder
+
+   integer :: highOrderHorizOn
+
+   public :: mpas_tracer_advection_std_tend, &amp;
+             mpas_tracer_advection_std_init
+
+   contains
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  routine mpas_tracer_advection_std_tend
+!
+!&gt; \brief MPAS non-monotonic tracer advection tendency
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This routine computes the tracer advection tendencity.
+!&gt;  Both horizontal and vertical.
+!
+!-----------------------------------------------------------------------
+   subroutine mpas_tracer_advection_std_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+   !
+   ! Input: s - current model state
+   !        grid - grid metadata
+   !
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+
+      real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !&lt; Input: current tracer values
+      real (kind=RKIND), dimension(:,:), intent(in) :: uh !&lt; Input: Thichness weighted velocitiy
+      real (kind=RKIND), dimension(:,:), intent(in) :: w !&lt; Input: Vertical velocitiy
+      real (kind=RKIND), dimension(:,:), intent(in) :: h !&lt; Input: Thickness
+      real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !&lt; Input: Distance between vertical interfaces of a cell
+      real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !&lt; Input: Tendency for thickness field
+      real (kind=RKIND), intent(in) :: dt !&lt; Input: Timestep
+      type (mesh_type), intent(in) :: grid !&lt; Input: Grid information
+      real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !&lt; Input/Output: Tracer tendency
+      integer, dimension(:), pointer, optional :: maxLevelCell_in !&lt; Input - optional: Index to max level at cell center
+      integer, dimension(:), pointer, optional :: maxLevelEdgeTop_in !&lt; Input - optional: Index to max level at edge with non-land cells on both sides
+      integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !&lt; Input- optional: Mask for high order advection
+
+      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, highOrderAdvectionMask, lowOrderAdvectionMask
+
+      real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
+      real (kind=RKIND) :: verticalWeightK, verticalWeightKm1
+      real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+      real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
+      real (kind=RKIND), dimension(:,:), pointer :: high_order_horiz_flux, high_order_vert_flux, tracer_cur
+
+      real (kind=RKIND), parameter :: eps = 1.e-10
+
+      type (field2dReal), pointer :: high_order_horiz_flux_field
+
+      ! Get dimensions
+      nCells = grid % nCells
+      nCellsSolve = grid % nCellsSolve
+      nEdges = grid % nEdges
+      nVertLevels = grid % nVertLevels
+      num_tracers = size(tracers,dim=1)
+
+      ! Initialize pointers
+      dvEdge      =&gt; grid % dvEdge % array
+      cellsOnEdge =&gt; grid % cellsOnEdge % array
+      cellsOnCell =&gt; grid % cellsOnCell % array
+      areaCell    =&gt; grid % areaCell % array
+
+      nEdgesOnCell =&gt; grid % nEdgesOnCell % array
+      nAdvCellsForEdge =&gt; grid % nAdvCellsForEdge % array
+      advCellsForEdge =&gt; grid % advCellsForEdge % array
+      adv_coefs =&gt; grid % adv_coefs % array
+      adv_coefs_3rd =&gt; grid % adv_coefs_3rd % array
+      
+      ! Setup arrays from optional arguments
+      if(present(maxLevelCell_in)) then
+        maxLevelCell =&gt; maxLevelCell_in
+      else
+        allocate(maxLevelCell(nCells+1))
+        maxLevelCell(:) = nVertLevels
+      end if
+
+      if(present(maxLevelEdgeTop_in)) then
+        maxLevelEdgeTop =&gt; maxLevelEdgeTop_in
+      else
+        allocate(maxLevelEdgeTop(nEdges+1))
+        maxLevelEdgeTop(:) = nVertLevels
+      end if
+
+      if(present(highOrderAdvectionMask_in)) then
+        highOrderAdvectionMask =&gt; highOrderAdvectionMask_in
+      else
+        allocate(highOrderAdvectionMask(nVertLevels, nEdges+1))
+        if(config_horiz_tracer_adv_order == 2) then
+          highOrderAdvectionMask = 0
+        else
+          highOrderAdvectionMask = 1
+        end if
+
+      end if
+
+      ! Setup high order horizontal flux field
+      allocate(high_order_horiz_flux_field)
+      nullify(high_order_horiz_flux_field % next)
+      high_order_horiz_flux_field % block =&gt; grid % block
+      high_order_horiz_flux_field % sendList =&gt; grid % xEdge % sendList
+      high_order_horiz_flux_field % recvList =&gt; grid % xEdge % recvList
+      high_order_horiz_flux_field % copyList =&gt; grid % xEdge % copyList
+      high_order_horiz_flux_field % dimSizes(1) = nVertLevels
+      high_order_horiz_flux_field % dimSizes(2) = nEdges+1
+      allocate(high_order_horiz_flux_field % array(high_order_horiz_flux_field % dimSizes(1), high_order_horiz_flux_field % dimSizes(2)))
+      high_order_horiz_flux =&gt; high_order_horiz_flux_field % array
+
+      ! allocate nVertLevels+1 and nCells arrays
+      allocate(high_order_vert_flux(nVertLevels+1, nCells))
+
+      !allocate nVertLevels and nCells+1 arrays
+      allocate(tracer_cur(nVertLevels, nCells+1))
+
+      ! 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
+        high_order_vert_flux = 0.0
+        high_order_horiz_flux = 0.0
+        tracer_cur(:,:) = tracers(iTracer, :, :)
+
+        !  Compute the high order vertical flux. Also determine bounds on tracer_cur. 
+        do iCell = 1, nCells
+          k = max(1, min(maxLevelCell(iCell), 2))
+          verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
+             
+          do k=3,maxLevelCell(iCell)-1
+             if(vert4thOrder) then
+               high_order_vert_flux(k, iCell) = mpas_tracer_advection_vflux4( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell),  &amp;
+                                      tracer_cur(k  ,iCell),tracer_cur(k+1,iCell), w(k,iCell))
+             else if(vert3rdOrder) then
+               high_order_vert_flux(k, iCell) = mpas_tracer_advection_vflux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell),  &amp;
+                                      tracer_cur(k  ,iCell),tracer_cur(k+1,iCell), w(k,iCell), coef_3rd_order )
+             else if (vert2ndOrder) then
+               verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+               verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+               high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
+             end if
+          end do 
+ 
+          k = max(1, maxLevelCell(iCell))
+          verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+          high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
+        end do ! iCell Loop
+
+        !  Compute the high order horizontal flux
+        do iEdge = 1, nEdges
+          cell1 = cellsOnEdge(1, iEdge)
+          cell2 = cellsOnEdge(2, iEdge)
+
+          ! Compute 2nd order fluxes where needed.
+          do k = 1, maxLevelEdgeTop(iEdge)
+            tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5) * uh(k, iEdge)
+
+            high_order_horiz_flux(k, iEdge) = high_order_horiz_flux(k, iedge) + tracer_weight * (tracer_cur(k, cell1) + tracer_cur(k, cell2))
+          end do ! k loop
+
+          ! Compute 3rd or 4th fluxes where requested.
+          do i = 1, nAdvCellsForEdge(iEdge) * highOrderHorizOn
+            iCell = advCellsForEdge(i,iEdge)
+            do k = 1, maxLevelCell(iCell)
+              tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+
+              tracer_weight = uh(k,iEdge)*tracer_weight
+              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
+
+        ! Accumulate the scaled high order horizontal tendencies
+        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) - 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) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell))
+          end do ! k loop
+        end do ! iCell loop
+
+      end do ! iTracer loop
+
+      deallocate(high_order_horiz_flux)
+      deallocate(high_order_vert_flux)
+      deallocate(high_order_horiz_flux_field)
+      deallocate(tracer_cur)
+
+      if(.not.present(maxLevelCell_in)) then
+        deallocate(maxLevelCell)
+      end if
+
+      if(.not.present(maxLevelEdgeTop_in)) then
+        deallocate(maxLevelEdgeTop)
+      end if
+
+      if(.not.present(highOrderAdvectionMask_in)) then
+        deallocate(highOrderAdvectionMask)
+      end if
+   end subroutine mpas_tracer_advection_std_tend!}}}
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  routine mpas_tracer_advection_std_init
+!
+!&gt; \brief MPAS initialize tracer advection tendency.
+!&gt; \author Doug Jacobsen
+!&gt; \date   03/09/12
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This routine initializes the tracer advection tendencity module.
+!
+!-----------------------------------------------------------------------
+   subroutine mpas_tracer_advection_std_init(err)!{{{
+      integer, intent(inout) :: err !&lt; Input: Error Flags
+
+      integer :: err_tmp
+
+      err = 0
+
+      vert2ndOrder = .false.
+      vert3rdOrder = .false.
+      vert4thOrder = .false.
+      highOrderHorizOn = 0
+
+      if ( config_horiz_tracer_adv_order == 3 .or. config_horiz_tracer_adv_order == 4) then
+        highOrderHorizOn = 1
+      end if
+
+      if ( config_horiz_tracer_adv_order == 3) then
+          coef_3rd_order = config_coef_3rd_order
+      else if(config_horiz_tracer_adv_order == 2 .or. config_horiz_tracer_adv_order == 4) then
+          coef_3rd_order = 0.0
+      end if
+
+      if (config_vert_tracer_adv_order == 3) then
+          vert3rdOrder = .true.
+      else if (config_vert_tracer_adv_order == 4) then
+          vert4thOrder = .true.
+      else
+          vert2ndOrder = .true.
+          if(config_vert_tracer_adv_order /= 2) then
+            write(6,*) 'Invalid value for config_vert_tracer_adv_order, defaulting to 2nd order'
+          end if
+      end if
+
+
+   end subroutine mpas_tracer_advection_std_init!}}}
+
+end module mpas_tracer_advection_std

</font>
</pre>