<p><b>dwj07@fsu.edu</b> 2012-02-10 22:04:46 -0700 (Fri, 10 Feb 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding shared tracer advection routines from previous advection_routines branch.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/monotonic_advection/namelist.input.ocean
===================================================================
--- branches/ocean_projects/monotonic_advection/namelist.input.ocean 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/namelist.input.ocean 2012-02-11 05:04:46 UTC (rev 1500)
@@ -80,6 +80,7 @@
config_vert_tracer_adv = 'stencil'
config_vert_tracer_adv_order = 2
config_tracer_adv_order = 2
+ config_tracer_vadv_order = 2
config_thickness_adv_order = 2
config_positive_definite = .false.
config_monotonic = .false.
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-02-11 05:04:46 UTC (rev 1500)
@@ -69,8 +69,10 @@
namelist character eos config_eos_type linear
namelist character advection config_vert_tracer_adv stencil
namelist integer advection config_vert_tracer_adv_order 4
+namelist integer advection config_tracer_vadv_order 2
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
+namelist real advection config_coef_3rd_order 1.0
namelist logical advection config_positive_definite false
namelist logical advection config_monotonic false
namelist logical restore config_restoreTS false
@@ -154,6 +156,12 @@
var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 - deriv_two mesh - -
var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
+% Added for monotonic advection scheme
+var persistent real adv_coefs ( FIFTEEN nEdges ) 0 - adv_coefs mesh - -
+var persistent real adv_coefs_3rd ( FIFTEEN nEdges ) 0 - adv_coefs_3rd mesh - -
+var persistent integer advCellsForEdge ( FIFTEEN nEdges ) 0 - advCellsForEdge mesh - -
+var persistent integer nAdvCellsForEdge ( nEdges ) 0 - nAdvCellsForEdge mesh - -
+
% !! NOTE: the following arrays are needed to allow the use
% !! of the module_advection.F w/o alteration
% Space needed for deformation calculation weights
@@ -172,6 +180,9 @@
var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
+var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
+var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
+var persistent real hRatioZLevelKm1 ( nVertLevels ) 0 - hRatioZLevelKm1 mesh - -
% Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -45,6 +45,7 @@
subroutine mpas_core_init(domain, startTimeStamp)!{{{
use mpas_grid_types
+ use mpas_tracer_advection
implicit none
@@ -88,6 +89,9 @@
call ocn_tendency_init(err_tmp)
err = ior(err,err_tmp)
+ call mpas_tracer_advection_init(err_tmp)
+ err = ior(err,err_tmp)
+
call mpas_timer_init(domain)
if(err.eq.1) then
@@ -218,6 +222,7 @@
use mpas_grid_types
use mpas_rbf_interpolation
use mpas_vector_reconstruction
+ use mpas_tracer_advection
implicit none
@@ -226,6 +231,7 @@
real (kind=RKIND), intent(in) :: dt
integer :: i, iEdge, iCell, k
+ call mpas_tracer_advection_coefficients(mesh)
call ocn_time_average_init(block % state % time_levs(1) % state)
call mpas_timer_start("diagnostic solve", .false., initDiagSolveTimer)
@@ -508,7 +514,7 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -519,6 +525,9 @@
h => block % state % time_levs(1) % state % h % array
referenceBottomDepth => block % mesh % referenceBottomDepth % array
nVertLevels = block % mesh % nVertLevels
+ hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
+ hRatioZLevelK => block % mesh % hRatioZLevelK % array
+ hRatioZLevelKm1 => block % mesh % hRatioZLevelKm1 % array
! mrp 120208 right now hZLevel is in the grid.nc file.
! We would like to transition to using referenceBottomDepth
@@ -526,10 +535,17 @@
! When the transition is done, hZLevel can be removed from
! registry and the following four lines deleted.
referenceBottomDepth(1) = block % mesh % hZLevel % array(1)
+ hMeanTopZLevel(1) = 0.0
+ hRatioZLevelK(1) = 0.0
+ hRatioZLevelKm1(1) = 0.0
do k = 2,nVertLevels
referenceBottomDepth(k) = referenceBottomDepth(k-1) + block % mesh % hZLevel % array(k)
+ hMeanTopZLevel(k) = 0.5*(block % mesh % hZLevel % array(k-1) + block % mesh % hZLevel % array(k))
+ hRatioZLevelK(k) = 0.5*block % mesh % hZLevel % array(k)/hMeanTopZLevel(k)
+ hRatioZLevelKm1(k) = 0.5*block % mesh % hZLevel % array(k-1)/hMeanTopZLevel(k)
end do
+
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
if (trim(config_time_integration) == 'unsplit_explicit') then
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -21,6 +21,8 @@
use mpas_constants
use mpas_timer
+ use mpas_tracer_advection
+
use ocn_thick_hadv
use ocn_thick_vadv
use ocn_gm
@@ -289,7 +291,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_scalar(tend, s, d, grid)!{{{
+ subroutine ocn_tend_scalar(tend, s, d, grid, dt)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -306,26 +308,41 @@
type (state_type), intent(in) :: s
type (diagnostics_type), intent(in) :: d
type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), intent(in) :: dt
+ real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1, zTopZLevel, hMeanTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- u, uTransport, h,wTop, h_edge, vertDiffTopOfCell
+ u, uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
- integer :: err
+ integer :: err, iEdge, k
call mpas_timer_start("ocn_tend_scalar")
u => s % u % array
- uTransport => s % uTransport % array
+ uTransport => s % uTransport % array
h => s % h % array
wTop => s % wTop % array
tracers => s % tracers % array
h_edge => s % h_edge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
+ hRatioZLevelK => grid % hRatioZLevelK % array
+ hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
+ hMeanTopZLevel => grid % hMeanTopZLevel % array
+
tend_tr => tend % tracers % array
+ tend_h => tend % h % array
+ allocate(uh(grid % nVertLevels, grid % nEdges+1))
+
+ do iEdge = 1, grid % nEdges
+ do k = 1, grid % nVertLevels
+ uh(k, iEdge) = uTransport(k, iEdge) * h_edge(k, iEdge)
+ end do
+ end do
+
!
! initialize tracer tendency (RHS of tracer equation) to zero.
!
@@ -339,18 +356,15 @@
! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
! tracer_edge at the boundary will also need to be defined for flux boundaries.
- call mpas_timer_start("hadv", .false., tracerHadvTimer)
- call ocn_tracer_hadv_tend(grid, uTransport, h_edge, tracers, tend_tr, err)
- call mpas_timer_stop("hadv", tracerHadvTimer)
+ call mpas_timer_start("adv", .false., tracerHadvTimer)
+ ! Monotonoic Advection, or standard advection
+ call mpas_tracer_advection_tend(tracers, uh, wTop, h, dt, hMeanTopZLevel, hRatioZLevelKm1, hRatioZLevelK, grid, tend_h, tend_tr)
- !
- ! tracer tendency: vertical advection term -d/dz( h \phi w)
- !
+ ! Only standadrd advection
+! call ocn_tracer_hadv_tend(grid, uTransport, h_edge, tracers, tend_tr, err)
+! call ocn_tracer_vadv_tend(grid, h, wtop, tracers, tend_tr, err)
+ call mpas_timer_stop("adv", tracerHadvTimer)
- call mpas_timer_start("vadv", .false., tracerVadvTimer)
- call ocn_tracer_vadv_tend(grid, h, wtop, tracers, tend_tr, err)
- call mpas_timer_stop("vadv", tracerVadvTimer)
-
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -182,7 +182,7 @@
call filter_btr_mode_tend_u(block % tend, provis, block % diagnostics, block % mesh)
endif
- call ocn_tend_scalar(block % tend, provis, block % diagnostics, block % mesh)
+ call ocn_tend_scalar(block % tend, provis, block % diagnostics, block % mesh, dt)
block => block % next
end do
call mpas_timer_stop("RK4-tendency computations")
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -721,7 +721,7 @@
call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
call ocn_tend_h (block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
- call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh, dt)
block => block % next
end do
Modified: branches/ocean_projects/monotonic_advection/src/operators/Makefile
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/Makefile 2012-02-11 04:41:50 UTC (rev 1499)
+++ branches/ocean_projects/monotonic_advection/src/operators/Makefile 2012-02-11 05:04:46 UTC (rev 1500)
@@ -1,6 +1,17 @@
.SUFFIXES: .F .o
-OBJS = mpas_rbf_interpolation.o mpas_vector_reconstruction.o mpas_spline_interpolation.o
+OBJS = mpas_rbf_interpolation.o \
+ mpas_vector_reconstruction.o \
+ mpas_spline_interpolation.o \
+ mpas_tracer_advection.o \
+ mpas_tracer_advection_mono.o \
+ mpas_tracer_advection_std.o \
+ mpas_tracer_advection_std_hadv.o \
+ mpas_tracer_advection_std_vadv.o \
+ mpas_tracer_advection_std_vadv2.o \
+ mpas_tracer_advection_std_vadv3.o \
+ mpas_tracer_advection_std_vadv4.o \
+ mpas_tracer_advection_helpers.o
all: operators
@@ -8,9 +19,29 @@
ar -ru libops.a $(OBJS)
mpas_vector_reconstruction.o: mpas_rbf_interpolation.o
+
mpas_rbf_interpolation.o:
+
mpas_spline_interpolation:
+mpas_tracer_advection.o: mpas_tracer_advection_std.o mpas_tracer_advection_mono.o
+
+mpas_tracer_advection_mono.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_std.o: mpas_tracer_advection_std_hadv.o mpas_tracer_advection_std_vadv.o
+
+mpas_tracer_advection_std_hadv.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_std_vadv.o: mpas_tracer_advection_std_vadv2.o mpas_tracer_advection_std_vadv3.o mpas_tracer_advection_std_vadv4.o
+
+mpas_tracer_advection_std_vadv2.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_std_vadv3.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_std_vadv4.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_helpers.o:
+
clean:
$(RM) *.o *.mod *.f90 libops.a
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,216 @@
+module mpas_tracer_advection
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+
+ use mpas_tracer_advection_std
+ use mpas_tracer_advection_mono
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_init, &
+ mpas_tracer_advection_coefficients, &
+ mpas_tracer_advection_tend
+
+ logical :: monotonicOn
+
+ contains
+
+ subroutine mpas_tracer_advection_coefficients( grid )!{{{
+
+ implicit none
+ type (mesh_type) :: grid
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
+ integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge
+
+ integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
+ integer :: cell_list(20), ordered_cell_list(20)
+ logical :: addcell
+
+ deriv_two => grid % deriv_two % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % array
+ cellsOnCell => grid % cellsOnCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+
+ do iEdge = 1, grid % nEdges
+ nAdvCellsForEdge(iEdge) = 0
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ !
+ ! do only if this edge flux is needed to update owned cells
+ !
+ if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then
+
+ cell_list(1) = cell1
+ cell_list(2) = cell2
+ n = 2
+
+ ! add cells surrounding cell 1. n is number of cells currently in list
+ do i = 1, nEdgesOnCell(cell1)
+ if(cellsOnCell(i,cell1) /= cell2) then
+ n = n + 1
+ cell_list(n) = cellsOnCell(i,cell1)
+ end if
+ end do
+
+ ! add cells surrounding cell 2 (brute force approach)
+ do iCell = 1, nEdgesOnCell(cell2)
+ addcell = .true.
+ do i=1,n
+ if(cell_list(i) == cellsOnCell(iCell,cell2)) addcell = .false.
+ end do
+ if(addcell) then
+ n = n+1
+ cell_list(n) = cellsOnCell(iCell,cell2)
+ end if
+ end do
+
+ ! order the list by increasing cell number (brute force approach)
+
+ do i=1,n
+ ordered_cell_list(i) = grid % nCells + 2
+ j_in = 1
+ do j=1,n
+ if(ordered_cell_list(i) > cell_list(j) ) then
+ j_in = j
+ ordered_cell_list(i) = cell_list(j)
+ end if
+ end do
+! ordered_cell_list(i) = cell_list(j_in)
+ cell_list(j_in) = grid % nCells + 3
+ end do
+
+ nAdvCellsForEdge(iEdge) = n
+ do iCell = 1, nAdvCellsForEdge(iEdge)
+ advCellsForEdge(iCell,iEdge) = ordered_cell_list(iCell)
+ end do
+
+ ! we have the ordered list, now construct coefficients
+
+ adv_coefs(:,iEdge) = 0.
+ adv_coefs_3rd(:,iEdge) = 0.
+
+ ! pull together third and fourth order contributions to the flux
+ ! first from cell1
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell1 ) j_in = j
+ end do
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,1,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(1,1,iEdge)
+
+ do iCell = 1, nEdgesOnCell(cell1)
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cellsOnCell(iCell,cell1) ) j_in = j
+ end do
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
+ end do
+
+ ! pull together third and fourth order contributions to the flux
+ ! now from cell2
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell2 ) j_in = j
+ enddo
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,2,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(1,2,iEdge)
+
+ do iCell = 1, nEdgesOnCell(cell2)
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cellsOnCell(iCell,cell2) ) j_in = j
+ enddo
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,2,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(iCell+1,2,iEdge)
+ end do
+
+ do j = 1,n
+ adv_coefs (j,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs (j,iEdge) / 12.
+ adv_coefs_3rd(j,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs_3rd(j,iEdge) / 12.
+ end do
+
+ ! 2nd order centered contribution - place this in the main flux weights
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell1 ) j_in = j
+ enddo
+ adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell2 ) j_in = j
+ enddo
+ adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+
+ ! multiply by edge length - thus the flux is just dt*ru times the results of the vector-vector multiply
+
+ do j=1,n
+ adv_coefs (j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs (j,iEdge)
+ adv_coefs_3rd(j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_3rd(j,iEdge)
+ end do
+
+ end if ! only do for edges of owned-cells
+
+ end do ! end loop over edges
+
+ end subroutine mpas_tracer_advection_coefficients!}}}
+
+! subroutine mpas_tracer_advection_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, dt, dminfo, cellsToSend, cellsToRecv, tend)
+ subroutine mpas_tracer_advection_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)!{{{
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh, w, h
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), intent(in) :: dt
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), dimension(:,:), intent(in) :: tend_h
+! real (kind=RKIND) :: dt
+
+! type (dm_info), intent(in) :: dminfo
+! type (exchange_list), pointer :: cellsToSend, cellsToRecv
+
+ if(monotonicOn) then
+! call mpas_tracer_advection_mono_tend( tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)
+ call mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)
+ else
+ call mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ endif
+ end subroutine mpas_tracer_advection_tend!}}}
+
+ subroutine mpas_tracer_advection_init(err)!{{{
+
+ integer, intent(inout) :: err
+
+ integer :: err_tmp
+
+ err = 0
+
+ call mpas_tracer_advection_std_init(err_tmp)
+
+ err = ior(err, err_tmp)
+
+ monotonicOn = .false.
+
+ if(config_monotonic .and. config_positive_definite) then
+ monotonicOn = .true.
+ endif
+
+ end subroutine mpas_tracer_advection_init!}}}
+
+end module mpas_tracer_advection
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_helpers.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_helpers.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_helpers.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,21 @@
+module mpas_tracer_advection_helpers
+
+ use mpas_kind_types
+
+ implicit none
+ save
+
+ contains
+
+ real function flux4(q_im2, q_im1, q_i, q_ip1, u)!{{{
+ real (kind=RKIND), intent(in) :: q_im2, q_im1, q_i, q_ip1, u
+ flux4 = u*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+ end function!}}}
+
+ real function flux3( q_im2, q_im1, q_i, q_ip1, u, coef)!{{{
+ real (kind=RKIND), intent(in) :: q_im2, q_im1, q_i, q_ip1, u, coef
+
+ flux3 = (u * (7.0 * (q_i + q_im1) - (q_ip1 + q_im2)) + coef * abs(u) * ((q_ip1 - q_im2) - 3.0*(q_i-q_im1)))/12.0
+ end function!}}}
+
+end module mpas_tracer_advection_helpers
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,293 @@
+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
+
+ public :: mpas_tracer_advection_mono_tend
+
+ contains
+
+ subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ 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, 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, 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
+ areaCell => grid % areaCell % array
+
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % array
+ 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(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(high_order_horiz_flux(nVertLevels, nEdges))
+
+ ! 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_min(k, iCell) = tracer_cur(k, iCell)
+! tracer_max(k, iCell) = tracer_cur(k, iCell)
+ 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
+! 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
+ 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
+ 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 )
+ 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)
+ 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))
+
+! high_order_vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
+
+ ! 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
+ 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))
+ 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 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)
+ 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 Ocean and Atmosphere are different in vertical
+! 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)))
+
+ 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 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)) &
+ + 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 Ocean and Atmosphere are different in vertical.
+! 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
+
+ ! 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)) + upwind_tendency(k,iCell)
+ 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 ! k loop
+! end do ! iCell loop
+ end do ! iTracer loop
+
+ 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)
+ end subroutine mpas_tracer_advection_mono_tend!}}}
+
+end module mpas_tracer_advection_mono
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,59 @@
+module mpas_tracer_advection_std
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+ use mpas_timer
+
+ use mpas_tracer_advection_std_hadv
+ use mpas_tracer_advection_std_vadv
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_tend, &
+ mpas_tracer_advection_std_init
+
+ contains
+
+ subroutine mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed scalar tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh, w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ type (mesh_type), intent(in) :: grid
+
+ call mpas_timer_start("tracer-hadv", .false.)
+ call mpas_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)
+ call mpas_timer_stop("tracer-hadv")
+ call mpas_timer_start("tracer-vadv", .false.)
+ call mpas_tracer_advection_std_vadv_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ call mpas_timer_stop("tracer-vadv")
+
+ end subroutine mpas_tracer_advection_std_tend!}}}
+
+ subroutine mpas_tracer_advection_std_init(err)!{{{
+ integer, intent(inout) :: err
+
+ integer :: err_tmp
+
+ err = 0
+
+ call mpas_tracer_advection_std_hadv_init(err_tmp)
+ err = ior(err, err_tmp)
+ call mpas_tracer_advection_std_vadv_init(err_tmp)
+ err = ior(err, err_tmp)
+
+ end subroutine mpas_tracer_advection_std_init!}}}
+
+end module mpas_tracer_advection_std
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_hadv.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_hadv.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,105 @@
+module mpas_tracer_advection_std_hadv
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_hadv_tend, &
+ mpas_tracer_advection_std_hadv_init
+
+ real (kind=RKIND) :: coef_3rd_order
+
+ logical :: hadvOn
+
+ contains
+
+ subroutine mpas_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_weight
+
+ real (kind=RKIND), dimension(:), pointer :: 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( s_old % num_tracers, grid % nVertLevels ) :: flux_arr
+ real (kind=RKIND), dimension(:,:), allocatable :: flux_arr
+ integer :: nVertLevels, num_tracers
+
+ if (.not. hadvOn) return
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ areaCell => grid % areaCell % array
+
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % array
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+
+ allocate(flux_arr(num_tracers, grid % nVertLevels))
+
+ ! horizontal flux divergence, accumulate in tracer_tend
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
+ flux_arr(:,:) = 0.
+ do i=1,nAdvCellsForEdge(iEdge)
+ iCell = advCellsForEdge(i,iEdge)
+ do k=1,grid % nVertLevels
+ tracer_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0, uh(k,iEdge))*adv_coefs_3rd(i,iEdge)
+ do iTracer=1,num_tracers
+ flux_arr(iTracer,k) = flux_arr(iTracer,k) + tracer_weight* tracers(iTracer,k,iCell)
+ end do
+ end do
+ end do
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - uh(k,iEdge)*flux_arr(iTracer,k)/areaCell(cell1)
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + uh(k,iEdge)*flux_arr(iTracer,k)/areaCell(cell2)
+ end do
+ end do
+ end if
+ end do
+
+ deallocate(flux_arr)
+
+ end subroutine mpas_tracer_advection_std_hadv_tend!}}}
+
+ subroutine mpas_tracer_advection_std_hadv_init(err)!{{{
+ integer, intent(inout) :: err
+
+ err = 0
+
+ hadvOn =.true.
+
+ coef_3rd_order = config_coef_3rd_order
+ end subroutine mpas_tracer_advection_std_hadv_init!}}}
+
+end module mpas_tracer_advection_std_hadv
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,66 @@
+module mpas_tracer_advection_std_vadv
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_std_vadv2
+ use mpas_tracer_advection_std_vadv3
+ use mpas_tracer_advection_std_vadv4
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv_tend, &
+ mpas_tracer_advection_std_vadv_init
+
+ logical :: order2, order3, order4
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ if(order2) then
+ call mpas_tracer_advection_std_vadv2_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ else if(order3) then
+ call mpas_tracer_advection_std_vadv3_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ else if(order4) then
+ call mpas_tracer_advection_std_vadv4_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ endif
+
+ end subroutine mpas_tracer_advection_std_vadv_tend!}}}
+
+ subroutine mpas_tracer_advection_std_vadv_init(err)!{{{
+ integer, intent(inout) :: err
+
+ err = 0
+
+ order2 = .false.
+ order3 = .false.
+ order4 = .false.
+
+ if (config_tracer_vadv_order == 2) then
+ order2 = .true.
+ else if (config_tracer_vadv_order == 3) then
+ order3 = .true.
+ else if (config_tracer_vadv_order == 4) then
+ order4 = .true.
+ else
+ print *, 'invalid value for config_tracer_vadv_order'
+ print *, 'options are 2, 3, or 4'
+ err = 1
+ endif
+
+ end subroutine mpas_tracer_advection_std_vadv_init!}}}
+
+end module mpas_tracer_advection_std_vadv
+
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,81 @@
+module mpas_tracer_advection_std_vadv2
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv2_tend
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv2_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_edge, tracer_weight
+ real (kind=RKIND) :: tracer_weight_cell1, tracer_weight_cell2
+
+
+! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
+ real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ integer :: nVertLevels, num_tracers
+ integer, dimension(:), pointer :: maxLevelCell
+
+ integer, parameter :: hadv_opt = 2
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+ maxLevelCell => grid % maxLevelCell % array
+
+ allocate(vert_flux(num_tracers, nVertLevels+1))
+
+ !
+ ! vertical flux divergence
+ !
+
+ ! zero fluxes at top and bottom
+
+ vert_flux(:,1) = 0.
+
+ do iCell=1,grid % nCellsSolve
+ do k = 2, maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ end do
+ end do
+
+ vert_flux(:,maxLevelCell(iCell)+1) = 0
+
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + ( vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
+! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
+! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+ end do
+
+ deallocate(vert_flux)
+
+ end subroutine mpas_tracer_advection_std_vadv2_tend!}}}
+
+end module mpas_tracer_advection_std_vadv2
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,93 @@
+module mpas_tracer_advection_std_vadv3
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv3_tend
+
+ real (kind=RKIND) :: coef_3rd_order
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv3_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+
+! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
+ real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ integer :: nVertLevels, num_tracers
+ integer, dimension(:), pointer :: maxLevelCell
+
+ integer, parameter :: hadv_opt = 2
+
+ coef_3rd_order = config_coef_3rd_order
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+ maxLevelCell => grid % maxLevelCell % array
+
+ allocate(vert_flux(num_tracers, nVertLevels+1))
+
+ vert_flux(:,1) = 0.
+! vert_flux(:,grid % nVertLevels+1) = 0.
+
+ do iCell=1,grid % nCellsSolve
+
+ k = 2
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ enddo
+
+ do k=3,maxLevelCell(iCell)-1
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = flux3( tracers(iTracer,k-2,iCell),tracers(iTracer,k-1,iCell), &
+ tracers(iTracer,k ,iCell),tracers(iTracer,k+1,iCell), &
+ w(k,iCell), coef_3rd_order )
+ end do
+ end do
+
+ k = maxLevelCell(iCell)
+
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ enddo
+
+ vert_Flux(:, maxLevelCell(iCell)+1) = 0.0
+
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
+
+! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
+! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+ end do
+
+ deallocate(vert_flux)
+
+ end subroutine mpas_tracer_advection_std_vadv3_tend!}}}
+
+end module mpas_tracer_advection_std_vadv3
Added: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F (rev 0)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F 2012-02-11 05:04:46 UTC (rev 1500)
@@ -0,0 +1,88 @@
+module mpas_tracer_advection_std_vadv4
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv4_tend
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv4_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+
+! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
+ real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ integer :: nVertLevels, num_tracers
+ integer, dimension(:), pointer :: maxLevelCell
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+ maxLevelCell => grid % maxLevelCell % array
+
+ allocate(vert_flux(num_tracers, nVertLevels+1))
+
+ ! vertical flux divergence
+ !
+
+ ! zero fluxes at top and bottom
+
+ vert_flux(:,1) = 0.
+
+ do iCell=1,grid % nCellsSolve
+
+ k = 2
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ enddo
+ do k=3,nVertLevels-1
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = flux4( tracers(iTracer,k-2,iCell),tracers(iTracer,k-1,iCell), &
+ tracers(iTracer,k ,iCell),tracers(iTracer,k+1,iCell), w(k,iCell) )
+ end do
+ end do
+
+ k = maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ enddo
+
+ vert_flux(:,maxLevelCell(iCell)+1) = 0.0
+
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
+! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
+! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+
+ end do
+
+ deallocate(vert_flux)
+
+ end subroutine mpas_tracer_advection_std_vadv4_tend!}}}
+
+end module mpas_tracer_advection_std_vadv4
</font>
</pre>